Three-dimensional polarimetric integral imaging in photon-starved conditions: performance comparison between visible and long wave infrared imaging.

Three-dimensional (3D) polarimetric integral imaging (InIm) to extract the 3D polarimetric information of objects in photon-starved conditions is investigated using a low noise visible range camera and a long wave infrared (LWIR) range camera, and the performance between the two sensors is compared. Stokes polarization parameters and degree of polarization (DoP) are calculated to extract the polarimetric information of the 3D scene while integral imaging reconstruction provides depth information and improves the performance of low-light imaging tasks. An LWIR wire grid polarizer and a linear polarizer film are used as polarimetric objects for the LWIR range and visible range cameras, respectively. To account for a limited number of photons per pixel using the visible range camera in low light conditions, we apply a mathematical restoration model at each elemental image of visible camera to enhance the signal. We show that the low noise visible range camera may outperform the LWIR camera in detection of polarimetric objects under low illumination conditions. Our experiments indicate that for 3D polarimetric measurements under photon-starved conditions, visible range sensing may produce a signal-to-noise ratio (SNR) that is not lower than the LWIR range sensing. We derive the probability density function (PDF) of the 2D and 3D degree of polarization (DoP) images and show that the theoretical model demonstrates agreement to that of the experimentally obtained results. To the best of our knowledge, this is the first report comparing the polarimetric imaging performance between visible range and infrared (IR) range sensors under photon-starved conditions and the relevant statistical models of 3D polarimetric integral imaging.

reflected beam can be described in terms of Stokes parameters and degree of polarization (DoP) [4]. Four quantities of Stokes parameters (S 0 , S 1 , S 2 and S 3 ) are required to measure the DoP. These quantities can be recorded by placing polarizing elements in front of an image sensor, and recording images with different orientations of the polarizing elements. In low light conditions, long wave infrared (LWIR) imaging sensors are commonly employed to obtain the polarimetric information of objects [5][6][7]. Employment of infrared (IR) imaging in low light is suitable for DoP calculation; however, the IR imaging systems can be bulky and expensive in comparison with visible range imaging systems. The measurement of DoP using visible cameras could be an effective technique in photon-starved conditions but the dominance of camera noise in photonstarved conditions can degrade the accuracy of DoP calculations. Passive three-dimensional (3D) polarimetric integral imaging (InIm) [8][9][10][11] is one of the prominent techniques to overcome the aforementioned issues with visible cameras operating in photon-starved conditions. 3D InIm [12][13][14][15][16][17][18][19][20][21][22] is performed by recording two-dimensional (2D) images from different perspectives, then reconstructing the three-dimensional (3D) scene through either optical or computational reconstruction methods. Different perspective images known as elemental images can be captured by a single image sensor with a lenslet or camera array, or by a single moving image sensor [12][13][14][15][16][17][18][19][20][21][22]. The computational reconstruction process is performed by the back projection of the optical rays through a virtual pinhole to the desired reconstruction distance [23].
In this paper, we compare a visible range polarimetric 3D integral imaging system with an LWIR polarimetric 3D integral imaging system in low light conditions in terms of signal-to-noise ratio (SNR) and polarimetric object detection. Passive 3D polarimetric InIm [8][9][10][11] is used to extract the polarimetric information of visible and LWIR polarimetric objects in low light conditions. Since the captured visible range images are significantly degraded by read noise in the photon-starved conditions, it becomes challenging to calculate the Stokes parameters and the corresponding DoP. In order to extract the polarimetric information for a visible range camera in these conditions, a mathematical model for signal restoration is applied to each elemental image [24][25][26]. Moreover, the camera read noise is reduced during 3D InIm reconstruction of polarimetric images, which is optimal in the maximum likelihood sense [27]. Furthermore, the total variation (TV) denoising [28] algorithm is used for noise reduction of DoP images for the visible range InIm systems. The mathematical model for signal restoration and TV denoising are not applied in case of the LWIR imaging system because enough thermal photons are present in the scene.
Our experiments for polarimetric 3D integral imaging under low light conditions indicate that the visible range imaging system produces an SNR that may be higher than the LWIR imaging system. In our experiments, we also find that a small polarimetric object may be detectable using the visible range 3D integral imaging system but is barely detectable by the LWIR 3D integral imaging system. Furthermore, we derive the theoretical probability density functions (PDF) for the 2D and 3D DoP images and show strong similarity between the theoretically derived distributions and the probability distribution functions found through experimental measurements.

Theoretical background of degree of polarization (DoP)
The Stokes parameters provide quantitative measurements of polarimetric imaging [4]. Let us define the expression for quasi-monochromatic light propagating along the z direction as where x and y are orthogonal unit vectors, E 0x and E 0y are the amplitude of electromagnetic components in the x and y directions, respectively, and φ x and φ y are the corresponding phase terms. The Stokes parameters for quasi-monochromatic light are defined as: where . represents time averaging and φ = φ y -φ x is the phase difference between the two orthogonal components of the electric field. S i (i=0,1,2,3) denotes the Stokes parameters. I θ is the intensity of polarized light recorded when the linear polarizer in front of the imaging sensor is placed at an angle of θ with respect to the x-axis, and I θ, π/2 is the intensity recorded after inserting the quarter wave plate (QWP) in addition to the linear polarizer. The degree of polarization (DoP) and degree of linear polarization (DoLP) in terms of Stokes parameters are defined as: The value of DoP ranges from 0 to 1, with 1 representing completely polarized light and 0 representing completely unpolarized light. Generally, in passive polarimetric imaging applications, circularly polarized light S 3 is very small and rarely measurable. As such, in these experiments, the value of S 3 is taken to be zero. When the value of S 3 is assumed as zero, the equations for DoLP and DoP become interchangeable. Therefore, in the following sections, all analysis is carried out for the DoP values with S 3 set to zero, and only four polarimetric images [I 0°, I 45°, I 90°, I 135°] are required for each perspective.

Experimental methods
The 3D scene consists of a linear polarizer film, a mannequin, and test tubes filled with hot water to provide an IR source, as shown by a reference image taken in high illumination conditions using the low noise visible range camera (Hamamatsu C11440-42U) by Fig. 1(a). Figure 1(b) shows the same 3D scene imaged using the visible range camera in low light illumination conditions. The estimated number of photons per pixel under low illumination conditions was calculated as 1.23 photons/pixel [10,29]. For the LWIR imaging, the same 3D scene is used, however the linear polarizer is replaced by an IR wire grid polarizer in the wavelength range of 7-15µm. The 3D scene is imaged using the LWIR camera (Tamarisk 320 LWIR camera 60 Hz) as shown in Fig. 1(c). Note, the ambient illumination in Fig. 1(c) is the same as that of Fig. 1(b) but the LWIR camera is not affected by the scene illumination. In this paper, the polarimetric visible and LWIR images were captured using synthetic aperture integral imaging (SAII) [23] with the visible range camera and the LWIR range camera in photon starved conditions. SAII consists of a camera on a moving platform [9][10] as shown in Fig. 2(a). Also shown in Fig. 2(a), a polarizer is placed in front of the imaging sensor for polarimetric imaging. The SAII pickup process of a 3D scene is shown in Fig. 2(b) and reconstruction of the 3D scene using the SAII algorithm is shown in Fig. 2(c). To provide a fair comparison of DoP between the visible and LWIR cameras, the pixel size of the visible camera is binned with 3 by 3 binning to 19.5×19.5 µm in order to be comparable with the LWIR camera which has a pixel size 17×17 µm. Both cameras use nearly the same f-number to capture the polarimetric object.  In 3D polarimetric imaging reconstruction process, the 3D image I z θ (x, y) at the reconstruction plane z from the camera system is expressed as [23]: In this equation, M and N are the number of polarimetric elemental images in the x-and ydirections and O(x, y) is the overlapping pixel number at (x, y). p x and p y are the camera pitch size in horizontal (H) and vertical (V) directions respectively. L x and L y are the total number of pixels in each column and row on images, c x × c y is the sensor size of camera. The focal length of camera lens is f, and ε is the additive camera noise. I θ m,n is the ideal 2D polarimetric elemental image, θ represents the different polarization states of the image while the subscripts m, n represent the location of elemental image. For 3D integral imaging in these experiments, a total of 25 perspective images [5(H) × 5(V)] with a camera pitch of 30 mm in both directions are recorded. Table 1 summarizes the values of different integral imaging and camera parameters used in the experiments. When using a visible range sensor in low light illumination conditions, very few photons are reflected from the scene and may be attenuated by the atmospheric particles in the environment, resulting in captured images that are dominated by camera read noise. Therefore, the captured Elemental images images using a visible camera in photon-starved conditions become degraded and sparse. In photon-starved conditions, the electron counts are very small due to the low number of photons reflected from the object scene. To alleviate the problem of camera noise, we apply a mathematical signal restoration model. Before applying the signal restoration model, first we must subtract the camera offset from the captured elemental images. The camera offset is required to prevent the clipping of very small signals during image capture that would otherwise occur because some signals, due to noise, may be less than zero when digitized [30]. This offset can be measured by recording of large number of single camera bias reference frames, which are obtained by setting the exposure time of camera to be a minimum (3 ms) and closing the aperture to avoid the atmospheric photons [10,30]. After subtraction of the camera offset image, a mathematical model for signal restoration in low light conditions is applied to restore the visibility of the 2D elemental images in photon-starved conditions. Dark channel prior-based mathematical models for haze removal are widely used for image restoration and enhancement [24]. Additionally, it has been reported that the dark channel prior method is very useful for low light image enhancement [25,26]. The dark channel prior method was derived from the statistics of outdoor images. From the set of outdoor images, it was shown that in most of sky-free regions of the images, some pixels (defined as dark pixels) have very low intensity in at least one of the color (RGB) channels. The process of selecting the minimum intensity pixels from an image serves as a prior, and the method of image enhancement is referred to as the dark channel prior method [24]. The prior information of the dark pixels can provide an accurate estimation of transmission function caused by the haze. This dehazing model has been used to enhance the image in low light illumination conditions since inversion of the low light image results in an image that appears very similar to a hazy image [25,26]. Therefore, in this work, we apply the dark channel prior method to the inverted low light images to enhance the image in photon starved conditions.
Since we are using the Hamamatsu C11440-42U sensor, which is a monochromatic image sensor, the monochromatic version of dark channel prior method is used to recover the image. Instead of performing the optimization over all the channels, in the monochromatic version of dark channel prior method, the minimization of global cost function is performed over the single channel. For the mathematical model of signal restoration, let us consider I(ξ) as the image to be recorded in photon starved condition where ξ is the vector representation of the pixel location (x, y). To facilitate the signal recovery using the monochromatic version of dark channel prior mathematical restoration model, we invert the observed image intensity I(ξ) using R(ξ) = 255-I(ξ). The inverted image R(ξ) is similar to the hazy image and fed into the widely used signal restoration model [24][25][26]: where R(ξ) is the recorded intensity of the inverted input image, J(ξ) is the ideal image to be recovered, A is the atmospheric light, and t(ξ) is portion of light that is not scattered and reaches the camera. Using the monochromatic version of dark channel prior, A and t(ξ) can be estimated. The transmission t(ξ) value ranges from 0 to 1, for complete dark regions, t(ξ) becomes 0 and for brightest regions, t(ξ) becomes 1. We have chosen the top 0.1% of brightest pixels in R dark . The highest pixel value among the brightest pixels is selected as estimated atmospheric light A . The final image can be recovered with improved quality in comparison to the original inverted image R(ξ) by estimating J(ξ) using the above method. Recovered image J(ξ) is then inverted back to the final enhanced image S(ξ) by S(ξ) = 255-J(ξ).
The noise in the processed images are further reduced using TV denoising algorithm [28]. The TV denoising algorithm is applied to the four polarimetric 2D images [I 0°, I 45°, I 90°, I 135°] in the 2D case, and after 3D reconstruction of polarimetric images [I 0 z°, I 45 z°, I 90 z°, I 135 z°] , before calculating the 3D DoP image, in the 3D case.

DoP calculation in low light using visible and LWIR range cameras
For the experiments conducted using the visible range sensor, an average of 1.23 photons per pixel were estimated prior to pixel binning [10,29]. The captured visible 2D elemental images are used to calculate the 2D Stokes parameters and 2D DoP using Eq. (2) and Eq. (4). When the number of photons per pixel is low, the noise present in the image is enhanced during the DoP calculation due to the nonlinear combination of the Stokes parameters. In this case, the polarimetric information in 2D image is very noisy and produces erroneous results [see Fig. 3(a)]. The reconstructed 3D integral image is obtained from the 2D elemental images using Eq. (5) and the 3D Stokes parameters and the 3D DoP image are calculated using Eq. (2) and Eq. (4), respectively. The amount of noise in the 3D DoP image [see Fig. 3(b)] is reduced in comparison to the 2D DoP image; but it is still noisy. Further reduction of noise in the DoP images can be achieved using dark channel modelling along with total variation (TV) denoising algorithm. To further demonstrate the proposed method of dark channel modeling along with TV denoising in comparison to prior arts, we compare the 2D DoP and 3D DoP images with no processing, using only dark channel modeling, using only TV denoising and using the combination of dark channel modeling and TV denoising (i.e. the proposed method). Figure  For quantitative comparison of each part of the image processing pipeline, the signal-to-noise ratio (SNR) of polarimetric (signal) and non-polarimetric (background) areas is measured. The

SNR is defined as SNR
where the mean of the polarimetric signal is µ s , and the mean of the non-polarimetric area (background) is µ b . σ s 2 and σ b 2 are the variances of the polarimetric and non-polarimetric areas, respectively [10,29]. For fair comparison of the visible range DoP images using different image processing pipelines, the same image regions were considered within the red (signal) and yellow (background) windows selected from DoP images as shown in Fig. 3(b). The dark channel modeling alone provides a small improvement in the SNR of the visible DoP image. The SNR of the visible DoP images is significantly improved when using the TV denoising algorithm. Further improvement of noise reduction can be achieved by using dark channel modeling along with TV denoising. Therefore, all the experimental results of visible DoP images in photon-starved condition are processed using both dark channel modeling and TV denoising. The results of the SNR comparison for different processing pipelines are shown in Table 2.  The process of polarimetric imaging with the LWIR range camera is the same as that of the visible range camera except that the mathematical restoration model and TV denoising are not needed for the LWIR sensor because there are sufficient thermal photons in the scene. Visible polarimetric images are corrupted by camera read noise which has Gaussian distribution in low light conditions. However, the polarimetric images in case of LWIR range were recorded in sufficient thermal photon conditions. Therefore, TV denoising in addition of dark channel were applied only on the visible data to remove the Gaussian noise in photon-starved conditions. Moreover, after applying the denoising algorithms in the case of the LWIR images, there may be a reduction in the SNR of the DoP images. The 2D elemental images recorded with an LWIR camera are used to reconstruct 3D Integral Images using Eq. (5). The Stokes parameters and DoP of LWIR range camera are calculated using Eq. (2) and Eq. (4), respectively. For the scene in Fig. 1(a), the 2D and 3D DoP images in the LWIR range are illustrated in Fig. 4  . For quantitative comparison of each part of the image-processing pipeline using the LWIR sensor, the SNR between the polarimetric (signal) and non-polarimetric (background) areas is measured. From the results shown in Fig. 4 and Table 3, it is evident that there is no benefit of applying the dark channel prior method and TV denoising to the LWIR data. The highest SNR is achieved for the 3D DOP image with no additional processing.

Summary of comparison between 3D visible and 3D LWIR polarimetric imaging systems
For an imaging system with a fixed f -number and a fixed object distance, the resolution of imaging systems is determined by the pixel size and wavelength of light source [31]. Since the pixel size and f -number of our imaging systems were nearly the same, the resolutions of our systems are expected to differ only based on the wavelength of light source. The LWIR range camera shows a low polarimetric signature of the object as shown by the 3D DoP image [see Fig. 4(b)]. However, the visible range camera displays a comparatively high polarimetric signature of polarimetric object as shown by the 3D DoP image [see Fig. 3(h)] using the proposed method of both dark channel modeling and total variation denoising. Moreover, the polarimetric information in the visible range 3D DoP image is comparatively enhanced, and the detailed information is recovered. We quantitatively compared the DoP images from the LWIR and visible range cameras by measuring the SNR of 3D DoP images. The maximum SNR achieved in the visible range is 23.25 of 3D DoP image [see Fig. 3(h)] using the dark channel modeling along with TV denoising. However, the maximum SNR achieved in the LWIR range is 8.24 of 3D DoP image [see Fig. 4(b)]. Therefore, in our experiments, an improvement of 182.16% was achieved in SNR using the 3D DoP image in the visible range in comparison to the 3D DoP image in the LWIR range.
In order to provide further comparison of polarization capabilities for both the visible and LWIR range cameras, a one-dimensional line plot of the DoP intensity is examined. An aperture with 3 mm diameter was placed in front of the polarimetric object as shown by a reference image taken in high illumination conditions using low noise visible range camera by Fig. 5(a). The experiments were repeated to record the polarimetric information in the visible range and the LWIR range as shown in Figs. 5(b) and 5(d), respectively. The visible range camera provides the ability to identify the polarimetric object of size 3 mm and has a higher DoP value as shown by the 1D intensity plot in Fig. 5(c). However, using the LWIR range camera, the polarimetric object of size 3 mm is barely detectable, and 1D intensity plot of the DoP value [ Fig. 5(e)] is nearly constant around zero.

Statistical analysis of 2D and 3D DoP images
Lastly, we propose a statistical approach to derive the statistical distribution of the 2D and 3D DoP images in low light conditions. Under low light conditions, the captured visible 2D image is dominated with camera read noise, which has a Gaussian distribution. To verify the distribution of 2D image, a Kolmogorov-Smirnov test [32] was applied to the 2D elemental image in Fig. 1(b). The polarimetric region was selected from 2D elemental image as a region of interest (ROI). The result shows that the test accepted the null hypothesis at significance level of 5%. Furthermore, in the following derivations for the degree of linear polarization, we assume S 3 =0 (i.e. no circular polarization). Therefore, only S i [i = 0, 1, 2] will be considered and the probability density function of polarimetric region in the DoP image for both the 2D and 3D cases can be written as [33]: where P represents the DoP of the image and the non-centrality parameters δ = 2 i=1 (µ i /σ) 2 and γ = (µ 0 /σ) 2 , where i is the index of Stokes parameter. B(.) is the beta function, j and k are the variables for summation operator. The details of this derivation and approach are described in Appendix A. In the 3D case, the mean and variance of the 3D polarimetric image are derived by the following equations: The intensity value µ z θ (x,y) at pixel (x,y) is averaged from the 2D elemental images [10]. All the parameters in Eq. (8) and Eq. (9) are the same as those in Eq. (5).
For experimental measurement of histograms, the regions of interest containing the polarizing element in the visible range 2D DoP image shown in Fig. 3(g), and the 3D DoP image shown in Fig. 3(h). The experimental histograms and computed theoretical PDFs are compared in Fig. 6. In Fig. 6(a), the red curve shows the theoretically derived PDF of the 2D DoP using Eq. (7), where µ 0 = 0.35, µ 1 = 0.158, µ 2 = 0.026 and σ 2 = 1.5×10 −3 . These values were computed from the experimental 2D polarimetric image. For the theoretically derived PDF of the 3D DoP [ Fig. 6(b)], the mean and variance of 3D reconstructed images were computed using Eq. (8) and Eq. (9) respectively, then the non-centrality parameters were computed using the mean and the square root of the variance of 3D reconstructed image. Finally, the computed theoretical PDF of the 3D DOP is compared with the experimental histogram [ Fig. 6(b)]. The experimental results agree well with the theoretically derived model. The differences may be due to the finite number of samples used to compute the histograms, and the assumption of statistical independence of elemental images in obtaining the theoretical results. The smaller standard deviation of PDF of 3D image confirms that 3D DoP has less noise compared to the 2D DoP image. Therefore, 3D polarimetric InIm can improve the DoP image quality over 2D DoP. This analysis agrees with the prior experimental results in section 4.1, which showed that the 3D DoP image may outperform the 2D DoP image in terms of SNR.

Conclusion
In this paper, passive 3D polarimetric integral imaging in the visible range and LWIR range have been compared in low light conditions. In order to extract the polarimetric information using a visible imaging sensor in photon-starved conditions, a mathematical model for signal restoration of elemental images is used. The total variation (TV) denoising is used for noise reduction of the 2D and 3D DoP images in the visible range. Quantitative comparison between polarimetric 3D integral imaging using a low noise visible band sCMOS camera and an LWIR band camera was presented. The polarimetric signature of a 3 mm-sized object was detected with the visible range camera and was barely detectable by the LWIR camera in low light conditions. The visible range camera similarly outperforms the LWIR range camera in detection of a polarimetric signature in low light conditions based on the calculated SNR for each camera. Finally, theoretically derived and experimentally measured probability distribution functions of the 2D and 3D DoP images were presented and compared. This comparison shows agreement between the experimental results and the theoretically derived models. To the best of our knowledge, this is the first report comparing the polarimetric imaging performance between visible range and IR range sensors under photon-starved conditions. Future work may examine the integration of the preprocessing steps using a joint numerical procedure rather than a piece-wise strategy with the integral imaging reconstruction algorithm for optimized performance.

Appendix A: theoretical statistical distribution of DoP image
In low light illumination conditions, the images recorded at different polarizing angles are corrupted with camera read noise, which has the Gaussian distribution. The Stokes parameters are the addition and subtraction of two orthogonal polarimetric images. Therefore, the calculated Stokes parameters are read noise dominated with Gaussian distribution with same variance and different means. The experimentally measured distribution (histogram) of Stokes parameters (S 0 , S 1 and S 2 ) in photon starved condition without any processing is shown in Fig. 7. The histograms in Fig. 7 show the Stokes parameters follow the Gaussian distribution with nearly same variance and different means.
Assuming I 0°a nd I 90°a re independent of each other, the variances of S 0 and S 1 are defined as, where σ 2 S 0 and σ 2 S 1 are the variances of Stokes parameters S 0 and S 1 , respectively and σ 2 I 0 • and σ 2 I 90 • are the variances of two orthogonal polarimetric 2D elemental images. Experimentally, the variances of Stokes parameters are found to be approximately equal (see Fig. 7).
The Stokes parameters for the polarimetric area of a scene follow the Gaussian distribution with non-zero mean S i ∼N (µ i , σ 2 ) [i = 0, 1,2].
Let Q = The probability density function of P can be obtained by substituting the value of Q in above equation and using the transformation of variable,