Reconstructing Stokes parameters from non-uniform division-of-focal-plane modulation
Zhaoxiang Jiang Qingchuan Zhang Shangquan Wu Tan Xu Yong Su Chuanbiao Bai
physics.optics eess.SP




mymainaddress]Zhaoxiang Jiang

mymainaddress]Qingchuan Zhangmycorrespondingauthor1
[mycorrespondingauthor1]Corresponding author
zhangqc@ustc.edu.cn

mymainaddress]Shangquan Wumycorrespondingauthor
[mycorrespondingauthor]Corresponding author
wushq@ustc.edu.cn

mymainaddress]Tan Xu

mymainaddress]Yong Su

mymainaddress]Chuanbiao Bai




[mymainaddress]CAS Key Laboratory of Mechanical Behavior and Design of Materials, Department of Modern Mechanics, University of Science and Technology of China, Hefei 230027, China




Division-of-focal-plane modulation is a powerful technique for real-time polarization imaging. This technique, however, suffers from the non-uniformity of the performance of linear polarization filters and photodetectors. We study the Stokes parameters reconstruction from the division-of-focal-plane modulation in the presence of the non-uniformity. Two reconstruction methods, named as ordinary least-squares and smoothing regularization methods, are proposed. The performance of the proposed methods are evaluated through Fourier analysis, numerical simulations, and experiments. The results indicate that the proposed methods can effectively mitigate the reconstruction errors and artifacts caused by the non-uniformity, and therefore may further facilitate the practical application of the division-of-focal-plane technique.



	Division-of-focal-plane Polarization imaging Fringe analysis Non-uniformity
	[2010] 00-0199-00









§ INTRODUCTION

Polarization imaging aims to measure the polarization information described by Stokes parameters and their derivatives such as degree of linear polarization (DoLP) and angle of polarization (AoP). It has a wide range of applications in various fields such as remote sensing <cit.>, biomedical diagnosis <cit.>, and interferometry <cit.>. The modulation technique in polarization imaging can be classified into four categories: division-of-time (DoT) <cit.>, division-of-amplitude <cit.>,  division-of-aperture <cit.>, and division-of-focal-plane (DoFP) <cit.>. The DoFP technique achieves the polarization modulation by integrating a polarization filter array in front of a photodetector array. The most commonly used polarization filter array arrangement is shown in Fig. <ref>, which is composed by 2×2 periodically patterned 0^∘, 45^∘, 90^∘, and 135^∘ linear polarization filters. The DoFP polarimeters have advantages of compact structure and high temporal resolution, and therefore are especially suitable for real-time polarization imaging.



A fundamental problem in the DoFP technique is the Stokes parameters reconstruction. Since each pixel can only capture the polarization information in one orientation, the measurement of the Stokes parameters is incomplete. Mathematically, reconstructing the Stokes parameters from the DoFP polarization modulation is an ill-posed inverse problem. Many methods have been proposed to reconstruct the Stokes parameters <cit.>. Most of these methods are spatial domain interpolation-based methods <cit.> and frequency domain filtering-based methods <cit.>. In the interpolation-based methods, the DoFP image is split into 0^∘, 45^∘, 90^∘, and 135^∘ polarization images. After interpolating the missing pixel values in these polarization images, the interpolation-based methods determine the Stokes parameters through the ordinary least-squares criterion. The interpolation algorithms are usually convolutional, including nearest-neighbor, bilinear, bicubic, and natural bicubic spline interpolation algorithms <cit.>. Some edge-preserved interpolation algorithms have also been proposed recently <cit.>. Based on the characteristics of the spectrum of the DoFP image, the filtering-based methods use the filter transfer functions constructed by window functions to reconstruct the Stokes parameters. The window functions used in previous studies including Hamming <cit.>, Gaussian <cit.>, and Planck-taper <cit.> window functions. However, the interpolation-based and filtering-based methods are only suitable for the theoretical case that the Stokes parameters are periodically modulated. In practice, the manufacturing imperfects of the DoFP polarimeters cause the non-uniformity of the performance of the linear polarization filters and photodetectors. The non-uniformity is characterized as the differences of the major and minor principal transmittances of the linear polarization filters, the differences of the gains and dark offsets of the photodetectors, and the deviations between the actual and designed orientations of the linear polarization filters <cit.>. The non-uniformity destroys the periodicity of the polarization modulation. Consequently, the interpolation-based and filtering-based methods will fail in practical applications since these methods are unable to tackle the reconstruction errors and artifacts caused by the non-uniformity. 


In this paper, we study the Stokes parameters reconstruction from the DoFP modulation in the presence of the non-uniformity. We propose two reconstruction methods that can tackle the reconstruction errors and artifacts caused by the non-uniformity. The proposed methods are inspired by the classical Lucas-Kanade method <cit.> and Horn-Schunck method <cit.> in optical flow estimation. One is the ordinary least-squares method (OLSM), which reconstructs the Stokes parameters under the local constant assumption that the Stokes parameters are constant functions in 2×2 subsets. The basic idea of this method has been reported in Ref. <cit.> and our previous study <cit.>. Here, we further add a four-subset averaging strategy, present in-depth theoretical analyses, and explain the relationship between the OLSM and interpolation-based methods. The other is the smoothing regularization method (SRM), which reconstructs the Stokes parameters under the global smoothing assumption that the Stokes parameters are spatially smooth. This method has more similarities to the filtering-based methods.

This paper is organized as follows: in Section <ref>, a linear pixel model is introduced to characterize the non-uniformity; in Section <ref>, the OLSM and SRM are presented to reconstruct the Stokes parameters; in Section <ref>, Fourier analysis and numerical simulations are used to evaluate the reconstruction errors of the OLSM and SRM; in Section <ref>, the choice of the regularization parameters in the SRM is discussed; in Section <ref>, the performance of the OLSM, SRM, and interpolation-based and filtering-based methods is evaluated and compared through two experiments; in Section <ref>, the performance of the OLSM and SRM is summarized.


§ LINEAR PIXEL MODEL

After integrating with the array composed by the linear polarization filters, the photodetector array can be regarded as being sensitive to the first three Stokes parameters. Assuming the performance of the linear polarization filters and photodetectors is linear, for each pixel, we have <cit.>

    i(x,y)=m_0(x,y)s_0(x,y)+m_1(x,y)s_1(x,y)+m_2(x,y)s_2(x,y)+d(x,y).

Here, x and y are the horizontal and vertical pixel coordinates, respectively, i represents the intensity of the DoFP image, s_0, s_1, and s_2 represent the first three Stokes parameters, m_0, m_1, and m_2 are the modulation parameters of s_0, s_1, and s_2, respectively, and d represents the dark offset of the photodetector. Specifically, m_0, m_1, and m_2 are expressed as

    {[           m_0(x,y)=1/2g(x,y)[k_1(x,y)+k_2(x,y)]; m_1(x,y)=1/2g(x,y)[k_1(x,y)-k_2(x,y)]cos2θ(x,y); m_2(x,y)=1/2g(x,y)[k_1(x,y)-k_2(x,y)]sin2θ(x,y) ]..

Here, g represents the gain of the photodetector, k_1 and k_2 represent the major and minor principal transmittances of the linear polarization filter, respectively, and θ represents the orientation of the linear polarization filter.

In the interpolation-based and filtering-based methods, the modulation parameters are regarded as periodically ideal values, expressed as <cit.>

    {[                      m_0(x,y)=1; m_1(x,y)=1/2[cos(π x)+cos(π y)]; m_2(x,y)=1/2[cos(π x)-cos(π y)] ],.

corresponding to g(x,y)=2, k_1(x,y)=1, k_2(x,y)=0, and θ(x,y) takes 0, π/4, π/2, or 3π/4 according to the arrangement shown in Fig. <ref>. In practice, due to the existence of the non-uniformity, the major and minor principal transmittances of the linear polarization filters and the gains and dark offsets of the photodetectors generally show individual differences, and the actual orientations of the linear polarization filters are also deviated from the designed orientations. To tackle the reconstruction errors and artifacts caused by the non-uniformity, these parameters need to be calibrated. The calibration of these parameters has been well discussed in Refs. <cit.>, therefore we consider m_0, m_1, m_2, and d as known. Notice that here the term calibration refers to experimental processes for measuring the modulation and offset parameters, instead of numerical methods for correcting the non-uniform intensity responses in the DoFP image to ideal <cit.>. To normalize the modulation parameters and choose an orientation as the reference of 0^∘, we apply the scaling and rotation transformation given in <ref> to the modulation parameters. And since the influence caused by the non-uniformity of d can be mitigated by simply subtracting d from i, we will omit d for brevity.


To generalize our discussions, we consider two sets of modulation parameters. One is the ideal modulation parameters given in Eq. (<ref>). The scatterplots of the ideal modulation parameters are shown in the first row of Fig. <ref>. We use the ideal modulation parameters to give the theoretical evaluations of the performance of the proposed methods, and explain the relationship between the proposed methods and the interpolation-based and filtering-based methods. The other is the non-uniform modulation parameters obtained from the calibration of our self-developed DoFP polarimeters. The second row of Fig. <ref> shows the scatterplots of the non-uniform modulation parameters. And Table <ref> gives the means and standard deviations of g·(k_1+k_2), g·(k_1-k_2), and θ in the non-uniform modulation parameters according to the designed polarization orientations. It can be seen that the non-uniform modulation parameters show strong non-uniformity deviations. We use the non-uniform modulation parameters to illustrate the abilities of the OLSM and SRM for mitigating the reconstruction errors and artifacts caused by the non-uniformity.






§ STOKES PARAMETERS RECONSTRUCTION

Reconstructing the Stokes parameters requires dealing with the underdetermined system of linear equations composed by Eq. (<ref>). This problem has infinite solutions. Generally, some prior constraints need to be introduced to find a desired solution. In this section, the OLSM and SRM are presented to solve the underdetermined system of linear equations by applying the local constant assumption and global smoothing assumption, respectively.


 §.§ Ordinary least-squares method

For each 2×2 subset in the DoFP image, there are four equality constraints and twelve unknown Stokes parameters. By applying the local constant assumption that the Stokes parameters are constant functions in 2×2 subsets, the number of the unknowns in each 2×2 subset is reduced to three. Then the Stokes parameters can be determined according to the ordinary least-squares criterion, expressed as

    [ ŝ_0(x + 12,y + 12); ŝ_1(x + 12,y + 12); ŝ_2(x + 12,y + 12) ] 
    = 
    [         m_0(x,y)         m_1(x,y)         m_2(x,y);     m_0(x,y + 1)     m_1(x,y + 1)     m_2(x,y + 1);     m_0(x + 1,y)     m_1(x + 1,y)     m_2(x + 1,y); m_0(x + 1,y + 1) m_1(x + 1,y + 1) m_2(x + 1,y + 1) ]
    ^†[         i(x,y);     i(x,y + 1);     i(x + 1,y); i(x + 1,y + 1) ].

Here, ŝ_0, ŝ_1, and ŝ_2 represent the reconstructed Stokes parameters, and [∙]^† represents the pseudo-inverse of a matrix. The feasibility of using Eq. (<ref>) to reconstruct the Stokes parameters has been demonstrated in Ref. <cit.>. However, when substituting the ideal modulation parameters into Eq. (<ref>), it can be found that the reconstruction results of Eq. (<ref>) are equivalent to that of the nearest-neighbor interpolation-based method. This indicates that when facing the spatial variations of the Stokes parameters, Eq. (<ref>) performs as poorly as the nearest-neighbor interpolation-based reconstruction.

Considering that each pixel is contained in four different 2×2 subsets, we further apply a four-subset averaging strategy to obtain the final reconstruction results, expressed as


    ŝ_k(x,y)  = 14 [ŝ_k(x - 12,y - 12) + ŝ_k(x - 12,y + 12) + ŝ_k(x + 12,y - 12) + ŝ_k(x + 12,y + 12)].

Here, the subscript k takes 0, 1, and 2. The four-subset averaging strategy is also an optimization in the sense of the ordinary least-squares. With this strategy, when substituting the ideal modulation parameters, it can be found that the reconstruction results of the OLSM (Eqs. (<ref>) and (<ref>)) are  equivalent to that of the bilinear interpolation-based method.


 §.§ Smoothing regularization method

Regularization is a commonly used optimization technique for solving underdetermined problems <cit.>. We assume that the Stokes parameters are spatially smooth and apply the regularization technique to find the globally smoothing solution of the system of linear equations composed by Eq. (<ref>). The Stokes parameters are determined by minimizing the objective function L defined as

    L(ŝ_0,ŝ_1,ŝ_2)    = ∑_x,y[ m_0(x,y)ŝ_0(x,y) + m_1(x,y)ŝ_1(x,y) + m_2(x,y)ŝ_2(x,y) - i(x,y)]^2
       + λ _0R(ŝ_0) + λ _1R(κ _1ŝ_1 + ŝ_22) + λ _2R(κ _2ŝ_1 - ŝ_22).

Here, the first term on the right side of Eq. <ref> is the fidelity term used to penalize the deviation between the reconstructed Stokes parameters and the constraint of Eq. (<ref>), R(ŝ_0), R(κ _1ŝ_1+ŝ_22), and R(κ _2ŝ_1-ŝ_22) are the regularization terms, R is the discrete thin-plate energy functional used to introduce the constraint of the spatial smoothness, defined as

    R(f)    = ∑_x,y[ f(x - 1,y) - 2f(x,y) + f(x + 1,y)]^2
       + ∑_x,y[ f(x,y - 1) - 2f(x,y) + f(x,y + 1)]^2
       + 2∑_x,y[ f(x,y) - f(x + 1,y) - f(x,y + 1) + f(x + 1,y + 1)]^2,

λ_0, λ_1, and λ_2 are the regularization parameters used to control the weights of the regularization terms, and κ_1 and κ_2 are two parameters used to compensate the change of the relative weights between the fidelity term and regularization terms caused by the non-uniformity, defined as

    {[ κ_1 = ∑_x,ym_1(x,y)cos (π x)  + ∑_x,ym_2(x,y)cos (π x)/∑_x,ym_0(x,y); κ_2 = ∑_x,ym_1(x,y)cos (π y)  - ∑_x,ym_2(x,y)cos (π y)/∑_x,ym_0(x,y) ]..


The desired Stokes parameters should satisfy the Euler-Lagrange equation of the objective function, expressed as

    {[                                                                   ∂ L/∂ŝ_0 = m_0^2ŝ_0 + m_0m_1ŝ_1 + m_0m_2ŝ_2 - m_0i + λ _0∇ ^4ŝ_0 = 0; ∂ L/∂ŝ_1 = m_0m_1ŝ_0 + m_1^2ŝ_1 + m_1m_2ŝ_2 - m_1i + 1/4( λ _1κ _1^2 + λ _2κ _2^2 )∇ ^4ŝ_1 + 1/4( λ _1κ _1^2 - λ _2κ _2^2 )∇ ^4ŝ_2 = 0; ∂ L/∂ŝ_2 = m_0m_2ŝ_0 + m_1m_2ŝ_1 + m_2^2ŝ_2 - m_2i + 1/4( λ _1κ _1^2 - λ _2κ _2^2 )∇ ^4ŝ_1 + 1/4( λ _1κ _1^2 + λ _2κ _2^2 )∇ ^4ŝ_2 = 0 ]..

Here, ∇ ^4 represents the discrete biharmonic operator, corresponding to the variation of the discrete thin-plate functional. The discrete biharmonic operator is implemented by the convolution operation

    ∇ ^4f = f * 
    [   0   0   1   0   0;   0   2 - 8   2   0;   1 - 8  20 - 8   1;   0   2 - 8   2   0;   0   0   1   0   0 ].

Here, * represents the convolution operation.

Equation (<ref>) can be solved by gradient descent algorithms. We gave a MATLAB code implementation of the SRM in Ref. <cit.>. Since the objective function is a convex function, the calculations can well converge to the global optimal solution that minimizes the objective function.

It is worth pointing out that the proposed methods and numerical calibration methods tackle the non-uniformity based on the same pixel model but have different purposes. The proposed methods aims to directly reconstruct the Stokes parameters, while the numerical calibration methods aims to correct the non-uniform intensity responses in the DoFP image. The relationship between them is further explained in <ref>.


§ RECONSTRUCTION ERROR EVALUATIONS

In the DoFP modulation, the measurement information from different pixel coordinates is combined to reconstruct the Stoke parameters. Consequently, the spatial variations of the Stokes parameters will introduce reconstruction errors, which is significantly different from the other modulation techniques. 
In this section, the reconstruction errors of the OLSM and SRM for the different frequency components of s_0, (s_1+s_2)/2, and (s_1-s_2)/2 are evaluated. Firstly, Fourier analysis is applied to evaluate the reconstruction errors in the ideal case that the Stokes parameters are modulated by the ideal modulation parameters. Secondly, numerical simulations are applied to evaluate the reconstruction errors in the non-uniform case that the Stokes parameters are modulated by the non-uniform modulation parameters.


 §.§ Fourier analysis

When the Stokes parameters are modulated by the ideal modulation parameters, substituting Eq. (<ref>) into Eq. (<ref>), the discrete Fourier transform (DFT) of i is expressed as

    I(u,v) = S_0(u,v) + 12[S_1(u + 12,v) + S_2(u + 12,v)] + 12[S_1(u,v + 12) - S_2(u,v + 12)].

Here, u and v represent the horizontal and vertical frequency coordinates, respectively, and the uppercase symbols represent the DFT of the corresponding lowercase symbols. Equation. (<ref>) indicates that the frequency components of s_0 are located in the center region of the spectrum of i, while the frequency components of (s_1+s_2)/2 and (s_1-s_2)/2 are shifted into the horizontal and vertical border regions, respectively. Figure. <ref>(a) gives an example of the spectrum of i where the Stokes parameters are modulated by the ideal modulation parameters (see Section <ref>). It can be seen that there are noticeable carrier peaks at the frequency coordinates (0,0), (±1/2,0), and (0,±1/2).



The distribution of the frequency components of s_0, (s_1+s_2)/2, and (s_1-s_2)/2 allows us to reconstruct the Stokes parameters through frequency domain filtering, expressed as

    {[                           Ŝ_0(u,v) = I(u,v) ·H_0(u,v); 1/2[Ŝ_1(u + 12,v) + Ŝ_2(u + 12,v)] = I(u,v) ·H_1(u,v); 1/2[Ŝ_1(u,v + 12) - Ŝ_2(u,v + 12)] = I(u,v) ·H_2(u,v) ]..

Here, H_0, H_1, and H_2 represent the filter transfer functions which aim to retrieve the frequency components of s_0, (s_1+s_2)/2, and (s_1-s_2)/2, respectively. In the filtering-based methods, the filter transfer functions are constructed by window functions such as Hamming, Gaussian, and Planck-taper window functions <cit.>. In fact, ignoring the differences around the image boundaries of the reconstruction results, we find that the convolutional interpolation-based methods, OLSM, and SRM have equivalent implementations in the frequency domain. The filter transfer functions of the mainstream convolutional interpolation-based methods are given in <ref>, including the nearest-neighbor, bilinear, bicubic, and natural bicubic spline interpolation-based methods. For the OLSM, substituting Eq. (<ref>) into Eqs. (<ref>) and (<ref>), and performing the DFT, after some simplifications, we have

    {[      H_0(u,v) = 1/4[ 1 + cos (2π u) ][ 1 + cos (2π v) ]; H_1(u,v) = 1/4[ 1 + cos (2π u + π ) ][ 1 + cos (2π v) ]; H_2(u,v) = 1/4[ 1 + cos (2π u) ][ 1 + cos (2π v + π ) ] ]..

The first row of Fig. <ref> shows the filter transfer functions of the OLSM and their full-widths-at-half-maximum (FWHMs). Notice that the OLSM is equivalent to the bilinear interpolation-based method in this case. For the SRM, substituting Eq. (<ref>) into Eq. (<ref>), and performing the DFT, after the same simplifications, we have

    {[           H_0(u,v) = [ 1 + λ _0G(u,v) + λ _0G(u,v)/λ _1G(u + 12,v) + λ _0G(u,v)/λ _2G(u,v + 12)]^ - 1; H_1(u,v) = [ 1 + λ _1G(u + 12,v) + λ _1G(u + 12,v)/λ _0G(u,v) + λ _1G(u + 12,v)/λ _2G(u,v + 12)]^ - 1; H_2(u,v) = [ 1 + λ _2G(u,v + 12) + λ _2G(u,v + 12)/λ _0G(u,v) + λ _2G(u,v + 12)/λ _1G(u + 12,v)]^ - 1 ]..

Here, G is the DFT of the discrete biharmonic operator, expressed as

    G(u,v)    = 20 - 16cos (2π u) - 16cos (2π v) + 4cos (2π u + 2π v) 
        + 4cos (2π u - 2π v) + 2cos (4π u) + 2cos (4π v).

The filter transfer functions of the SRM are adjustable by changing the values of the regularization parameters. For arbitrary λ_1/λ_0 and λ_2/λ_0, when λ_0 tends to 0, these filter transfer functions satisfy

    {[ H_0( ±1/πarctan( (λ _1/λ _0)^1/4),0) = H_1( ±1/πarctan( (λ _1/λ _0)^1/4),0) = 0.5; H_0(0, ±1/πarctan( (λ _2/λ _0)^1/4)) = H_2(0, ±1/πarctan( (λ _2/λ _0)^1/4)) = 0.5 ]..

The second, third, and fourth rows of Fig. <ref> show the filter transfer functions of the SRM and their FWHMs with the regularization parameters chosen as "λ_0=0.001, λ_1=0.001, λ_2=0.001",  "λ_0=0.001, λ_1=0.0407, λ_2=0.0204", and "λ_0=0.001, λ_1=0.0407, λ_2=0.0285", respectively. It can be seen that the FWHMs along u=0 and v=0 in these examples are well consistent with Eq. (<ref>).


In the perspective of sampling and reconstruction, the reconstruction errors come from the pre-aliasing occurring in the sampling process and the post-aliasing occurring in the reconstruction process. The pre-aliasing refers to the overlap of the frequency components of s_0, (s_1+s_2)/2, and (s_1-s_2)/2. Sufficient sampling is a necessary prerequisite to apply the DoFP technique, otherwise the strong pre-aliasing will always cause large reconstruction errors, which means the failure of the measurement. It has been show in Ref. <cit.> that if the sampling can make the frequency components of s_0, (s_1+s_2)/2, and (s_1-s_2)/2 satisfy a band-limit condition, theoretically the Stokes parameters can be perfectly reconstructed without errors. Practically, the pre-aliasing is hard to avoid, while for most polarization imaging targets, the main energy of s_0, (s_1+s_2)/2, and (s_1-s_2)/2 to be measured are concentrated in their low-frequency components, and thus the pre-aliasing reconstruction errors are usually small.

The post-aliasing refers to the misscategorization of the frequency components of s_0, (s_1+s_2)/2, and (s_1-s_2)/2. For the nearest-neighbor interpolation-based method, its filter transfer functions cannot well attenuate the undesired frequency components, and therefore usually result in large post-aliasing reconstruction errors. For the OLSM and the bilinear, bicubic, and natural bicubic spline interpolation-based methods, their filter transfer functions indicate that only when the bandwidths occupied by the frequency components of s_0, (s_1+s_2)/2, and (s_1-s_2)/2 do not exceed 0.5 cycles per pixel, can the Stokes parameters be reconstructed with less post-aliasing reconstruction errors. For the SRM and filtering-based methods, their filter transfer functions are adjustable by changing regularization parameters and window function parameters, respectively. Considering that the practical bandwidths occupied by the frequency components of s_0, (s_1+s_2)/2, and (s_1-s_2)/2 are different for different polarization imaging targets, the adjustable filter transfer functions make these methods have more powerful abilities to reduce the post-aliasing reconstruction errors. Comparing the filter transfer functions of the SLM and the filter transfer functions constructed by the window functions used in previous studies, a noticeable difference is that the frequency responses of the filter transfer functions of the SRM for the frequency components of s_0, (s_1+s_2)/2, and (s_1-s_2)/2 are anisotropic. It can be seen from Fig. <ref>(a) that the probabilities that the pre-aliasing occurs along different directions are different. This indicates that the anisotropic frequency responses can better reduce the reconstruction errors.


 §.§ Numerical simulations

Figure. <ref>(b) gives an example of the spectrum of i where the Stokes parameters are modulated by the non-uniform modulation parameters (see Section <ref>). There are still noticeable carrier peaks at the frequency coordinates (0,0), (±1/2,0), and (0,±1/2). However, due to the existence of the non-uniformity, the spectrum is disturbed. In this case, if the ideal modulation parameters are still used in the reconstruction, according to the Fourier analysis, the reconstruction results are inevitably influenced by the non-uniformity. However, combining with the non-uniform modulation parameters, the OLSM and SRM are able to construct the filters with space-varying filter kernels to deal with the non-uniformity.


Since the non-uniform modulation parameters are no longer periodic, we used numerical simulations to evaluate the reconstruction errors of the OLSM and SRM for the different frequency components of s_0, (s_1+s_2)/2, and (s_1-s_2)/2 in the non-uniform case. Firstly, we generated the cosine patterns of s_0, (s_1+s_2)/2, and (s_1-s_2)/2 by setting "α=1, β=0, γ=0", "α=0, β=√(2)/2, γ=0", and "α=0, β=0, γ=√(2)/2", respectively, in the following equation

    {[             s_0(x,y) = 1 + αcos (2π u'x + 2π v'y); 1/2[s_1(x,y) + s_2(x,y)] = βcos (2π u'x + 2π v'y); 1/2[s_1(x,y) - s_2(x,y)] = γcos (2π u'x + 2π v'y) ]..

Here, u' and v' are the horizontal and vertical frequencies varying from 0 to 0.5 cycles per pixel. The choices of α, β, and γ make the generated cosine patterns satisfy the physical constraint 0 ≤√(s_1^2 + s_2^2)/s_0≤ 1. Secondly, we substituted the generated cosine patterns and the non-uniform modulation parameters into Eq. (<ref>) to generate the DoFP images. Lastly, we used the OLSM and SRM combined with the ideal and non-uniform modulation parameters to reconstruct the Stokes parameters from the generated DoFP images. The regularization parameters used in the numerical simulations were chosen as "λ_0=0.001, λ_1=0.0407, λ_2=0.0285" as example. The first and second rows of Fig. <ref> show the root mean square errors (RMSEs) of the cosine patterns reconstructed by the OLSM combined with the ideal and non-uniform modulation parameters, respectively, the third and fourth rows of Fig. <ref> show the RMSEs of the cosine patterns reconstructed by the SRM combined with the ideal and non-uniform modulation parameters, respectively, and the last row of Fig. <ref> shows the comparisons of the above RMSEs along v'=0.



For the both methods, it can be seen that the high-frequency components of s_0, (s_1+s_2)/2, and (s_1-s_2)/2 tend to have large reconstruction errors. Therefore, it is still necessary to ensure sufficient sampling to reduce the energy in the high-frequency components. When the ideal modulation parameters are used, the low-frequency components of s_0, (s_1+s_2)/2, and (s_1-s_2)/2 have noticeable reconstruction errors, but when the non-uniform modulation parameters are used, the reconstruction errors for the low-frequency components of s_0, (s_1+s_2)/2, and (s_1-s_2)/2 are reduced, and the RMSEs reach zero for the zero-frequency components of s_0, (s_1+s_2)/2, and (s_1-s_2)/2. This indicates that when the non-uniformity is correctly characterized, the space-varying filter kernels constructed by the OLSM and SRM can reduce the reconstruction errors caused by the non-uniformity. It is worth pointing out that when the non-uniform modulation parameters are used, the RMSEs in Fig. <ref> show exactly opposite trends compared with the frequency responses of the filter transfer functions in Fig. <ref>. This phenomenon indicates that the overall frequency responses of the space-varying filter kernels for the frequency components of s_0, (s_1+s_2)/2, and (s_1-s_2)/2 are consistent with that of the filter transfer functions. The reason for this phenomenon is that the non-uniform modulation parameters are still close to the ideal modulation parameters in average sense.



§ CHOICE OF THE REGULARIZATION PARAMETERS

In this section, we give some rules of thumb for choosing suitable regularization parameters in the SRM. We represente the regularization parameters as λ_0, λ_1/λ_0, and λ_2/λ_0 to discuss their choices.

Keeping λ_1/λ_0 and λ_2/λ_0 invariant, λ_0 determines the overall weight of the regularization terms. The major factor influencing the choice of λ_0 is the level of the additional Gaussian noise. However, in this paper, the DoFP images are considered to be noiseless, so λ_0 needs to be set to a small value so that the reconstruction results can meet the constraint of Eq. (<ref>). Empirically, λ_0 is chosen as 0.001, smaller values will not cause significant changes of the reconstruction results.




λ_1/λ_0 determines the relative smoothness between the reconstructed s_0 and (s_1+s_2)/2. According to the evaluations in Section <ref>, we need to choose a suitable value of λ_1/λ_0 so that the filters constructed by the SRM can well categorize the frequency components of s_0 and (s_1+s_2)/2. Our strategy for choosing λ_1/λ_0 is illustrated in Fig. <ref>. As shown in Fig. <ref>(a), to reduce the non-uniformity disturbances, the spectrum of i/m_0 is used for choosing λ_1/λ_0. By sliding the rectangular window whose center is located at the horizontal line v=0 and averaging the spectrum in the rectangular window, a spectrum curve can be obtained. Empirically, we choose the size of the rectangular window as 0.02×0.4. The frequency coordinate of the lowest point u_p in the spectrum curve represents our estimate of the relative bandwidth occupied by the frequency components of s_0 and (s_1+s_2)/2 in the horizontal direction. As illustrated in Fig. <ref>(b), we determine λ_1/λ_0 through the formula

    λ _1/λ _0 = [ tan (πu_p)]^4.

so that the overall frequency responses of the filters for reconstructing s_0 and (s_1+s_2)/2 can be approximately equal at the frequency coordinate (± u_p,0).
Likewise, λ_2/λ_0 can be chosen by a similar strategy. In practice, inappropriate choices of λ_1/λ_0 and λ_2/λ_0 usually result in significantly serrated artifacts in the reconstructed results. If the above strategy fails, the regularization parameters can still be suitably chosen through visual inspection.





§ EXPERIMENTS

Two experiments were performed to evaluate and compare the performance of the OLSM, SRM, bicubic interpolation-based, Newton's polynomial interpolation-based method <cit.>, and Planck-Taper window filtering-based methods. The polarization imaging target in the experiments is a model car. A DoT polarimeter and the self-developed DoFP polarimeter were used to capture the polarization information. Three measures were implemented to ensure sufficient sampling, including making the target occupy a large field of view, defocusing the target slightly, and reducing the lens aperture. 
The images recorded by the polarimeters were captured 100 times and averaged to decouple the effects of noise.



The DoT polarimeter is composed by a monochrome CCD camera (FLIR GS3-U3-15S5M-C) and a linear polarization filter (THORLABS WP25M-VIS) placed on a motorized rotating platform (THORLABS PRM1/MZ8). The DoT polarimeter captured 0^∘, 45^∘, 90^∘, and 135^∘ polarization images. These polarization images were used to calculate the full-resolution Stokes parameters, DoLP, and AoP. The results are regarded as the true values. By re-sampling the four polarization images, a DoFP image was synthesized. Benefiting from the high positioning accuracy of the rotating platform and the uniform performance of the camera and linear polarization filter, the Stokes parameters in the synthesized DoFP image can be regarded as being modulated by the ideal modulation parameters. The spectrum of the synthesized DoFP image is given in Fig. <ref>(a). It can be seen that the frequency components of s_0 occupy more sized bandwidth than these of (s_1+s_2)/2 and (s_1-s_2)/2. As a matter of fact, this is a common phenomenon which has been shown in Ref. <cit.>. The OLSM, bicubic interpolation-based method, Newton's polynomial interpolation-based method, Planck-Taper window filtering-based method, and SRM were used to reconstruct the Stokes parameters, DoLP, and AoP from the synthesized DoFP image. The window function parameters of the Planck-Taper window functions are given in <ref>. Since there is no good method for choosing suitable window function parameters, the best choice minimizing the RMSEs of ŝ_0, (ŝ_1+ ŝ_2)/2, and (ŝ_1-ŝ_2)/2 was used. The regularization parameters used in the SRM were chosen as "λ_0=0.001, λ_1=0.0407, λ_2=0.0204" according to the discussions given in Section <ref>. Figure. <ref> shows the true and reconstructed s_0, DoLP, and AoP. s_0 ranges from 0 to 255, while the range of the grayscale bar in Fig. <ref> is set to [0,160] to increase the image contrast. The most noticeable differences are marked in the white boxes. s_0 reconstructed by the OLSM and bicubic interpolation-based method are blurred, and the reconstructed DoLP and AoP show obviously serrated artifacts. The main reason for these phenomena is that the frequency components of s_0 that more than ± 0.25 cycles per pixel are incorrectly categorized. The Newton's polynomial interpolation-based method is an edge-preserved reconstruction method. It can be seen that the reconstructed s_0 shows more details, and the reconstructed DoLP image avoids the serrated artifacts. However, due to the wrong edge discriminations, the reconstructed AoP image has "×" artifacts. The reconstructed results of the Planck-Taper window filtering-based method and SRM all show good visual effects. This benefits from the suitably choices of the window function parameters and regularization parameters.
Table <ref> gives the RMSEs and normalized RMSEs of the reconstructed results. The normalized RMSEs are defined as the ratio of the RMSE to the difference between the maximum and minimum value of the true value. All the reconstruction methods can achieve relatively small reconstruction errors, which illustrates the effectiveness of the DoFP technique. Comparing the relative reconstruction errors of the different reconstruction items, it can be seen that the relative reconstruction errors of the AoP images are larger than that of the others. This is consistent with the results given in Ref. <cit.>. The SRM shows minimal reconstruction errors for most reconstruction results. Meanwhile, due to the anisotropic frequency responses, the SRM shows smaller reconstruction errors than the Planck-Taper window filtering-based method for all the reconstruction results.


     



The modulation parameters of the self-developed DoFP polarimeter are shown in Fig. <ref>. The spectrum of the DoFP image captured by this polarimeter is shown in Fig. <ref>(b), which is disturbed by the non-uniformity. The first three rows of Fig. <ref> show the s_0, DoLP, and AoP reconstructed by the bicubic interpolation-based method, Newton's polynomial interpolation-based method, Planck-Taper window filtering-based method. The window function parameters used in the reconstruction were the same as the previous. The reconstruction results of these three methods contain strong non-uniformity artifacts, which indicates the interpolation-based and filtering-based methods fail in the presence of the non-uniformity. The fourth and fifth rows of Fig. <ref> show the s_0, DoLP, and AoP reconstructed by the OLSM and SRM, respectively, with the non-uniform modulation parameters used in the reconstructions. The regularization parameters were chosen as "λ_0=0.001, λ_1=0.0407, λ_2=0.0285" according to the discussions given in Section <ref>. Since the non-uniformity has been effectively characterized by the non-uniform modulation, the non-uniformity artifacts in the reconstruction results of these two methods are well mitigated. This proves the effectiveness of our methods for reconstructing the Stokes parameters in the presence of the non-uniformity. Comparing the reconstruction results of the OLSM and SRM in Fig. <ref> and Fig. <ref>, it can be seen that the performance of the OLSM and SRM in  these two experiments is basically consistent. In Fig. <ref>, s_0 reconstructed by the OLSM is still blurred, and the reconstructed DoLP and AoP also show obviously serrated artifacts. While the s_0 reconstructed by the SRM shows more details, and the reconstructed DoLP and AoP show better visual effects. The reason for these phenomena is that when the non-uniformity has been considered in the reconstruction, the frequency responses of the two methods for the different frequency components of the Stokes parameters under the non-uniform modulation is basically the same as that under the ideal modulation. Notice that although the OLSM and SRM tackle the non-uniformity based on the same pixel model, the results are different at defect pixels. This phenomenon is marked by the white ellipses in Fig. <ref>. It can be seen that the reconstruction results of the OLSM show more steep changes in flat background compared with that of the SRM. The performance of the defects pixels far deviates from the design performance. Since two adjacent defect pixels can easily cause the degeneration of the 4×3 matrix in Eq.(<ref>), the OLSM is more likely to fail at defect pixels. In contrast, the SRM is more robust to the defect pixels, this is due to the use of the global measurement information. However, empirically, we found that it is still difficult for the SRM to obtain satisfactory reconstruction results when the size of the defect area exceeds that of the biharmonic operator. 




§ CONCLUSION

We proposed the OLSM and SRM to reconstruct the Stokes parameters. The performance of the OLSM and SRM was investigated through Fourier analysis, numerical simulations, and experiments. The proposed methods can effectively mitigate the reconstruction errors and artifacts caused by the non-uniformity. The OLSM reconstructs the Stokes parameters under the local constant assumption. This method can be regraded as a generalization of the bilinear interpolation-based methods. Since the pseudo-inverse of the matrix in Eq. <ref> can be pre-calculated, the OLSM consumes less computational costs, and is therefore more suitable for hardware implementation. The performance of the OLSM can be further improved if an optical low-passed filter is applied to restrict the bandwidths occupied by the frequency components of the Stokes parameters. The SRM reconstructs the Stokes parameters under the global smoothing assumption. This method has more flexible filtering performance. With a suitable choice of the regularization parameters, the SRM can reconstruct the Stokes parameters with good visual effects and low reconstruction errors. Some edge-preserved constraints can be applied to further improve the performance of the SRM.




§ SCALING AND ROTATION TRANSFORMATION

Let m_0', m_1', and m_2' be the original modulation parameters obtained from the calibration, the modulation parameters used in this paper are obtained by applying a scaling and rotation transformation to m_0', m_1', and m_2', expressed as

    [ m_0(x,y); m_1(x,y); m_2(x,y) ]
    = 1/m[           1           0           0;           0   cos (2δ ) - sin (2δ );           0   sin (2δ )   cos (2δ ) ][ m_0'(x,y); m_1'(x,y); m_2'(x,y) ].

Here, m is the mean of m_0', and δ is the mean of the differences between the designed and the actual orientations of the linear polarization filters. Notice that the differences between the designed and the actual orientations should be wrapped between -π/2 and π/2. The purpose of this transformation is to normalize the modulation parameters and choose a relative orientation in the linear polarization filter array itself as the reference of 0^∘. The changes caused by this transformation can be compensated by applying the inverse of this transformation to the reconstructed Stokes parameters.





§ EXPLANATION OF NUMERICAL CALIBRATION

Numerical calibration methods are proposed to correct the non-uniform intensity responses in the DoFP image to ideal. A commonly used numerical calibration method is Gruev et. al.'s super-pixel calibration method, which can be expressed as

    [ [         î(2p,2q);     î(2p,2q + 1);     î(2p + 1,2q); î(2p + 1,2q + 1) ] ]
    = 
    
    [ [  1  1  0;  1  0  1;  1 -1  0;  1  0 -1 ] ]
    [ [         m_0(2p,2q)         m_1(2p,2q)         m_2(2p,2q);     m_0(2p,2q + 1)     m_1(2p,2q + 1)     m_2(2p,2q + 1);     m_0(2p + 1,2q)     m_1(2p + 1,2q)     m_2(2p + 1,2q); m_0(2p + 1,2q + 1) m_1(2p + 1,2q + 1) m_2(2p + 1,2q + 1) ] ]
    	^†[ [         i(2p,2q);     i(2p,2q + 1);     i(2p + 1,2q); i(2p + 1,2q + 1) ] ]_Stokes parameters reconstruction

Here, î represents the estimate of the ideal DoFP image, p and q are the coordinates which take half the ranges of x and y, respectively. The super-pixel calibration method can be divided into two steps. The first step is also to reconstruct the Stokes parameters in the presence of the non-uniformity. The second step is to multiply the reconstructed results by a matrix defined by the ideal modulation parameters to generate î. The Stokes parameters reconstruction in the super-pixel calibration method is similar to Eq. <ref>. This indicates that the super-pixel calibration methods will definitely cause additional correction errors due to the spatial variations of the Stokes parameters.

Some studies suggest to first applying the super-pixel calibration methods to correct the non-uniformity, and then use the interpolation-based or filtering-based methods to reconstruct the Stokes parameters. However, we believe that this suggestion has potential problems. Firstly, the whole reconstruction process will be redundant, since the Stokes parameters are actually reconstructed twice. Secondly, with the advancement of the manufacturing technique, the actual modulation parameters of the DoFP polarimeter will be close to the ideal modulation parameters, it is foreseeable that the correction errors introduced by the super-pixel calibration method will exceed the non-uniformity errors mitigated by this method. In fact, the OLSM and SRM can also be used as the numerical calibration methods, achieved by substituting the ideal modulation parameters and the reconstructed Stokes parameters into Eq. <ref> to generate î. But based on the above concerns, we do not suggest to use the OLSM and SRM as the numerical calibration methods.





§ CONVOLUTIONAL INTERPOLATION-BASED METHODS





The Stokes parameters reconstructed by the convolutional interpolation-based methods can be expressed as

    {[        ŝ_0 = i * 1/4h; ŝ_1 = (i ∘m_1) * 1/2h; ŝ_2 = (i ∘m_2) * 1/2h ]..

Here, ∘ represents the Hadamard product, m_1 and m_2 take the ideal values given in Eq. (<ref>), h represents the interpolation kernel. The interpolation kernels of the nearest-neighbor, bilinear, bicubic, and natural bicubic spline interpolation algorithms are given in Table. <ref>. The natural bicubic spline interpolation algorithm is a piece-wise interpolation algorithm, therefore the interpolation kernel can have different sizes. The 7×7 and 11×11 sized interpolation kernels are given as example.



Performing the DFT of Eq. (<ref>), after some simplifications, we have

    {[                               Ŝ_0(u,v) = I(u,v) ·1/4H(u,v); 12[Ŝ_1(u + 12,v) + Ŝ_2(u + 12,v)] = I(u,v) ·1/4H(u + 12,v); 12[Ŝ_1(u,v + 12) - Ŝ_2(u,v + 12)] = I(u,v) ·1/4H(u,v + 12) ]..

The DFTs of the interpolation kernels of the nearest-neighbor, bilinear, bicubic, and natural bicubic spline interpolation algorithms and their amplitudes along v=0 are plotted in Fig. <ref>. It can be found that the interpolation kernels of the bilinear, bicubic, and natural bicubic spline interpolation algorithms are the approximations of sinc(x/2)sinc(y/2) in the spatial domain, and the DFTs are the approximations of 4rect(2u)rect(2v) in the frequency domain.





§ PLANCK-TAPER WINDOW FUNCTION



The Planck-Taper window function is defined as <cit.>

    H(r) = {[                                            1,                                   r ≤ l - w/2; 1/1 + exp(  - w/r - l + w/2 - w/r - l - w/2),                         l - w/2 < r < l + w/2;                                            0,                                   r ≥ l + w/2 ] ..

Here, r=√(u^2+v^2), l is the radius of the FWHM, w is the falloff range. The filter transfer functions constructed by the Planck-Taper window function used in Section <ref> are illustrated in Fig. <ref>.






§ ACKNOWLEDGEMENTS

This work was supported by the National Natural Science Foundation of China (Grant Nos. 11627803, 11872354, 11872355) and the Strategic Priority Research Program of the Chinese Academy of Sciences (XDB22040502).


