PCA-based denoising method for division of focal plane polarimeters

Division of focal plane (DoFP) polarimeters are composed of interlaced linear polarizers overlaid upon a focal plane array sensor. The interpolation is essential to reconstruct polarization information. However, current interpolation methods are based on the unrealistic assumption of noise-free images. Thus, it is advantageous to carry out denoising before interpolation. In this paper, we propose a principle component analysis (PCA) based denoising method, which works directly on DoFP images. Both simulated and real DoFP images are used to evaluate the denoising performance. Experimental results show that the proposed method can effectively suppress noise while preserving edges. c © 2017 Optical Society of America OCIS codes: (100.3020) Image reconstruction-restoration; (110.5405) Polarimetric imaging; (120.5410) Polarimetry; (260.5430) Polarization; (100.3190) Inverse problems. References and links 1. T. Krishna, C. Creusere, and D. Voelz,“Passive polarimetric imagery-based material classification robust to illumination source position and viewpoint,” IEEE Trans. Image Process. 20(1), 288–292 (2011). 2. S. Fang, X. Xia, X. Huo, and C. Chen, “Image dehazing using polarization effects of objects and airlight,” Opt. Express 22(16),19523–19537 (2014). 3. Z. Chen, X. Wang, and R. Liang, “Snapshot phase shift fringe projection 3D surface measurement,” Opt. Express 23(2), 667–673 (2015). 4. A. Pierangelo, A. Benali, M. R. Antonelli, T. Novikova, P. Validire, B. Gayet, and A. De Martino, “Ex-vivo characterization of human colon cancer by Mueller polarimetric imaging,” Opt. Express 19(2), 1582–1593 (2011). 5. S. Alali and A. Vitkin, “Polarized light imaging in biomedicine: emerging Mueller matrix methodologies for bulk tissue assessment,” J. Biomed. Opt. 20(6), 061104 (2015). 6. B. Kunnen, C. Macdonald, A. Doronin, S. Jacques, M. Eccles, and I. Meglinski, “Application of circularly polarized light for non-invasive diagnosis of cancerous tissues and turbid tissue-like scattering media,” J. Biophotonics 8(4), 317–323 (2015). 7. G. Myhre, W. Hsu, A. Peinado, C. LaCasse, N. Brock, R. A. Chipman, and S. Pau, “Liquid crystal polymer fullStokes division of focal plane polarimeter,” Opt. Express 20(25), 27393–27409 (2012). 8. J. S. Tyo, “Hybrid division of aperture/division of a focal-plane polarimeter for real-time polarization imagery without an instantaneous field-of-view error,” Opt. Lett. 31(20), 2984–2986 (2006). 9. 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). 10. J. S. Tyo, D. L. Goldstein, D. B. Chenault, and J. A. Shaw, “Review of passive imaging polarimetry for remote sensing applications,” Appl. Opt. 45(22), 5453–5469 (2006). 11. B. Ratliff, C. LaCasse, and S. Tyo, “Interpolation strategies for reducing IFoV artifacts in microgrid polarimeter imagery,” Opt. Express 17(11), 9112–9125 (2009). 12. S. Gao and V. Gruev, “Bilinear and bicubic interpolation methods for division of focal plane polarimeters,” Opt. Express 19(27), 26161–26173 (2011). 13. S. Gao and V. Gruev, “Image interpolation methods evaluation for division of focal plane polarimeters,” Proc. SPIE 8012, 80120N (2011). Vol. 25, No. 3 | 6 Feb 2017 | OPTICS EXPRESS 2391 #281508 https://doi.org/10.1364/OE.25.002391 Journal © 2017 Received 24 Nov 2016; revised 9 Jan 2017; accepted 9 Jan 2017; published 30 Jan 2017 14. S. Gao and V. Gruev, “Gradient-based interpolation method for division of focal plane polarimeters,” Opt. Express 21(1), 1137–1151 (2013). 15. 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–20807 (2016). 16. E. Gilboa, J. P. Cunningham, A. Nehorai, and V. Gruev, “Image interpolation and denoising for division of focal plane sensors using Gaussian processes,” Opt. Express 22(12), 15277–15291 (2014). 17. S. Chang, B. Yu, and M. Vetterli, “Spatially adaptive wavelet thresholding with context modeling for image denoising,” IEEE Trans. Image Process. 9(9), 1522–1531 (2000). 18. M. Elad and M. Aharon, “Image denoising via sparse and redundant representations over learned dictionaries,” IEEE Trans. Image Process. 15(20), 3736–3745 (2006). 19. A. Beck and M. Teboulle, “Fast Gradient-Based Algorithms for Constrained Total Variation Image Denoising and Deblurring Problems,” IEEE Trans. Image Process. 18(11), 2419–2434 (2009). 20. G. Liu, Z. Lin, S. Yan, J. Sun, Y. Yu, and Y. Ma, “Robust Recovery of Subspace Structures by Low-Rank Representation,” IEEE Trans. Pattern Anal. Machine Intel. 35(1), 171–184 (2013). 21. L. Zhang, W. Dong, D. Zhang, and G. Shi, “Two-stage image denoising by principal component analysis with local pixel grouping,” Pattern Recogn. 43(4), 1531–1549 (2010). 22. L. Zhang, R. Lukac, X. Wu, and D. Zhang, “PCA-Based Spatially Adaptive Denoising of CFA Images for SingleSensor Digital Cameras,” IEEE Trans. Image Process. 18(4), 2419–2434 (2009). 23. Philip R. Bevington and D. Keith Robinson, Data Reduction and Error Analysis for the Physical Sciences (McGrawHill, New York, 1992). 24. D. L. Donoho and I. M. Johnstone, “Ideal spatial adaptation by wavelet shrinkage,” Biometrika, 81(3),425–455 (1994). 25. J. Zhang, H. Luo, B. Hui, and Z. Chang, “Non-uniformity correction for division of focal plane polarimeters with a calibration method,” Appl. Opt. 55(26), 7236–7240 (2016).


Introduction
The four primary physical characteristics of light are intensity, wavelength, coherence, and polarization.Conventional imaging sensors can only record the intensity and wavelength as brightness and color respectively.However, polarization imaging is an emerging technology and it can provide extra useful information.Thus, it has been widely applied in the fields of material classification [1], fog de-hazing [2], 3D shape reconstruction [3], and biomedical imaging [4][5][6].
Current polarization imaging sensors are mainly categorized into division of time (DoT), division of amplitude, division of aperture, and division of focal plane (DoFP) [7][8][9][10].DoFP polarimeters are composed of a collection of pixelated polarizers aligned in four different orientations for 2-by-2 pixels.DoFP polarimeters are inherently temporally synchronized and all four polarized image are well-aligned.DoFP polarimeters have attracted increasing attention in recent years because of these advantages.
However, each pixel only records one out of four necessary intensity measurements, so the other three missing intensity measurements need to be interpolated from their neighbors for each pixel.Many interpolation methods [11][12][13][14][15][16] proposed in the past are based on the assumption of noise-free DoFP images.However, this assumption is unrealistic because an image is usually corrupted by noise during the processes of being captured, recorded or transmitted.The presence of noise in the DoFP images not only degrades the image quality, but also leads to reconstruction of false polarization information.Many denoising methods [17][18][19][20][21] are designed for monochromatic images and not directly applied to DoFP images due to the interlaced arrangement of polarizers.Even though we can carry out denoising by using some denoising methods for gray-scale images after interpolation, this strategy will distort the properties of noise and make it more difficult to remove.A similar problem exists in color images, using the RGB color filter array.Zhang et al. [22] proposed a PCA-based denoising method that is directly applicable to color filter array images.However, this method can't effectively suppress noise in the degree of linear polarization (DoLP) images because the residual noises are magnified during reconstructing the Stokes parameters.Meanwhile, denoising methods that are directly applied to DoFP images have not been previously investigated to the best of the authors' knowledge.
In this paper, we propose a new denoising method for DoFP polarimeters by using the dimen-sionality reduction property of PCA.Firstly, we model a variable block composed of intensity measurements orientated in four different orientations.Secondly, we select training samples by grouping the similar spatial structures to the variable block in the local window.Then, denoising is implemented in the PCA domain by using dimensionality reduction and linear minimum mean square error (LMMSE) estimation technology.Finally, we estimate residual noise in the Stokes parameters according to the standard error propagation techniques and carry out denoising for Stokes parameters.Both simulated and real DoFP images are used to evaluate the performance of the denoising method.Experimental results show that the proposed method can effectively suppress noise while preserving edges.

Denoising for DoFP images
Assume that an observed noisy DoFP image I n θ is generated by adding Gaussian noise n θ to the noise-free image I θ , as formulated by Eq. ( 1): where n θ ∼ N (0, σ 2 θ ) is a Gaussian noise with zero mean and the standard deviation σ θ , for θ = 0 • , 45 In order to fully exploit the correlation of the four intensity measurements at different orientations, the variable block should include all four intensity measurements equally, i.e. the number of these four intensity measurements should be same.For the sake of simplicity, we choose a variable block shown in Fig. 1 in the size of 2-by-2, but this block can be bigger in practice.We describe the variable block in the form of a column vector, as following: T is the corresponding noise-free variable block and the noise To suppress noise from the variable block, we need to obtain abundant training samples from the training block as shown in Fig. 1.The size of training block should be larger than that of variable block in order to obtain adequate training samples.Then, we just select some training samples that are similar to the center variable block among all the candidates for the training sample, which will ensure that the estimation of the covariance matrix is more accurate.The similarity measure is formulated by Eq. (3): ). ( where x n 0 , x n i are the variable block and the candidate for the training sample respectively, d(x n i , x n 0 ) is the distance between x n i and x n 0 , and σ d and σ r provide a tradeoff between the distance similarity and the intensity similarity.The two parameters (σ d and σ r ) should be determined according to the requirements in practice.We set the two parameters to be the same value because the distance similarity and the intensity similarity are equally important.When a candidate for the training sample is more similar to the variable block, the corresponding value of the similarity measure is larger.Thus, all the candidates for the training samples will be sorted by the similarity measure in descending order.We choose first M candidates as the training samples to estimate the covariance matrix.The item M should be large enough to guarantee a reasonable estimation.
We denote these training samples as The corresponding noise-free training samples are denoted as x 0 , x 1 , x 2 , . . ., x M −1 .The whole training samples can be described as following: where , we obtain Eq. ( 5): Next we calculate the covariance matrix of X n , denoted by C X n .
Since X and N are not correlated, items XN T and NX T are nearly zero matrix.Thus, Eq. ( 6) can be modified as following: where C N = diag{σ 2 0 , σ 2 45 , σ 2 135 , σ 2 90 }.Then, we can estimate C X by subtracting C N from C X n .The eigenvalue decomposition of C X can be written as: where Q is the eigenvector matrix and S = diag{λ 1 , λ 2 , . . ., λ m } is the diagonal eigenvalue matrix with λ 1 ≥ λ 2 ≥ . . .≥ λ m .Then, the orthonormal PCA transformation matrix can be calculated according to Eq. ( 9).
Next, we transform the training samples into PCA domain as following: In order to remove noise, we need to exploit the dimensionality reduction property of PCA.We reset the last several rows of Y n to be zeros.Meanwhile, the noise is further suppressed by using LMMSE estimation technology.
where C Y and C Y n are the covariance matrices of Y and Y n respectively.The de-noised variable block can be calculated by transforming Y n back to the time domain as following: where the subscript center represents the center 2-by-2 block of the variable block.The final de-noised single-pixel value is calculated by using the average of four results produced by four adjacent center blocks.Then, we can reconstruct Stokes parameters and DoLP images by using interpolation methods.
In the proposed denoising algorithm, the number of training samples and the sizes of variable block and training block should be determined in practice.We hereby give the relationship among these three parameters.The sizes of variable block and training block are denoted by k and L respectively.Because the number of training samples should be less than or equal to the number of all the candidates for the training samples, we can get the relationship: Meanwhile, the number of training samples is generally larger than the dimension, i.e.M ≥ k 2 .Thus, we can get the relationship between L and k as: L ≥ 3k − 2. In order to guarantee the uniqueness of center variable block, the constraint condition that L minus k is a multiple of 4 should be satisfied.Thus, once we determine the size of variable block, we can get the size of training block and the number of training samples.

Denoising for the Stokes parameters
As shown in Fig. 2, there are some residual noises in the DoLP image.These residual noises are generated during denoising the DoFP image and magnified in the following process of reconstructing the Stokes parameters.In order to remove these residual noises, we need to accurately estimate the variance of residual noises.The variance of the residual noise in DoFP image can be estimated as following: where I n θ and I θ are original DoFP image and de-noised DoFP image for θ channel respectively.The item n r represents the residual noise.Thus, we can estimate the level of residual noise δ 2 θ by subtracting E[( T ] from σ 2 θ .Next, we estimate the standard variance of the noise in the Stokes parameters.The first three Stokes parameters can be calculated as following: Thus, the standard variance of the Stokes parameters can be calculated according to Eq. ( 14) and the standard error propagation techniques [23].
Then, we apply the proposed denoising method to the Stokes parameters according to the estimation variance of residual noise.Finally, we reconstruct DoLP images by using the de-noised Stokes parameters.

Experimental setup
In order to get four true images in different polarization orientations, we apply DoT polarization imaging system to record a static scene.

Experimental results
In this section, visual comparison and PSNR analysis are applied to the simulated DoFP images.Furthermore, real DoFP images are used to visually compare.The proposed method is implemented in Matlab R2012b and runs on a desktop computer with Intel(R) Core(TM) i3-2130 @3.40GHz(CPU) and 4GB of RAM.Assume that the number of image pixels is n and the computational complexity of similarity measure is O(T ).It takes most of the computational cost in PCA transformation.It requires k 2 (2M − 1) + k 4 (M − 1) additions, k 4 M + k 2 multiplications and an eigenvalue decomposition for each block.Thus, the time complexity of the proposed method is O([T + k 4 M + k 6 ]n).In order to reduce the computational complexity, the size of variable block is set to 6 × 6 and that of training block is set to 34 × 34.We choose first 50 candidates for the training samples to estimate the covariance matrix.

Results on simulated DoFP images
Denoising results for simulated DoFP image are shown in Fig. 2. The first column presents intensity images, i.e. S 0 .The second column presents DoLP images.The original images in the first row are produced by using four true images which are recorded by the DoT polarization imaging system.The original images are used to visually compare the denoising results and as the noise-free images to calculate the PSNR.The simulated noise-free DoFP image is interpolated to generate four interpolated images which are in different polarization orientations and the four interpolated images are of the same size as the simulated noise-free DoFP image.The interpolation method proposed in [15] is used.Then, these four images are used to reconstruct the intensity image and DoLP image as shown in the second row.The interpolation error is inevitably introduced during the interpolation process.The third row presents noisy intensity and DoLP images.They are reconstructed as following: we firstly add Gaussian noise to four channels separately with σ 0 • = σ 90 • = 1, σ 45 • = σ 135 • = 2 before the interpolation.Next, four interpolated images are generated by using the interpolation method.Finally, intensity and DoLP images are reconstructed by using these four interpolated images.The intensity image is slightly corrupted by noise, but the DoLP image is full of noise.The fourth row shows the de-noised images by using the proposed method without implementing residual noise-removed.Although the noise is obviously mitigated, some residual noise (RN) is strong in the DoLP image because the RN is generated during denoising the DoFP image and magnified in the pro- The PSNR is used to as an objective evaluation criterion to analyze the denoising performance.The standard variance of the noise is increased by every 0.5 from 1 to 25.The corresponding PSNR curves are drawn in Fig. 3.The PSNR of noisy images decreases quickly when the noise increases.Experimental results show that the PSNR is improved significantly after denoising.Especially, the image quality is further improved after the residual noises are removed.
The denoising performance for different illuminations should be evaluated because the DoFP polarimeters work with various dynamic ranges.The scene shown in Fig. 3 is used.The ratios of original illumination range from 0.05 to 1.0 and the step is 0.05.In order to obtain the true  images for RMSE calculation, the image at each illumination is sampled 300 times and averaged to generate essentially true image.Then, the single frame images and de-noised images are used to calculate the root mean square error (RMSE).The de-noised images are obtained by using the proposed method.The RMSE comparison results are shown in Fig. 4. The RMSE of intensity image increases with the illumination increasing.The RMSE of DoLP is approximately to a constant because DoLP is independent of intensity measurements.

Results on real DoFP images
The standard variance of noise for each channel should be estimated before implementing denoising.The H × W DoFP image is firstly down-sampled into four H/2 × W/2 sub-images.Then, a one-stage wavelet transformation is applied to each sub-image.According to the noise estimation theory [24], the standard variance of noise can be estimated as formulated in Eq. ( 16).
where CD is the diagonal coefficient matrix obtained by wavelet decomposition.Real DoFP images of cars on a snowy day are recorded by the DoFP camera as shown in Fig. 5.The first and second columns present intensity and DoLP images respectively.The first row presents raw images that are full of random noises and non-uniformity errors.Previous research has shown that these non-uniformity errors can be removed by using correction methods [25].The images in the second row are obtained by using non-uniformity correction (NUC) technology proposed in [25].However, NUC technology does not remove random noise effectively.Thus, we focus on denoising random noises in this paper which can be implemented after finishing NUC.The third row shows de-noised images without removing residual noise.Slight residual noises exist in the DoLP image.The last row presents RN-removed images.The noise is removed effectively and the best visual effects are achieved after removing residual noises.

Conclusion
In this paper, we propose a PCA-based denoising method that is directly applied to DoFP images.The spatial correlation of different polarization orientations is fully exploited during the denoising process.The dimensionality reduction and linear minimum mean squared-error estimation in the transformed domain are used to suppress noise.Meanwhile, we estimate residual noise accurately according to the standard error propagation theory, which is the key to remove residual noises in the following process.Both simulated and real DoFP images are used to evaluate the denoising performance.Experimental results show that the proposed denoising method for DoFP polarimeters can effectively mitigate noise while preserving edges and details.However, the proposed method does not work well for non-Gaussian noise or Gaussian noise with large magnitude.Therefore, we will investigate more general denoising methods for DoFP images in the future.
A Newport 10LP-VIS-B linear polarizer is mounted on a rotation stage (Thorlabs PRM1Z8E) to generate four different polarization orientations from 0 • to 135 • offset by 45 • .The linear polarizer is placed in front of a gray-scale camera (DMK 22AUC03, 640 × 480 pixels) that records four true images in different polarization orientations.Each image is sampled 300 times and averaged to generate essentially noise-free images.Next, the simulated noise-free DoFP image is generated from these four true images as the sampling pattern shown in Fig.1, e.g. the pixels in the odd rows and columns of the simulated noisefree DoFP image are copied from the pixels in the corresponding rows and columns of the true image in 0 • orientation.Thus, the simulated noise-free DoFP image and four true images are the same size.Then, we add Gaussian white noise to the simulated noise-free DoFP image for four channels respectively.Finally, we evaluate the denoising performance by visual comparison and peak signal to noise ratio (PSNR) comparison.Real DoFP images are recorded by our self-developed DoFP camera in the visible waveband (752 × 582 pixels, extinction ratio: 9.58, transmittance: 71%).
Figure 4 also shows the smaller RMSE is obtained after denoising.