Modified Newton-residual interpolation for division of focal plane polarization image demosaicking

Xinlong Liu; Linjie Yang; Luping Wang

doi:10.1364/OE.460495

1. Introduction

1.1 Background

Light is an important intermediary of information transfer, and it has physical properties such as wavelength, amplitude, and polarization. The human visual system can process amplitude and wavelength information, but cannot detect any polarization characteristics. However, polarized light expresses richer information in natural scenes by signifying other characteristics such as surface roughness, medium properties, or biochemistry of the measured object [1]. In some vision tasks, it is often difficult to find image features through information with light intensity and wavelength, such as metal material recognition [2], low-light road recognition [3], three-dimensional reconstruction [4], biomedical imaging [5], target tracking [6,7], and astronomical detection [8,9]. Using polarization properties, we can handle these tasks more efficiently in polarized images. Therefore, people are paying more and more attention to the acquisition of polarization information and obtain polarization images with excellent imaging quality.

Currently, there are many studies of polarization imaging methods. The typical acquisition methods of polarization images are: (1) division of time (DoT); (2) division of amplitude (DoAM); (3) division of aperture (DoAP) and (4) division of focal plane polarimeters (DoFP) [10,11]. The DoFP method simultaneously acquires the polarization of incident light at different angles by combining pixelated microstructures with image element arrays, without considering the dynamic state of the or needing complex auxiliary equipment [12]. Typically, a nanofabrication-grade polarization filter array (PFA) is assembled on a conventional sensor to form a polarization sensor [13]. The common PFA is composed of four angles of 0, 45, 90, and 135 into a 2×2 periodic structure [14,15]. When the PFA is combined with the sensor directly to form a monochrome polarization filter array (MPFA) and with the sensor that has color filter array (CFA) to form a color polarization filter array (CPFA), the structure is shown in Fig. 1. But the pixels of the DoFP structure only obtain the polarization information for one angle using micro-polarizer array (MPA), and the missing pixel components need to be restored by interpolation of neighboring pixels [16].This interpolation process is called polarization demosaicking (PDM). Bayer CFA demosaicking has been developed and studied for many years [17]. However, due to the different imaging modes in the color domain and the polarization domain, many excellent color algorithms cannot be transplanted to the polarization domain. Because in the color domain, the green channel accounts for half of all color channels, so people reconstruct the image of the other color channels after restoring the full green channel image [18]. For polarized arrays, the 4 angles have the same pixel ratio, which makes reconstruction more difficult. The goal of PDM is to recover the full-resolution images of the four polarization channels, calculate the degree of linear polarization (DoLP) and angle of line polarization (AoLP) parameters more accurately, and avoid the effects of errors.

Fig. 1. The schematic diagram of DoFP structure with MPFA and CPFA.

Download Full Size | PDF

The current PDM algorithms are roughly divided into three categories [19], including: interpolation algorithm, sparse representation and deep learning. Interpolation algorithms are the most classical and common algorithms. Methods such as quadratic interpolation, cubic interpolation, and weight interpolation are algorithms that use the relationship between pixels in the spaced field for direct interpolation. The gradient-based method proposed by Gao [20] extends the direct interpolation method, using strategies such as bicubic convolution to consider the gradient relationship between pixels. Li [21] proposed a Newton interpolation method from mathematics and image edge detection method to improve the edge preservation effect. Wu [22] proposed an interpolation method for edge preservation, through the detection of the edge for different intensity interpolation in the orthogonal direction. This type of interpolation gives good results in images with small pixel gradients or neat edges. However, in high-frequency regions such as edges or sharp gradient changes, some direct interpolation methods over-smooth edge pixels to the point of causing artifacts, which seriously affects image quality and is difficult to suppress. Therefore, distinguishing edges during interpolation is increasingly important, and the interpolation algorithm still has great potential.

The sparse representation algorithm is a common image reconstruction algorithm, which is used in many image processing tasks [23]. Zhang [24] proposed a PDM algorithm based on sparse representation. This algorithm constructs a dictionary to represent the numerical relationship between pixel blocks, and then restores the final demosaic image from the dictionary. Although the dictionary ensures the correlation of different polarization channels, sparse expression takes a lot of time, which leads to poor real-time performance.

Along with the development of deep learning, some PDM algorithms based on deep learning have started to be proposed [25,26]. Zeng et al. [27] proposed an end-to-end convolutional neural network (CNN) model called Fork-NET, which does not follow the traditional interpolation process to build the network structure during the learning process, but learns the polarization properties in the mosaic images directly, thus outputting DoLP and AoLP directly. Their method achieved good performance on the test set. However, deep learning methods do not extensively test performance in all scenes, the optimum test set is roughly the same as the scenes captured in the training set and the algorithm results are not tested with realistic unprocessed images. In addition, in some scenarios, it is critical to not use dedicated devices such as a GPU to compute the processing. So, there are still limitations to using deep learning methods for PDM.

Compared with other interpolation algorithms, residual interpolation (RI) is a newer interpolation algorithm and has shown good performance [28–30]. The residual is the difference between the estimated pixel value and the known pixel value. The current RI has some key points that notably affect the result. First, obtaining an ideal guide image is challenging, and an excellent guide image improves the edge retention effect during demosaicking; secondly, ensuring the continuity of high-frequency information during the guide process is arduous, which makes the edge degenerate easily. These problems make the edges of the resulting image blurry, prone to halo and zipper effects, which shows that the edge of the image presents a zipper-like burr due to the wrong interpolation calculation or the wrong interpolation direction. In this paper, we propose a modified Newton-residual interpolation (MNRI) method with better edge performance based on RI. In the interpolation process, we continue to pay attention to the judgment and processing of the edge during interpolation. First, we use the modified Newton interpolation method for up-sampling to generate the guide image by bidirectional and weighted Newton interpolation. Next, when the residual image is used for the guiding process, we increase the flexibility of the guiding coefficients with different weights by judging the local features of the image edges. This guiding method can better deal with the recovery at the edge of the image. Finally, we employ MNRI to get missing pixel values of 4 polarization angles and extensive experimental comparisons with existing algorithms. The structure of MNRI is shown in Fig. 2.

Fig. 2. The structure of MNRI demosaicking algorithm.

Download Full Size | PDF

1.2 Linear polarization computation

The Stokes parameter is a tool for mathematically expressing the polarization state of light. Since the proportion of circularly polarized components is extremely small in nature, and the existing focal plane polarization sensors can only measure the linear polarization state, so the Stokes vector is calculated as:

(1)$$\overrightarrow S = \left[ {\begin{array}{{c}} {{S_0}}\\ {{S_1}}\\ {{S_2}} \end{array}} \right] = \left[ {\begin{array}{{c}} {\frac{1}{2}({I_0} + {I_{45}} + {I_{90}} + {I_{135}})}\\ {{I_0} - {I_{90}}}\\ {{I_{45}} - {I_{135}}} \end{array}} \right]$$

where ${I_0}$, ${I_{45}}$, ${I_{90}}$ and ${I_{135}}$ are the intensity images obtained at the four polarization angles of 0°, 45°, 90° and 135°. Among them, ${S_0}$ is the total intensity of incident light, while ${S_1}$ and ${S_2}$ describe the state of linear polarization. The two polarization characteristics that are most interesting for observing polarization are DoLP and AoLP. The definition for the DoLP and the AoLP are:

(2)$$\begin{array}{l} \textrm{DoLP} = \sqrt {S_1^2 + S_2^2} /{S_0}\\ \textrm{AoLP} = \frac{1}{2}\arctan ({S_2}/{S_1}) \end{array}$$

To obtain high-quality DoLP and AoLP images, accurate images in 4 directions are required. Thus, the goal of the PDM algorithm is to restore the real image as much as possible [31].

The rest of the article is organized as follows. The second Section introduces RI and its improved algorithm. Section 3 introduces the details of the proposed MNRI algorithm. In Section 4, we carry out experiments and comparisons with other methods. The final conclusion and future research directions are expressed in Section 5.

2. Related work

The RI framework first appeared to calculate the RGB Bayer array demosaicking algorithm [32], which mainly includes the following five steps (here we assume that the array arrangement is GRBG, and solve the R band):

1. obtain the guiding image G through linear interpolation;
2. through guided up-sampling of the observed R value, a tentative estimate of the R-band (represented as $\mathrm{\hat{R}}$) is generated from the interpolated G-band;
3. the gradient-based residual (represented as R-$\mathrm{\hat{R}}$) is calculated by using the difference between the tentative estimate and the observed R value;
4. calculate the residual image interpolation;
5. add the tentative estimate and the interpolation residual to obtain the R band of the interpolation (represented as $\mathrm{\tilde{R}}$).

Figure 3 illustrates the interpolation process of the R band in the horizontal direction. We use the same residual interpolation method in the vertical direction of R band to get a complete R image.

Fig. 3. Schematic diagram of RI interpolation in the horizontal direction of the R band.

Download Full Size | PDF

RI uses a guide filter (GF) to perform an up-sampling process to reconstruct the resulting image from the sparse pixel image. GF is a widely used algorithm in computer vision and image processing [33]. The original RI uses GF when generating the preliminary estimation map, and uses the G map generated by interpolation as the guide map. The preliminary estimate of the window ${\omega _r}$ at pixel $(p,q)$ is expressed as:

(3)$${\hat{R}_{i,j}} = {a_{p,q}}{G_{i,j}} + {b_{p,q}},{\forall _{i,j}} \in {\omega _r}$$

where ${\hat{R}_{i,j}}$ represents the tentative estimate of the pixel $(i,j)$ in the window ${\omega _r}$, and (${a_{p,q}}$, ${b_{p,q}}$) is a pair of guide coefficients assumed to be constant. The minimized cost function for linear coefficient (${a_{p,q}}$, ${b_{p,q}}$) is expressed as:

(4)$$\begin{aligned} E({a_{p,q}},{b_{p,q}}) & =\sum {{{({M_{i,j}}({R_{i,j}} - {{\hat{R}}_{i,j}}))}^2}} \\ & =\sum {{{({M_{i,j}}({R_{i,j}} - {a_{p,q}}{G_{i,j}} - {b_{p,q}}))}^2}} \end{aligned}$$

where the binary image ${M_{i,j}}$ is the binary mask at pixel $(i,j)$. The value of the observed R pixel is 1 and the value of other pixels is 0. There is an adjustment parameter in the original GF, which is used to adjust the smoothness of the result. Kiku et al. completely eliminated this parameter but added a tiny coefficient to avoid dividing by 0. Therefore, there are no adjustment parameters in the above formula. The above cost function derivation shows that the original RI minimizes the residuals (i.e., ${R_{i,j}} - {\hat{R}_{i,j}}$) themselves. In the final calculation of the linear parameter term, since the window slides evenly across the image, the overlapping position of the window is calculated with the weighted result. Thus, the weight calculation method based on residuals is:

(5)$${W_{p,q}} = 1/\frac{1}{{|{{\omega_{p,q}}} |}}\sum\limits_{i,j \in {\omega _{p,q}}} {{{({{R_{i,j}} - {a_{p,q}}G_{i,j}^M - {b_{p,q}}} )}^2}}$$

where $|{{\omega_{p,q}}} |$ refers to the number of ${R_{i,j}}$ pixels subsampled in the window. Then the result of linear parameters ${\bar{a}_{i,j}}$ and ${\bar{b}_{i,j}}$ is improved and the calculation method is:

(6)$${\bar{a}_{i,j}} = \frac{{\sum\limits_{p,q \in {\omega _{i,j}}} {{W_{p,q}}} {a_{p,q}}}}{{\sum\limits_{p,q \in {\omega _{i,j}}} {{W_{p,q}}} }},{\bar{b}_{i,j}} = \frac{{\sum\limits_{p,q \in {\omega _{i,j}}} {{W_{p,q}}} {b_{p,q}}}}{{\sum\limits_{p,q \in {\omega _{i,j}}} {{W_{p,q}}} }}.$$

Kiku et al. proposed an improved calculation method for the linear coefficient pair (${a_{p,q}}$, ${b_{p,q}}$): Minimized-Laplacian RI (MLRI) [28]. In the guided filtering of MLRI, the minimized Laplacian energy residual is introduced into the GF for each local window. The cost function of ${a_{p,q}}$ is:

(7)$$\begin{aligned} E({a_{p,q}}) & =\sum {{{({M_{i,j}}{\nabla ^2}({R_{i,j}} - {{\hat{R}}_{i,j}}))}^2}} \\ & =\sum {{{({M_{i,j}}{\nabla ^2}({R_{i,j}} - {a_{p,q}}{G_{i,j}} - {b_{p,q}}))}^2}} \\ & =\sum {{{({\nabla ^2}(R_{i,j}^M) - {a_{p,q}}{\nabla ^2}(G_{i,j}^M))}^2}} \end{aligned}$$

Kiku et al. proves that RI with Laplace energy has a better probability of estimating the residual image. Since the observed pixels are down-sampled, a sparse Laplacian filter is needed to calculate the Laplacian energy for convolution to calculate the approximate Laplacian map of the obscured G image and the observed R image. The sparse Laplacian filter is expressed as:

(8)$$\left[ {\begin{array}{{ccccc}} 0&0&{ - 1}&0&0\\ 0&0&0&0&0\\ { - 1}&0&4&0&{ - 1}\\ 0&0&0&0&0\\ 0&0&{ - 1}&0&0 \end{array}} \right]$$

Among them, the bias coefficient ${b_{p,q}}$ does not affect the Laplace energy, it can be calculated by formula (3) when ${a_{p,q}}$ is known:

(9)$$E({{b_{p,q}}} )= \sum\limits_{i,j \in {\omega _{p,q}}} {{{({{R_{i,j}} - {a_{p,q}}G_{i,j}^M - {b_{p,q}}} )}^2}}$$

After MLRI [28], Manno et al. proposed an iterative RI [34]. By using the R band result as a guide image, the residual calculation process was iterated again to obtain a new R band, and the iteration was stopped when the difference between the obtained estimated value and the observed value met an appropriate threshold. In this way, the results can be continuously corrected to improve the signal-to-noise ratio (SNR). Manno et al. proposed the ARI algorithm, which can flexibly combine different RI algorithms and iteration times to optimize the demosaicking results.

In the latest polarization edge-aware RI (EARI) algorithm [29], Morimatsu et al. apply RI in the field of PDM. EARI divides the interpolated direction polarization pixels into four direction blocks, and then calculates the weight of each direction block and combines it into a half intensity image as a guide image, which is called edge-aware. Then the captured intensity guide map and the sub-sampled image are subjected to the residual interpolation process to obtain the reconstructed interpolation image. The guide image generation method perceives the position of the edge and improves the accuracy of the edge region interpolation. The subsequent residual process follows the MLRI process, and Laplace energy is introduced in the guide process. This method has achieved State-of-the-art (SOTA) results. The EARI algorithm structure is shown in Fig. 4. Although EARI calculates the edge distribution of the pixels to be interpolated in 4 directions during interpolation, it does not consider the distinction between flat and edge regions in the image. MNRI follows the processing of different regions of the image during the up-sampling to generate the guide image and the guide process. Thus, MNRI achieves better results than EARI in PDM.

Fig. 4. The structure of ERI demosaicking algorithm.

Download Full Size | PDF

3. Proposed method

3.1 Modified Newton interpolation

According to the RI framework, we use a modified second-order Newton interpolation algorithm to generate effective guide images, which has better edge retention than EARI.

In the field of polarization, using the pixel value of intensity we use ${{S0} / 2}$ as the guide image just like EARI, there are two calculation methods of $S0$ as follows [22,29]:

(10)$$S{0_1} = {I_0} + {I_{90}} = {I_{45}} + {I_{135}}$$

(11)$$S{0_2} = \frac{1}{2}({I_0} + {I_{90}} + {I_{45}} + {I_{135}})$$

Through experiments, we found that the calculation method of S0 affects the final interpolation result. We show the difference in results in the next section.

In order to obtain an image with sufficient accuracy, it is necessary to obtain the most realistic possible 4-angle polarization images. In the interpolation process, for the non-edge area to be interpolated in the image, we select the two-direction weighting calculation of the interpolation point to obtain the target pixel value. For edge pixels, the target pixel value is obtained by calculating weight distribution in multiple directions. This region differentiation processing method has low computational complexity and also handles the boundary and flat interpolation distinction. The traditional Newton interpolation uses the difference quotient to calculate the interpolation point value. So the different interpolation points are linearly recursive, which leads to the non-smooth interpolation curve at the node [35]. We improve it according to local features in two-dimensional images in order to improve the accuracy of interpolation.

It is known that the function value of the function $f(x)$ at the equidistant node ${x_i}$ is expressed as

(12)$$f({x_i}) = {f_i}(i = 0,1,\ldots ,n)$$

Then the n-order Newton interpolation polynomial is expressed as:

(13)$$\begin{aligned} {N_n}(x) & =f({{x_0}} )+ f[{{x_0},{x_1}} ]({x - {x_0}} )+ \ldots \\& + f[{{x_0},{x_1}, \ldots ,{x_n}} ]({x - {x_0}} )({x - {x_1}} )({x - {x_{n - 1}}} )\end{aligned}$$

The interpolation method is illustrated in Fig. 5. Assuming that there are 4 known adjacent points $P1,P2,P3$ and $P4$ in the horizontal direction, the interval between each point is 1, and I is the position of the point to be calculated. The $\Delta x$ is the position between I and $P2$, with a value range of 0≤$\Delta x$<1, and $t = 1 + \Delta x$.

Fig. 5. Schematic diagram of the horizontal interpolation of the target pixel.

Download Full Size | PDF

We consider that the pixel is correlated with three neighboring pixels, then calculate the pixel to be interpolated by the second-order Newton interpolation formula. We choose $P1$, $P2$ and $P3$ three original pixels as known data, then the point I calculation method is:

(14)$${F_1}(I) =P1 + \Delta P1 + \frac{{{\Delta ^2}P1}}{{2!}}(t - 1)$$

When three original pixels of $P2$, $P3$, and $P4$ are known data, the calculation method is:

(15)$${F_2}(I) =P2 + \Delta P2(t - 1) + \frac{{{\Delta ^2}P2}}{{2!}}(t - 1)(t - 2)$$

where $\Delta P1$ and $\Delta P2$ are the first-order differences of points $P1$ and $P2$, as well as ${\Delta ^2}P1$ and ${\Delta ^2}P2$ are the second-order differences of $P1$ and $P2$ separately. Bringing $t = 1 + \Delta x$ into Eq. (14) and Eq. (15), this can be simplified to get:

(16)$${F_1}(I) = \frac{1}{2}(\Delta {x^2} - \Delta x)P1 + ( - \Delta {x^2} + 1)P2 + \frac{1}{2}(\Delta {x^2} + \Delta x)P3$$

(17)$${F_2}(I) = \frac{1}{2}(\Delta {x^2} - \Delta x)P4 + \frac{1}{2}(\Delta {x^2} - \frac{3}{2}\Delta {x^2} + 1)P2 + ( - \Delta {x^2} + 2\Delta x)P3$$

When point I is in a flat area, the results of the two calculation methods are almost the same. When the point I is in the edge area, the calculation results of the two methods will be different, and the pixel distribution around the pixel position needs to be considered. It can be seen from the above formula that the weights of the far points $P1$ and $P4$ in the left and right directions of point I have the same effect on the result. According to interpolation theory, if point I is close to $P2$, the calculated value is also close to $P2$, and vice versa. Therefore, in order to take into account the pixel information in both directions, the weights of the pixels on both sides are considered for interpolation. Then the pixel value of point I is:

(18)$$f(I) = L{F_1}(I) + R{F_2}(I)$$

In Eq. (18), L and R respectively represent the correlation weighting coefficients in the left and right directions. The larger the coefficient, the greater the regional consistency at the pixel. We consider that:

(19)$${F_L} = |(P1 - P2) - (P3 - P4)|$$

(20)$${F_R} = |(P4 - P3) - (P2 - P1)|$$

Then the coefficients L and R are calculated as:

(21)$$L = \frac{{{F_L}}}{{{F_L} + {F_R}}},R = \frac{{{F_R}}}{{{F_L} + {F_R}}}$$

So far, the image interpolation in the horizontal direction has been completed, and then the interpolation in the other directions is performed to obtain the final guide image. Perform vertical interpolation in the image shown in Fig. 6 to obtain an interpolation result of 0° at point $(i,j)$. In the 3×3 area, we select the vertical direction of the point, and the 45° direction and the 135° direction to calculate the second-order difference results, which are expressed as: ${\Delta ^2}{f_1}$, ${\Delta ^2}{f_2}$, and ${\Delta ^2}{f_3}$.

Fig. 6. Schematic diagram of vertical interpolation of target pixel. The green line indicates the 45° direction, the red line indicates the vertical direction, and the blue line indicates the 135° direction.

Download Full Size | PDF

We use the inaccurate Newton interpolation to roughly interpolate the image. This method has higher accuracy in the flat areas of the image. Near the edge of images, only by introducing more directional values can we get closer to the edge. We have obtained the edge-preserved guide image through secondary processing, so as to make sufficient preparations for obtaining polarization interpolation pixels in the RI frame later. In this way, it is judged whether the second-order difference between the vertical direction, the 45° direction and the 135° direction are less than T. If so, it means that the image gradient changes less in that direction. When the condition is met for the first time, then the target pixel value is calculated by second-order Newton interpolation. When the three difference results are greater than the threshold T, it is considered that the point is at the boundary position. As the second-order difference value is smaller, the position of the point is more related to its direction, the results of the two calculations so the normalization weight calculations are as follows:

(22)$${e_1} = 1 - \frac{{|{{\Delta ^2}{f_1}} |}}{{|{{\Delta ^2}{f_1}} |+ |{{\Delta ^2}{f_2}} |+ |{{\Delta ^2}{f_3}} |}}$$

(23)$${e_2} = 1 - \frac{{|{{\Delta ^2}{f_2}} |}}{{|{{\Delta ^2}{f_1}} |+ |{{\Delta ^2}{f_2}} |+ |{{\Delta ^2}{f_3}} |}}$$

(24)$${e_3} = 1 - \frac{{|{{\Delta ^2}{f_3}} |}}{{|{{\Delta ^2}{f_1}} |+ |{{\Delta ^2}{f_2}} |+ |{{\Delta ^2}{f_3}} |}}$$

After the three weight coefficients are obtained, the calculation method of the weight pixel value of the predicted point is:

(25)$$\begin{aligned} {\mathop{I}\limits^\frown_0}(i,j) & =\frac{1}{2}[\frac{1}{2}{e_1}({I_0}(i,j - 1) + {I_0}(i,j + 1)) + \\ &\frac{1}{2}{e_2}({I_0}(i + 1,j - 1) + {I_0}(i - 1,j + 1)) + \\& \frac{1}{2}{e_3}({I_0}(i - 1,j - 1) + {I_0}(i\textrm{ + }1,j\textrm{ + }1))] \end{aligned}$$

Since the sum of the three directional weighting factors is 2, the final target pixel value is obtained by multiplying by 1/2. The threshold value T is used as the control parameter considering the edge position judgment. Through experiments, we set T = 25 to obtain better results. With the above steps, we obtain a guide image with good edge retention for subsequent steps.

3.2 Improved GF

In RI, the regular penalty term $\lambda$ of the original GF is set to a very small smoothing parameter to avoid dividing by 0. We believe that this parameter can achieve the effect of adjusting edge retention through flexible settings [36]. On this basis, we added a parameter to ensure that the problem of edge judgment is solved, making the GF more robust in determining whether the region is flat or an edge.

In the guide image filtering, the output image ${Q_i}$ is linearly transformed from the k pixels of the guide image ${I_i}$ in the window $\omega$:

(26)$${Q_i} = {a_k}{I_i} + {b_k}\quad\textrm{ }\forall i \in {\omega _k}$$

where $i$ and $k$ are pixel index, and ${a_k}$ and ${b_k}$ are constants of the window $\omega$. Their values are obtained by minimizing the cost function $E({a_k},{b_k})$:

(27)$$E({a_k},{b_k}) = \min \sum\limits_{i \in {\omega _k}} {(({a_k}{I_i} + {b_k} - {p_i}) + \lambda a_k^2)}$$

Among them, $\lambda$ is the regular penalty term for ${a_k}$, and ${p_i}$ is the filter input. The parameter $\lambda$ determines the flatness, and it determines the standard of “flat patch” or “high variance”. More specifically, blocks with variance ($a_k^2$) much smaller than λ are smoothed, while blocks with variance much larger than λ are kept [37].

We have added a parameter to indicate the flatness coefficient of the partial image. The mean value of the local variables of all the pixel values in the image is set to:

(28)$$\Gamma = {\bar{\sigma }^2} = \frac{1}{N}\sum\limits_{k = 1}^N {\sigma _k^2}$$

Among them, $\sigma _k^2$ is the local variance of the image I in the window ${\omega _k}$. Then the new cost function becomes:

(29)$$E({a_k},{b_k}) = \min \sum\limits_{i \in {\omega _k}} {(({a_k}{I_i} + {b_k} - {p_i}) + \lambda \Gamma a_k^2)}$$

Then the calculation method of the parameter ${a_k}$ and ${b_k}$ is:

(30)$${a_k} = \frac{{\frac{1}{{|\omega |}}\sum\limits_{i \in {\omega _k}} {{I_i}} {p_i} - {\mu _k}{{\bar{p}}_k}}}{{\sigma _k^2 + \lambda \Gamma }}$$

(31)$${b_k} = {\overline p _k} - {a_k}{\mu _k}$$

Among them, ${\mu _k}$ and $\sigma _k^2$ are the mean and variance of the I image in window ${\omega _k}$, $|\omega |$ is the number of pixels in ${\omega _k}$, and ${\bar{p}_k}$ is the mean of filter input ${p_i}$ in ${\omega _k}$.

In the application of GF to local contrast enhancement, I and ${p_i}$ are the same, and the following equation holds:

(32)$$\begin{aligned} \frac{1}{{|\omega |}}\sum\limits_{i \in {\omega _k}} {{I_i}} {p_i} - {\mu _k}{{\bar{p}}_k} & =\sigma _k^2\\ {\mu _k} & ={{\bar{p}}_k} \end{aligned}$$

Substituting the above formula Eq. (32) into Eq. (30) and Eq. (31), Then ${a_k}$ and ${b_k}$ can be calculated:

(33)$${a_k} = \frac{{\sigma _k^2}}{{\sigma _k^2 + \lambda \Gamma }}$$

(34)$${b_k} = ({1 - {a_k}} ){\mu _k}$$

From the above formula, ${a_k}$ and ${b_k}$ are obtained. If it is in a high difference area, the local mean should be greater than the mean, ${a_k}$ and ${b_k}$ are close to 1 and 0; if in a flat area, the local variance should be less than the average variance, while ${a_k}$ and ${b_k}$ are close to 0 and 1, and the average of the pixel values is near the average pixel. The improved GF can maintain a larger ${a_k}$ at the edge and near the edge, so that the resulting image has a better protective effect on the edge. The result with more complete and prominent edges is more in line with the real situation, according to the principle of guided filtering [37]. Figure 7 shows the processing results of the improved GF compared with the original GF. As shown in Fig. 7, the improved GF is sharper than the original GF at the transition edge, with no artifacts appearing.

Fig. 7. Comparison results of (b) original GF and (c) improved GF for (a) input image.

Download Full Size | PDF

Bringing this guided filtering parameter result into the RI framework, we obtain the new calculation of cost function $E({a_{p,q}})$:

(35)$$\begin{aligned} E({a_{p,q}}) & =\sum {{{({M_{i.j}}{\nabla ^2}({R_{i,j}} - {{\hat{R}}_{i,j}}))}^2}} \\ & =\sum {{{({\nabla ^2}(R_{i,j}^M) - \frac{{\sigma _{i,j}^2}}{{\sigma _{i,j}^2 + \lambda \Gamma }}{\nabla ^2}(G_{i,j}^M))}^2}} \end{aligned}$$

So, the line guide coefficient of residual image, ${\bar{a}_{i,j}}$ and ${\bar{b}_{i,j}}$, are obtained as:

(36)$${\bar{a}_{i,j}} = \frac{{\sum\limits_{p,q \in {\omega _{i,j}}} {{W_{p,q}}} \sigma _k^2}}{{\sum\limits_{p,q \in {\omega _{i,j}}} {{W_{p,q}}(\sigma _k^2 + \lambda \Gamma )} }},{\bar{b}_{i,j}} = \frac{{\sum\limits_{p,q \in {\omega _{i,j}}} {{W_{p,q}}} {b_{p,q}}}}{{\sum\limits_{p,q \in {\omega _{i,j}}} {{W_{p,q}}} }}$$

The size of the partial window determines the judgment range of the edge field of view, and the traditional window size is 3. However, the 3×3 field of view is very small, and sliding around the edge may confuse the edge with the flat area. We expanded the window size to have a better perception field of vision and obtain better results. The experimental results can be seen in Section 4.

3.3 CPFA demosaicking

For the CPFA demosaicking method, we implement it by extending MNIR in MPFA. The overall demosaicking framework of our proposed method is the same as EARI, by decomposing the known color images into monochromatic polarized images, which are then processed. First, we transform the original CPFA data into a 4-directional Bayer mode color image by down-sampling, and then perform a color demosaicking algorithm to obtain a 4-directional filled image. Next, the obtained 12-channel image is up-sampled to obtain a 3-color fully polarized image. Finally, the CPFA demosaicking process is completed by using the MPFA demosaicking method three times on the fully polarized image. Figure 8 illustrates the general flow of the CPFA framework.

Fig. 8. The outline of MNRI for CPFA.

Download Full Size | PDF

4. Experimental result

To quantitatively evaluate the performance of our algorithm, we select two public datasets for experiment. According to the review [38], we use the proposed MNRI algorithm to compare with 6 existing algorithms, including bilinear interpolation (Bilinear), Bicubic interpolation (Bicubic), intensity correlation among polarization channels (ICPC) [39], pseudo-panchromatic image difference (PPID) method [40], Newton's polynomial interpolation (NP) [21], and edge-aware residual interpolation (EARI) [29]. Bilinear and Bicubic are processed using MATLAB code. NP, ICPC and EARI come from the open source code [41–43]. The code of PPID is provided by the author. We use 4 sets of complete polarization images (0°, 45°, 90°, 135°) to verify our algorithm, and calculate the corresponding intensity and polarization images through the comparison algorithm. The results of the algorithm and calculation results are shown in Fig. 9. Finally, all algorithms of the GT images are used for evaluation.

Fig. 9. The 4 angles polarization images and display method of AoP-DoP.

Download Full Size | PDF

In order to accurately compare MPFA's demosaicking algorithm, we use green channel image as the GT image. We selected two objective indicators to evaluate image quality. Peak Signal-to-Noise Ratio (PSNR) is a classic indicator of the demosaicking algorithm, which measures the difference between two images by calculating the ratio of the GT image to the resulting image. A larger PSNR means better image recovery performance. For AoLP images, we use angle error as an indicator, and lower it means better performance.

(37)$$\textrm {PSNR} = 10{\log _{10}}\left( {\frac{{{{({\max {I_{GT}}} )}^2}}}{{\textrm{MSE} ({I_{GT}},{I_{IP}})}}} \right)$$

(38)$$\textrm{Angle error} = \sqrt {\frac{1}{{MN}}\sum\limits_{i = 1}^M {\sum\limits_{j = 1}^N {{{({{I_{GT\textrm{ }}}(i,j) - {I_{IP\textrm{ }}}(i,j)} )}^2}} } }$$

4.1 Experiment 1

The first monochrome and color data from Japan [29], which contains 40 images with 12 channels. The data of 12 channels consists of 4 polarization angles and RGB images. The experimental equipment is a JAI CV-M9GE 3-CCD camera, and SIGMAKOKI SPF-50C-32 linear polarizer attached to PH-50-ARS rotating polarizer mount. The scene is under non-polarized light during acquisition. For each polarization angle, the author collected 1000 images and averaged them to reduce noise as GT image [44]. All image sizes are 1024×768 with 10bit-encode. The dataset [45] contains the scenes as shown in Fig. 10.

Fig. 10. The collection of 40 images included in the dataset.

Download Full Size | PDF

4.1.1 MPFA demosaicking

In the MPFA experiment, the size of the GF window slightly affects the final result, and we test the impact of all parameters. At the same time, we choose a different S0 so that there is a variance in the results.

First, we experimentally demonstrate the effect of the size of the GF window on the results, which are shown in Table 1. The results in the table illustrate that an increase in the window within a certain range will improve the metrics, due to the fact that the perceptual field of view is expanded to better determine whether the window is a smooth or edge region. Based on the experimental results, we chose the window radius of 15 and used it as the default radius in our subsequent methods.

Table 1. Comparison results of GF windows with different radius (Average of 40 scenes)

View Table | View all tables in this article

According to the results in Table 1, it can be seen that in the single polarization angle image, the MNRI algorithm has an exceptional improvement. For polarized images, our results also have strong performance. From the experimental results, as the filter window increases, the indicators gradually increase, which is consistent with our previous analysis. For the two methods of calculating S0, the MNRI results are stronger than other algorithms in some indicators. In summary, it can be seen that the MNRI algorithm obtains better demosaicking results and reduces polarization image calculation errors by focusing on and calculating edges.

Next, we conduct a comparison experiment with other methods for objective indicators. In order to illustrate the effectiveness of our method, we performed the ablation experiments. The results of the comparison of different methods are shown in Table 2. The results show the better achieved by our algorithm in single-angle polarization image recovery as well as partial Stokes parameters. Although we did not get the first place in the PSNR metrics of S2 and DoLP, we were behind by a small margin.

Table 2. Comparison results of different methods for MPFA demosaicking (Average of 40 scenes)

View Table | View all tables in this article

To visualize the robustness of the proposed method, we show the previous SOTA EARI algorithm results to zoom in and compare, as shown in Fig. 11. The text area of EARI has a strong zipper effect, different colors are mixed up seriously, and the number “9” can no longer be distinguished. Our algorithm has a slight interpolation confusion problem, the zipper effect is lighter, and the “9” can be clearly recognized. At the boundary between the window and the door below, EARI has different degrees of edge aliasing. In comparison, the boundary of the S0₁ is more continuous than that of S0₂, and there is no point-like noise. The result of MNRI with R = 15 is that the extension of the boundary is more stable, and there is a small fault in the boundary of MNRI with R = 3.

Fig. 11. Subjective comparison of the intensity image and the AoP-DoP visualization for CPFA demosaicking.^aThe result is MNRI(S0₁) with GF window radius R = 3. ^bThe result is MNRI(S0₁) with GF window radius R = 15.

Download Full Size | PDF

4.1.2 CPFA demosaicking

For CPFA, we used 12 channels of fully polarized color images for the experiment, and the results are shown in Table 3.

Table 3. Comparison results of different methods for CPFA demosaicking (Average of 40 scenes)

View Table | View all tables in this article

According to the best results (in bold), the MNRI algorithm has achieved an advantage. Compared with the significant improvement of MPFA, the addition of color makes it difficult to guarantee the accuracy of the interpolation method. The separation of color and polarization information makes the relationship between pixels sparser, and it takes more effort to take into account the interpolation accuracy of the two. We still choose the enlarged comparison chart to compare the actual results. According to Fig. 12, it can be seen that there is minor difference between the MNRI and EARI algorithms. The MNRI algorithm reduces the error relative to EARI in the eyebrow area, and has slightly better results, but compared with the GT image, there is still a jagged defect.

Fig. 12. Subjective comparison of the intensity image and the AoP-DoP visualization for MPFA demosaicking.

Download Full Size | PDF

4.2 Experiment 2

We also tested our technique using data from Qiu [46]. This test set consists of 40 color polarized images, 30 in unpolarized light and 10 in polarized light. These images are averaged over 100 images as GT, and have a resolution of 1024×1024. We test all 40 images and use the G-channel image as the grayscale test image as in the previous experiment. The dataset contains the scene as shown in Fig. 13. The complete dataset is available in Ref. [47].

Fig. 13. The collection of 40 images included in the dataset.

Download Full Size | PDF

4.2.1 MPFA demosaicking

Table 4 shows the comparison of our proposed method with other methods. To ensure the robustness of the algorithm, we mosaicked 40 images together without distinguishing the effect of the lighting environment. This experimental result still demonstrates the advantage of our algorithm in single-angle polarization image recovery and AoLP angle error. In this test, the MNRI using the S0₂ calculation method obtained a better performance compared to S0₁, indicating that the calculation method is more effective. In the actual image performance, we still compare with EARI, as shown in Fig. 14.

Fig. 14. Subjective comparison of the intensity image and the AoP-DoP visualization for CPFA demosaicking.

Download Full Size | PDF

Table 4. Comparison results of different methods for MPFA demosaicking (Average of 40 scenes)

View Table | View all tables in this article

In the magnified image shown in the red box, the branches and the blackboard border in the upper right corner are the places where the difference is obvious. In our results, there are fewer burrs on the edges of the branches, and the interior of the frame is cleaner. In the blue box, the white stripes at the bottom left are separated in GT, and our results have less aliasing than EARI. However, in high-frequency regions such as dense fringes, edge protrusions appear in all algorithms, which are difficult to control. There is still much room for improvement in our method.

4.2.2 MPFA demosaicking

As before, we test 40 12-channel images and the results are shown in Table 5. In CPFA experiments, our metric wins by a small margin. In the actual image performance, the gap between the algorithms is not obvious. As shown in Fig. 15, our algorithm makes fewer errors than EARI on the right red edge and is purer in the light green beam on the right. However, in a flat background, all algorithms have an unavoidable increase in noise, which is quite obvious compared with GT images.

Fig. 15. Subjective comparison of the intensity image and the AoP-DoP visualization for MPFA demosaicking.

Download Full Size | PDF

Table 5. Comparison results of different methods for CPFA demosaicking (Average of 40 scenes)

View Table | View all tables in this article

In this dataset test, the results of S0₂ are more stable and outstanding. We believe this is due to the fact that some of the images in this set were imaged in an environment of polarized light, which caused better results for the single polarized images at different angles. The calculation of S0₂ ultimately causes it to give generally better results than S0₁ due to the way. Regardless of the calculation method used, our approach is superior in most metrics.

Overall, for both CPFA and MPFA, our algorithm maintains the leading PSNR value in a single polarized image, and the accuracy of the Stokes parameter also has good results. Moreover, the edge processing of the actual image is closer to that of the GT image, with fewer wrong pixel values. But, the demosaicing algorithm still has huge room for improvement, so that the output result can match the real environment as much as possible and present a clear image.

5. Conclusion

In this paper, we proposed a modified Newton residual interpolation demosaicking algorithm based on the EARI. For the features of the RI framework, we focus on the improvement of the guided filtered image and the guided process. We draw lessons from the previous article about the enhancement of edge processing to enhance the demosaicking effect. We first improved the Newton interpolation algorithm by judging and correcting the edges of the predicted image after the initial interpolation, resulting in a guide image with excellent edge retention. Then, we added a parameter about the variance value in the local window during the guide process, through this parameter, the flat area and the edge area have different weights, so as to get better guiding parameters, and finally get demosaicking with improved effect image. The algorithm can be extended to CPFA demosaicking framework. Through experimental comparison, we have a greater improvement in objective indicators than the latest demosaicking algorithm. As shown in comparing Fig. 11, Fig. 12, Fig. 14, and Fig. 15, our algorithm has less zipper effect and has fewer errors in the calculation results of Stokes parameters, and the demosaicking results are more accurate. In future work, we can improve this algorithm to include denoising methods and enhance the effect.

Funding

Science and Technology Planning Project of Guangdong Province (2019B121203006); Guangdong Key Laboratory of Advanced IntelliSense Technology.

Acknowledgment

We thank the authors who provided the open source code and the experimental data. Thanks to their work, we can get better results.

Disclosures

The authors declare no conflicts of interest.

Data availability

All experimental results are available in Ref. [48].

References

1. J. J. Foster, S. E. Temple, M. J. How, I. M. Daly, C. R. Sharkey, D. Wilby, and N. W. Roberts, “Polarisation vision: overcoming challenges of working with a property of light we barely see,” Naturwissenschaften (1913-2014) 105(3-4), 27 (2018). [CrossRef]

2. G. Courtier, P. J. Lapray, J. B. Thomas, and I. Farup, “Correlations in joint spectral and polarization imaging,” Sensors 21(1), 6–12 (2020). [CrossRef]

3. O. A. Elgendy, A. Gnanasambandam, S. H. Chan, and J. Ma, “Low-Light Demosaicking and Denoising for Small Pixels Using Learned Frequency Selection,” IEEE Trans. Comput. Imaging 7, 137–150 (2021). [CrossRef]

4. M. Reda, Y. Zhao, and J. C. W. Chan, “Polarization Guided Autoregressive Model for Depth Recovery,” IEEE Photonics J. 9(3), 1–16 (2017). [CrossRef]

5. C. He, H. He, J. Chang, Y. Dong, S. Liu, N. Zeng, Y. He, and H. Ma, “Characterizing microstructures of cancerous tissues using multispectral transformed Mueller matrix polarization parameters,” Biomed. Opt. Express 6(8), 2934 (2015). [CrossRef]

6. Y. W. Zhou, Z. F. Li, J. Zhou, N. Li, X. H. Zhou, P. P. Chen, Y. L. Zheng, X. S. Chen, and W. Lu, “High extinction ratio super pixel for long wavelength infrared polarization imaging detection based on plasmonic microcavity quantum well infrared photodetectors,” Sci. Rep. 8(1), 15070 (2018). [CrossRef]

7. K. Gurton, M. Felton, R. Mack, D. LeMaster, C. Farlow, M. Kudenov, and L. Pezzaniti, “MidIR and LWIR polarimetric sensor comparison study,” Proc. SPIE 7664, 76640L (2010). [CrossRef]

8. A. J. Brown, T. I. Michaels, S. Byrne, W. Sun, T. N. Titus, A. Colaprete, M. J. Wolff, G. Videen, and C. J. Grund, “The case for a modern multiwavelength, polarization-sensitive LIDAR in orbit around Mars,” J. Quant. Spectrosc. Radiat. Transfer 153, 131–143 (2015). [CrossRef]

9. A. J. Brown, “Equivalence relations and symmetries for laboratory, LIDAR, and planetary Müeller matrix scattering geometries,” J. Opt. Soc. Am. A 31(12), 2789 (2014). [CrossRef]

10. X. Zhao, A. Bermak, and F. Boussaid, “A CMOS digital pixel sensor with photo-patterned micropolarizer array for real-time focal-plane polarization imaging,” in 2008 IEEE-BIOCAS Biomedical Circuits and Systems Conference (2008), pp. 145–148.

11. V. Gruev, R. Perkins, and T. York, “CCD polarization imaging sensor with aluminum nanowire optical filters,” Opt. Express 18(18), 19087 (2010). [CrossRef]

12. C. Xu, J. Ma, C. Ke, Y. Huang, Z. Zeng, and W. Weng, “Numerical study of a DoFP polarimeter based on the self-organized nanograting array,” Opt. Express 26(3), 2517 (2018). [CrossRef]

13. V. Gruev, J. Spiegel, N. B. T.-C. and, and S. Engheta, “Image sensor with focal plane extraction of polarimetric information,” in ISCAS 2006 Proceedings (2006).

14. X. Zhao, F. Boussaid, A. Bermak, and V. G. Chigrinov, “Thin photo-patterned micropolarizer array for CMOS image sensors,” IEEE Photonics Technol. Lett. 21(12), 805–807 (2009). [CrossRef]

15. W. Kang, J. Chu, X. Zeng, and Y. Fan, “Large-area flexible infrared nanowire grid polarizer fabricated using nanoimprint lithography,” Appl. Opt. 57(18), 5230 (2018). [CrossRef]

16. B. M. Ratliff, C. F. LaCasse, and J. Scott Tyo, “Interpolation strategies for reducing IFOV artifacts in microgrid polarimeter imagery,” Opt. Express 17(11), 9112 (2009). [CrossRef]

17. E. S. Kaur and V. K. Banga, “Enriched image demosaicing using illuminate normalization and content based color filter array,” in International Conference on Information Communication and Embedded Systems (2016).

18. Q. Jin, G. Facciolo, and J. M. Morel, “A review of an old dilemma: Demosaicking first, or denoising first?” in IEEE Comput. Soc. Conf. Comput. Vis. Pattern Recognit. Work. (2020), 2169–2179.

19. Y. Giménez, P.-J. Lapray, A. Foulonneau, and L. Bigué, “Calibration algorithms for polarization filter array camera: survey and evaluation,” J. Electron. Imag. 29(04), 1 (2020). [CrossRef]

20. S. Gao and V. Gruev, “Bilinear and bicubic interpolation methods for division of focal plane polarimeters,” Opt. Express 19(27), 26161 (2011). [CrossRef]

21. N. Li, Y. Zhao, Q. Pan, and S. G. Kong, “Demosaicking DoFP images using Newton’s polynomial interpolation and polarization difference model,” Opt. Express 27(2), 1376 (2019). [CrossRef]

22. R. Wu, Y. Zhao, N. Li, and S. G. Kong, “Polarization image demosaicking using polarization channel difference prior,” Opt. Express 29(14), 22066 (2021). [CrossRef]

23. R. Wang, W. Gao, F. Wang, and C. Shen, “Non-Local Sparse Representation Method for Demosaicing of Single DoFP Polarimetric Image,” 12th Int. Conf. Commun. Softw. Networks (2020), 259–263.

24. J. Zhang, H. Luo, R. Liang, A. Ahmed, X. Zhang, B. Hui, and Z. Chang, “Sparse representation-based demosaicing method for microgrid polarimeter imagery,” Opt. Lett. 43(14), 3265 (2018). [CrossRef]

25. J. Zhang, J. Shao, H. Luo, X. Zhang, B. Hui, Z. Chang, and R. Liang, “Learning a convolutional demosaicing network for microgrid polarimeter imagery,” Opt. Lett. 43(18), 4534 (2018). [CrossRef]

26. G. C. Sargent, B. M. Ratliff, and V. K. Asari, “Conditional generative adversarial network demosaicing strategy for division of focal plane polarimeters,” Opt. Express 28(25), 38419 (2020). [CrossRef]

27. X. Zeng, Y. Luo, X. Zhao, and W. Ye, “An end-to-end fully-convolutional neural network for division of focal plane sensors to reconstruct S 0, DoLP, and AoP,” Opt. Express 27(6), 8566 (2019). [CrossRef]

28. D. Kiku, Y. Monno, M. Tanaka, and M. Okutomi, “Minimized-Laplacian residual interpolation for color image demosaicking,” Digit. Photogr. X 9023, 90230L (2014). [CrossRef]

29. M. Morimatsu, Y. Monno, M. Tanaka, and M. Okutomi, “Monochrome and Color Polarization Demosaicking Using Edge-Aware Residual Interpolation,” in Proc. - Int. Conf. Image Process. (2020), 2571–2575.

30. T. Jiang, D. Wen, Z. Song, W. Zhang, Z. Li, X. Wei, and G. Liu, “Minimized Laplacian residual interpolation for DoFP polarization image demosaicking,” Appl. Opt. 58(27), 7367 (2019). [CrossRef]

31. T. Mu, C. Zhang, and R. Liang, “Demonstration of a snapshot full-Stokes division-of-aperture imaging polarimeter using Wollaston prism array,” J. Opt. 17(12), 125708 (2015). [CrossRef]

32. D. Kiku, Y. Monno, M. Tanaka, and M. Okutomi, “Residual interpolation for color image demosaicking,” in IEEE Int. Conf. Image Process. (2013), 2304–2308.

33. K. He, J. Sun, and X. Tang, “Guided Image Filtering BT - Computer Vision – ECCV 2010,” in K. Daniilidis, P. Maragos, and N. Paragios, eds. (Springer, 2010), pp. 1–14.

34. Y. Monno, D. Kiku, M. Tanaka, and M. Okutomi, “Adaptive residual interpolation for color and multispectral image demosaicking,” Sensors 17(12), 2787 (2017). [CrossRef]

35. S. Liu, J. Chen, Y. Xun, X. Zhao, and C. H. Chang, “A New Polarization Image Demosaicking Algorithm by Exploiting Inter-Channel Correlations with Guided Filtering,” IEEE Trans. on Image Process. 29, 7076–7089 (2020). [CrossRef]

36. Z. Lu, B. Long, K. Li, and F. Lu, “Effective Guided Image Filtering for Contrast Enhancement,” IEEE Signal Process. Lett. 25(10), 1585–1589 (2018). [CrossRef]

37. K. He, J. Sun, and X. Tang, “Guided image filtering,” IEEE Trans. Pattern Anal. Mach. Intell. 35(6), 1397–1409 (2013). [CrossRef]

38. S. Mihoubi, P. J. Lapray, and L. Bigué, “Survey of demosaicking methods for polarization filter array images,” Sensors 18(11), 3688 (2018). [CrossRef]

39. J. Zhang, H. Luo, B. Hui, and Z. Chang, “Image interpolation for division of focal plane polarimeters with intensity correlation,” Opt. Express 24(18), 20799 (2016). [CrossRef]

40. S. Mihoubi, O. Losson, B. Mathon, and L. Macaire, “Multispectral Demosaicing Using Pseudo-Panchromatic Image,” IEEE Trans. Comput. Imaging 3(4), 982–995 (2017). [CrossRef]

41. “EARI,” https://github.com/ymonno/EARI-Polarization-Demosaicking.

42. “ICPC,” https://github.com/Junchao2018/DoFP-Interpolation.

43. “NP,” https://github.com/polwork/Newton-Polynomial-Interpolation-1.0/.

44. X. Liu, M. Tanaka, and M. Okutomi, “Practical Signal-Dependent Noise Parameter Estimation from a Single Noisy Image,” IEEE Trans. on Image Process. 23(10), 4361–4371 (2014). [CrossRef]

45. “PolarDem,” http://www.ok.sc.e.titech.ac.jp/res/PolarDem/index.html.

46. S. Qiu, Q. Fu, C. Wang, and W. Heidrich, “Linear Polarization Demosaicking for Monochrome and Colour Polarization Focal Plane Arrays,” Computer Graphics Forum 40(6), 77–89 (2021). [CrossRef]

47. “PolarData,” doi:10.25781/KAUST-2VA2X.

48. “MNRI,” https://github.com/xxwudi508/MNRI/.

Method		$I_{0}$	$I_{45}$	$I_{90}$	$I_{135}$	$S 0$	$S 1$	$S 2$	DoLP	AoLP
EARI		47.39	44.91	47.84	44.63	49.62	51.48	47.83	36.79	17.13
MNRI	R = 3(S0₁)^a,^b	47.74	45.27	48.21	44.71	49.85	52.58	47.59	35.9	16.60
	R = 6(S0₁)	47.80	45.40	48.27	44.86	49.92	52.64	47.73	36.08	16.45
	R = 9(S0₁)	47.80	45.43	48.28	44.90	49.94	52.64	47.77	36.09	16.43
	R = 12(S0₁)	47.81	45.45	48.27	44.92	49.95	52.64	47.79	36.09	16.42
	R = 15(S0₁)	47.81	45.46	48.28	44.94	49.96	52.63	47.80	36.09	16.42
	R = 18(S0₁)	47.81	45.47	48.27	44.94	49.96	52.63	47.81	36.09	16.42

Method	PSNR								Angle error
Method	$I_{0}$	$I_{45}$	$I_{90}$	$I_{135}$	$S 0$	$S 1$	$S 2$	DoLP	AoLP
Bilinear	42.34	41.58	42.50	41.58	44.89	46.14	45.03	33.70	21.36
Bicubic	43.45	42.48	43.63	42.48	46.22	47.00	45.73	34.46	20.64
ICPC	43.10	42.22	43.23	42.22	45.78	47.01	45.73	34.75	20.50
PPID	46.52	44.52	46.91	44.34	48.94	50.56	47.59	36.96	17.65
Newton	47.47	44.42	47.96	44.20	49.91	51.92	46.88	35.38	17.47
EARI	47.39	44.91	47.84	44.63	49.62	51.48	47.83	36.79	17.13
NoMN^a	47.51	44.92	47.89	44.65	49.70	51.58	47.44	36.78	17.08
NoGF^b	47.71	45.15	48.18	44.72	49.75	52.54	47.59	35.89	16.62
MNRI(S02)	47.97	45.05	48.66	44.58	49.84	52.83	47.39	35.95	16.58
MNRI(S01)	47.81	45.46	48.28	44.94	49.96	52.63	47.80	36.09	16.42

Method (color\|polarization)		CPSNR								Angle error
		$I_{0}$	$I_{45}$	$I_{90}$	$I_{135}$	$S 0$	$S 1$	$S 2$	DoLP	AoLP
Bilinear		35.32	34.94	35.47	35.00	36.31	43.29	41.28	31.20	24.98
RI	Bilinear	38.34	37.79	38.50	37.86	40.03	44.09	42.58	31.96	24.15
	Bicubic	38.65	38.05	38.81	38.12	40.43	44.30	42.72	32.05	24.11
	ICPC	38.61	38.01	38.77	38.09	40.33	44.49	42.87	32.36	23.86
	PPID	39.37	38.68	39.57	38.71	40.73	46.34	44.04	33.98	22.40
	NP	39.45	38.63	39.64	38.65	40.49	47.88	44.45	33.79	21.27
	EARI	39.41	38.72	39.62	38.72	40.76	46.49	44.10	33.74	22.18
	MNRI(S02)	39.43	38.74	39.67	38.72	40.38	48.04	44.68	34.09	20.95
	MNRI(S01)	39.43	38.78	39.61	38.76	40.40	48.20	44.60	34.03	21.00

Method	PSNR								Angle error
Method	$I_{0}$	$I_{45}$	$I_{90}$	$I_{135}$	$S 0$	$S 1$	$S 2$	DoLP	AoLP
Bilinear	40.71	42.01	42.23	41.11	44.49	44.58	44.81	34.78	12.72
Bicubic	40.64	41.95	42.16	41.06	45.60	43.71	43.92	33.76	13.40
ICPC	44.52	45.87	46.16	44.82	48.51	48.60	48.80	38.97	11.16
PPID	45.98	47.23	46.92	46.53	48.93	49.12	50.32	40.68	10.01
NP	46.44	47.42	47.06	47.06	51.07	49.95	50.59	38.76	9.79
EARI	46.64	47.91	48.07	47.07	50.58	50.83	51.16	40.95	9.96
MNRI(S0₂)	47.19	48.34	48.38	47.67	51.18	50.95	51.41	40.15	9.55
MNRI(S0₁)	46.47	48.23	47.54	47.29	50.64	50.26	51.05	39.70	9.73

Method (color\|polarization)		CPSNR								Angle error
		$I_{0}$	$I_{45}$	$I_{90}$	$I_{135}$	$S 0$	$S 1$	$S 2$	DoLP	AoLP
Bilinear		37.71	38.80	38.98	37.98	39.57	45.23	45.79	35.53	12.93
RI	Bilinear	40.27	41.41	41.58	40.61	43.83	44.22	44.50	34.38	13.64
	Bicubic	40.19	41.34	41.50	40.54	44.57	43.45	43.70	33.47	14.22
	ICPC	42.92	44.02	44.23	43.18	46.36	47.37	47.64	37.55	12.47
	PPID	43.21	44.51	44.34	44.02	46.51	48.31	48.54	38.32	12.12
	NP	43.78	44.67	44.56	44.24	47.01	48.69	49.28	37.85	11.09
	EARI	43.84	44.90	45.05	44.14	46.99	48.74	49.10	38.79	11.63
	MNRI(S0₂)	44.11	45.08	45.20	44.44	46.92	49.23	49.68	38.77	10.99
	MNRI(S0₁)	43.78	45.03	44.82	44.19	46.75	48.89	49.22	38.44	11.16

Modified Newton-residual interpolation for division of focal plane polarization image demosaicking

Abstract

1. Introduction

1.1 Background

1.2 Linear polarization computation

2. Related work

3. Proposed method

3.1 Modified Newton interpolation

3.2 Improved GF

3.3 CPFA demosaicking

4. Experimental result

4.1 Experiment 1

4.1.1 MPFA demosaicking

4.1.2 CPFA demosaicking

4.2 Experiment 2

4.2.1 MPFA demosaicking

4.2.2 MPFA demosaicking

5. Conclusion

Funding

Acknowledgment

Disclosures

Data availability

References

Data availability

Cited By

Figures (15)

Tables (5)

Equations (38)

Optics Express