- Email: [email protected]
Haze removal based on multiple scattering model with superpixel algorithm
Author’s Accepted Manuscript Haze removal based on multiple scattering model with superpixel algorithm Rui Wang, Rui Li, Haiyan Sun
www.elsevier.com/locate/sigpro
PII: DOI: Reference:
S0165-1684(16)00041-4 http://dx.doi.org/10.1016/j.sigpro.2016.02.003 SIGPRO6050
To appear in: Signal Processing Received date: 14 September 2015 Revised date: 10 January 2016 Accepted date: 1 February 2016 Cite this article as: Rui Wang, Rui Li and Haiyan Sun, Haze removal based on multiple scattering model with superpixel algorithm, Signal Processing, http://dx.doi.org/10.1016/j.sigpro.2016.02.003 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Haze removal based on multiple scattering model with superpixel algorithm Rui Wanga, Rui Lia, Haiyan Sunb a. Key Laboratory of Precision Opto-mechatronics Technology, Ministry of Education, School of Instrumentation Science and Opto-electronics Engineering, Beihang University, Beijing, 100191, China b. School of Mathematics and System Science, Beihang University, Beijing, 100191 China Abstract: An efficient method to remove haze from single image based on dark channel prior and the multiple scattering described by atmosphere point spread function (APSF) is proposed in this paper. By means of statistical analysis of the quality assessments for more than one hundred images, we apply generalized Gaussian distribution to approximate APSF kernel in the image domain with parameter mapping from both shape similarity and numeric reasonability. Furthermore, a superpixel method is employed for estimating the transmission on the sky and non-sky region, in order to mitigate the halo artifact around the sharp edges and reduce color distortion in the sky region. Therefore, the haze-free images with abundant distinguished details and little halo artefacts can be finally restored by first applying deconvolution to the original hazy image. A series of experiments are additionally implemented to demonstrate the effectiveness and robustness of the proposed dehazing algorithm from both qualitative and quantitative comparison with the state-of-art.
Key words: Dehazing; Multiple scattering; Atmosphere point spread function; Dark channel prior; Superpixel Address all correspondence to: Rui Li, Beihang University, School of Instrumentation Science and Opto-electronics Engineering, Laboratory of Precision Opto-Mechatronics Technology, No. 37 Xueyuan Road, Haidian District, Beijing, China, 100191; Tel:+8601082338896; Email: [email protected]
1. Introduction Computer vision systems have enjoyed remarkable success in controlled and structured indoor environments, but applications is limited when it comes to be deployed outdoors [1]. Under bad weather conditions especially on a hazy or foggy day, amounts of sub-micron suspended particles (dust, smoke, airborne droplets etc.) interact with the light that passes through the air, for instance, the absorption and scattering. These interactions would severely cause degradation of images, like blurring effects, poor contrast, fade color and so on [1, 2]. Therefore, it is imperative to include some effective haze removal mechanisms in order to enable outdoor vision systems to function more reliably.
The common haze removal techniques can be broadly categorized in two types: methods based on perspective contrast enhancement [3, 4] and methods based on the degradation mechanism of hazy images. Comparatively, the first type method cannot completely overcome haze effects due to just using conventional image processing tools [1]. Hence, most of researchers prefer exploring methods rely on a physical image formation model that describes the hazy image as a convex combination between the scene radiance and the atmospheric light. The most reasonable model, which is widely adopted in dehazing algorithms, is the single scattering model (SSM) proposed by Narasimhan and Nayar. However, haze removal rely on this model is a challenging problem due to the unknown transmission coefficient that is dependent on the scene depth. Thus many methods are aimed at extracting the depth information from multiple images and extra information. For instance, polarization properties of different scattered light are used by Schechner et al. to restore the depth information through the polarized light in different directions [5-7]. Narasimhan et al. exploited two or more images captured in different weather conditions to obtain the 3D structure of the hazy scene for haze removal [8-10]. Kopf et al. introduced a kind of known 3D model [11] to obtain the depth. Although these methods can achieve reasonable transmission estimation and produce impressive results, they need to use at least two images or other extra information, which limits their applications. More recently, single image dehazing methods have been developed with significant progress. Researchers put forward different prior constraints to solve the model-based imaging formula to alleviate the requirements for scene depth. Fattal [12] provided an approach for transmission estimation by assuming the shading and transmission functions are locally independent. However, this method may be invalid when dealing with heavily hazed images. [2, 13]. Tan [14] argued that a haze-free image has more contrasts than the hazy one and gives an enhancement method for hazy images by maximizing their local contrast. This method could tackle very dense 2
haze conditions, yet, it also often produces oversaturated results. Tarel and Hautière [15] proposed a real time dehazing algorithm based on median filter and restored the visibility of the scene. Unfortunately it usually produces both oversaturated results and serious halo artefacts. He et al. [2] proposed Dark Channel Prior (DCP) technology for single image dehazing based on a statistic observation that most local patches in haze-free images often contain some low-intensity pixels. In most cases, DCP can recover high-quality haze-free images so that it receives widespread attention from researchers. However, it would fail when hazy images contain some scene objects similar to the atmospheric light especially in the sky region. As the consequence, color distortions appear in the restored results. Besides, based on the assumption that the depth is constant in local patch, DCP suffers halo artifacts in abrupt depth discontinuities unless using greatly time-consuming soft matting to refine the transmission. Although the guided filter [16] is suggested to improve the dark channel prior to substitute for the soft matting for time reduction, it has no effect on halo artefacts yet and even inevitably causes more serious color shift in sky region. Including the latest works of Nishino [17], Fattal [13] and Meng [18] that are based on statistics, the balance among robustness, recovery quality and complexity of these dehazing algorithms is still a tough issue. Moreover, the physical model these aforementioned methods based on is the single scattering model (SSM) and it ignores atmospheric multiple scattering, which would cause blurring effects on the image [19]. In this article, inspired by the robust patch-based dark channel prior, we propose a novel method for single-image dehazing that takes advantage of the SLIC superpixel [20], in which small image patches will be generated to share approximately the same depth for transmission estimation. The guided filter [16] is used to refine the transmission instead of soft matting, while the sky region is separated with region merging and the color distortion is rendered by properly adjusting the transmission on it. Moreover, an image degradation model based on multiple scattering which is described by atmosphere point spread function (APSF) is used. Our approach applies the generalized Gaussian distribution to approximate APSF kernel on the image domain by 3
parameter mapping from both shape similarity and numeric reasonability. In contrast to existing DCP [2] approaches that employs SSM and fixed size image patches, our algorithm validates its improvements and hence obtains high-quality haze-free images with abundant remarkable texture details and little halo artefacts as well as more less computational complexity. The remainder of this paper is organized as follows: Section 2 introduces the background of image formation in the presence of haze and gives the refined multiple scattering model. Section 3 presents the proposed image recovering procedure in detail. Comparative experiments as well as the analysis are undertaken in Section 4 to verify the feasibility and robustness of the proposed algorithm and conclusion are given in Section 5.
2. Image formation model The widespread dichromatic atmospheric scattering model, single scattering model (SSM), defines a simple and elegant expression for the radiance of an object at a distance d, i.e. Ls(d), in haze or fog condition [1, 2, 12]: d
Ls (d ) Lo e
dl 0
d
L e l dl Lo e d L (1 e d )
(1)
0
where Lo and L∞ are the radiance at the surface of the object and the airlight radiance at the horizon, respectively. β is the atmospheric attenuation coefficient that is assumed to be constant in homogenous atmosphere. Furthermore, the term Loe-βd and L∞(1-e-βd) on the right hand side describes respectively the attenuation mechanism (absorption and single scattering) of incident rays from the object and, airlight mechanism caused by the environmental illumination that are scattered into the field of view [2, 12, 14]. The SSM well explains the degradation causation of hazy images, but it is usually treated as a random phenomenon and assumes that the scattered light would not be dispersed by other suspended particles in the air. In fact, multiple scattering accounts for significant percentage in haze and thus can’t be neglected. Affected by multiple scattering, the extended light source of Ls that can be regarded as a combination of isotropic point 4
sources, spreads onto the neighboring zones and causes blurring effects or glows on the captured image, which could be simulated by using the range dependent Atmospheric Point Spread Function (APSF) as a convolution kernel [21, 22]:
Lm Ls APSF
(2)
where Lm is radiance captured by the camera, and ∗ stands for the convolution operator. Accordingly, a refined pattern in the multiple scattering model (MSM) for hazy image formation is proposed and illustrated in Fig. 1. Since the influence resulted from atmospheric multiple scattering is assumed to occur after the attenuation caused by atmospheric absorption and single scattering, Eq. (1) can be accepted as Ls in the MSM of Eq. (2). In another perspective, the airlight mechanism can be considered as the atmospheric interaction along the transmission path of length d, meanwhile radiance from the target object are attenuated through the atmospheric medium, both of which would be scattered multiple times by particles and reach the camera with different directions. Hence, the MSM can be derived by implementing the convolution on both sided of Eq. (1): d
Lm Lo e
dl 0
d
APSF(d ) L e l APSF(l )dl [Lo e d L (1 e d )] APSF(d )
(3)
0
In computer graphics, Eq. (3) is usually expressed in the form of intensity in RGB image:
I (x) [ J (x)t(x) A(1 t(x))] hAPSF
(4)
where t(x)=e-βd(x) indicates the medium transmission. I and J are the hazed and inherent haze-free image respectively. A is the global airlight color vector, x=(x,y)T denotes the pixel position on the image. hAPSF is the convolution matrix that can be obtained from APSF kernel hr(x) in image domain. According to Eq. (4), the goal of haze removal is to recover the true scene appearance J with the assumption that the transmission map t, the convolution matrix hAPSF and airlight A from the input hazy image are already estimated (the estimation of these parameters will be detailed later in subsections 3.1 and 3.2.2). Addressed to
5
the problem of imprecise estimation for t and A, we proposed a novel image recovery scheme based on DCP[2] with superpixel method and MSM. On another hand, with the research work on the radiation transfer theory, a good solution for APSF calculation by projecting which onto the 2D image domain hAPSF could be acquired [23].
Atmosphere
Image Plane
Field of View
α Lens
y
θ x d Attenuated direct light and scattered airlight
Particals Scene Point Assumed point source
Multiply scattered light
Fig. 1. Image formation based on multiple scattering model
3. Description of the proposed Image recovery
Airlight
Generate Superpixel
Restoring Equation
Sky separation
Transmission map Input hazy image
Haze-free image
Deconvolve with Wiener Filtering
Deconvolved image Convolution matrix APSF
Fig. 2. The overall procedure of the proposed image recovery scheme. The whole flowchart of the proposed dehazing approach is illustrated in Fig. 2. The scheme begins with a superpixel method for accurately representing depth information of the scene with image patch. Then the sky detection is implemented for estimating the global atmospheric light A and the medium transmission t(x). In the meantime, using the APSF kernel approximation that will be discussed in subsection 3.1, a set of convolution matrixes under different optical thickness T is determined. Based on the relationship between t and T, Wiener filtering is applied to deconvolve the input hazy image with suitable hAPSF for each pixel. Finally, the haze-free 6
image characterized by high visual quality and plentiful details can be restored with the deconvolved image.
3.1 Estimating hAPSF In order to quantify hAPSF, the APSF kernel hr(x), which is defined as the version of APSF in image domain, should be estimated. To the best of our knowledge, publications on APSF are confined to describe the relationship between APSF and the scattering angle θ (see Fig. 1)[1]. What’s more, it is impossible to establish the mapping about the image coordinate x in terms of θ just from single image in practice. An alternative solution, which we pursue here, is to approximate the analytic expression for hr(x) by observing the variation tendency of APSF versus θ. Although Metari and Deschênes have already proposed a good approximation for hr(x) with generalized Gaussian distribution (GGD) [24] through parameter mapping, they just consider qualitatively the peak shape similarity and ignore the quantitative numeric reasonability of hAPSF. In convolution operation, the values in hAPSF usually represent the contribution of the corresponding neighbor points, and unsuitable assignment for the weights of them might not result in satisfactory results. Actually, it is proved experimentally that a kind of negligible weight to the neighbors is produced by the mapping method of Metari and Deschênes so that the generated hAPSF does little improvement on the blur effects of multiple scattering. Consequently, aiming for better visual quality, we propose to refine the mapping between GGD parameters and APSF parameters by additionally taking the numeric reasonability into account. Taking the object point as an isotropic point light source, APSF can be calculated by solving the spherical radiation transfer equation (RTE) utilizing the Henyey-Greenstein phase function expanded with Legendre polynomials [21]:
7
I (T , ) ( g m (T ) g m 1 (T )) Pm ( ) m 0 g (T ) I e mT m log T 0 m m m 1 m 2m 1 (1 q m 1 ) m
(5)
where μ=cos θ. T=βd is the optical thickness, the value of which at each pixel can be calculated according to the transmission T (x) ln t (x) [22]. I0 is the radiant intensity of the isotropic source, and I(T, μ) is the captured radiant intensity for a given radial direction θ at the distance d. Pm and gm are the mth-order Legendre polynomial and the corresponding coefficient respectively. q is called the forward scattering parameter, whose value is assigned 0.2, 0.75, 0.9 corresponding to clear, hazy and foggy weather condition respectively. From Eq. (5), for a given value of T and q, the three-dimensional APSF as a function of θ can be generated by taking I0 as unit radiant intensity. Under various weather conditions (q varies) for mild (T=1.2) and strong dense atmosphere (T=4) [21, 24], normalized cross-sections of APSF, which are obtained from actual 3D APSFs, are displayed in Fig. 3 (a) and (b) respectively. From [21], we know that Eq. (5) is not convergent when T≥1. For hazy day ( q=0.75), APSF cross-sections in terms of θ under different T are illustrated in Fig.3(c). A conclusion can be drawn that the variance of APSF becomes less remarkable when T is larger than 2, especially when T≥4. For computational efficiency and without loss of generality, we focus on the estimation of hr(x) only under several representative values of T, i.e., T=1.2, 1.5, 2 and 4. 1
0.4 0.2
0.6 0.4 0.2
-100
0
=arccos
(a) T=1.2
100
200
0 -200
T=1.2 T=1.5 T=2 T=4 T=6
0.8
APSF
0.6
q=0.2 q=0.75 q=0.9
0.8
APSF
APSF
0.8
0 -200
1
1 q=0.2 q=0.75 q=0.9
0.6 0.4 0.2
-100
0
=arccos
(b) T=4 8
100
200
0 -200
-100
0
=arccos
100
200
(c) APSF for hazy day( q=0.75)
Fig. 3. Cross-sections of APSF normalized to [0-1] under different weather conditions and atmosphere types. On the other hand, the generalized Gaussian distribution (GGD) can be described as follows: |
e
GGD( x; , , p)
2(1
x p | A( p , )
1 ) A( p, ) p
(6)
, x
1 3 where A( p, ) ( ( ) / ( ))1/2 , and (∙) is the Gamma function. ξ is the mean that is assumed to be zero here. p p p and σ is the shape parameter and scale parameter respectively. p=0.6 p=0.8 p=1.5 p=2 p=2.5 p=3
GGD
0.8 0.6
1
0.6
0.4
0.4
0.2
0.2
0
-2
-1
0 x/pixel
1
p=0.6 p=0.8 p=1.5 p=2 p=2.5 p=3
0.8 GGD
1
0
2
(a) σ =0.5
-2
-1
0 x/pixel
1
2
(b) σ =0.7
Fig. 4. Sample lots of GGD normalized to [0-1] under different p and σ. Various GGD generated by varying p and σ are shown in Fig. 4. Comparing Fig. 3 with Fig. 4, we can see that T determines the peak shape of APSF in a similar manner to that of the parameter p dose on GGD, while the influence of q on the width of APSF is just opposite to that σ does on GGD. In consequence, a proper mapping relationship between p and T, σ and q could be respectively deduced from both peak shape similarity and numeric reasonability for approximating hr(x) with GGD. On one hand, it can be seen that the peak shape of GGD diverges a lot from APSF when the value of p is too large (p>2) or too small (p<1.5). So, for a good shape similarity, we suggest that the corresponding p at T=4 should be near 2, i.e., p (T )
T 4
2.
In addition, experimental results demonstrate that the quality of blurring effect removal during dehazing is affected by σ in some degree, which suggests that we should seek an admissible σ for q=0.75 through image 9
quality assessment. Since the commonly used criteria PSNR (Peak Signal to Noise Ratio), contrast and GMG (Gray Mean Grads) [25, 26] respectively evaluate the noise level, gray distribution and textural detail of an image, they are adopted to describe statistically the relationship between σ and the quality of the dehazed images for choosing the value of σ yielding outstanding results. Taking Fig. 5(a) as instance, the procedure for realizing the goal is carried out as follows: 1) A series of convolution matrixes are firstly generated from Eq.(6) by varying σ from 0.3 to 0.8 with a step of 0.05 and fixing p=2; 2) Then the three criteria (PSNR, contrast and GMG) correspond to each convolution matrix are calculated respectively on the dehazed image that is obtained with the proposed dehazing scheme shown in Fig. 2; 3) Finally the relationship between these criteria and σ are generated, as shown in Fig. 5(e). During experiment, the steps 1)-3) are repeated on about 110 hazy images that are randomly selected from Fattal’s database, the Internet and snapshots in real life by authors. Fig. 5 illustrates example images and the corresponding variation of these evaluators with respect to σ. Besides, the statistical result, derived through the averages of the three evaluators on the 110 images, is shown in Fig. 6.
(a)
(b)
(e)
(c)
(f)
(d)
(g)
(h)
Fig. 5. Example hazy images and the three criteria with respect to σ on each image. (a)-(d) hazy images, (e)-(h) the criteria PSNR, contrast and GMG in terms of σ for (a)-(h), respectively. 10
From Fig. 5 and Fig. 6, it is clear that the variation of all of the three criteria in terms of σ is monotone in spite of the difference in the change of their values. Specifically, PSNR declines while the others are in an uptrend as σ increases. Moreover, all of them vary slightly at the beginning (σ<0.5) and then have more and more sharp variations when σ becomes larger. In consideration of the fact that larger value of these evaluators means better perceptive quality for the recovered image after dehazing, we deduce that 0.5<σ<0.7 should be a better trade-off for a good evaluator from all of the three criteria. Then, we further determine the best value for σ from the averages of the three criteria in Fig. 6. Due to the variation of the average contrast is comparatively mild, the value of σ is derived mainly with the other two criteria. From Fig. 6, GMG has a dramatic increase near σ =0.6, while PSNR decays more and more severely when σ >0.6. Thus, the best balance for these criteria is suggested near σ =0.6. To make the expression concise, the mapping relationships between (p, σ), and (T, q) are ultimately established as follows:
p(T ) T , 0.5 / q .
(7)
In view of the rotation invariance of APSF, the analytical model of hr(x) approximated with GGD in two-dimension is given as follows in which the parameters p and σ can be obtained by Eq. (7): (
hr ( x, y , , p, )
x 2 y 2 0.5 p ) A( p , )
e , x ( x, y ) T 2(1 1 / p) A( p, )
2
(8)
Fig. 6. Average PSNR, contrast and GMG on the used 110 images. With Eq. (7) and Eq. (8), different hr(x) under T=1.2, 1.5, 2 and 4 for hazy day are obtained, from which the 11
corresponding convolution matrixes can be acquired to implement the deconvolution. In fact, owing to the similarity in ASPF between hazy day and foggy day (see Fig. 3), the mapping formula in Eq. (7) can also applied in foggy day (q=0.9). By substituting q=0.9 in Eq. (7)-(8), we can estimate the hr(x) for foggy day in the same way. Fig. 7 gives the 3D surface of the hr(x) in the case of mild (T=1.2) and strongly dense atmosphere (T=4) for both hazy (q=0.75) and foggy (q=0.9) weather conditions. Obviously, they match well with the cross-sections of the APSF in Fig. 3(a)(b) in terms of shape.
(a) T=1.2, q=0.75
(b) T=4, q=0.75
(c)T=1.2, q=0.9
(b) T=4, q=0.9
Fig. 7. Estimated 3D APSF kernel hr(x) normalized to [0-1] for both hazy and foggy weather conditions. Fig. 8 shows the representative convolution matrixes under T=4 and the deconvolved images generated by Metari’s and our methods. Intuitively, the proposed mapping approach yields more harmonious hAPSF and thus produces more distinguished detail information on the deconvolved image.
(a)
(b)
(c)
(d)
(e)
Fig. 8. Comparison with Metari’s method. (a)(b) 5×5 convolution matrixes obtained by Metari’s and our methods. (c) Hazy image, (d)(e) Deconvolved images obtained by Metari’s and our methods respectively.
3.2 Estimate t and A with the proposed method He et al. [2] summarize the dark channel prior (DCP) that the dark channel Jdark of the haze-free image J
12
tends to be zero: J dark min (min( J c (y))) 0 , based on which t(x) can be estimated through the dark channel c{r , g ,b} yx
of the normalized haze image Ic/Ac: t (x) 1 min (min( c{r , g ,b} yx
I c (y) Ac
))
(9)
where Ic, Ac is a color channel of I and A, respectively. Ωx represents a local patch of fixed size centered at x. According to [2], the parameter ω∈[0,1] is used for more natural dehazing results and usually set to 0.95. In addition, He et al. give the estimation for A: the top 0.1% brightest pixels in the dark channel of the original hazy image I are first picked, and then select the pixels with highest intensity in I as the airlight A [2]. With the assumption that the depth in the local patch Ωx is hardly variant, the transmission map obtained by Eq. (9) encounters severe block and halo artefacts. Although the pixel-wise soft matting algorithm [2] could well ameliorate these problems, it requires too much computing time due to the calculation of large sparse matrixes. In addition, DCP is invalid in sky region and easily causes serious color distortion, as shown in Fig. 12(e). Addressed to these issues, several solutions are already suggested [16, 27], however, robust performance has not been achieved. Inspired by the clue that adjacent pixels approximately sharing the same depth usually exhibit similar color and textural distribution in the image, we conduct a series of improvements for the conventional DCP with superpixel algorithm. Furthermore, the sky region detection is undertaken to render color distortion by identifying and then merging the possible superpixels which belong to the sky. In the detected sky region, A is calculated from the averages of three color channels, while the transmission t(x) is also adjusted following Eq. (9). The detailed procedure of A and t estimation is illustrated in Fig. 9.
13
Original hazy image I
Refine t in the skt region
Generate superpixels with SLIC
Estimate t with guided filter
Identify possible sky superpixel
Estimate A in the sky region
Sky superpixel Merging
Detect precisely the sky region
Fig. 9. Flowchart of the proposed transmission and airlight estimation method. 3.2.1 Superpixel generation and sky detection Since the superpixel is a small group of neighboring pixels with homogeneous color, it could represent the depth consistency in some degree and is more appropriate than the fixed size image patches. As for the choice of superpixel algorithm, the computational complexity should be concerned first. Then the segmentation performance, like boundary adherence, compactness of superpixels, is taken into account. According to the literatures, compared with some of the state-of-the-art superpixel algorithms such as TurboPixel, Normalized-cuts, Linear Spectral Clustering (LSC), et al [28], the Simple Linear Iterative Clustering (SLIC) [20] is less complex and requires less time. Even though the SEEDS[28] method runs a little faster than SLIC, it is difficult to control the superpixel number and to generate compact superpixels [28]. Thus, superpixel algorithm SLIC is adopted for image partition during dark channel calculation in our work. Contrast to the k-means clustering, SLIC dramatically speeds up execution with distance-limited clustering. Furthermore, a weighted distance measure combining color and spatial proximity is also applied in SLIC, which is described as follows:
D (ls lk )2 (as ak )2 (bs bk )2 r 2 (( xs xk )2 ( ys yk )2 ) / S 2
(10)
where (ls, as, bs)T and (xs, ys)T indicate the values of L, a, b in the CIELab color space and position vector of pixel s respectively. S is the interval between two adjacent superpixels, and r is a factor for adjusting the weight of position similarity. The effect of r is discussed in detail in [20] and here, we fix it to 20 for all results in this 14
paper. By iteratively clustering, a specific number of superpixels which well match object boundaries are ultimately constructed with low complexity, as it can be seen in Fig. 11(b)(e). Owing to the high intensity distribution in sky region, the possible sky superpixels, which are referred to as Sky Candidate Set (SCS) hereinafter, are loosely picked out with the condition that the mean values of R, G, B channel are larger 120. Then in SCS, we merge the most similar superpixels by 4 features: 3 color values in the CIELab color space and textural energy of the superpixel (the textural energy is calculated as: fi
( x, y )i
n 1,2,3
th
superpixel
| Y Fn | / 3 , where Fn , n=1,2,3, denotes three Laws’ masks shown in Fig.10, Y is the
intensity channel). The similarity of two adjacent superpixels in SCS, the ith and the jth superpixel, is measured on the basis of edge-focused strategy: Dij ij ( (li l j )2 (ai a j )2 (bi b j )2 | fi f j |) ij ( Dc Dt ) ,
(11)
where li, ai, bi are the mean color values of the ith superpixel, fi is the textural feature of the ith superpixel. τij is 1 if there is no edge between the two superpixels, otherwise ∞. The first and second term on the right hand side in Eq.(11) describe respectively the color and spacial texture measurement. ρ is an adaptive factor for adjusting the weight of textural distance. The statistical results in experiments show that, the Dt is usually 40 to100 times larger of Dc. Thus it is reasonable to set ρ between 40 to 300, then the two measurements could reach the same order of magnitude. Through the experiment and observation, it is found that despite the color measurmen is almost able to help us successfully pick out the sky superpixels from the input image, some bright objects whose intensity is similar to the sky might also be mixed in the sky superpixels by mistake. Meanwhile, the texture in the sky is not uniform especially near the boundary. According to the above statement, a small ρ is preferable as the auxiliary for more precisely distinguishing the sky superpixels. In our work,
the
nice sky identification result for most of all images is achieved with ρ fixed to 50. After merging the most similar superpixels by Eq.(11), the following criterion discribed in Eq.(12), in which ℤ and µ are empirical 15
integer values, is adopted to exactly identify whether a region k is the sky region or not:
min (mean( y c )) & & f k min fi .
c{r , g ,b}
(12)
yk
Fig. 10. Three Laws’ masks used to extract textural feature. Note that the haze concentration would affect the intensity distribution in the sky, the threshold ℤ and µ could be adjusted based on application as well as the scene structure information. Through repeatedly testing on more than 100 hazy images, we summarize that the values of ℤ and µ among 130~180 and 1~3 respectively could produce satisfying sky detection results. Fig.11 shows the results of superpixel generation and sky detection for images in thin and thick haze conditions respectively. As it can be seen, the superpixel algorithm is proved to yield good matches with object boundary information, and then a clearly demarcated sky region can be produced with the constraints after region merging.
(a)
(b)
(c)
(d)
(e)
(f)
Fig. 11. Results of superpixel generation and sky detection in thin and thick haze. (a)(d) Hazy image, (b)(e) Generated superpixels, (c)(f) Detected sky region which is painted black. Considering the specific secne structure and the density of haze, we set ℤ=160, µ=2.5 for (a), and ℤ=130, µ=2 for (d). 3.2.2 Estimate t and A After sky detection, a binary image like Fig. 11(c)(f) is generated by Eq. (12), guiding the estimation of the transmission t and airlight A under two situations. For hazy images with sky, the detected sky region is painted black while the non-sky is in white on the binary image. In this case, A is calculated by averaging the marked sky region. As it will be mentioned in the second paragraph in section 4, the small value of superpixel number N 16
would produce superpixels with non-uniform intensity distribution. Specifically, on ith suerpixel, the dark channel of its high intensity region is affected by its low intensity region, leading to a larger estimation for t than the true value. As a consequence, its high intensity region exhibits under-correction (i.e., incompletely removing the effect of haze) result. To avoid this, a correction factor λ is introduced in Eq. (9) as follows: t (x) 1 (1 ) min ( min (
I c (y)
c{r , g ,b} y,xsupi
Ac
))
(13)
where λ is experimentally set to 0.025 for non-sky region and -0.25 for sky region to reduce the color distortion caused by overcorrection due to the invalidation of dark channel prior. After that, the guided filter [16, 29] is
utilized to refine the transmission on the non-sky region instead of soft matting. If the image does not contain a sky region, a binary image is produced with no black pixel, meaning no sky is identified. Under such situation, A can be obtained with He’s method reported in [2], and the estimation of t is performed in the same way with the non-sky region. Therefore, the proposed dehazing algorithm can be applied to tackle almost any hazy image under various conditions.
3.3 Recover haze-free image After the transmission t is estimated, the optical thickness T can be calculated according to the relationship between T and t mentioned in subsection 3.1. Then we can deconvolve the input hazy image using the corresponding APSF kernel at each pixel, which is described as follows:
I (x) deconv( I (x), hr (x))
(14)
where I´ denotes the deconvolved hazy image, so the haze-free image is ultimately obtained with I´, t, and A: J
I A A. max(t , 0.1)
17
(15)
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
Fig. 12. Comparison with DCP and GDCP on Mountain image. (a) Hazy image, (b)-(d) Transmission map obtained by DCP, GDCP and our method respectively, (e)-(g) Dehazed results obtained with (b)-(d) respectively, (h) Enlarged results at the local patch marked with yellow in (e)-(g). Fig. 12 shows the transmission map as well as the dehazed images obtained by DCP, GDCP (DCP based on guided filter) and the proposed method. DCP produces excellent transmission with the most details, but it gives inaccurate transmission estimation in the sky region and outputs unnatural sky appearance. As the consequence of ambiguity object edges in the transmission obtained with GDCP, severe halo effects are encountered. Comparatively, our method avoids these unpleasant phenomena and demonstrates more detail information. As is reported, dense foggy images recovery is tougher due to the poorer visibility and serious loss of detail with gray appearance especially in the distant scene. Most current dehazing approaches always under- or over-correct and distort the color information in this region. Fortunately, by experiment, such gray region can be detected through similar color and texture strategy to the sky detection mentioned in subsection 3.2.1, which motivates us to control the degree of dehazing by adjusting the correction factor λ in Eq. (13) among -0.25~0.35. Based on this modification, we achieve better visual quality compared with GDCP and one of the state-of-the-art Gao’s [30] on both uniform and non-uniform illumination foggy images, as shown in Fig. 13. In addition, since 18
GMG can be used to distinguish the weather condition (haze or fog) from the perspective of textural details, it provides us the possibility to adaptively restore the degraded images in practical applications.
(a)
(d)
(b)
(c)
(e)
(f)
(g)
Fig. 13. Comparison with GDCP and Gao’s method on dense foggy images. (a) (d) Hazy images, (b)(e) GDCP, (f) Gao’s method and (c)(g) Our method.
4. Results and Analysis A series of hazy images that includes the benchmark images used by prestigious researchers are utilized to fully evaluate our methods. When restoring images in our way, Wiener Filtering is adopted to implement the deconvolution. The source codes written in MATLAB 2012a perform on a computer with Intel core E8400 3.0 GHz CPU and 3.0 GB RAM. In the proposed dehazing approach, the effect of superpixel number N should be considered. On one hand, the greater the N, the more similar the depth is on each superpixel, so the better the halo alleviation. On the other hand, the runtime for sky detection will increase with the increase of N. To balance the effect and efficiency of the haze removal, the value of N=200 for a 400×600 sized image, and N=400 for a large one about 768×1024 is adopted for most images in our experiments. 19
4.1 Qualitative Comparison
(a)
(b)
(c)
(d)
Fig. 14. Comparison with DCP and GDCP on Road image. (a) Hazy image, (b) DCP (c) GDCP and (d) Our method. The size of image patches for Eq.(9) in DCP and GDCP is 15×15. In Fig. 12-14, we compare our proposed approach with DCP and GDCP 1 under various circumstances. Even though DCP yields better dehazing performance than GDCP which suffers halo artefacts a lot, both of them fail to achieve an appropriate perceptual quality in the sky areas since they completely trust dark channel prior in all regions, as shown in Fig. 12(e)(f) and Fig. 14 (b)(c). In contrast, from Fig. 12(g)(h) and Fig.14(d), it is obvious that our results achieve dramatic visual quality with little halo effects and natural sky appearance caused by the proposed transmission estimation method based on superpixel algorithm and sky separation. Fig. 15 shows the comparison between our method and the three state-of-art haze removal algorithms: Fattal’s method (2008)2 [12], Tarel’s method (2009)3 [15] and He’s method (2009) (DCP)[2]. From Fig. 15(b), Fattal’s method (2008) could not reliably estimate the transmission sometimes, which can be attributed to the statistics-based hypothesis and the requirement for sufficient color information and variance. In dense haze or low variance region, therefore, this method under-corrects the distant scene in Landscape image and overcorrects the hut’s roof (marked with red rectangle) in Mountain image. Similar observations can also be found in Tarels’s results, which would bring out some halo artefacts at the same time, shown in the Mountain image in Fig. 15(c).
1
The matlab code of GDCP and DCP are obtained from https://github.com/akutta/Haze-Removal. The haze-free images of Fattal’s method and He’s method are downloaded from http://research.microsoft.com/en-us/um/people/kahe/cvpr09/comparisons.html. 3 The matlab code of Tarel’s method is obtained from http://perso.lcpc.fr/tarel.jean-philippe/publis/iccv09.html. 20 2
(a)
(b)
(c)
(d)
(e)
Fig. 15. Compared with Fattal’s, Tarel’s, and He’s method. From top to bottom: Landscape, Mountain and Cones image. (a)Hazy image, (b)Fattal’s method, (c)Tarel’s method, (d) DCP and (e)Our method. Compared with Fattal’s method and Tarel’s method, DCP and our method produce the most satisfactory results and demonstrate good robustness in most cases. Furthermore, with the multiple scattering based image formation model, our method exhibits more distinguished details and vivid contrasts in the restored results. Such superiority can be observed from the well refined edges of the objects in the regions marked with yellow rectangles in Fig. 12(g), Fig. 14(d) and Fig. 15(e). In Fig. 16 and Fig. 17, the proposed method is also compared with the latest dehazing works: methods of Nishino (2012) [31], Gibson (2013) [32] and Fattal (2014) [13]. In total, all the methods well remove the effect of haze with satisfactory visual results. However, Nishino’s and Gibson’s methods are prone to produce distorted sky color (see Fig. 16(b)-(c)), while the color information of the restored images obtained by Fattal’s method (2014) looks somewhat faint in Fig. 17(g). Once again, our results achieve the best performance in both contrast and the sharpness of the object edges, as shown in Fig. 16(f) and Fig. 17(d)(h).
21
(a)
(b)
(d)
(c)
(e)
(f)
Fig. 16. Compared with Nishino’s (2012), Gibson’s method (2013), Fattal’s (2014) and DCP. (a) Hazy image,(b) Nishino’s method, (c) Gibson’s method, (d) Fattal’s method, (e) DCP and (f) Our method.
(a)
(e)
(b)
(c)
(f)
(g)
(d)
(h)
Fig. 17. More comparisons with Nishino’s (2012), Gibson’s method (2013), Fattal’s (2014) on Pumpkins image and Cityscape image. (a)(e) Hazy images, (b) Nishino’s method, (c)(g) Fattal’s results, (f) Gibson’s method, (d)(h) Our method.
4.2 Quantitative Comparison In order to quantitatively assess the performance of the mentioned dehazing algorithms, two types of 22
assessment methods based on visible edges and image color are used: the ratio of new visible edges e, the ratio of the gradient on visible edges after and before restoration r , color naturalness index CNI and color colorfulness index CCI. e , r are proposed by Tarel and Hautière to evaluate the performance of a dehazing approach, and they measure the ability of the method to restore edges and contrast in gray-scale. CNI and CCI mainly assess the naturalness of a color image from the perspective of color saturation. Thus, the former two evaluators would give more reasonable assessment for dehazing methods. Usually, larger value of the four criteria means better image quality, and details of these four criteria can be found in publications [33, 34].
Table 1 Four criteria calculated on the images in Fig. 15. Landscape Method
Mountain
Cones
e
r
CNI
CCI
e
r
CNI
CCI
Ori
--
--
0.255
3.489
--
--
0.846
24.328
--
Fattal 08
2.640
2.277
0.510
8.004
0.101
1.717
0.904
33.937
Tarel 09
2.829
2.975
0.331
6.685
0.257
1.463
0.798
He 09
3.771
2.703
0.619
7.972
0.151
1.180
0.905
9.221
0.220
2.376
0.874
Our
4.672
4.651
0.686
r
CNI
CCI
--
0.579
7.990
0.239
2.324
0.798
13.099
28.763
0.279
2.317
0.944
17.009
29.521
0.294
1.577
0.906
11.689
2.477
0.924
15.749
40.096
e
0.334
Table 2 Four criteria calculated on the images in Fig. 16-17. Buildings Method
Pumpkins
Cityscape
e
r
CNI
CCI
e
r
CNI
CCI
e
r
CNI
CCI
Ori
--
--
0.292
3.520
--
--
0.685
21.010
--
--
0.493
5.497
Fattal 14
-0.0881
2.9334
0.570
13.166
0.1809
1.8978
0.910
31.974
1.0239
4.4375
0.659
11.327
Gibson 13
0.0165
1.6504
0.494
6.950
--
--
--
--
0.5859
2.3260
0.690
11.282
Nishinio 12
-0.0038
2.4656
0.745
11.250
0.3255
2.2853
0.631
42.129
--
--
--
--
Our
0.0895
2.0128
0.680
8.549
0.3659
3.1051
0.887
28.525
1.0727
4.9645
0.847
15.228
Table 1-2 calculates the four assessment criteria for the images in Fig.16-17. In most cases, our method generates the largest e and r , indicating that our method could restore the most visible edges and contrasts compared with others. On visual perception, our method almost yields the highest or the second highest CNI and CCI. These high values achieved by our method demonstrate the good capability of restoring high-quality haze-free images. 23
With the analysis above, our proposed approach performs well for most haze images. Nevertheless, it is inherently a combination of GDCP and sky detection, the iterative clustering in which would need extra runtime for the algorithm. Thus, our method could not be comparable to GDCP in the time consumption. More specifically, for an image sized 400×600 such as Fig. 16(a) and Fig. 17(a)(e), DCP, GDCP and our method require about 70 seconds, 2-3seconds and 6.5 seconds respectively. So, future works would focus on speeding up execution to make our method more efficient.
5. Conclusion We present a novel single-image dehazing method based on multiple scattering and dark channel prior (DCP). Unlike existing DCP approaches that employ SSM and fixed size image patches, the SLIC superpixel as a kind of classification tool for signal and pattern feature is first used to generate constant depth image patches. Moreover, the atmospheric multiple scattering described by atmosphere point spread function (APSF) is employed to refine the physics-model of hazy image formation. We derived a kind of reasonable parameter mapping relationships through image quality assessment to approximate the kernel function of APSF with the generalized Gaussian distribution. These techniques make the proposed method to achieve expressive results from both qualitative and quantitative assessment compared to state of the art algorithms as illustrated in the experiments. Since the more natural recovery effect can be developed with a kind of patch-wise self-adaption learning, other advanced models/features might be applied to strike a new balance between the efficiency and the effect. We leave this for future studies.
Acknowledgements The authors thank the anonymous reviewers for helping. This work was supported in part by a grant from National Natural Science Foundation of China (60974108).
24
References [1] S. G. Narasimhan, Models and algorithms for vision through the atmosphere, Doctor, Columbia University, 2003. [2] H. Kaiming and S. Jian and T. Xiaoou, Single image haze removal using dark channel prior, In: Proceesings of IEEE Conference on Computer Vision and Pattern Recognition, 2009, pp. 1956-1963. [3] H. Zhang and Q. Ye, Fog-degraded image clearness based on wavelet fusion, In: Proceesings of International Conference on Intelligent System Design and Engineering Application, 2011, pp. 759-761. [4] K. R. Joshi and R. S. Kamathe, Quantification of retinex in enhancement of weather degraded images, In: Proceesings of International Conference on Audio, Language and Image Processing, 2008, pp. 1229-1233. [5] S. Shwartz and E. Namer and Y. Y. Schechner, Blind haze separation, In: Proceesings of IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 2006, pp. 1984-1991. [6] Y. Y. Schechner and S. G. Narasimhan and S. K. Nayar, Instant dehazing of images using polarization, In: Proceesings of IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 2001, pp. I-325-I-332. [7] S. YY and N. SG and N. SK, Polarization-based vision through haze, Applied Optics, 42 (3) (2003) 511-525. [8] S. G. Narasimhan and S. K. Nayar, Chromatic framework for vision in bad weather, In: Proceesings of IEEE Conference on Computer Vision and Pattern Recognition, 2000, pp. 598-605. [9] S. G. Narasimhan and S. K. Nayar, Contrast restoration of weather degraded images, Pattern Analysis and Machine Intelligence, IEEE Transactions on, 25 (6) (2003) 713-724. [10] S. Narasimhan, Nayar, S., Interactive (de) weathering of an image using physical models, IEEE Workshop on Color and Photometric Methods in Computer Vision, 6 (4) (2003) 1-8. [11] J. Kopf, B. Neubert, B. Chen et al., Deep photo: model-based photograph enhancement and viewing, ACM 25
Transactions on Graphics, 27 (5) (2008) 1-10. [12] R. Fattal, Single image dehazing, In: Proceedings of ACM SIGGRAPH, 2008, pp. 1-9. [13] R. Fattal, Dehazing using color-lines, ACM Transactions on Graphics (TOG), 34 (1) (2014) 13:1-13:14. [14] R. T. Tan, Visibility in bad weather from a single image, In: Proceesings of IEEE Conference on Computer Vision and Pattern Recognition, , 2008, pp. 1-8. [15] J. P. Tarel and N. Hautiere, Fast visibility restoration from a single color or gray level image, In: Proceesings of 12th IEEE International Conference on Computer Vision, 2009, 2009, pp. 2201-2208. [16] X. Lan, L. Zhang, H. Shen et al., Single image haze removal considering sensor blur and noise, EURASIP Journal on Advances in Signal Processing, 2013 (1) (2013) 1-13. [17] L. Kratz and K. Nishino, Factorizing Scene Albedo and Depth from a Single Foggy Image, In: Proceesings of 12th IEEE International Conference on Computer Vision, 2009, pp. 1701-1708. [18] G. Meng, Y. Wang, J. Duan et al., Efficient image dehazing with boundary constraint and contextual regularization, In: Proceesings of IEEE International Conference on Computer Vision, 2013, pp. 617-624. [19] P. Wei, Y. Liu, Y. Liu et al., Dehazing model based on multiple scattering, In: Proceesings of 3rd International Congress on Image and Signal Processing (CISP), 2010, pp. 249-252. [20] R. Achanta, A. Shaji, K. Smith et al., SLIC superpixels compared to state-of-the-art superpixel methods, IEEE Transactions on Pattern Analysis and Machine Intelligence, 34 (11) (2012) 2274-2282. [21] S. G. Narasimhan and S. K. Nayar, Shedding light on the weather, In: Proceesings of IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 2003, pp. I/665-I/672. [22] S. Tao, H. Feng, Z. Xu et al., Image degradation and recovery based on multiple scattering in remote sensing and bad weather condition, Optics Express, 20 (15) (2012) 16584-16595. [23] R. N. Mahalati and J. M. Kahn, Effect of fog on free-space optical links employing imaging receivers, Optics 26
Express, 20 (2) (2012) 1649-1661. [24] S. Metari and F. Deschenes, A new convolution kernel for atmospheric point spread function applied to computer vision, In: Proceesings of 11th IEEE International Conference on Computer Vision, 2007, pp. 1-8 [25] J. Xiao-qiang, Research on fast image defogging and visibility restoration, Doctor, Changchun Institute of Optics, Fine Mechanics and Physics, Chinese Academy of Science, Changchun, Chinese Academy of Science, Changchun, 2012. [26] R. Wang and Y.-z. Cui and Y. Yuan, Image quality assessment using full-parameter singular value decomposition, Optical Engineering, 50 (5) (2011) 057005-1-057005-9. [27] T. M. Bui, H. N. Tran, W. Kim et al., Segmenting dark channel prior in single image dehazing, Electronics Letters, 50 (7) (2014) 516-518. [28] L. Zhengqin and C. Jiansheng, Superpixel segmentation using Linear Spectral Clustering, In: Proceesings of IEEE Conference on Computer Vision and Pattern Recognition, 2015, pp. 1356-1363. [29] S. Hao, D. Pan, Y. Guo et al., Image detail enhancement with spatially guided filters, Signal Processi ng, 120 (3) (2015) 789-796. [30] Y. Gao, H.-M. Hu, S. Wang et al., A fast image dehazing algorithm based on negative correction, Signal Processing, 103 (2014) 380-398. [31] K. Nishino and L. Kratz and S. Lombardi, Bayesian Defogging, International Journal of Computer Vision, 98 (3) (2012) 263-278. [32] K. B. Gibson and T. Q. Nguyen, An analysis of single image defogging methods using a color ellipsoid framework, EURASIP Journal on Image and Video Processing, 2013 (1) (2013) 1-14. [33] N. HAUTIèRE, J.-P. TAREL, D. AUBERT et al., Blind contrast restoration assessment by gradient rationing at visible edges, Image Analysis and Stereology, 27 (2) (2008) 87-95. 27
[34] K.-Q. Huang and Q. Wang and Z.-Y. Wu, Natural color image enhancement and evaluation algorithm based on human visual system, Computer Vision and Image Understanding, 103 (1) (2006) 52-63.
Highlights
A multiple scattering model for hazy image formation is proposed.
The atmosphere point spread function is estimated with shape-numeric reasonability.
Sky detection with superpixel merging to reduce the color distortion.
The superpixel method is employed to estimate the transmission for halo prevention.
28
www.elsevier.com/locate/sigpro
PII: DOI: Reference:
S0165-1684(16)00041-4 http://dx.doi.org/10.1016/j.sigpro.2016.02.003 SIGPRO6050
To appear in: Signal Processing Received date: 14 September 2015 Revised date: 10 January 2016 Accepted date: 1 February 2016 Cite this article as: Rui Wang, Rui Li and Haiyan Sun, Haze removal based on multiple scattering model with superpixel algorithm, Signal Processing, http://dx.doi.org/10.1016/j.sigpro.2016.02.003 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Haze removal based on multiple scattering model with superpixel algorithm Rui Wanga, Rui Lia, Haiyan Sunb a. Key Laboratory of Precision Opto-mechatronics Technology, Ministry of Education, School of Instrumentation Science and Opto-electronics Engineering, Beihang University, Beijing, 100191, China b. School of Mathematics and System Science, Beihang University, Beijing, 100191 China Abstract: An efficient method to remove haze from single image based on dark channel prior and the multiple scattering described by atmosphere point spread function (APSF) is proposed in this paper. By means of statistical analysis of the quality assessments for more than one hundred images, we apply generalized Gaussian distribution to approximate APSF kernel in the image domain with parameter mapping from both shape similarity and numeric reasonability. Furthermore, a superpixel method is employed for estimating the transmission on the sky and non-sky region, in order to mitigate the halo artifact around the sharp edges and reduce color distortion in the sky region. Therefore, the haze-free images with abundant distinguished details and little halo artefacts can be finally restored by first applying deconvolution to the original hazy image. A series of experiments are additionally implemented to demonstrate the effectiveness and robustness of the proposed dehazing algorithm from both qualitative and quantitative comparison with the state-of-art.
Key words: Dehazing; Multiple scattering; Atmosphere point spread function; Dark channel prior; Superpixel Address all correspondence to: Rui Li, Beihang University, School of Instrumentation Science and Opto-electronics Engineering, Laboratory of Precision Opto-Mechatronics Technology, No. 37 Xueyuan Road, Haidian District, Beijing, China, 100191; Tel:+8601082338896; Email: [email protected]
1. Introduction Computer vision systems have enjoyed remarkable success in controlled and structured indoor environments, but applications is limited when it comes to be deployed outdoors [1]. Under bad weather conditions especially on a hazy or foggy day, amounts of sub-micron suspended particles (dust, smoke, airborne droplets etc.) interact with the light that passes through the air, for instance, the absorption and scattering. These interactions would severely cause degradation of images, like blurring effects, poor contrast, fade color and so on [1, 2]. Therefore, it is imperative to include some effective haze removal mechanisms in order to enable outdoor vision systems to function more reliably.
The common haze removal techniques can be broadly categorized in two types: methods based on perspective contrast enhancement [3, 4] and methods based on the degradation mechanism of hazy images. Comparatively, the first type method cannot completely overcome haze effects due to just using conventional image processing tools [1]. Hence, most of researchers prefer exploring methods rely on a physical image formation model that describes the hazy image as a convex combination between the scene radiance and the atmospheric light. The most reasonable model, which is widely adopted in dehazing algorithms, is the single scattering model (SSM) proposed by Narasimhan and Nayar. However, haze removal rely on this model is a challenging problem due to the unknown transmission coefficient that is dependent on the scene depth. Thus many methods are aimed at extracting the depth information from multiple images and extra information. For instance, polarization properties of different scattered light are used by Schechner et al. to restore the depth information through the polarized light in different directions [5-7]. Narasimhan et al. exploited two or more images captured in different weather conditions to obtain the 3D structure of the hazy scene for haze removal [8-10]. Kopf et al. introduced a kind of known 3D model [11] to obtain the depth. Although these methods can achieve reasonable transmission estimation and produce impressive results, they need to use at least two images or other extra information, which limits their applications. More recently, single image dehazing methods have been developed with significant progress. Researchers put forward different prior constraints to solve the model-based imaging formula to alleviate the requirements for scene depth. Fattal [12] provided an approach for transmission estimation by assuming the shading and transmission functions are locally independent. However, this method may be invalid when dealing with heavily hazed images. [2, 13]. Tan [14] argued that a haze-free image has more contrasts than the hazy one and gives an enhancement method for hazy images by maximizing their local contrast. This method could tackle very dense 2
haze conditions, yet, it also often produces oversaturated results. Tarel and Hautière [15] proposed a real time dehazing algorithm based on median filter and restored the visibility of the scene. Unfortunately it usually produces both oversaturated results and serious halo artefacts. He et al. [2] proposed Dark Channel Prior (DCP) technology for single image dehazing based on a statistic observation that most local patches in haze-free images often contain some low-intensity pixels. In most cases, DCP can recover high-quality haze-free images so that it receives widespread attention from researchers. However, it would fail when hazy images contain some scene objects similar to the atmospheric light especially in the sky region. As the consequence, color distortions appear in the restored results. Besides, based on the assumption that the depth is constant in local patch, DCP suffers halo artifacts in abrupt depth discontinuities unless using greatly time-consuming soft matting to refine the transmission. Although the guided filter [16] is suggested to improve the dark channel prior to substitute for the soft matting for time reduction, it has no effect on halo artefacts yet and even inevitably causes more serious color shift in sky region. Including the latest works of Nishino [17], Fattal [13] and Meng [18] that are based on statistics, the balance among robustness, recovery quality and complexity of these dehazing algorithms is still a tough issue. Moreover, the physical model these aforementioned methods based on is the single scattering model (SSM) and it ignores atmospheric multiple scattering, which would cause blurring effects on the image [19]. In this article, inspired by the robust patch-based dark channel prior, we propose a novel method for single-image dehazing that takes advantage of the SLIC superpixel [20], in which small image patches will be generated to share approximately the same depth for transmission estimation. The guided filter [16] is used to refine the transmission instead of soft matting, while the sky region is separated with region merging and the color distortion is rendered by properly adjusting the transmission on it. Moreover, an image degradation model based on multiple scattering which is described by atmosphere point spread function (APSF) is used. Our approach applies the generalized Gaussian distribution to approximate APSF kernel on the image domain by 3
parameter mapping from both shape similarity and numeric reasonability. In contrast to existing DCP [2] approaches that employs SSM and fixed size image patches, our algorithm validates its improvements and hence obtains high-quality haze-free images with abundant remarkable texture details and little halo artefacts as well as more less computational complexity. The remainder of this paper is organized as follows: Section 2 introduces the background of image formation in the presence of haze and gives the refined multiple scattering model. Section 3 presents the proposed image recovering procedure in detail. Comparative experiments as well as the analysis are undertaken in Section 4 to verify the feasibility and robustness of the proposed algorithm and conclusion are given in Section 5.
2. Image formation model The widespread dichromatic atmospheric scattering model, single scattering model (SSM), defines a simple and elegant expression for the radiance of an object at a distance d, i.e. Ls(d), in haze or fog condition [1, 2, 12]: d
Ls (d ) Lo e
dl 0
d
L e l dl Lo e d L (1 e d )
(1)
0
where Lo and L∞ are the radiance at the surface of the object and the airlight radiance at the horizon, respectively. β is the atmospheric attenuation coefficient that is assumed to be constant in homogenous atmosphere. Furthermore, the term Loe-βd and L∞(1-e-βd) on the right hand side describes respectively the attenuation mechanism (absorption and single scattering) of incident rays from the object and, airlight mechanism caused by the environmental illumination that are scattered into the field of view [2, 12, 14]. The SSM well explains the degradation causation of hazy images, but it is usually treated as a random phenomenon and assumes that the scattered light would not be dispersed by other suspended particles in the air. In fact, multiple scattering accounts for significant percentage in haze and thus can’t be neglected. Affected by multiple scattering, the extended light source of Ls that can be regarded as a combination of isotropic point 4
sources, spreads onto the neighboring zones and causes blurring effects or glows on the captured image, which could be simulated by using the range dependent Atmospheric Point Spread Function (APSF) as a convolution kernel [21, 22]:
Lm Ls APSF
(2)
where Lm is radiance captured by the camera, and ∗ stands for the convolution operator. Accordingly, a refined pattern in the multiple scattering model (MSM) for hazy image formation is proposed and illustrated in Fig. 1. Since the influence resulted from atmospheric multiple scattering is assumed to occur after the attenuation caused by atmospheric absorption and single scattering, Eq. (1) can be accepted as Ls in the MSM of Eq. (2). In another perspective, the airlight mechanism can be considered as the atmospheric interaction along the transmission path of length d, meanwhile radiance from the target object are attenuated through the atmospheric medium, both of which would be scattered multiple times by particles and reach the camera with different directions. Hence, the MSM can be derived by implementing the convolution on both sided of Eq. (1): d
Lm Lo e
dl 0
d
APSF(d ) L e l APSF(l )dl [Lo e d L (1 e d )] APSF(d )
(3)
0
In computer graphics, Eq. (3) is usually expressed in the form of intensity in RGB image:
I (x) [ J (x)t(x) A(1 t(x))] hAPSF
(4)
where t(x)=e-βd(x) indicates the medium transmission. I and J are the hazed and inherent haze-free image respectively. A is the global airlight color vector, x=(x,y)T denotes the pixel position on the image. hAPSF is the convolution matrix that can be obtained from APSF kernel hr(x) in image domain. According to Eq. (4), the goal of haze removal is to recover the true scene appearance J with the assumption that the transmission map t, the convolution matrix hAPSF and airlight A from the input hazy image are already estimated (the estimation of these parameters will be detailed later in subsections 3.1 and 3.2.2). Addressed to
5
the problem of imprecise estimation for t and A, we proposed a novel image recovery scheme based on DCP[2] with superpixel method and MSM. On another hand, with the research work on the radiation transfer theory, a good solution for APSF calculation by projecting which onto the 2D image domain hAPSF could be acquired [23].
Atmosphere
Image Plane
Field of View
α Lens
y
θ x d Attenuated direct light and scattered airlight
Particals Scene Point Assumed point source
Multiply scattered light
Fig. 1. Image formation based on multiple scattering model
3. Description of the proposed Image recovery
Airlight
Generate Superpixel
Restoring Equation
Sky separation
Transmission map Input hazy image
Haze-free image
Deconvolve with Wiener Filtering
Deconvolved image Convolution matrix APSF
Fig. 2. The overall procedure of the proposed image recovery scheme. The whole flowchart of the proposed dehazing approach is illustrated in Fig. 2. The scheme begins with a superpixel method for accurately representing depth information of the scene with image patch. Then the sky detection is implemented for estimating the global atmospheric light A and the medium transmission t(x). In the meantime, using the APSF kernel approximation that will be discussed in subsection 3.1, a set of convolution matrixes under different optical thickness T is determined. Based on the relationship between t and T, Wiener filtering is applied to deconvolve the input hazy image with suitable hAPSF for each pixel. Finally, the haze-free 6
image characterized by high visual quality and plentiful details can be restored with the deconvolved image.
3.1 Estimating hAPSF In order to quantify hAPSF, the APSF kernel hr(x), which is defined as the version of APSF in image domain, should be estimated. To the best of our knowledge, publications on APSF are confined to describe the relationship between APSF and the scattering angle θ (see Fig. 1)[1]. What’s more, it is impossible to establish the mapping about the image coordinate x in terms of θ just from single image in practice. An alternative solution, which we pursue here, is to approximate the analytic expression for hr(x) by observing the variation tendency of APSF versus θ. Although Metari and Deschênes have already proposed a good approximation for hr(x) with generalized Gaussian distribution (GGD) [24] through parameter mapping, they just consider qualitatively the peak shape similarity and ignore the quantitative numeric reasonability of hAPSF. In convolution operation, the values in hAPSF usually represent the contribution of the corresponding neighbor points, and unsuitable assignment for the weights of them might not result in satisfactory results. Actually, it is proved experimentally that a kind of negligible weight to the neighbors is produced by the mapping method of Metari and Deschênes so that the generated hAPSF does little improvement on the blur effects of multiple scattering. Consequently, aiming for better visual quality, we propose to refine the mapping between GGD parameters and APSF parameters by additionally taking the numeric reasonability into account. Taking the object point as an isotropic point light source, APSF can be calculated by solving the spherical radiation transfer equation (RTE) utilizing the Henyey-Greenstein phase function expanded with Legendre polynomials [21]:
7
I (T , ) ( g m (T ) g m 1 (T )) Pm ( ) m 0 g (T ) I e mT m log T 0 m m m 1 m 2m 1 (1 q m 1 ) m
(5)
where μ=cos θ. T=βd is the optical thickness, the value of which at each pixel can be calculated according to the transmission T (x) ln t (x) [22]. I0 is the radiant intensity of the isotropic source, and I(T, μ) is the captured radiant intensity for a given radial direction θ at the distance d. Pm and gm are the mth-order Legendre polynomial and the corresponding coefficient respectively. q is called the forward scattering parameter, whose value is assigned 0.2, 0.75, 0.9 corresponding to clear, hazy and foggy weather condition respectively. From Eq. (5), for a given value of T and q, the three-dimensional APSF as a function of θ can be generated by taking I0 as unit radiant intensity. Under various weather conditions (q varies) for mild (T=1.2) and strong dense atmosphere (T=4) [21, 24], normalized cross-sections of APSF, which are obtained from actual 3D APSFs, are displayed in Fig. 3 (a) and (b) respectively. From [21], we know that Eq. (5) is not convergent when T≥1. For hazy day ( q=0.75), APSF cross-sections in terms of θ under different T are illustrated in Fig.3(c). A conclusion can be drawn that the variance of APSF becomes less remarkable when T is larger than 2, especially when T≥4. For computational efficiency and without loss of generality, we focus on the estimation of hr(x) only under several representative values of T, i.e., T=1.2, 1.5, 2 and 4. 1
0.4 0.2
0.6 0.4 0.2
-100
0
=arccos
(a) T=1.2
100
200
0 -200
T=1.2 T=1.5 T=2 T=4 T=6
0.8
APSF
0.6
q=0.2 q=0.75 q=0.9
0.8
APSF
APSF
0.8
0 -200
1
1 q=0.2 q=0.75 q=0.9
0.6 0.4 0.2
-100
0
=arccos
(b) T=4 8
100
200
0 -200
-100
0
=arccos
100
200
(c) APSF for hazy day( q=0.75)
Fig. 3. Cross-sections of APSF normalized to [0-1] under different weather conditions and atmosphere types. On the other hand, the generalized Gaussian distribution (GGD) can be described as follows: |
e
GGD( x; , , p)
2(1
x p | A( p , )
1 ) A( p, ) p
(6)
, x
1 3 where A( p, ) ( ( ) / ( ))1/2 , and (∙) is the Gamma function. ξ is the mean that is assumed to be zero here. p p p and σ is the shape parameter and scale parameter respectively. p=0.6 p=0.8 p=1.5 p=2 p=2.5 p=3
GGD
0.8 0.6
1
0.6
0.4
0.4
0.2
0.2
0
-2
-1
0 x/pixel
1
p=0.6 p=0.8 p=1.5 p=2 p=2.5 p=3
0.8 GGD
1
0
2
(a) σ =0.5
-2
-1
0 x/pixel
1
2
(b) σ =0.7
Fig. 4. Sample lots of GGD normalized to [0-1] under different p and σ. Various GGD generated by varying p and σ are shown in Fig. 4. Comparing Fig. 3 with Fig. 4, we can see that T determines the peak shape of APSF in a similar manner to that of the parameter p dose on GGD, while the influence of q on the width of APSF is just opposite to that σ does on GGD. In consequence, a proper mapping relationship between p and T, σ and q could be respectively deduced from both peak shape similarity and numeric reasonability for approximating hr(x) with GGD. On one hand, it can be seen that the peak shape of GGD diverges a lot from APSF when the value of p is too large (p>2) or too small (p<1.5). So, for a good shape similarity, we suggest that the corresponding p at T=4 should be near 2, i.e., p (T )
T 4
2.
In addition, experimental results demonstrate that the quality of blurring effect removal during dehazing is affected by σ in some degree, which suggests that we should seek an admissible σ for q=0.75 through image 9
quality assessment. Since the commonly used criteria PSNR (Peak Signal to Noise Ratio), contrast and GMG (Gray Mean Grads) [25, 26] respectively evaluate the noise level, gray distribution and textural detail of an image, they are adopted to describe statistically the relationship between σ and the quality of the dehazed images for choosing the value of σ yielding outstanding results. Taking Fig. 5(a) as instance, the procedure for realizing the goal is carried out as follows: 1) A series of convolution matrixes are firstly generated from Eq.(6) by varying σ from 0.3 to 0.8 with a step of 0.05 and fixing p=2; 2) Then the three criteria (PSNR, contrast and GMG) correspond to each convolution matrix are calculated respectively on the dehazed image that is obtained with the proposed dehazing scheme shown in Fig. 2; 3) Finally the relationship between these criteria and σ are generated, as shown in Fig. 5(e). During experiment, the steps 1)-3) are repeated on about 110 hazy images that are randomly selected from Fattal’s database, the Internet and snapshots in real life by authors. Fig. 5 illustrates example images and the corresponding variation of these evaluators with respect to σ. Besides, the statistical result, derived through the averages of the three evaluators on the 110 images, is shown in Fig. 6.
(a)
(b)
(e)
(c)
(f)
(d)
(g)
(h)
Fig. 5. Example hazy images and the three criteria with respect to σ on each image. (a)-(d) hazy images, (e)-(h) the criteria PSNR, contrast and GMG in terms of σ for (a)-(h), respectively. 10
From Fig. 5 and Fig. 6, it is clear that the variation of all of the three criteria in terms of σ is monotone in spite of the difference in the change of their values. Specifically, PSNR declines while the others are in an uptrend as σ increases. Moreover, all of them vary slightly at the beginning (σ<0.5) and then have more and more sharp variations when σ becomes larger. In consideration of the fact that larger value of these evaluators means better perceptive quality for the recovered image after dehazing, we deduce that 0.5<σ<0.7 should be a better trade-off for a good evaluator from all of the three criteria. Then, we further determine the best value for σ from the averages of the three criteria in Fig. 6. Due to the variation of the average contrast is comparatively mild, the value of σ is derived mainly with the other two criteria. From Fig. 6, GMG has a dramatic increase near σ =0.6, while PSNR decays more and more severely when σ >0.6. Thus, the best balance for these criteria is suggested near σ =0.6. To make the expression concise, the mapping relationships between (p, σ), and (T, q) are ultimately established as follows:
p(T ) T , 0.5 / q .
(7)
In view of the rotation invariance of APSF, the analytical model of hr(x) approximated with GGD in two-dimension is given as follows in which the parameters p and σ can be obtained by Eq. (7): (
hr ( x, y , , p, )
x 2 y 2 0.5 p ) A( p , )
e , x ( x, y ) T 2(1 1 / p) A( p, )
2
(8)
Fig. 6. Average PSNR, contrast and GMG on the used 110 images. With Eq. (7) and Eq. (8), different hr(x) under T=1.2, 1.5, 2 and 4 for hazy day are obtained, from which the 11
corresponding convolution matrixes can be acquired to implement the deconvolution. In fact, owing to the similarity in ASPF between hazy day and foggy day (see Fig. 3), the mapping formula in Eq. (7) can also applied in foggy day (q=0.9). By substituting q=0.9 in Eq. (7)-(8), we can estimate the hr(x) for foggy day in the same way. Fig. 7 gives the 3D surface of the hr(x) in the case of mild (T=1.2) and strongly dense atmosphere (T=4) for both hazy (q=0.75) and foggy (q=0.9) weather conditions. Obviously, they match well with the cross-sections of the APSF in Fig. 3(a)(b) in terms of shape.
(a) T=1.2, q=0.75
(b) T=4, q=0.75
(c)T=1.2, q=0.9
(b) T=4, q=0.9
Fig. 7. Estimated 3D APSF kernel hr(x) normalized to [0-1] for both hazy and foggy weather conditions. Fig. 8 shows the representative convolution matrixes under T=4 and the deconvolved images generated by Metari’s and our methods. Intuitively, the proposed mapping approach yields more harmonious hAPSF and thus produces more distinguished detail information on the deconvolved image.
(a)
(b)
(c)
(d)
(e)
Fig. 8. Comparison with Metari’s method. (a)(b) 5×5 convolution matrixes obtained by Metari’s and our methods. (c) Hazy image, (d)(e) Deconvolved images obtained by Metari’s and our methods respectively.
3.2 Estimate t and A with the proposed method He et al. [2] summarize the dark channel prior (DCP) that the dark channel Jdark of the haze-free image J
12
tends to be zero: J dark min (min( J c (y))) 0 , based on which t(x) can be estimated through the dark channel c{r , g ,b} yx
of the normalized haze image Ic/Ac: t (x) 1 min (min( c{r , g ,b} yx
I c (y) Ac
))
(9)
where Ic, Ac is a color channel of I and A, respectively. Ωx represents a local patch of fixed size centered at x. According to [2], the parameter ω∈[0,1] is used for more natural dehazing results and usually set to 0.95. In addition, He et al. give the estimation for A: the top 0.1% brightest pixels in the dark channel of the original hazy image I are first picked, and then select the pixels with highest intensity in I as the airlight A [2]. With the assumption that the depth in the local patch Ωx is hardly variant, the transmission map obtained by Eq. (9) encounters severe block and halo artefacts. Although the pixel-wise soft matting algorithm [2] could well ameliorate these problems, it requires too much computing time due to the calculation of large sparse matrixes. In addition, DCP is invalid in sky region and easily causes serious color distortion, as shown in Fig. 12(e). Addressed to these issues, several solutions are already suggested [16, 27], however, robust performance has not been achieved. Inspired by the clue that adjacent pixels approximately sharing the same depth usually exhibit similar color and textural distribution in the image, we conduct a series of improvements for the conventional DCP with superpixel algorithm. Furthermore, the sky region detection is undertaken to render color distortion by identifying and then merging the possible superpixels which belong to the sky. In the detected sky region, A is calculated from the averages of three color channels, while the transmission t(x) is also adjusted following Eq. (9). The detailed procedure of A and t estimation is illustrated in Fig. 9.
13
Original hazy image I
Refine t in the skt region
Generate superpixels with SLIC
Estimate t with guided filter
Identify possible sky superpixel
Estimate A in the sky region
Sky superpixel Merging
Detect precisely the sky region
Fig. 9. Flowchart of the proposed transmission and airlight estimation method. 3.2.1 Superpixel generation and sky detection Since the superpixel is a small group of neighboring pixels with homogeneous color, it could represent the depth consistency in some degree and is more appropriate than the fixed size image patches. As for the choice of superpixel algorithm, the computational complexity should be concerned first. Then the segmentation performance, like boundary adherence, compactness of superpixels, is taken into account. According to the literatures, compared with some of the state-of-the-art superpixel algorithms such as TurboPixel, Normalized-cuts, Linear Spectral Clustering (LSC), et al [28], the Simple Linear Iterative Clustering (SLIC) [20] is less complex and requires less time. Even though the SEEDS[28] method runs a little faster than SLIC, it is difficult to control the superpixel number and to generate compact superpixels [28]. Thus, superpixel algorithm SLIC is adopted for image partition during dark channel calculation in our work. Contrast to the k-means clustering, SLIC dramatically speeds up execution with distance-limited clustering. Furthermore, a weighted distance measure combining color and spatial proximity is also applied in SLIC, which is described as follows:
D (ls lk )2 (as ak )2 (bs bk )2 r 2 (( xs xk )2 ( ys yk )2 ) / S 2
(10)
where (ls, as, bs)T and (xs, ys)T indicate the values of L, a, b in the CIELab color space and position vector of pixel s respectively. S is the interval between two adjacent superpixels, and r is a factor for adjusting the weight of position similarity. The effect of r is discussed in detail in [20] and here, we fix it to 20 for all results in this 14
paper. By iteratively clustering, a specific number of superpixels which well match object boundaries are ultimately constructed with low complexity, as it can be seen in Fig. 11(b)(e). Owing to the high intensity distribution in sky region, the possible sky superpixels, which are referred to as Sky Candidate Set (SCS) hereinafter, are loosely picked out with the condition that the mean values of R, G, B channel are larger 120. Then in SCS, we merge the most similar superpixels by 4 features: 3 color values in the CIELab color space and textural energy of the superpixel (the textural energy is calculated as: fi
( x, y )i
n 1,2,3
th
superpixel
| Y Fn | / 3 , where Fn , n=1,2,3, denotes three Laws’ masks shown in Fig.10, Y is the
intensity channel). The similarity of two adjacent superpixels in SCS, the ith and the jth superpixel, is measured on the basis of edge-focused strategy: Dij ij ( (li l j )2 (ai a j )2 (bi b j )2 | fi f j |) ij ( Dc Dt ) ,
(11)
where li, ai, bi are the mean color values of the ith superpixel, fi is the textural feature of the ith superpixel. τij is 1 if there is no edge between the two superpixels, otherwise ∞. The first and second term on the right hand side in Eq.(11) describe respectively the color and spacial texture measurement. ρ is an adaptive factor for adjusting the weight of textural distance. The statistical results in experiments show that, the Dt is usually 40 to100 times larger of Dc. Thus it is reasonable to set ρ between 40 to 300, then the two measurements could reach the same order of magnitude. Through the experiment and observation, it is found that despite the color measurmen is almost able to help us successfully pick out the sky superpixels from the input image, some bright objects whose intensity is similar to the sky might also be mixed in the sky superpixels by mistake. Meanwhile, the texture in the sky is not uniform especially near the boundary. According to the above statement, a small ρ is preferable as the auxiliary for more precisely distinguishing the sky superpixels. In our work,
the
nice sky identification result for most of all images is achieved with ρ fixed to 50. After merging the most similar superpixels by Eq.(11), the following criterion discribed in Eq.(12), in which ℤ and µ are empirical 15
integer values, is adopted to exactly identify whether a region k is the sky region or not:
min (mean( y c )) & & f k min fi .
c{r , g ,b}
(12)
yk
Fig. 10. Three Laws’ masks used to extract textural feature. Note that the haze concentration would affect the intensity distribution in the sky, the threshold ℤ and µ could be adjusted based on application as well as the scene structure information. Through repeatedly testing on more than 100 hazy images, we summarize that the values of ℤ and µ among 130~180 and 1~3 respectively could produce satisfying sky detection results. Fig.11 shows the results of superpixel generation and sky detection for images in thin and thick haze conditions respectively. As it can be seen, the superpixel algorithm is proved to yield good matches with object boundary information, and then a clearly demarcated sky region can be produced with the constraints after region merging.
(a)
(b)
(c)
(d)
(e)
(f)
Fig. 11. Results of superpixel generation and sky detection in thin and thick haze. (a)(d) Hazy image, (b)(e) Generated superpixels, (c)(f) Detected sky region which is painted black. Considering the specific secne structure and the density of haze, we set ℤ=160, µ=2.5 for (a), and ℤ=130, µ=2 for (d). 3.2.2 Estimate t and A After sky detection, a binary image like Fig. 11(c)(f) is generated by Eq. (12), guiding the estimation of the transmission t and airlight A under two situations. For hazy images with sky, the detected sky region is painted black while the non-sky is in white on the binary image. In this case, A is calculated by averaging the marked sky region. As it will be mentioned in the second paragraph in section 4, the small value of superpixel number N 16
would produce superpixels with non-uniform intensity distribution. Specifically, on ith suerpixel, the dark channel of its high intensity region is affected by its low intensity region, leading to a larger estimation for t than the true value. As a consequence, its high intensity region exhibits under-correction (i.e., incompletely removing the effect of haze) result. To avoid this, a correction factor λ is introduced in Eq. (9) as follows: t (x) 1 (1 ) min ( min (
I c (y)
c{r , g ,b} y,xsupi
Ac
))
(13)
where λ is experimentally set to 0.025 for non-sky region and -0.25 for sky region to reduce the color distortion caused by overcorrection due to the invalidation of dark channel prior. After that, the guided filter [16, 29] is
utilized to refine the transmission on the non-sky region instead of soft matting. If the image does not contain a sky region, a binary image is produced with no black pixel, meaning no sky is identified. Under such situation, A can be obtained with He’s method reported in [2], and the estimation of t is performed in the same way with the non-sky region. Therefore, the proposed dehazing algorithm can be applied to tackle almost any hazy image under various conditions.
3.3 Recover haze-free image After the transmission t is estimated, the optical thickness T can be calculated according to the relationship between T and t mentioned in subsection 3.1. Then we can deconvolve the input hazy image using the corresponding APSF kernel at each pixel, which is described as follows:
I (x) deconv( I (x), hr (x))
(14)
where I´ denotes the deconvolved hazy image, so the haze-free image is ultimately obtained with I´, t, and A: J
I A A. max(t , 0.1)
17
(15)
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
Fig. 12. Comparison with DCP and GDCP on Mountain image. (a) Hazy image, (b)-(d) Transmission map obtained by DCP, GDCP and our method respectively, (e)-(g) Dehazed results obtained with (b)-(d) respectively, (h) Enlarged results at the local patch marked with yellow in (e)-(g). Fig. 12 shows the transmission map as well as the dehazed images obtained by DCP, GDCP (DCP based on guided filter) and the proposed method. DCP produces excellent transmission with the most details, but it gives inaccurate transmission estimation in the sky region and outputs unnatural sky appearance. As the consequence of ambiguity object edges in the transmission obtained with GDCP, severe halo effects are encountered. Comparatively, our method avoids these unpleasant phenomena and demonstrates more detail information. As is reported, dense foggy images recovery is tougher due to the poorer visibility and serious loss of detail with gray appearance especially in the distant scene. Most current dehazing approaches always under- or over-correct and distort the color information in this region. Fortunately, by experiment, such gray region can be detected through similar color and texture strategy to the sky detection mentioned in subsection 3.2.1, which motivates us to control the degree of dehazing by adjusting the correction factor λ in Eq. (13) among -0.25~0.35. Based on this modification, we achieve better visual quality compared with GDCP and one of the state-of-the-art Gao’s [30] on both uniform and non-uniform illumination foggy images, as shown in Fig. 13. In addition, since 18
GMG can be used to distinguish the weather condition (haze or fog) from the perspective of textural details, it provides us the possibility to adaptively restore the degraded images in practical applications.
(a)
(d)
(b)
(c)
(e)
(f)
(g)
Fig. 13. Comparison with GDCP and Gao’s method on dense foggy images. (a) (d) Hazy images, (b)(e) GDCP, (f) Gao’s method and (c)(g) Our method.
4. Results and Analysis A series of hazy images that includes the benchmark images used by prestigious researchers are utilized to fully evaluate our methods. When restoring images in our way, Wiener Filtering is adopted to implement the deconvolution. The source codes written in MATLAB 2012a perform on a computer with Intel core E8400 3.0 GHz CPU and 3.0 GB RAM. In the proposed dehazing approach, the effect of superpixel number N should be considered. On one hand, the greater the N, the more similar the depth is on each superpixel, so the better the halo alleviation. On the other hand, the runtime for sky detection will increase with the increase of N. To balance the effect and efficiency of the haze removal, the value of N=200 for a 400×600 sized image, and N=400 for a large one about 768×1024 is adopted for most images in our experiments. 19
4.1 Qualitative Comparison
(a)
(b)
(c)
(d)
Fig. 14. Comparison with DCP and GDCP on Road image. (a) Hazy image, (b) DCP (c) GDCP and (d) Our method. The size of image patches for Eq.(9) in DCP and GDCP is 15×15. In Fig. 12-14, we compare our proposed approach with DCP and GDCP 1 under various circumstances. Even though DCP yields better dehazing performance than GDCP which suffers halo artefacts a lot, both of them fail to achieve an appropriate perceptual quality in the sky areas since they completely trust dark channel prior in all regions, as shown in Fig. 12(e)(f) and Fig. 14 (b)(c). In contrast, from Fig. 12(g)(h) and Fig.14(d), it is obvious that our results achieve dramatic visual quality with little halo effects and natural sky appearance caused by the proposed transmission estimation method based on superpixel algorithm and sky separation. Fig. 15 shows the comparison between our method and the three state-of-art haze removal algorithms: Fattal’s method (2008)2 [12], Tarel’s method (2009)3 [15] and He’s method (2009) (DCP)[2]. From Fig. 15(b), Fattal’s method (2008) could not reliably estimate the transmission sometimes, which can be attributed to the statistics-based hypothesis and the requirement for sufficient color information and variance. In dense haze or low variance region, therefore, this method under-corrects the distant scene in Landscape image and overcorrects the hut’s roof (marked with red rectangle) in Mountain image. Similar observations can also be found in Tarels’s results, which would bring out some halo artefacts at the same time, shown in the Mountain image in Fig. 15(c).
1
The matlab code of GDCP and DCP are obtained from https://github.com/akutta/Haze-Removal. The haze-free images of Fattal’s method and He’s method are downloaded from http://research.microsoft.com/en-us/um/people/kahe/cvpr09/comparisons.html. 3 The matlab code of Tarel’s method is obtained from http://perso.lcpc.fr/tarel.jean-philippe/publis/iccv09.html. 20 2
(a)
(b)
(c)
(d)
(e)
Fig. 15. Compared with Fattal’s, Tarel’s, and He’s method. From top to bottom: Landscape, Mountain and Cones image. (a)Hazy image, (b)Fattal’s method, (c)Tarel’s method, (d) DCP and (e)Our method. Compared with Fattal’s method and Tarel’s method, DCP and our method produce the most satisfactory results and demonstrate good robustness in most cases. Furthermore, with the multiple scattering based image formation model, our method exhibits more distinguished details and vivid contrasts in the restored results. Such superiority can be observed from the well refined edges of the objects in the regions marked with yellow rectangles in Fig. 12(g), Fig. 14(d) and Fig. 15(e). In Fig. 16 and Fig. 17, the proposed method is also compared with the latest dehazing works: methods of Nishino (2012) [31], Gibson (2013) [32] and Fattal (2014) [13]. In total, all the methods well remove the effect of haze with satisfactory visual results. However, Nishino’s and Gibson’s methods are prone to produce distorted sky color (see Fig. 16(b)-(c)), while the color information of the restored images obtained by Fattal’s method (2014) looks somewhat faint in Fig. 17(g). Once again, our results achieve the best performance in both contrast and the sharpness of the object edges, as shown in Fig. 16(f) and Fig. 17(d)(h).
21
(a)
(b)
(d)
(c)
(e)
(f)
Fig. 16. Compared with Nishino’s (2012), Gibson’s method (2013), Fattal’s (2014) and DCP. (a) Hazy image,(b) Nishino’s method, (c) Gibson’s method, (d) Fattal’s method, (e) DCP and (f) Our method.
(a)
(e)
(b)
(c)
(f)
(g)
(d)
(h)
Fig. 17. More comparisons with Nishino’s (2012), Gibson’s method (2013), Fattal’s (2014) on Pumpkins image and Cityscape image. (a)(e) Hazy images, (b) Nishino’s method, (c)(g) Fattal’s results, (f) Gibson’s method, (d)(h) Our method.
4.2 Quantitative Comparison In order to quantitatively assess the performance of the mentioned dehazing algorithms, two types of 22
assessment methods based on visible edges and image color are used: the ratio of new visible edges e, the ratio of the gradient on visible edges after and before restoration r , color naturalness index CNI and color colorfulness index CCI. e , r are proposed by Tarel and Hautière to evaluate the performance of a dehazing approach, and they measure the ability of the method to restore edges and contrast in gray-scale. CNI and CCI mainly assess the naturalness of a color image from the perspective of color saturation. Thus, the former two evaluators would give more reasonable assessment for dehazing methods. Usually, larger value of the four criteria means better image quality, and details of these four criteria can be found in publications [33, 34].
Table 1 Four criteria calculated on the images in Fig. 15. Landscape Method
Mountain
Cones
e
r
CNI
CCI
e
r
CNI
CCI
Ori
--
--
0.255
3.489
--
--
0.846
24.328
--
Fattal 08
2.640
2.277
0.510
8.004
0.101
1.717
0.904
33.937
Tarel 09
2.829
2.975
0.331
6.685
0.257
1.463
0.798
He 09
3.771
2.703
0.619
7.972
0.151
1.180
0.905
9.221
0.220
2.376
0.874
Our
4.672
4.651
0.686
r
CNI
CCI
--
0.579
7.990
0.239
2.324
0.798
13.099
28.763
0.279
2.317
0.944
17.009
29.521
0.294
1.577
0.906
11.689
2.477
0.924
15.749
40.096
e
0.334
Table 2 Four criteria calculated on the images in Fig. 16-17. Buildings Method
Pumpkins
Cityscape
e
r
CNI
CCI
e
r
CNI
CCI
e
r
CNI
CCI
Ori
--
--
0.292
3.520
--
--
0.685
21.010
--
--
0.493
5.497
Fattal 14
-0.0881
2.9334
0.570
13.166
0.1809
1.8978
0.910
31.974
1.0239
4.4375
0.659
11.327
Gibson 13
0.0165
1.6504
0.494
6.950
--
--
--
--
0.5859
2.3260
0.690
11.282
Nishinio 12
-0.0038
2.4656
0.745
11.250
0.3255
2.2853
0.631
42.129
--
--
--
--
Our
0.0895
2.0128
0.680
8.549
0.3659
3.1051
0.887
28.525
1.0727
4.9645
0.847
15.228
Table 1-2 calculates the four assessment criteria for the images in Fig.16-17. In most cases, our method generates the largest e and r , indicating that our method could restore the most visible edges and contrasts compared with others. On visual perception, our method almost yields the highest or the second highest CNI and CCI. These high values achieved by our method demonstrate the good capability of restoring high-quality haze-free images. 23
With the analysis above, our proposed approach performs well for most haze images. Nevertheless, it is inherently a combination of GDCP and sky detection, the iterative clustering in which would need extra runtime for the algorithm. Thus, our method could not be comparable to GDCP in the time consumption. More specifically, for an image sized 400×600 such as Fig. 16(a) and Fig. 17(a)(e), DCP, GDCP and our method require about 70 seconds, 2-3seconds and 6.5 seconds respectively. So, future works would focus on speeding up execution to make our method more efficient.
5. Conclusion We present a novel single-image dehazing method based on multiple scattering and dark channel prior (DCP). Unlike existing DCP approaches that employ SSM and fixed size image patches, the SLIC superpixel as a kind of classification tool for signal and pattern feature is first used to generate constant depth image patches. Moreover, the atmospheric multiple scattering described by atmosphere point spread function (APSF) is employed to refine the physics-model of hazy image formation. We derived a kind of reasonable parameter mapping relationships through image quality assessment to approximate the kernel function of APSF with the generalized Gaussian distribution. These techniques make the proposed method to achieve expressive results from both qualitative and quantitative assessment compared to state of the art algorithms as illustrated in the experiments. Since the more natural recovery effect can be developed with a kind of patch-wise self-adaption learning, other advanced models/features might be applied to strike a new balance between the efficiency and the effect. We leave this for future studies.
Acknowledgements The authors thank the anonymous reviewers for helping. This work was supported in part by a grant from National Natural Science Foundation of China (60974108).
24
References [1] S. G. Narasimhan, Models and algorithms for vision through the atmosphere, Doctor, Columbia University, 2003. [2] H. Kaiming and S. Jian and T. Xiaoou, Single image haze removal using dark channel prior, In: Proceesings of IEEE Conference on Computer Vision and Pattern Recognition, 2009, pp. 1956-1963. [3] H. Zhang and Q. Ye, Fog-degraded image clearness based on wavelet fusion, In: Proceesings of International Conference on Intelligent System Design and Engineering Application, 2011, pp. 759-761. [4] K. R. Joshi and R. S. Kamathe, Quantification of retinex in enhancement of weather degraded images, In: Proceesings of International Conference on Audio, Language and Image Processing, 2008, pp. 1229-1233. [5] S. Shwartz and E. Namer and Y. Y. Schechner, Blind haze separation, In: Proceesings of IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 2006, pp. 1984-1991. [6] Y. Y. Schechner and S. G. Narasimhan and S. K. Nayar, Instant dehazing of images using polarization, In: Proceesings of IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 2001, pp. I-325-I-332. [7] S. YY and N. SG and N. SK, Polarization-based vision through haze, Applied Optics, 42 (3) (2003) 511-525. [8] S. G. Narasimhan and S. K. Nayar, Chromatic framework for vision in bad weather, In: Proceesings of IEEE Conference on Computer Vision and Pattern Recognition, 2000, pp. 598-605. [9] S. G. Narasimhan and S. K. Nayar, Contrast restoration of weather degraded images, Pattern Analysis and Machine Intelligence, IEEE Transactions on, 25 (6) (2003) 713-724. [10] S. Narasimhan, Nayar, S., Interactive (de) weathering of an image using physical models, IEEE Workshop on Color and Photometric Methods in Computer Vision, 6 (4) (2003) 1-8. [11] J. Kopf, B. Neubert, B. Chen et al., Deep photo: model-based photograph enhancement and viewing, ACM 25
Transactions on Graphics, 27 (5) (2008) 1-10. [12] R. Fattal, Single image dehazing, In: Proceedings of ACM SIGGRAPH, 2008, pp. 1-9. [13] R. Fattal, Dehazing using color-lines, ACM Transactions on Graphics (TOG), 34 (1) (2014) 13:1-13:14. [14] R. T. Tan, Visibility in bad weather from a single image, In: Proceesings of IEEE Conference on Computer Vision and Pattern Recognition, , 2008, pp. 1-8. [15] J. P. Tarel and N. Hautiere, Fast visibility restoration from a single color or gray level image, In: Proceesings of 12th IEEE International Conference on Computer Vision, 2009, 2009, pp. 2201-2208. [16] X. Lan, L. Zhang, H. Shen et al., Single image haze removal considering sensor blur and noise, EURASIP Journal on Advances in Signal Processing, 2013 (1) (2013) 1-13. [17] L. Kratz and K. Nishino, Factorizing Scene Albedo and Depth from a Single Foggy Image, In: Proceesings of 12th IEEE International Conference on Computer Vision, 2009, pp. 1701-1708. [18] G. Meng, Y. Wang, J. Duan et al., Efficient image dehazing with boundary constraint and contextual regularization, In: Proceesings of IEEE International Conference on Computer Vision, 2013, pp. 617-624. [19] P. Wei, Y. Liu, Y. Liu et al., Dehazing model based on multiple scattering, In: Proceesings of 3rd International Congress on Image and Signal Processing (CISP), 2010, pp. 249-252. [20] R. Achanta, A. Shaji, K. Smith et al., SLIC superpixels compared to state-of-the-art superpixel methods, IEEE Transactions on Pattern Analysis and Machine Intelligence, 34 (11) (2012) 2274-2282. [21] S. G. Narasimhan and S. K. Nayar, Shedding light on the weather, In: Proceesings of IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 2003, pp. I/665-I/672. [22] S. Tao, H. Feng, Z. Xu et al., Image degradation and recovery based on multiple scattering in remote sensing and bad weather condition, Optics Express, 20 (15) (2012) 16584-16595. [23] R. N. Mahalati and J. M. Kahn, Effect of fog on free-space optical links employing imaging receivers, Optics 26
Express, 20 (2) (2012) 1649-1661. [24] S. Metari and F. Deschenes, A new convolution kernel for atmospheric point spread function applied to computer vision, In: Proceesings of 11th IEEE International Conference on Computer Vision, 2007, pp. 1-8 [25] J. Xiao-qiang, Research on fast image defogging and visibility restoration, Doctor, Changchun Institute of Optics, Fine Mechanics and Physics, Chinese Academy of Science, Changchun, Chinese Academy of Science, Changchun, 2012. [26] R. Wang and Y.-z. Cui and Y. Yuan, Image quality assessment using full-parameter singular value decomposition, Optical Engineering, 50 (5) (2011) 057005-1-057005-9. [27] T. M. Bui, H. N. Tran, W. Kim et al., Segmenting dark channel prior in single image dehazing, Electronics Letters, 50 (7) (2014) 516-518. [28] L. Zhengqin and C. Jiansheng, Superpixel segmentation using Linear Spectral Clustering, In: Proceesings of IEEE Conference on Computer Vision and Pattern Recognition, 2015, pp. 1356-1363. [29] S. Hao, D. Pan, Y. Guo et al., Image detail enhancement with spatially guided filters, Signal Processi ng, 120 (3) (2015) 789-796. [30] Y. Gao, H.-M. Hu, S. Wang et al., A fast image dehazing algorithm based on negative correction, Signal Processing, 103 (2014) 380-398. [31] K. Nishino and L. Kratz and S. Lombardi, Bayesian Defogging, International Journal of Computer Vision, 98 (3) (2012) 263-278. [32] K. B. Gibson and T. Q. Nguyen, An analysis of single image defogging methods using a color ellipsoid framework, EURASIP Journal on Image and Video Processing, 2013 (1) (2013) 1-14. [33] N. HAUTIèRE, J.-P. TAREL, D. AUBERT et al., Blind contrast restoration assessment by gradient rationing at visible edges, Image Analysis and Stereology, 27 (2) (2008) 87-95. 27
[34] K.-Q. Huang and Q. Wang and Z.-Y. Wu, Natural color image enhancement and evaluation algorithm based on human visual system, Computer Vision and Image Understanding, 103 (1) (2006) 52-63.
Highlights
A multiple scattering model for hazy image formation is proposed.
The atmosphere point spread function is estimated with shape-numeric reasonability.
Sky detection with superpixel merging to reduce the color distortion.
The superpixel method is employed to estimate the transmission for halo prevention.
28