Active Polarization Imaging for Cross-Linear Image Histogram Equalization and Noise Suppression in Highly Turbid Water

Abstract

The absorption and scattering of impurity particles in turbid water cause the target signal light to be attenuated and to produce backscattered light, resulting in the reduced quality of underwater polarimetric imaging. As water turbidity increases, the effect of backscattered light becomes greater, making polarization imaging in highly turbid water a challenge. Theory and experiment show that the increase in the intensity of backscattered light leads to high noise gain in the underwater active polarization imaging model. In order to enhance image contrast and suppress noise gain in highly turbid water, we propose an underwater imaging enhancement method that appropriately combines the non-physical and physical models. The method uses contrast limited adaptive histogram equalization (CLAHE) for a certain number of cross-linear images (I_min) before calculating their polarization enhancement images, and it constructs joint filtering (multi-frame averaging and bilateral filtering) to suppress the high noise gain introduced by the imaging model and CLAHE. The experimental results in highly turbid water validate the rationality and feasibility of the proposed method, and the comparative processing results (52.7~98.6 NTU) outperform those of the conventional non-physical and physical model methods. The method maintains the complexity of the system and facilitates the application of conventional polarimetric imaging in harsher underwater environments.

1. Introduction

The underwater environment contains a large number of impurity particles, and light will be scattered and absorbed by the particles when propagating underwater, resulting in a reduction in the quality of underwater imaging, mainly in the form of reduced image brightness and loss of contrast. The turbidity in natural waters typically varies in the range of a few dozen to a few tens of nephelometric turbidity units (NTU) [1,2], and water is considered highly turbid when it is above 50 NTU [1]. Some underwater tasks [3,4,5] (e.g., submarine microbial imaging, pipeline maintenance, and shipwreck rescue) can rapidly increase the turbidity of the water and maintain it for a period of time due to the particles that are stirred up when operating close to the bottom area. This significantly reduces the effective distance of underwater imaging and severely hampers the execution of underwater missions. In recent years, in order to improve the quality of underwater imaging by reducing the effects of underwater backscattered light, different de-scattering imaging methods have been pioneered and are mainly divided into non-physical and physical model-based methods [6,7,8,9]. They provide strategies to address the different needs of underwater imaging accordingly.

Among the non-physical methods, the contrast limited adaptive histogram equalization (CLAHE) method is a common operation for underwater de-scattering [3,10]. Compared to using conventional histogram equalization (HE) to improve the operation of the dynamic range of the greyscales in the raw image, CLAHE operates on small regions of the image that are tiled, and it uses bilinear differences to reduce the effect of boundaries between the small regions, resulting in a more stable histogram equalization [11]. Thus, CLAHE can effectively enhance the contrast of the target signal light in the raw image. Among the physical model-based methods, the underwater polarimetric imaging model proposed by Schechner et al. is widely used for underwater imaging due to its compactness, ease of piggybacking, and operation [5,6,12]. It also belongs to the category of a time-sharing rotating polarization imaging system [3,13]. In order to obtain a better de-scattering effect, many optimization measures have been developed around the Schechner model, using, for example, optical correlation [14,15,16] and extrapolation fitting [17] and finding feasibility intervals [18] to optimize the selection of relevant parameters, as well as combining Mueller matrices [19,20] and computer algorithms [21,22,23] to suppress the backscattered light. In highly turbid water, Liu et al. found that red light was more effective than other wavelengths in suppressing scattering, and they combined this with active polarization imaging to make the leap from “undetectable” to “detectable” [1]. At the same time, Hu et al. cleverly applied the memory effect of circularly polarized light to underwater polarization imaging and realized target detection in highly turbid water through the combined action of linearly and circularly polarized light [24]. In contrast to Liu’s and Hu’s methods of modulating the light source and increasing the complexity of the imaging system, Li et al. used a histogram stretching operation with greyscale normalized for a cross-linear image (I_min) and combined it with polarimetric imaging to achieve experimental results in highly turbid water that were superior to those of the traditional methods [25]. Considering the high backscattered light intensity [26] and high noise gain characteristics in highly turbid water [6], the raw image will have overexposed regions due to backscattered light enhancement, and CLAHE will be more limited in its processing. In addition, normalized histogram stretching can weaken the contrast enhancement effect of the target signal light under high noise conditions. Therefore, it is necessary to design a method for polarization enhancement in highly turbid water without increasing the complexity of the system.

In this paper, the problem of strong backscattering and the high noise gain of the imaging model caused by it under high-turbidity conditions is discussed theoretically and experimentally. Based on the presence of two influencing factors, we appropriately combine the non-physical and the physical models and propose an active polarization imaging method based on the CLAHE histogram equalization of the cross-linear image and the joint noise suppression. CLAHE is used to equalize the intensity distribution of the target signal light in I_min. Meanwhile, the polarization enhancement results are processed using a combination of multi-frame averaging [27,28] and bilateral filtering [29] to counteract the high noise gain generated by the imaging model and CLAHE. Only a certain number of I_min need to be acquired in our method to achieve the multi-frame averaging effect, and experimentally, we discuss the dependence of the noise suppression effect on the number of averaged frames. Finally, real-world experiments in different water samples with high turbidity are set up to verify the effectiveness and advantages of our method.

2. Methodology

2.1. Underwater Polarization Imaging Model

Turbid water bodies contain a large number of impurities, including soluble substances as well as insoluble particulate matter and microorganisms. The light shining on the object in the turbid water is absorbed and scattered by the water molecules and impurities, resulting in the received image degradation. The attenuation of the target signal light intensity and the backscattered light covering the object result in a loss of contrast in the image. According to the underwater imaging model, the received image is the incoherent sum of the target signal light and the backscattered light [6].

$$I\left(x,y\right)=S\left(x,y\right)+B\left(x,y\right)$$

(1)

where $S\left(x,y\right)$ denotes the target signal light, derived from the irradiance of the target object. $B\left(x,y\right)$ denotes the backscattered light. When $S\left(x,y\right)$ propagates in the water, it is affected by two interfering factors; one is the absorption of water molecules and impurities, resulting in its intensity attenuation; the other is caused by the scattering of impurities by the backscattered light $B\left(x,y\right)$ superimposed on $S\left(x,y\right)$, resulting in a reduction in the image contrast, ultimately leading to a decline in imaging quality. Backscattered light has the characteristics of partial polarization; so, rotating the analyzer before the receiver will obtain the darkest (cross-linear) and brightest (co-linear) two orthogonal polarization images, noted as I_min and I_max. According to Equation (1), $I_{\min }$ and $I_{\max }$ can be expressed as:

$$\left\{\begin{array}{c}{I}_{\min }\left(x,y\right)=S_{\min }\left(x,y\right)+B_{\min }\left(x,y\right)\\ I_{\max }\left(x,y\right)=S_{\max }\left(x,y\right)+B_{\max }\left(x,y\right)\end{array}\right.$$

(2)

$S_{\min }\left(x,y\right)$ and $S_{\max }\left(x,y\right)$ denote the darkest and brightest target signal light in the orthogonal polarization images, and $B_{\min }\left(x,y\right)$ and $B_{\max }\left(x,y\right)$ denote the darkest and brightest backscattered light in the orthogonal polarization images. As partial polarization light can be decomposed in two orthogonal directions, the total light intensity $I\left(x,y\right)$ can be expressed as:

$$I\left(x,y\right)=I_{\min }\left(x,y\right)+I_{\max }\left(x,y\right)$$

(3)

According to the definition of the degree of polarization (DOP), the DOP of the target signal light $p_{\mathrm{targ}}$ and the backscattered light $p_{\mathrm{scat}}$ can be calculated from the corresponding components in $I_{\min }\left(x,y\right)$ and $I_{\max }\left(x,y\right)$,

$$\left\{\begin{array}{c}{p}_{\mathrm{targ}}=\frac{S_{\max }-S_{\min }}{S_{\max }+S_{\min }}\\ p_{\mathrm{scat}}=\frac{B_{\max }-B_{\min }}{B_{\max }+B_{\min }}\end{array}\right.$$

(4)

The expressions for the target signal light and the backscattered light can be obtained by combining Equations (1)–(4),

$$\left\{\begin{array}{c}S\left(x,y\right)=\frac{1}{p_{\mathrm{scat}}-p_{\mathrm{targ}}}\left[I_{\min }\left(1+p_{\mathrm{scat}}\right)-I_{\max }\left(1-p_{\mathrm{scat}}\right)\right]\\ B\left(x,y\right)=\frac{1}{p_{\mathrm{scat}}-p_{\mathrm{targ}}}\left[I_{\max }\left(1-p_{\mathrm{targ}}\right)-I_{\min }\left(1+p_{\mathrm{targ}}\right)\right]\end{array}\right.$$

(5)

2.2. Noise Analysis of Polarization Imaging in High-Turbidity Water

There are three main components of the noise sources in the imaging process; these are readout noise, dark current noise, and photon noise. Define $g_{\mathrm{electr}}$ as the number of photogenerated electrons in the imaging system that cause a change per unit of gray level. Therefore, the noise variance of the readout noise can be represented by ${\rho }^2/g_{\mathrm{electr}}^2$, where $\rho$ is the standard deviation (STD) of the electronic readout noise of the imaging system. Dark current noise is related to the exposure time of the imaging system, and its noise variance is expressed by $Dt/g_{\mathrm{electr}}^2$, where $D$ represents the magnitude of the dark current, and $t$ is the exposure time. Photon noise can be denoted as $I\left(x,y\right)/g_{\mathrm{electr}}$, which is the signal-related noise independent of the parameters and quality of the camera. Thus, for an image, the noise variance of its pixel gray level is [6,30],

$${\sigma }_I^2={\rho }^2/g_{\mathrm{electr}}^2+Dt/g_{\mathrm{electr}}^2+I\left(x,y\right)/g_{\mathrm{electr}}$$

(6)

The underwater polarization imaging model relies on two polarization orthogonal sub-images, $I_{\min }$ and $I_{\max }$, which are considered as two statistically independent measurements and are linearly related to $S\left(x,y\right)$ and $B\left(x,y\right)$ in Equation (5). Assuming that the noise variances of $I_{\min }$ and $I_{\max }$ are ${\sigma }_{\min }^2$ and ${\sigma }_{\max }^2$, respectively, and combining Equations (5) and (6), we can obtain the following relationship between the noise variance of the target signal light and the two polarization orthogonal sub-images,

$$\begin{array}{cc}{\sigma }_S^2& ={\left(\frac{1+p_{\mathrm{scat}}}{p_{\mathrm{targ}}-p_{\mathrm{scat}}}\right)}^2\left({\rho }^2/g_{\mathrm{electr}}^2+Dt/g_{\mathrm{electr}}^2+I_{\min }\left(x,y\right)/g_{\mathrm{electr}}\right)\hfill \\ & +{\left(\frac{1-p_{\mathrm{scat}}}{p_{\mathrm{targ}}-p_{\mathrm{scat}}}\right)}^2\left({\rho }^2/g_{\mathrm{electr}}^2+Dt/g_{\mathrm{electr}}^2+I_{\max }\left(x,y\right)/g_{\mathrm{electr}}\right)\hfill \\ & ={\left(\frac{1+p_{\mathrm{scat}}}{p_{\mathrm{targ}}-p_{\mathrm{scat}}}\right)}^2{\sigma }_{\min }^2+{\left(\frac{1-p_{\mathrm{scat}}}{p_{\mathrm{targ}}-p_{\mathrm{scat}}}\right)}^2{\sigma }_{\max }^2\hfill \end{array}$$

(7)

For Equation (7), we consider it in two cases; first, we assume that the non-signal-related noise is dominant in the imaging process; that is, when ${\sigma }_{\min }^2={\sigma }_{\max }^2={\sigma }_0^2$, the noise variance of the target signal light becomes of the form,

$$\begin{array}{cc}{\sigma }_S^2& ={\left(\frac{1+p_{\mathrm{scat}}}{p_{\mathrm{targ}}-p_{scat}}\right)}^2{\sigma }_0^2+{\left(\frac{1-p_{\mathrm{scat}}}{p_{\mathrm{targ}}-p_{\mathrm{scat}}}\right)}^2{\sigma }_0^2\hfill \\ & =2{\sigma }_0^2\left[\frac{1+p_{\mathrm{scat}}^2}{{\left(p_{\mathrm{targ}}-p_{\mathrm{scat}}\right)}^2}\right]\hfill \end{array}$$

(8)

From Equation (8), it can be seen that the noise in the target signal light gains constantly, and the magnitude of the noise variance is related to $p_{\mathrm{scat}}$ and $\Delta p$ ($\Delta p=p_{\mathrm{scat}}-p_{\mathrm{targ}}$, the difference between $p_{\mathrm{scat}}$ and $p_{\mathrm{targ}}$). The variation of the noise amplification ${\sigma }_S/{\sigma }_0$ in the target signal light with $p_{\mathrm{scat}}$ and $p_{\mathrm{targ}}$ is shown in Figure 1a. It can be seen that when the DOP of the backscattered light is $p_{\mathrm{scat}}=0$ and the DOP of the target signal light is $p_{\mathrm{targ}}=1$, the noise amplification of the target signal light reaches the minimum value of $\sqrt{2}$. However, when the two DOPs gradually approach, that is, $\Delta p\to 0$, the noise is continuously increased. Figure 1b,c show the changes in $p_{\mathrm{scat}}$, $p_{\mathrm{targ}}$, and $\Delta p$ measured in the experiment with the increasing water turbidity and the increasing imaging distance, respectively. $p_{\mathrm{scat}}$ and $p_{\mathrm{targ}}$ are calculated according to the DOP definition after pre-labeling the target and the background regions in the experiment. It can be seen that in the imaging process of the increasing turbidity and distance of the water, the value of $p_{\mathrm{targ}}$ has a tendency to gradually approach $p_{\mathrm{scat}}$; so, $\Delta p$ gradually decreases, which leads to the continuous increase in the noise amplification rate in the recovered target signal light.

Next, we consider the practical case where the photon noise is not negligible, i.e., ${\sigma }_{\min }\ne {\sigma }_{\max }$. From the expression for $S\left(x,y\right)$ in Equation (5), it is clear that the contribution of $p_{\mathrm{targ}}$ to $S\left(x,y\right)$ is a scaling factor that adjusts the global intensity. Here, we assume that $p_{\mathrm{targ}}=0$ [6], in which case Equation (5) can be simplified to,

$$\left\{\begin{array}{l}S\left(x,y\right)=\frac{1}{p_{\mathrm{scat}}}\left[I_{\min }\left(x,y\right)\left(1+p_{\mathrm{scat}}\right)-I_{\max }\left(x,y\right)\left(1-p_{\mathrm{scat}}\right)\right]\hfill \\ B\left(x,y\right)=\frac{1}{p_{\mathrm{scat}}}\left[I_{\max }\left(x,y\right)-I_{\min }\left(x,y\right)\right]\hfill \end{array}\right.$$

(9)

According to Equation (9), it is obtained that

$$\left\{\begin{array}{c}{I}_{\min }\left(x,y\right)=\frac{1}{2}\left[B\left(x,y\right)\left(1-p_{\mathrm{scat}}\right)+S\left(x,y\right)\right]\\ I_{\max }\left(x,y\right)=\frac{1}{2}\left[B\left(x,y\right)\left(1+p_{\mathrm{scat}}\right)+S\left(x,y\right)\right]\end{array}\right.$$

(10)

Through Equations (10) and (7), the noise variance expression of $S\left(x,y\right)$ with respect to the target signal light and backscattered light can be obtained,

$${\sigma }_S^2=\frac{1}{p_{\mathrm{scat}}^2}\left\{\left[2\left({\rho }^2/g_{\mathrm{electr}}^2+Dt/g_{\mathrm{electr}}^2\right)+\frac{S\left(x,y\right)}{g_{\mathrm{electr}}}\right]\right.\left(1+p_{\mathrm{scat}}^2\right)\left.+\frac{B\left(x,y\right)}{g_{\mathrm{electr}}}\left(1-p_{\mathrm{scat}}^2\right)\right\}$$

(11)

From Equation (11), we can see that the noise variance of the target signal light $S\left(x,y\right)$ is related to the intensity of the backscattered light $B\left(x,y\right)$. Therefore, it is beneficial to reduce the intensity of the backscattered light appropriately in the actual scene to improve the signal-to-noise ratio (SNR) of the target signal light. Figure 2 shows the raw images, $I_{\min }$ and $I_{\max }$, acquired under different high-turbidity environments and the DOP distribution maps calculated using Equation (4). In fact, as the turbidity of the water increases, the content of impurity particles in the medium continues to rise, resulting in an increased probability of the photons interacting with the impurity particles. It can be seen from the raw images and the $I_{\max }$ in Figure 2 that the backscattered light saturates the receiver response to light intensity when a certain turbidity is reached. In addition, as the turbidity increases, the large-sized particles in the water will also increase, resulting in an increase in the scattering coefficient of the medium, and the photons gradually lose their initial polarization state [26]. Consequently, the intensity of the backscattered light in $I_{\min }$ increases, as shown in $I_{\min }$ in Figure 2. At the same time, the DOP of the backscattered light decreases, as shown in the DOP distribution maps in Figure 2.

The noise analysis of the underwater polarization imaging model shows that there are two main sources of noise in highly turbid water; one is the change in the DOP, and the other is the increase in the intensity of the backscattered light. Together, they cause a reduction in the SNR of the recovered image in a highly turbid environment. In addition, it can be seen from I_min in Figure 2 that the depolarization phenomenon of the backscattered light in the highly turbid water makes the target signal in I_min gradually submerge in the backscattered light. Therefore, in order to improve the quality of polarization differential imaging in highly turbid water, two aspects need to be addressed; one is contrast enhancement, and the other is noise suppression.

2.3. CLAHE-Based Cross-Linear Image Histogram Equalization and Joint Noise Suppression

According to the analysis in the previous section, the polarization enhancement effect of the target signal light in highly turbid water is affected by two aspects; one is the increased backscattered light intensity, and the other is the high noise gain with the increase in turbidity. Aiming at these two effects, we propose an underwater polarization differential imaging enhancement method based on CLAHE-based cross-linear image histogram equalization and joint noise suppression to improve the imaging quality. Figure 3 shows the processing flowchart of the entire method.

First, adjust the analyzer in front of the receiver to a co-linear state and capture a co-linear image $I_{\max }$. Then, adjust the analyzer to the cross-linear state and continuously capture $N$ cross-linear images $I_{\min ,i}\left(i=1,2,\dots ,N\right)$. Using the captured $I_{\max }$ and $I_{\min ,i=1}$, and according to the definition of the DOP, the DOP distribution map can be obtained,

$$p\left(x,y\right)=\frac{I_{\max }\left(x,y\right)-I_{\min ,i=1}\left(x,y\right)}{I_{\max }\left(x,y\right)+I_{\min ,i=1}\left(x,y\right)}$$

(12)

Meanwhile, the $p_{\mathrm{scat}}$ and $p_{\mathrm{targ}}$ required for the polarization differential processing can be obtained according to Equation (4).

As the grayscale of the image in the turbid water is compressed to a very narrow range, especially when the water turbidity increases, the intensity of the backscattered light on the target signal light increases. This results in a grayscale overexposed area in the $I_{\max }$ image and the attenuation of the target signal light intensity in the $I_{\min }$ image, ultimately causing further compression of the gray space of the captured image. To mitigate the extent to which the image grayscale is compressed, a grayscale histogram equalization of $I_{\min ,i}\left(i=1,2,\dots ,N\right)$ is performed by using CLAHE to obtain the adjusted grayscale $I_{\min ,i}^{\mathrm{pro}}\left(i=1,2,\dots ,N\right)$,

$$I_{\min ,i}^{\mathrm{pro}}\left(i=1,2,\dots ,N\right)=\mathrm{CLAHE}\left[I_{\min ,i}\left(i=1,2,\dots ,N\right)\right]$$

(13)

It is worth noting that the $I_{\max }$ image acquired in highly turbid water contains a large number of grayscale overexposed areas, which is not conducive to the grayscale histogram equalization operation. Therefore, only $I_{\min ,i}\left(i=1,2,\dots ,N\right)$ needs to be processed during the actual processing. The DOP distribution map $p\left(x,y\right)$ reflects the polarization relationship between two polarization orthogonal images and contains the information on the target and the backscattered light. To preserve the polarization relationship [25], the processed $I_{\max }^{\mathrm{pro}}$ can be obtained using Equation (12),

$$I_{\max }^{\mathrm{pro}}\left(x,y\right)=\frac{1+p\left(x,y\right)}{1-p\left(x,y\right)}{I}_{\min ,i=1}^{\mathrm{pro}}\left(x,y\right)$$

(14)

Based on the already obtained $I_{\min ,1}^{\mathrm{pro}}$, $I_{\max }^{\mathrm{pro}}$, $p_{\mathrm{scat}}$, and $p_{\mathrm{targ}}$, the recovered target signal light can be obtained by using Equation (5). However, due to the high-turbidity conditions, the difference between the two DOPs gradually decreases, which reduces the effect of the common-mode suppression and leads to the poor recovery of the target signal light. To improve the image enhancement result, the values of $p_{\mathrm{scat}}$ and $p_{\mathrm{targ}}$ are adjusted by setting ${\epsilon }_1\in \left[0,\frac{1}{p_{\mathrm{scat}}}\right]$ and ${\epsilon }_2\in \left[0,\frac{1}{p_{\mathrm{targ}}}\right]$. The values of the two DOPs are optimal at this time when the values of ${\epsilon }_1$ and ${\epsilon }_2$ are taken so that the grayscale histogram of the recovered image $S\left(x,y\right)$ is the most uniform, i.e., $p_{\mathrm{scat}}^{\mathrm{optimum}}={\epsilon }_{\mathrm{scat}}^{\mathrm{optimum}}{p}_{\mathrm{scat}}$ and $p_{\mathrm{targ}}^{\mathrm{optimum}}={\epsilon }_{\mathrm{targ}}^{\mathrm{optimum}}{p}_{\mathrm{targ}}$.

With the known results given above, Equation (5) can be used to calculate the $N$ processed target signal light results obtained from $I_{\min ,i}^{\mathrm{pro}}\left(i=1,2,\dots ,N\right)$, as follows:

$$\begin{array}{cc}{S}_i\left(i=1,2,\dots ,N\right)& =\frac{1}{p_{\mathrm{scat}}^{\mathrm{optimum}}-p_{\mathrm{targ}}^{\mathrm{optimum}}}\left[I_{\min ,i}^{\mathrm{pro}}\left(i=1,2,\dots ,N\right)\left(1+p_{\mathrm{scat}}^{\mathrm{optimum}}\right)\right.\hfill \\ & \left.-I_{\max }^{\mathrm{pro}}\left(1-p_{\mathrm{scat}}^{\mathrm{optimum}}\right)\right]\hfill \end{array}$$

(15)

The $N$ $S_i\left(i=1,2,\dots ,N\right)$ obtained by using Equation (15) are then used with CLAHE to perform histogram equalization; so, there are,

$$S_i^{\mathrm{pro}}\left(i=1,2,\dots ,N\right)=\mathrm{CLAHE}\left[S_i\left(i=1,2,\dots ,N\right)\right]$$

(16)

Next, in order to reduce the gained noise in the target signal light, the $N$ histogram equalized $S_i^{\mathrm{pro}}\left(i=1,2,\dots ,N\right)$ is subjected to multi-frame averaging processing, and then, bilateral filtering is used,

$$\mathrm{BF}\left(S_{\mathrm{average}}^{\mathrm{pro}}\right)=\mathrm{BF}\left(\raisebox{1ex}{${\displaystyle \sum _{i=1}^N{S}_i^{\mathrm{pro}}}$}\!\left/ \!\raisebox{-1ex}{$N$}\right.\right)$$

(17)

Here, $\mathrm{BF}\left(\cdot \right)$ denotes the bilateral filtering operation, and $\mathrm{BF}\left(S_{\mathrm{average}}^{\mathrm{pro}}\right)$ is the result of polarization enhancement by equalizing the histogram of the cross-linear polarization image and combining the noise suppression processing.

3. Real-World Experiment and Results

3.1. Experimental Setup

The experimental setup is shown in Figure 4, showing the schematic diagram and a photograph of the experimental setup, respectively. The emitting end of the imaging system is composed of a white LED light source (Supfire, X90, China, wavelength in the visible range) and a rotatable linear polarizer (Sigma, USP-30C0.4–38, Japan). The illumination beam is transformed into a stable output of linearly polarized light after passing through the polarizer. The experimental environment is a glass water tank of 60 cm (length) × 35 cm (width) × 40 cm (height), and the liquid level during the experiment is 25 cm. Milk is mainly divided into skimmed milk, semi-skimmed milk, and whole milk, which contain casein molecules (size $0.04\sim 0.3 \mathsf{\mu }\mathrm{m}$) and fat globules (size $1\sim 20 \mathsf{\mu }\mathrm{m}$) of different sizes [14,31]. The scattering types of the milk are mainly divided into Mie scattering and Rayleigh scattering, which is similar to the scattering type of seawater. As the scattering coefficient of the same volume of skimmed milk is smaller [31], which facilitates the preparation of water with different turbidities, the experiment uses skimmed milk to simulate different underwater turbid environments. The illumination beam passes through the medium to the target object (the target is 20 cm away from the wall of the water tank at the transmitting end) and is reflected and finally reaches the receiving end of the imaging system. The receiving end consists of a rotatable analyzer (Sigma, USP-30C0.4–38, Japan) and a monochrome charge-coupled device (CCD) camera (Daheng Optics, R-125–30UM, China, wavelength response range $400\sim 1000 \mathrm{nm}$). During the acquisition process, we rotate the analyzer to acquire a co-linear polarization image and multiple cross-linear images, respectively. The acquisition time interval of the multi-frame images is 10 ms, and the image size is 964 (width) × 1292 (height) pixels.

3.2. Rationality and Feasibility Analysis of the Processing Flow

To verify the rationality and feasibility of the processing process in the proposed method, we added 200 mL of skimmed milk to the water tank; the turbidity of the water was 84.5 NTU. The analysis in Section 2.2 shows that as the turbidity of the water increases the backscattered light is continuously enhanced, which makes a large overexposure area in $I_{\max }$. With the existing histogram equalization operation, it would be difficult to achieve a better processing effect. Therefore, our method chooses to perform the histogram equalization operation on $I_{\min }$. At the same time, the continuous acquisition of $I_{\min }$ in a short period of time ensures the statistical characteristics of the noise distribution. Additionally, the DOP distribution map $p\left(x,y\right)$ in the whole scene can be considered as stable; so, we acquired 1 frame of $I_{\max }$ and 200 frames of $I_{\min }$ (with an acquisition interval of 10 ms).

According to the processing flow, the CLAHE operation is first performed on the captured $I_{\min }$. CLAHE itself belongs to the classical approach to dealing with de-scattering in the non-physical model [3,10], and here, we combine it appropriately with the classical physical de-scattering model. CLAEH is a variant of the adaptive histogram equalization (AHE) algorithm. Compared with AHE, its contrast amplification is limited, which can better preserve the overall contrast distribution of the raw image and effectively reduce the noise amplification problem [11]. However, the additional noise introduced by CLAHE still needs to be considered. Figure 5a,d show the raw intensity image and the normalized grayscale histogram of $I_{\min }$. Figure 5b,e show the results and the normalized grayscale histogram of $I_{\min }$ after processing using Li’s method [25], and Figure 5c,f show the results of the $I_{\min }$ processing using CLAHE [11]. From Figure 5, it can be seen that both Li’s method and the CLAHE method improve the contrast of the target signal in the image compared to the raw intensity of $I_{\min }$. The grayscale distribution ranges of the two methods are basically the same, but the grayscale distribution after CLAHE processing is more uniform.

It is worth noting that the image acquired in the highly turbid water contains a large amount of noise, which may cause some grayscale extreme points in the image. Therefore, the histogram stretching operation for the grayscale normalization of I_min using Li’s method will be disturbed. Figure 5g shows the row indexes of the grayscale values in Figure 5b that are through the minimum (green line) and the maximum (red line), and the presence of these extreme value points will interfere with the grayscale normalization effect of Li’s method. Therefore, we remove the extreme value points and then normalize to obtain Figure 5h, and the grayscale histogram is shown in Figure 5i. It can be seen that the grayscale distribution of the region to the right of the black dashed line in Figure 5i is higher than that in Figure 5d, which proves that removing the extreme points is effective for enhancing the image brightness. However, the removal of the extreme points is relatively troublesome and uncertain in the actual processing, which is not conducive to the automation of the processing process. Therefore, we use CLAHE to improve the contrast of the target signal light in $I_{\min }$ in the highly turbid water scene. It should be emphasized that the original scene of the DOP distribution map $p\left(x,y\right)$, which contains the polarization information of the water, was calculated before processing $I_{\min }$; so, it will not change the polarization relationship between $I_{\min }^{\mathrm{pro}}$ and $I_{\max }^{\mathrm{pro}}$ [25].

The reason for using the multi-frame averaging method as the first step of joint noise reduction is twofold: the first reason is that after the CLAHE process not only is the contrast of the target signal light improved, but also the noise in the image is gained, causing the overall peak signal-to-noise ratio (PSNR) of the image to decrease. Figure 6b,d show the results obtained using Li’s method and the $S_1^{\mathrm{pro}}$ obtained using our method, respectively, where Figure 6a shows the target image in clear water. Figure 6c,e show the spatial frequency spectrum of Figure 6b,d, respectively. It can be seen from the PSNR and the spectrogram that the noise amplification in the results containing the CLAHE processing step is significantly higher than that of Li’s method.

The second reason is that the polarization differential imaging itself belongs to the time-sharing rotating polarization imaging system [3,13], and the polarization images in the different directions are obtained by rotating the analyzer placed in front of the receiver. Hence, the imaging process itself does not exclude the noise suppression operation of the multi-frame averaging. At the same time, the method of multi-frame averaging has been widely used in reducing imaging noise [27,28]. In addition, multi-frame averaging preserves the high-frequency details of the image well and maintains the sharpness of the image compared to the methods that directly use low-pass filtering (such as Gaussian low-pass filtering [32]). It is worth noting that in our method, only a certain number of images $I_{\min }$ need to be continuously acquired to achieve the multi-frame averaging effect. We therefore choose multi-frame averaging as the first step in joint noise suppression. However, multi-frame averaging has certain limitations in the improvement of the image SNR. The SNR of $N$ images after multi-frame averaging is,

$$I_{\mathrm{SNR}}^N=\sqrt{N}\frac{I\cdot {\eta }_{\mathrm{QE}}}{\sqrt{I\cdot {\eta }_{\mathrm{QE}}+I_{\mathrm{RON}}{}^2+D\cdot t}}$$

(18)

Here, $I$ denotes the image signal, ${\eta }_{\mathrm{QE}}$ denotes the quantum efficiency of the optical signal, $I_{\mathrm{RON}}$ denotes the readout noise, and $D$ is the dark current, and the dark current noise is time-dependent. According to Equation (18), it can be known that the SNR can be improved by $\sqrt{N}$ times when the multi-frame averaging operation is performed on $N$ images. However, as $N$ increases, the SNR increases more and more slowly. Figure 7a,b show the results of the images acquired in clear water and with different noises artificially added (a) Gaussian random noise with mean $\mu =0$ and variance ${\sigma }^2=0.02$ and (b) Gaussian random noise with mean $\mu =0$ and variance ${\sigma }^2=0.2$). Figure 7c shows the PSNR curves after averaging over different numbers of frames, with the maximum number of frames averaged at 200. It can be seen that the PSNR rises slowly after $N$ reaches the value of 50. Figure 7d,e show the multi-frame averaging results of two images with different noise intensities added for $N=50$. Figure 7f–h show the results obtained from the actual experiment, which are consistent with the simulation results, and the improvement of the PSNR tends to be stable after the value of $N$ reaches 50. Comparing the enlarged part (green dashed area) in Figure 7f,h, it can be seen that the noise of the image is significantly suppressed. Therefore, the result of $N=50$ is used in the multi-frame averaging operation in the following experiments.

For the SNR analysis of multi-frame averaging, it is known that the ability of multi-frame averaging to suppress noise is limited. Bilateral filtering has the effect of double filtering, which can suppress noise on the one hand and preserve the edge strength of the object well on the other hand [29,33]. Therefore, the bilateral filtering operation is used as the last step of joint noise suppression. Figure 8 shows the polarization enhancement results after comparing the different processing methods; Figure 8a is the raw image intensity; Figure 8b is the processing result of Li’s method [25]; and Figure 8c is the processing result using Liu’s method [32]. Figure 8d–f are the results of $S_1^{\mathrm{pro}}$, 50 frames of $S_i^{\mathrm{pro}}$ multi-frame averaging $S_{\mathrm{average}}^{\mathrm{pro}}$, and the result $\mathrm{BF}\left(S_{\mathrm{average}}^{\mathrm{pro}}\right)$ after bilateral filtering on $S_{\mathrm{average}}^{\mathrm{pro}}$, respectively. It can be seen from the enlarged part in Figure 8f that the joint noise suppression method can effectively reduce the interference of noise and improve the clarity of the target signal. Figure 8g shows the pixel grayscale row index curves of the raw image intensity and the different processing methods, and the indexed pixel range is the yellow dashed line in Figure 8b. It can be seen from Figure 8g that our method better maintains the intensity distribution of the target signal and restores the details of the target signal compared to other methods.

3.3. Results in Different High-Turbidity Water Samples

In the experiment, the raw target image in clear water was first acquired as a reference image for calculating the PSNR. Then, 120 mL of skimmed milk was added to the water tank, at which time the initial turbidity of the medium in the water tank was 52.7 NTU, and the images were acquired. Then, the volume of skimmed milk was continuously increased, and the corresponding images were acquired until 240 mL of skimmed milk was added and the water turbidity was 98.6 NTU (120 mL, 160 mL, 200 mL, 220 mL, and 240 mL, corresponding to the turbidity of 52.7 NTU, 68.9 NTU, 84.5 NTU, 92.5NTU, and 98.6 NTU, respectively). The first column in Figure 9 shows the intensity of the raw image, and it can be found that as the turbidity of the water increases the intensity of the backscattered light will rise rapidly and generate a large area of grayscale overexposure area. The target object information in the raw image gradually decreases with the increase in turbidity until it is completely lost. This also reflects the necessity of improving the imaging quality of the target objects in high-turbidity underwater scenes. The second and fourth columns in Figure 9 show the processing results using Schechner’s method, Li’s method, and the CLAHE method, respectively, and the fifth column is the processing result of our method. Under the high-turbidity conditions, the signal light in $I_{\min }$ decreases rapidly, and the results of Schechner’s method are low in grayscale value, and the image is dark overall. Li’s method and the CLAHE method have adjusted the image grayscale histogram so that the grayscale is more uniform. However, as the intensity of the backscattered light increases and the noise of the imaging process continues to gain, the intensity of the target signal recovered by Li’s method decreases significantly with the increase in turbidity. CLAHE is a histogram equalization operation for the raw intensity image, and the CLAHE method will gradually fail when the backscattered light intensity exceeds the dynamic range of the CCD. Our method takes into account the high backscattered light intensity and the low SNR characteristics of the imaging process in highly turbid water and uses CLAHE for cross-linear image processing combined with joint noise suppression for polarization enhancement of the target signal. In contrast, as can be seen from the results in the fifth column of Figure 9, although the recovery of the target signal decreases with the increasing turbidity, our processing results outperform the conventional method under the same turbidity conditions.

To quantify the effect of the proposed method on the improvement of the recovered image contrast and noise suppression, the pixel coordinates of the yellow dashed line in Figure 9 are used as the indexing range to compare the differences in grayscale values between the raw image intensity and the results of the different methods of processing, as shown in the last column of Figure 9. The results of the image grayscale indexing under the different turbidity conditions are also shown. It can be seen from the grayscale indexing results that the backscattered light has a significant effect on the intensity of the raw image (purple line) and the CLAHE processing method (orange line). In addition, the lower signal light intensity makes the results recovered by Schechner’s method (green line) have a low contrast. Where the blue line indicates the results of Li’s method and the red line indicates our results, it can be seen that the red line has the best depiction effect on the details of the object under different turbidity conditions. Moreover, it can be verified from the degree of fluctuation of the curves that noise in the imaging results of the highly turbid water has a significant impact on the imaging quality. In order to objectively evaluate the effect of the image quality improvement, we use image entropy and PSNR to evaluate the effect of the different methods of processing, and the corresponding results are shown in

Table 1. Quantitative comparisons of the recovered images in Figure 9.

		Raw	Schechner’s	Li’s	CLAHE	Ours
52.7NTU	PSNR(dB)	8.90	6.49	13.73 a	10.43	17.33
	Entropy	1.69	1.62	2.04	2.05	2.27
68.9NTU	PSNR(dB)	6.61	5.95	11.41	8.87	16.14
	Entropy	1.68	1.24	1.34	1.96	2.21
84.5NTU	PSNR(dB)	5.05	5.64	11.18	5.63	13.05
	Entropy	1.09	1.05	1.18	1.17	2.00
92.5NTU	PSNR(dB)	4.73	5.58	10.66	4.99	11.12
	Entropy	0.48	1.05	1.21	0.66	1.87
98.6NTU	PSNR(dB)	4.63	5.56	10.71	4.68	10.22
	Entropy	0.12	0.99	1.24	0.19	1.75

^a For each evaluation metric, the best results are marked in bold red, and the second best results are marked in blue with underlining.

.

It can be seen in Table 1 that the image entropy value obtained by the proposed method is the highest under the different turbidity conditions; the PSNR is not as good as Li’s method under the condition of 98.6 NTU, but it is also at the second best value and the difference is not significant. The effectiveness of our method is verified by the quantized results; the PSNR can be improved by a maximum of about 2.6 times (84.5 NTU), and the image entropy can be improved by a maximum of about 14.6 times (98.6 NTU) compared to the raw image.

The limitations of image histogram equalization operating on a large number of overexposed regions can be seen from the CLAHE processing results in Figure 9. Our method applies it to the cross-linear image $I_{\min }$ and combines it with joint noise suppression to obtain better results than using CLAHE directly. However, the interaction probability between the photons and particles increases due to the rise in large-size particles in the highly turbid water, causing the photons to gradually lose their initial polarization state and become randomly polarized [26]. This eventually manifests itself as an increase in the intensity of the backscattered light in $I_{\min }$, which is also detrimental to the operation of CLAHE. To illustrate the applicability of the proposed method, peak-to-correlation energy (PCE) was used to evaluate the correlation between the polarization orthogonal sub-images as the turbidity increased [15,16,34]. Figure 10 shows the PCE curves calculated by $I_{\min }$ and $I_{\max }$ under the different turbidity conditions. It can be found that the intensity of the backscattered light in $I_{\min }$ and $I_{\max }$ gradually increases as the turbidity increases, leading to a gradual decrease in the correlation between them. The effect of the histogram equalization of $I_{\min }$ using CLAHE is thus weakened because the target signal light is gradually submerged in the background. In addition, the effect of CLAHE is also related to its setting of the number of tiles (NT) and the clip limit (CL). Therefore, in order to improve the applicability of this method in more complex environments, the optimization of NT and CL [35] and the means of signal enhancement under low-light conditions [36] can be combined to increase the signal proportion in $I_{\min }$ so as to achieve better imaging quality.

4. Discussion

Underwater imaging enhancement methods based on non-physical and physical models have been widely used in dealing with turbid media. However, as the turbidity increases, especially in highly turbid water, the processing effect of the traditional methods becomes limited. Therefore, it is necessary to optimize the imaging system [1,24] or combine different algorithms [5,25,37] to improve the imaging quality. We theoretically and experimentally analyzed the influence of highly turbid water on the underwater active polarization imaging model and showed that the increase in the backscattered light intensity and the high noise gain of the imaging model mainly affect the imaging quality. Inspired by the normalized histogram stretching of $I_{\min }$ by Li’s method, we proposed an active polarization imaging processing method combining non-physical and physical models. Our method uses CLAHE to process a certain number of $I_{\min }$ and then performs polarization enhancement and uses joint filtering noise suppression (multi-frame averaging and bilateral filtering) on the polarization enhancement results to reduce the high noise gain generated by the imaging model and CLAHE. The quantitative results of the experiments verify that our method is suitable for polarization imaging in highly turbid water.

Underwater active polarization imaging takes into account the DOP of target signal light $p_{\mathrm{targ}}$; so, we compared the experimental results of a high-DOP (HDOP) object. Figure 11 shows the experimental results of two metal blades glued to a metal plate. From the results, it is clear that our method is capable of handling the HDOP object. However, the HDOP object and the backscattered light both have a certain degree of polarization-maintaining ability; so, the overexposed area produced by the backscattered light leads to a distortion of the restoration result of the corresponding area, as shown by the red dotted line in Figure 11.

Our method is implemented on the basis of the non-physical model (CLAHE) and the physical model (underwater active polarization imaging) and borrows the processing idea from Li’s method. Therefore, we compare their processing results in Figure 9. To compare our result more fully, we show the results of the classic and the effective dark channel prior (DCP) [7,38] and retinex [39,40] methods in computer vision in Figure 12. Due to the active illumination and with the increasing turbidity, the backscattered light produces inhomogeneity and causes a large number of overexposed areas, which affects the processing effectiveness of the DCP and retinex.

While our method uses only classical processing methods from non-physical and physical models, it is effective to cross them reasonably to form a new processing method to achieve target imaging in harsh underwater environments. In future work, the development from “detectable” to “better detectable” can be achieved by setting up more suitable light sources (e.g., red light [1] or circularly polarized light illumination [24], etc.) and combining more advanced computational methods (e.g., computer vision [7,38,39,40] or machine learning [41,42], etc.).

5. Conclusions

In this paper, it is recognized through theoretical and experimental analysis that the influence of backscattered light on underwater active polarization imaging becomes more and more obvious as the turbidity increases, especially in the case of highly turbid water (>50 NTU). It is mainly reflected in two aspects: one is the increase in backscattered light intensity, and the other is the decrease in $\Delta p$. Together, they produce a processing result with a significantly lower contrast and high noise gain. We proposed an active polarization imaging method based on the CLAHE of cross-linear images and joint noise suppression (multi-frame averaging and bilateral filtering). The experimental results in the different high-turbidity water samples ($52.7\sim 98.6$ NTU) show that our recovery is better than the conventional methods. The PSNR and entropy of our method are improvements of up to 2.6 times and 14.6 times compared to the raw image. This method, which appropriately combines non-physical and physical methods [25], facilitates the application of underwater active polarization imaging in harsher environments without changing the complexity of the system.