Scattering of light caused by water particles will reduce quality of images captured in foggy weather. Hitherto, there are kinds of methods to attempt to remove haze in the images, including methods based on image processing [1,2], methods based on physical models [3–5], and methods based on deep learning [6–8]. Among these methods, the polarimetric dehazing method [9–11], as a physical-model-based method, has attracted the attention of researchers benefiting from development of polarimetric imaging techniques recently. In contrast to image processing methods and deep learning methods, the polarimetric method is based on a simple, reasonable, and general physical model, i.e., optical degradation model [12,13], meaning that the polarimetric method has a strong interpretability and generalization ability. On the other hand, essentially, the polarimetric method sufficiently utilizes multimodal information of light (intensity and polarization) to reconstruct a clear image, so the method inherently enables highly informative reconstruction [14]. Therefore, the polarimetric dehazing method can effectively improve the quality of hazy images and generate good results in general [15,16].

Although the polarimetric dehazing method is effective, there still are some problems that make the method unstable and unrobust. This study mainly focuses on two key problems. At first, the polarimetric dehazing method is vulnerable to noise. Note that in the polarimetric method, polarization parameters measured by polarization imaging, such as degree of linear polarization (DoLP) and angle of polarization (AoP), are used to estimate physical characteristics of haze for further dehazing. The quality of the result of the polarimetric dehazing algorithm heavily depends on the accuracy of these parameters because some mathematical operation executed in the algorithm may amplify noises. However, these parameters are sensitive to noise [17]. Besides, currently popular polarization imaging solutions, i.e., division-of-aperture polarization cameras and division-of-focal-plane polarization cameras, have unavoidable defects that lead to more measurement errors (registration errors, instantaneous field of view errors, etc.) compared to traditional intensity imaging, which means there are more uncertainties in the polarization parameters. So, how to improve the error resistance of the polarimetric dehazing method is the first problem. Second, results of the polarimetric dehazing method are images with a high dynamic range (HDR) in general, meaning that the contrast of such images may be extremely high. In some cases, a high contrast means more details in the image and is conducive to improving the quality of the image. However, most ordinary monitors are only supporting the standard dynamic range (SDR). If the dynamic range of an image is so high that it far exceeds the dynamic range of the monitors, the amount of information that can be transferred by the monitors will decrease on the contrary. Therefore, how to keep more details of dehazed images within SDR is the second problem.

To solve the two problems of the polarimetric dehazing method, this paper proposes a novel dehazing algorithm. For the first problem, the proposed algorithm uses a low-rank approximation method to estimate the DoLP, instead of the conventional Stokes parameter method. It should be emphasized that although the low-rank approximation method as a mathematic method has been used for general image denoising [18,19], we find that the method is especially appropriate to denoise DoLP images due to the physical and mathematical nature of such a polarization parameter. Therefore, the method is not a pure mathematical algorithm, but a method utilizing physical properties of polarized images to reduce noises. Compared to our previously proposed noise suppressing method for polarization parameters [9,17], the proposed method has stronger physical interpretability. Besides, the low-rank approximation method can be realized by a deterministic algorithm without iteration, meaning less cost on calculation.

For the second problem, i.e., the exaggeratedly high dynamic range of the dehazed image usually obtained through polarimetric dehazing, we propose a multiple virtual-exposure fusion (MVEF) scheme to covert the HDR image to an SDR image, meanwhile, keeping more details as much as possible.

In the following content of this paper, the principle and the complete flow path of the proposed method are introduced in Section 2 with necessary details. Some experimental results of the proposed method are shown in Section 3, where evaluation of such results and comparison to other dehazing methods are provided for reference. In Section 4, we discuss the superiority of the proposed method compared to conventional methods and our previous method.

The polarimetric dehazing method is based on a physical model, i.e., the optical degradation model. The optical degradation model is shown in Fig. 1 schematically, where a scenario of imaging with scattering media is modeled. The light emitted or reflected by the object is called object light in the figure. Before the incidence of object light into the camera, the light will be scattered by the scattering media. The scattering will lead to the attenuation of the object light and generate the scattering light. The attenuated object light is the degraded object light.

According to the model described in Fig. 1, the total intensity at (x,y) of the image plane of the camera can be expressed as I(x,y)=D(x,y)+A(x,y),(1)where D(x,y) and A(x,y) are intensities of the degraded object light and the scattering light. D(x,y) can be further expressed as a product of the intensity of object light and an attenuation factor, as shown in Eq. (2): D(x,y)=t(z)L(x,y),(2)where L(x,y) is the object light on the imaging plane (assuming that there is no scattering), and t(z) is the attenuation factor that is related to z, i.e., distance between the imaging plane and the object plane. In general, t(z) is modeled as a negative exponential function [17] for homogeneous scattering media. In the model, the scattering media are assumed homogeneous, so t(∞)=0 and t(0)=1. Retrieval of L(x,y) is the goal of the polarimetric dehazing method.

If the distance z is infinity, then the object light will be fully converted to scattering light, because t(∞)=0. Denoting the scattering light from the infinity as A∞(x,y), there is a relation between A(x,y) and A∞(x,y) [9]: A(x,y)=[1−t(z)]A∞(x,y).(3)

In Eq. (3), we assumed that the two-dimensional distributions of A(x,y) and A∞(x,y) are similar. It should be emphasized that, if the scattering media are inhomogeneous, then t(∞) may not equal zero so that Eq. (3) will not hold. So, the model is limited to deal with homogeneous scattering media.

According to Eqs. (1)–(3), the following equation is obtained: L(x,y)=I(x,y)−A(x,y)1−A(x,y)/A∞(x,y).(4)

Equation (4) shows that if I(x,y), A(x,y), and A∞(x,y) are known, then the object light L(x,y) will be reconstructed. The total incidence on the camera, i.e., I(x,y), is easy to obtain. So, the problem is how to estimate A(x,y) and A∞(x,y).

According to the Mie’s scattering theory, the scattering light is partially polarized [20]. A beam of partially polarized light can always be regarded as a superposition of beams of fully linearly polarized light, fully circularly polarized light, and unpolarized light. Assuming that the fully linearly polarized part in A(x,y) is Ap(x,y), then the DoLP distribution of the scattering light, noted as pA(x,y), can be expressed as pA(x,y)=Ap(x,y)A(x,y).(5)

Further, if the object light is fully unpolarized, then Ap(x,y) will also be the fully linearly polarized part in I(x,y), meaning that the DoLP of the total incidence I(x,y) satisfies p(x,y)=Ap(x,y)I(x,y),(6)where p(x,y) is the distribution of DoLP of the total incidence, which can be directly measured via polarization imaging, meaning that Ap(x,y) can be obtained from Eq. (6) with measured I(x,y). Further, according to Eq. (5), if A(x,y) is known, then pA(x,y) can be calculated from the equation. However, A(x,y) is unknown, so there is not an efficient method to precisely estimate pA(x,y). Therefore, to simplify the model, we assume that pA(x,y) is constant on the imaging plane. It is obvious that pA(x,y)≥p(x,y), because Eq. (1) leads to I(x,y)≥A(x,y). So, a reasonable way to estimate pA will be pA=max[p(x,y)].(7)

Once pA is estimated, then one can obtain A(x,y) from Eq. (5) with the Ap(x,y). The last problem is to estimate A∞(x,y). In general, A∞(x,y) can also be assumed as a constant, but one needs the sky region in the hazy image to estimate such a constant accurately [15], meaning that the generality of the algorithm is limited. Therefore, to avoid the problem, the A∞ is regarded as a value that is larger than the largest value of A(x,y) in this paper: A∞=1.02×max[A(x,y)].(8)

According to Eq. (3), A∞(x,y) is always larger than A(x,y) because 1−t(z)∈(0,1) when z∈(0,∞). Therefore, we assume that A∞ is larger than the largest value of A(x,y). The factor of 1.02 is used to prevent the denominator of Eq. (4) from being zero: if A∞=max[A(x,y)], then the denominator will be zero at (x,y)=argmax[A(x,y)]. On the other hand, the value of 1.02 is an empirical value. According to our experiments conducted before, if the factor is relatively large (for example, 1.1), the contrast of a dehazed image will decrease; if the factor is relatively small (for example, 1.001), the denominator of Eq. (4) will be close to zero at a certain (x,y), which will lead to numerical instability. Therefore, we choose the factor to be 1.02 to balance the image quality and the numerical stability.

Let us give a brief summary. At first, one needs to obtain p(x,y) and I(x,y) via polarization imaging. Second, the fully linearly polarized part of the scattering light, i.e., Ap(x,y), can be obtained from Eq. (6), and the DoLP of the scattering light, i.e., pA, can be estimated from Eq. (7). Third, A(x,y) and A∞ can be estimated from Eq. (5) and Eq. (8), respectively. Finally, one can reconstruct the object light, i.e., L(x,y), from Eq. (4).

In the end, as a supplement, I(x,y) and p(x,y) can be obtained by the Stokes method from polarization images captured under different angles of polarization of 0°, 45°, 90°, and 135°: [S0S1S2]=[I0°+I90°I0°−I90°I45°−I135°],(9)where (S0,S1,S2) is a linear Stokes vector of light; I0°, I45°, I90°, and I135° are intensities of light captured at polarizer angles of 0°, 45°, 90°, and 135°, respectively. The total incidence on camera I(x,y) and the DoLP p(x,y) can be calculated from the following equations: {I(x,y)=S0(x,y),p(x,y)=S12(x,y)+S22(x,y)S0(x,y).(10)

It should be noted that, if a time-divided polarization camera is used to capture the polarization images at polarization angles of 0°, 45°, 90°, and 135°, then the time interval of each of snap will need to be as fast as possible, because the intensity changes during the snaps usually lead to an abnormally high DoLP [according to Eqs. (9) and (10), if the differences among I0°, I45°, I90°, and I135° are high, then the p(x,y) will be high], which may cause a low-quality dehazed image.

The DoLP is a key parameter in the polarimetric dehazing algorithm. Normally, the DoLP is calculated by the Stokes method. However, the DoLP is sensitive to noise. Therefore, in order to address this problem for improving the robustness of the algorithm, we propose a novel way to estimate the DoLP, i.e., estimating the DoLP via low-rank approximation.

Similar to the conventional way to calculate the DoLP, the proposed method also needs four polarized images that are captured at polarizer angles of 0°, 45°, 90°, and 135°. In mathematics, the images are matrices. In this paper, the matrices corresponding to the polarized images with size of m×n are noted as I0°, I45°, I90°, and I135°, respectively. Based on the matrices, a matrix M is defined as M=[vec(I0°)vec(I45°)vec(I90°)vec(I135°)],(11)where vec(·) is the vectorization operator, which reshapes a matrix to a column vector in column-major order. It is obvious that the size of M is mn×4. Figure 2(a) shows the construction of M figuratively. It should be emphasized that, if there are no noises in the images, and the captured lights are linearly unpolarized (i.e., DoLP is zero), the rank of M will be one, since I0°, I45°, I90°, and I135° are the same in such a condition. However, in practice, M is a full column rank in general because there always are noises and the captured lights are partially linearly polarized. Therefore, we can assume that the M satisfies M=M′+P+N,(12)where M′ is the linearly unpolarized component in M; P is the fully linearly polarized component; N is the noise component. We assumed that there is no noise in M′ and P. As mentioned before, the matrix of the linearly unpolarized component without noises, i.e., M′, should be rank 1. On the other hand, in general, the DoLP in hazy scene is lower than 0.01 [17], and the intensities of noises are far lower than the intensity of captured light if the illumination was sufficient. So, in the general case, ‖M′‖F≫‖P+N‖F, or ‖M′‖F≈‖M‖F equivalently. These properties of M and M′ provide us with an idea to estimate M′ by finding a rank-1 approximation of M through M′≈argminX∈Rmn×4,rank(X)=1‖M−X‖F,(13)where mn is the total number of pixels of the polarized image. Such an M′ obtained from Eq. (13) is the rank-1 approximation of M under the Frobenius norm. As mentioned above, we assumed that there is no noise in M′, so a low-noise image of the linearly unpolarized part, noted as Iu, can be obtained from M′ via averaging columns of M′ and reshaping the result, as shown in Fig. 2(b).

Once the intensity of the linearly unpolarized part (i.e., Iu) is obtained, the DoLP can be calculated from the following equation: pE=I(x,y)−2Iu(x,y)I(x,y),(14)where pE is the estimated DoLP, and Iu(x,y) is the intensity of the linearly unpolarized part at (x,y) in the imaging plane. In practice, a value of Iu(x,y) is corresponding to an element in Iu. However, the estimated DoLP, i.e., pE, may be different from the real DoLP on scale. Therefore, to keep the scales consistent, we use the following equation to rescale the estimated DoLP: p(x,y)=p¯Sp¯EpE(x,y),(15)where p¯S is the mean value of the DoLP calculated from the traditional Stokes method, and p¯E is the mean value of pE(x,y). In the proposed dehazing algorithm, the rescaled estimated DoLP, as an estimation of the real DoLP, will substitute for the DoLP calculated from the traditional Stokes method. It should be emphasized that, although there is noise in the DoLP estimated by the conventional Stokes method, the mean value of such a DoLP, i.e., p¯S, is still reliable because the averaging operation can suppress noise. Therefore, although p(x,y) is obtained based on the DoLP estimated by the low-rank approximation and the conventional Stokes method simultaneously, such a p(x,y) can also be regarded as low-noise.

In the end, as a supplement, it is worth to note that M′ can be calculated stably via the compact singular value decomposition of M: M=UΣVT,(16)where U and V are semi-orthogonal and orthogonal matrices with sizes of mn×4 and 4×4, respectively; Σ is a 4×4 diagonal matrix that satisfies Σ=diag(σ1,σ2,σ3,σ4), where σ1 to σ4 are four non-zero singular values of M from large to small. It can be proved that the M′ can be calculated from the following equation [21]: M′=σ1uvT,(17)where u and v are singular vectors corresponding to the maximal singular value σ1. The main computational cost for calculating M′ includes the cost for calculating Eq. (17) and the cost for the singular value decomposition. We conducted a simple experiment to measure the time cost with MATLAB on a processor of 13th Gen Intel(R) Core(TM) i9-13900H with 32 GB RAM, where M is generated randomly, and a built-in function svds(·) is used to calculate the singular value decomposition of M. The result shows that the time complicity for calculating M′ is approximately equal to O(m) when mn∈[106,107]. When the size of a single polarization image is 1000×1000 (mn=106), the time cost is about 0.043 s.

In this paper, the SDR and HDR images are defined in a specific way: all the grayscale values of the SDR image should be in 0–255, while the grayscale value of the HDR image can be higher than 255. Because of the bigger grayscale range, an HDR image can show higher contrast. In general, dehazing results generated from the polarimetric dehazing algorithm are such HDR images with high contrasts; this is due to the fact that the image arithmetic involved in the algorithm usually leads to grayscales exceeding 255. However, a high contrast does not represent a high quality. In fact, an HDR image with an extremely high contrast is low-quality if the image is shown on a monitor that does not support the HDR display. Therefore, it is necessary to compress the dynamic range of a dehazed image to be within an appropriate range.

In order to generate a high-quality SDR dehazed image, we proposed an MVEF scheme to compress the dynamic range of the HDR dehazed image. In principle, such an MVEF scheme is inspired by the idea of the exposure fusion. Exposure fusion is an imaging technique to obtain a high-quality SDR image for direct display, which fuses a bracketed exposure sequence that consists of a series of images for the same scene but captured with different exposure times. As usual, the image via exposure fusion shows a higher dynamic range than the imaging dynamic range of the camera itself. Note that in our MVEF scheme, based on an HDR dehazed image we creatively construct multiple SDR dehazed images, which, equivalent to a bracketed exposure sequence, will be used for the subsequent exposure fusion. This is the reason why the proposed scheme is called the multiple virtual-exposure fusion (MVEF).

A flow chart of the exposure fusion is shown in Fig. 3. First, we capture the objective scene several times with different exposure levels to form a bracketed exposure sequence [22,23]; second, we calculate weight maps for each image referred to as qualities of the image, which include the contrast, the saturation, and the good exposedness; finally, we fuse images in the bracketed exposure sequence into one image with the guidance of weight maps. The detailed algorithm used in these procedures can be referred to in the original work [24].

The proposed MVEF scheme uses the exposure fusion to generate a high-quality SDR dehazed image from an HDR dehazed image. In the scheme, the bracketed exposure sequence used to generate the SDR dehazed image is generated from the HDR dehazed image, instead of a direct physical acquisition. Principally, we generate the bracketed exposure sequence by saturating bright areas of the HDR dehazed image.

At first, we calculate four saturation threshold values based on the grayscale histogram of the HDR dehazed image. The four values are 40th, 60th, 85th, and 100th percentiles of grayscale values. Figure 4 shows an example, where the four percentiles of grayscale values of the HDR dehazed image are 162, 218, 358, and 598, respectively. Further, as shown in Fig. 5, for each saturation threshold value, we saturate bright pixels in the HDR dehazed image: if the grayscale value of a pixel is larger than the saturation threshold, we will set the grayscale value to the saturation threshold value. After the saturation operation, the grayscale range of a saturated image will be linearly rescaled to 0–255. Therefore, the four saturation threshold values will generate four saturated SDR images. Finally, the four saturated SDR images, as the generated bracketed exposure sequence, will be inputted into the exposure fusion algorithm to obtain the SDR dehazed image.

In fact, the four saturated dehazed images are similar to images captured with different exposure times, and the corresponding four saturation threshold values are similar to the saturated intensities of a sensor under different exposure times. We use such a digital multiple saturation operation to simulate the multiple physical exposure. The exposure fusion algorithm can extract well-exposed parts of the four images and fuse these parts to form an SDR dehazed image with high quality. Therefore, the MVEF scheme is a way to fully utilize the information in the HDR dehazed image.

Figure 6 shows the flow chart of the proposed polarimetric dehazing algorithm. Inputs of the algorithm are polarized images captured at polarizer angles of 0°, 45°, 90°, and 135°. At first, we estimate DoLP of the scene via the low-rank approximation method introduced in Section 2.B. In the process, we divide the images into m×n blocks with equal size and estimate the DoLP by blocks. In our experiment, m and n are set to 100. The estimated DoLP of the scene can be obtained by splicing all the DoLP blocks. Further, based on the estimated DoLP and the polarization images, the object light L(x,y) can be calculated from equations introduced in Section 2.A. At this time, the image of the object light is an HDR dehazed image. As the final procedure of the proposed dehazing algorithm, the HDR dehazed image will be converted to an SDR dehazed image via the MVEF introduced in Section 2.C. The SDR dehazed image is the final output of the proposed dehazing algorithm.

In this section, dehazing experiments are conducted to show the effectiveness of the proposed method. In the experiment, some images captured in dense foggy weather are dehazed by the proposed method, and the quality of dehazed images is quantitatively evaluated to provide objective comparative references. Besides, to show the superiority of the proposed method, we also use other dehazing methods to dehaze the images. The dehazing methods used for comparison include the proposed method without the MVEF scheme, a dark channel prior method [25], and the Schechner’s method [11]. In the proposed method without the MVEF scheme, we still use low-rank approximation to estimate the DoLP, but the HDR dehazed image is converted to the SDR dehazed image by simple linear rescaling (i.e., rescale the range of grayscales to 0–255), instead of using the MVEF scheme. On the other hand, the dark channel prior method is a classical computer-vision-based dehazing method, and the Schechner’s method is a canonical polarimetric dehazing method.

In the experiment, some referenceless evaluation indictors are used to evaluate the quality of dehazed images quantitatively. The indictors are standard variance, information entropy, score of the natural image quality evaluator (NIQE) [26], and score of the perception-based image quality evaluator (PIQE) [27]. The standard variance is a reference of contrast of an image; a larger standard variance means a higher contrast. The information entropy is a reference of the amount of information of an image; a higher information entropy means more information. Both NIQE and PIQE are machine-learning-based image quality evaluators; each of them calculates an evaluative score based on the distortion and the perceptual quality of an input image. A lower score of NIQE or PIQE represents a higher quality of an image.

Figure 7 shows the seven hazy images for experiment and the corresponding dehazing results, where the hazy image, the dehazed image by the proposed method, the dehazed image by the proposed method without MVEF, the dehazed image by the dark channel prior method, and the dehazed image by the Schechner’s method form an experimental group. Original polarization images of experimental groups from (a) to (e) are captured by a general industrial camera with an external rotatable polarizer. The polarizer is separated with the camera, and the camera is well-fixed, so the rotation of the polarizer does not cause a significant instability on imaging. On the other hand, the aperture of the polarizer is larger than the aperture of the camera lens, so a minor displacement of the polarizer during the rotation also does not cause a significant instability. Original polarization images of experimental groups (f) and (g) are provided by Ref. [28], which are captured by a polarization camera developed by Hefei Shizhan Photoelectric Technology Company.

It should be clarified that the polarization images provided by Ref. [28] are captured at polarization angles of 0°, 60°, and 120°. However, the proposed method needs four polarization images at polarization angles of 0°, 45°, 90°, and 135° to dehaze. Therefore, we convert the three polarization images to the needed four polarization images according to relative principles in polarization optics (note that the linear Stokes vector can be solved from the three images at 0°, 60°, and 120°, and then we can predict the needed four images based on the solved linear Stokes vector). Theoretically, the key parameter of polarization used in the proposed method, i.e., the DoLP, is invariant after such a conversion.

Tables 1–4 show the standard variance, the information entropy, the score of the NIQE, and the score of the PIQE of images shown in Fig. 7, respectively, where optimal values of each indicator in each experimental group are in bold. Objectively, according to the tables, the proposed method with or without the MVEF is better than the other two methods in most cases. Besides, the better performance on experimental groups (g) and (h) shows a good transferability of the proposed method, i.e., the proposed method can be applied on different configurations of polarization angles.Table 1.

Standard Variances of Dehazing Results

At first, let us compare the proposed method with and without the MVEF scheme. The dehazed image by the proposed method without the MVEF can be regarded as a direct show of the HDR dehazed image, because we only rescale the grayscale range of the HDR dehazed image to 0–255 in this case. According to Fig. 7, qualitatively, if the dynamic range of the original HDR dehazed image is within an acceptable range, such as those shown in experimental groups (a), (b), and (c), the MVEF scheme cannot improve the image quality significantly because the visibility of objects in such images is sufficient. But, in experimental groups (b) and (c), the visibility of buildings at close range is still enhanced by the MVEF. On the other hand, if the dynamic range of the HDR dehazed image is too high, such as those shown in experimental groups (e), (f), and (g), the MVEF scheme will improve the visibility of dark areas significantly. Therefore, the MVEF scheme plays an indispensable role for improving the robustness of the dehazing algorithm.

Quantitatively, according to Tables 1 and 2, the MVEF scheme can improve the contrast and the amount of information of dehazed images in general. It is reasonable because the MVEF can enhance the visibility of dark areas. However, the cost of such a visibility enhancement is more noises: noises in the dark areas are enhanced simultaneously. It can be seen in Tables 3 and 4 that the two noise-related indictors, i.e., scores of NIQE and PIQE, may be degraded after the application of the MVEF. Although there may be more noises in the dehazed image with the MVEF, the MVEF scheme is still necessary because the scheme can make an almost invisible target be significantly visible; in addition, the scores of NIQE and PIQE are improved after the application of the MVEF in most cases due to the enhancement of the visibility. Therefore, overall, the MVEF scheme plays an indispensable role for improving the robustness of the dehazing algorithm.

After the discussion about the MVEF scheme, let us compare the proposed method and the other two classical methods, i.e., dark channel prior method and Schechner’s method. Note that, as for the experimental group (c) shown in Table 1, the standard variance of the dehazed image by Schechner’s method is significantly higher than that of the dehazed image by the proposed method, meaning that the contrast of the image dehazed by the Schechner’s method is higher. However, such a high contrast is reached at the cost of image details. Figure 8(a) shows the sky areas of images in group (c). It can be seen that the sky area of the dehazed image by Schechner’s method seems over-exposed. Therefore, as we mentioned in Section 2.C, a high contrast does not represent a high quality. It is necessary to make the contrast be an appropriate range for higher image quality. Besides, sky area of the dehazed image by the dark channel prior method suffered from severe noise and color distortion; meanwhile there is less noise and color distortion in the dehazed image by the proposed method. Meanwhile, there also is not a significant color distortion in the dehazed images by the Schechner’s method. In fact, polarization-based dehazing methods inherently have the ability to reduce color distortion. The color polarization images used to dehaze include extra physical information of the three color channels, which ensures us being able to precisely retrieve the object lights with different wavelengths.

Another noteworthy case is the pylon in images of the experimental group (a), which are shown in Fig. 8(b). It can be seen that the pylon in the dehazed image by the dark channel prior method is glowing, which is unnatural. A similar case is Fig. 8(c), where the pylon is still glowing in the dehazed image by the dark channel prior method. In fact, how to handle the glow effect is a scientific problem in the field of haze removal [29,30]. The glow effect is caused by a fuzzy estimation of intensity of air light in general. The polarimetric dehazing method can suppress such an effect because the intensity of air light can be more accurately estimated from polarimetric information. Therefore, from Figs. 8(b) and 8(c), it can be seen that there is no glowing effect in the dehazed image by the proposed method or Schechner’s method. However, the dehazed image by the proposed method is significantly clearer than by the Schechner’s method, which verified the superiority of the proposed method.

In Section 2.B, we mentioned that, compared to traditional Stokes methods, the DoLP estimated by the low-rank approximation method can reduce noise and improve the robustness of the algorithm. Now let us discuss how the method works and how superior the method is.

As an intuitive comparison, Figs. 9(b) and 9(c) show estimated DoLP images of the same scene by the low-rank approximation method and the Stokes method, respectively. The hazy scene corresponding to the DoLP images is shown in Fig. 9(a). Dehazed images based on the estimated DoLP images by the low-rank approximation method and the Stokes method are shown in Fig. 9(d) and Fig. 9(e), respectively. It can be seen from Figs. 9(b) and 9(c) that the estimated DoLP by the low-rank approximation method is smoother, because the range of DoLP in Fig. 9(b) is about 0.029 to 0.033 while that in Fig. 9(c) is 0 to 0.242. In the optical degradation model, the object light is assumed to be fully unpolarized, meaning that such a high DoLP of 0.2421 violates the basic assumption of the model. In fact, high DoLPs in Fig. 9(c) are caused by registration errors in the inputted four polarization images, i.e., positions of the same object in the four polarization images are not the same. Therefore, from Fig. 9(c) we can see that the DoLP at edges of buildings is high, just like what is shown in Fig. 9(g). The abnormally high DoLPs finally lead to errors in the dehazed image, which can be seen in Fig. 9(i), where the edge is abnormally bright. However, because of the smoother DoLP in Fig. 9(b), such errors are not shown in Fig. 9(d), which can be seen in Fig. 9(h).

On the other hand, in the polarimetric dehazing algorithm, DoLP is used to estimate the intensity of scattering light. Figures 10(a) and 10(b) show the estimated intensity map of scattering light based on DoLP shown in Fig. 9(b) and Fig. 9(c), respectively. Obviously, in Fig. 10(a), the intensity of scattering light is positively correlated with the depth of the scene and the intensity of object light, which accords with basic assumptions in the optical degradation model. However, this point cannot be seen in Fig. 10(b). There are heavy noises in Fig. 10(b) because of the rough DoLP distribution in Fig. 9(c). Besides, because of the abnormally high DoLP shown in Fig. 9(g), the intensity of scattering light at the edge in Fig. 10(b) is also abnormally high, which can be seen in Fig. 10(d). In the end, Fig. 10(e) shows intensity curves of pixels on the line of y=400 in Figs. 10(a) and 10(b). It can be seen that the curve of Fig. 10(a) is smoother than the curve of Fig. 10(b), meaning that there are less noises in Fig. 10(a).

As a summary of the above discussion, noises and errors in the inputted four polarization images significantly reduce the accuracy of the estimation of DoLP and intensity of scattering light by the Stokes method. However, such a negative impact is suppressed by the low-rank approximation method efficiently, because we can obtain a smoother estimation of DoLP and intensity of scattering light by the method. Therefore, the low-rank approximation method can reduce noises and improve robustness of the dehazing algorithm.

Gamma correction is a common dynamic range compression method [31,32], which applies a power function on the grayscale value of each pixel to adjust the dynamic range. In this section, a comparative experiment is conducted to illustrate the effectiveness and superiority of the MVEF scheme, where the gamma correction is the baseline.

The formula of gamma correction is shown as follows: Iout(x,y)=Iinγ(x,y),(18)where Iin(x,y) is the normalized pixel value of the input image; Iout(x,y) is the pixel value of the output image of gamma correction; γ is an adjustable parameter. In practice, the input image is the HDR dehazed image. In order to compress the dynamic range of the input image, γ should be set between zero and one. Obviously, in this case, the smaller the γ, the narrower the dynamic range.

Figure 11 shows a hazy image, an HDR dehazed image, an SDR dehazed image obtained by the MVEF, and SDR dehazed images obtained by gamma correction with different values of γ. In the HDR dehazed image, grayscale values on the roof area are so high that other areas appear dark visually. Therefore, in order to reduce the differences of intensity between the roof and other areas, and to improve the visual brightness of other areas, we need to compress the dynamic range of the HDR dehazed image. It can be seen from Figs. 8(d)–8(g) that if the value of γ is too high (for example, γ=0.8), the non-roof areas will still be visually dark, meaning that the dynamic range will be not compressed sufficiently. Conversely, if the value of γ is too low (for example, γ=0.2), the image will become hazy again, meaning that the contrast of the image will be too low. However, medium values of gamma (γ=0.6 or 0.4) still do not balance the degree of dynamic range compression and the contrast: compared to Fig. 11(c), the non-roof areas in Fig. 11(e) are relatively dark; meanwhile, Fig. 11(f) is relatively hazy.

For quantitative comparison, Fig. 12 shows standard variances of images shown in Fig. 11, where more values of γ are considered. It can be seen that the standard variance of the SDR dehazed image by the MVEF is slightly higher than that of the HDR image, meaning that using the MVEF to compress the dynamic range can preserve image details efficiently. However, it also can be seen that the standard variances of SDR dehazed images by the gamma correction are always lower than the HDR image (or equal to the HDR image), meaning that image details are lost with the compression of dynamic range. Therefore, regardless of the value of γ, the gamma correction cannot reach the same performance as the MVEF.

In this section, to show the novelty and superiority, the proposed method is compared to our recent work [17] in the polarimetric dehazing field with comparative experiments.

A main difference between the proposed method and our previous method is the introduction of the MVEF scheme. In our previous method, to ensure the brightness of objects, we use a simple strategy to convert the HDR dehazed image to the SDR dehazed image: set the average grayscale value to 120, and then set grayscale values higher than 255 to 255. So, there is not a dynamic mechanism (similar to the MVEF) to adjust the dynamic range of dehazed images in our previous work. Figure 13 shows dehazed images by the proposed method and the previous method, where original hazy images are from experimental groups shown in Fig. 7. It is obvious that sky regions in the dehazed images by the previous method seem over-exposed, such as shown in experimental groups (a), (b), (c), and (e), where the other regions are relatively dark. The reason for such a phenomenon is the excessive dynamic range of the HDR dehazed image, where the intensity difference between the sky area and other areas is too high. The simple strategy used in the previous method to convert the HDR to the SDR cannot dynamically adjust the dynamic range, meaning that the previous method cannot handle various illuminating conditions. On the other hand, because of the MVEF scheme, the proposed method can appropriately adjust the brightness of different areas in dehazed images, which ensures the visibility of all objects shown in the image, no matter what the dynamic ranges are. Therefore, the introduction of the MVEF scheme is an important improvement.

Another main difference between the proposed method and the previous method is the denoising method. In the previous method, we use a low-pass filter to suppress noises, instead of the low-rank approximation method used in the proposed method. However, a well-known fact is that the low-pass filter usually leads to image blurring. Figure 14 shows a hazy image, intensity maps of scattering light obtained from the low-rank approximation method and the low-pass filtering method, and corresponding dehazed images. It can be seen in Fig. 14(c) that the intensity maps of scattering light obtained from the low-pass filtering method are severely blurred. Although noises in the intensity map of scattering light are well suppressed by the low-pass filtering method [the noisy intensity map of scattering light is shown in Fig. 10(b)], the subsequent image blurring leads to the glow effect. It can be seen in Fig. 14(g) that the edge of the building is glowing. As we mentioned in Section 3, the glow effect is caused by a fuzzy estimation of the intensity of air light (scattering light) in general. Therefore, the glow effect is inevitable for the previous method. Figure 15 shows local areas in dehazed images shown in Fig. 13, where edges of buildings are glowing in the dehazed image by the previous method. On the other hand, because of a more acute estimation of the intensity of scattering light, the proposed method inherently has the ability to prevent the glowing effect. Therefore, it can also be seen in Fig. 14(f) and Fig. 15 that there is not such a glowing effect in the dehazed images by the proposed method.

In this paper, by using the low-rank approximation method to suppress the noise inherently existing in the calculated DoLP image and using the MVEF scheme to compress the dynamic range of the calculated dehazed image in which usually a large number of pixels have grayscale values beyond 255, a novel and robust polarimetric dehazing algorithm is proposed to enhance the quality of hazy images. Comparative experiments show that the low-rank approximation method can suppress noise efficiently, the MVEF scheme method can keep more image details after the compression of the dynamic range, and in most cases the qualities of dehazed images by the proposed algorithm are better than that of the other two classical algorithms (i.e., the dark channel prior method and the Schechner’s method) and our previous method.