Polarimetric data reduction : a Bayesian approach

In this paper, we introduce a general Bayesian approach to estimate polarization parameters in the Stokes imaging framework. We demonstrate that this new approach yields a neat solution to the polarimetric data reduction problem that preserves the physical admissibility constraints and provides a robust clustering of Stokes images in regard to image noises. The proposed approach is extensively evaluated by using synthetic simulated data and applied to cluster and retrieves the Stokes image issuing from a set of real measurements. ©2007 Optical Society of America OCIS codes: (120.5410) Polarimetry; (110.2960) Image analysis; (100.3190) Inverse problems; (000.3860) Mathematical methods in physics. References and Links 1. D. Miyazaki, M. Saito, Y. Sato, and K. Ikeuchi, “Determining surface orientations of transparent objects based on polarization degrees in visible and infrared wavelengths,” J. Opt. Soc. Am. A 19, 687-694 (2002). 2. D. Miyazaki, M. Kagesawa, and K. Ikeuchi, “Transparent surface modeling from a pair of polarization images,” IEEE Trans. PAMI 26,920-932 (2004). 3. J. M. Bueno and P. Artal, “Double-pass imaging polarimetry in the human eye,” Opt. Letters. 24 64-66 (1999). 4. S. D. Giattina, et al., “Assessment of coronary plaque collagen with polarization sensitive optical coherence tomography (PS-OCT),” Int. J. Cardiol. 107, 400-409 (2006). 5. D. H. Goldstein, D. B. Chenault, and Society of Photo-optical Instrumentation Engineers, Polarization: measurement, analysis, and remote sensing II, 19-21 July, 1999, Denver, Colorado. 1999, Bellingham, Washington: SPIE. ix, 426 p. 6. M. H. Smith, “Optimizing a dual-rotating-retarder Mueller matrix polarimeter,” in Polarization Analysis and Measurements IV, SPIE (2001). 7. J. S. Tyo, “Design of optimal polarimeters: maximization of signal-to-noise ratio and minimization of systematic error,” Appl. Opt. 41 619-630 (2002). 8. S. Ainouz, J. Zallat, A. de Martino, and C. Collet., “Physical interpretation of polarization-encoded images by color preview,” Opt. Express 14, 5916-5927 (2006). 9. S. N. Savenkov, “Optimization and structuring of the instrument matrix for polarimetric measurements,” Opt. Eng. 41, 965-972 (2002). 10. J. Bernardo and A. Smith, Bayesian Theory, (Wiley, 2000). 11. A. Gelman, J. Carlin, H. Stern, and D. Rubin, Bayesian data analysis, Second ed., (CRC Press, 2003). 12. S. Z. Li, Markov random field modeling in image analysis, Second ed., (Springer, 2001). 13. A. Gray, J. Kay, and D. Titterington, “An empirical study of the simulation of various models used for images,” IEEE Trans. PAMI, 16, 507-513 (1994). 14. A. Dunmur and D. Titterington, “Computational Bayesian analysis of hidden Markov mesh models,” IEEE PAMI. 19, 1296-1300 (1997). 15. S. Richardson and P. Green, “On Bayesian analysis of mixtures with an unknown number of components (with discussion),” J. R. Stat. Soc. Ser. B. 59, 731-792 (1997). 16. J. Besag, “On the statistical analysis of dirty pictures (with discussion),” J. R. Stat. Soc. B 48, 259-302 (1986).


Introduction
The main interest of the Stokes-Mueller formalism in optical imaging is mainly due to the definition of light polarization parameters in terms of real quadratic observables (intensities) which are directly sensed by CCD detectors.This allows extending classical intensity-wise imaging systems to acquire Stokes images through the use of Polarization State Analyzers (PSA) in front of the camera.Many such systems that use different polarization modulation techniques, have been built in recent years for many application areas ranging from metrology [1,2] to biomedical imaging [3,4], and remote sensing [5].
Whatever the modulation principle used in the PSA configuration, the measurement procedure remains the same.Indeed, the PSA analyses the incoming Stokes vector by measuring its projections onto N ( 4 N ≥ ) independent basis states.The complete set of the N measurements yields a matrix equation, which relates the out-coming Stokes vector S from the sample to the measured raw intensity data vector I for each pixel.This matrix relation has to be inverted properly.This implies that an efficient calibration procedure of polarimetric imaging systems must be employed in order to extract the desired polarization-encoded images effectively.Many interesting studies in regard to this problem have been published in the recent literature, see for example [6][7][8][9].All of these papers focused on the calibration strategies to adopt and on the suited optical configurations in order to obtain the best conditioning of the aforementioned matrix equation to yield the polarization-encoded images.
We note further that the problem of error equalization and its propagation to the estimated polarization parameters channels was also analyzed in details in Refs.[7] and [8].Indeed, systematic errors related to misalignments and image noises inherently related to the use of photonic detectors that propagate to the polarization channels alter the analyzing algorithms performances and may prevent the physical admissibility of the estimated quantities.Even though the proposed solutions are of great interest for polarimetric data reduction and provide the optimal pixel-wise approach in the mean least square sense they do not benefit from the spatial sampling naturally provided by the bi-dimensional nature of the measurements.
In this paper, we reformulate the polarization data reduction problem in the general framework of Bayesian inference theory and demonstrate that this new approach allows a neat solution to the polarimetric data reduction problem that preserves the physical admissibility constraints and provides a robust clustering of Stokes images in regard to image noises.

Observation model in the Stokes imaging polarimetry framework
In such a system, the polarization state of the light coming from the scene is acquired.This can be done by inserting a complete state analyzer in front of the camera.At least, four independent states of the analyzer are needed to acquire all elements of the Stokes vector for each pixel location ( ) , y x j j .This can be summarized by the following equation: where I m is the acquired image, P is the "polarization measurement matrix" (PMM) that depends on a parameter vector (η), while S in is the Stokes vector to be estimated.We note further that ( ) , m y x j j I will denote the intensity vector measured at the ( ) In the following, the image I m , and the Stokes image will be arranged in data cube structures such that: I j j j j J j J j N S j j j j J j J j Classical data reduction employs the pseudoinverse # P of the PMM to obtain a minimumnorm, least-square estimate ˆin S for the Stokes vector at each pixel location: If one uses four probing states to measure the Stokes vector, # P is identical to the inverse of P provided that η is chosen such that P is invertible.The reason for using N intensity measurements where 4 N > that yield an over-determined system for Eq. ( 1) is to reduce the system sensitivity to systematic errors as well as image noises.The drawback of such a strategy is an increase in acquisition, and data processing and handling time.
Systematic errors and the use of CCD detectors lead to noises that contaminate the measured raw data.This includes misalignment errors of optical elements, readout noise, thermal noise, dark current noise and photon noise among others.These perturbations can be accounted for by replacing Eq. ( 1) by: ( ) ( ) ( ) We note further that the number of detected photons depends mainly on the quantum efficiency ( E Q ) of the used photosensitive device (e.g.CCD camera) which in turn depends on the wavelength.This dependency can be expressed formally as: ( ) where ph i N is the number of photons impinging onto the detector and λ is the considered wavelength.
If now we consider the case of N intensity measurements where each one contains nearly the same amount of statistical noise, we can assume that the same number of photons ( ) contributes to each one.Finally, p δ I can be written as: We note further that in the case of coherent illumination, the ρ factor includes also the contribution of the speckle noise.
The final observation model can be written down: where ( ) The Stokes images estimation problem will be tackled in a Bayesian framework [10,11], where all variables involved will be considered as random variables.

Underlying segmentation model
The main feature of the estimation procedure is to explicitly model and estimate the unknown underlying segmentation image and to account for the physical admissibility constraints.The segmentation image will be denoted where K is the number of classes.
Considering such segmentation is justified since each Stokes vector is taken from a set of reduced cardinality: a Stokes vector is assumed to be invariant for a given type of class and the cardinality of the set is equal to the number of classes types in the analyzed sample.Even though the model can handle any number of probing states, without loss of generality, we consider here the case of a four probing states polarimeter, i.e., four intensity images are acquired to estimate the Stokes image.Let Φ being the matrix defined as: where Since there is a high probability that two neighboring pixels in a Stokes channel image have the same φ value, 0 I will be modeled as a Markov random field (MRF) [12].Having a closed form partition function is a desirable feature to ease estimation of the MRF hyperparameter.We consider here the Markov mesh model studied by Gray et al. [13] and by Dunmur and Titterington [14].Only the main lines of the model are summarized here, the interested reader may refer to the cited references for more details.The MRF is governed by the hyperparameter 0 β , which controls the size of the clusters: values of 0 β close to 1 yield monochrome images, whereas 0 1/ K β = corresponds to an independent pixel wise prior.The proposed probabilistic model reads: are the four neighbors labels values.This third-order Markov mesh model, where a given pixel admits three causal neighbors, is equivalent to a second order MRF where a pixel admits eight neighbors.This model yields ( ) [ ]

I
, where the integers n i depend only on I 0 .

Noise model
The noise will be supposed to be independent and Gaussian distributed according to , μ σ are the noise parameters (mean and standard deviation) that affect the k th class in the j n channel image.
Let us define Θ as the set , where All these parameters are unknown and will have to be estimated.

General probabilistic model
The overall probabilistic model may be depicted using a directed acyclic graph (DAG) (see Fig. 1).All variables in square boxes but I m represent fixed hyperparameters.Hyperparameters values are set so that the corresponding priors are weakly informative in the region of interest (see for example p. 735 of Ref. [15] and the appendix A).All prior distributions are chosen to be conjugate, as is classically done, to allow closed form computation in the optimization procedure which will be detailed in the sequel.Conjugacy may loosely be defined as follows.A prior ( ) p θ is conjugate to the likelihood ( ) θ is in the same class of distributions as the prior ( ) p θ .Conjugate priors are often good approximations and they simplify computations.In the present case, conjugacy allows parameters to be updated using closed form formulas, whereas non conjugate priors would have imposed resorting to iterative algorithms.Interested reader may refer to [10,11] for extensive details on this topic.Let us emphasize that all probability laws derive directly from the Gaussian model used for the noise and from the MRF model used for I 0 .Considering these assumptions and the conjugacy property, the other laws are fixed.This explains for example why the ( ) n j k μ 's have a Gaussian distribution and why the ( ) 's follow a Gamma distribution.The corresponding joint probability distribution π can be written as: where dependencies upon the hyperparameters were dropped for sake of clarity.
The stochastic dependencies are formulated as follows: ( ( ) Analogous choices for the distribution of Q were made in Ref. [15].U, Ga, and Be define respectively the uniform, Gamma, and Beta laws, as:

Estimation procedure
All information of this estimation problem is contained in a stochastic point of view in the posterior distribution ( ) We retain here the maximum a posteriori (MAP) estimate, which is given by the maximization of the posterior or equivalently by the joint distribution.There is no analytical expression for the solution to this optimization problem.We have to resort to an iterative algorithm, where each variable is optimized in turn.This yields a local maximum of the posterior distribution.Details are omitted since each elementary optimization step is straightforward.Let us merely mention that I 0 is updated using a single raster scan (Iterated Conditional Modes, see Ref. [16]) and that a local optimization method is used to update The proposed model can handle different practical situations that may arise in the framework of Stokes imaging measurements.Accounting for a number of probing states greater than four or considering different noise characteristics is straightforward.

Simulation results
To illustrate the relevance and the efficiency of the Bayesian approach in polarimetric data reduction, we synthesized a 64-pixel × 64-pixel Stokes image by using the four class label map shown in Fig. 2(a The PMM of an optimal rotating-retarder polarimeter as defined in Ref. [7] was used to generate the corresponding intensity images by mean of Eq. ( 1).The parameter vector η that controls the PMM is given by ( 132 ,{ 15.12 , 51.7 }) i η δ θ = = °= ± °± ° where δ stands for the retardance of the birifringent wave-plate while θ i represent its angular positions.Practically, acquired intensity images are not instantaneous snapshots but result from photons accumulation due to exposure time settings.Moreover, each intensity image combines many sources of errors and noises that make reasonable to assume that ε follows a centred Gaussian distribution.Hence, a zero mean white Gaussian noise with variance of 2 n σ was added to each intensity image for this simulation providing the necessary inputs to our model, e.g.intensity images as well as the PMM P. The aim of this section is to compare the estimated Stokes image ˆLS S resulting from Eq. ( 3) with the one resulting from our Bayesian model ( ) Figure 3 shows the four channels images corresponding to ˆLS S for a noise variance 2 0.01 . We apply now the Bayesian approach introduced in the preceding sections to the same estimation problem.We note that our model provides one Stokes estimate per class.The noise estimated variances are 0.0103 for 0.01 and 0.1015 for 2 0.1 . Table 1  × − s s s as well as the ratio of misclassified pixels of each class for the two considered variances.As expected the estimated ρ factor is practically zero.Figure 5 shows the true polarization locations over the Poincaré sphere as well as the estimated ones for the two considered noise variances.Furthermore, we point out that contrary to the min-norm solution, our model ensures always the physical admissibility of the estimated Stokes parameters.This can be seen clearly by comparing Figs. ( 4) and (5).   Figure (6) presents the estimated noise variance versus the true one for different noise variance levels.Circles represent the estimated variance by using the pixels that belong to the largest class.Horizontal bars correspond to the variations of the estimates over the four classes.We conclude that the proposed Bayesian approach has an excellent behavior with regard to different SNR.Table 1 and Fig. 5 emphasize the accuracy and the performance of the proposed model in estimating the Stokes image as well as the noise parameters.
At this stage, we relax the hypothesis of a uniform noise field over the whole image and consider the case of different noises that contaminate independently each class.One can imagine easily situations where the polarization state of the light emergent from a particular class of pixels to be nearly orthogonal to the analyzer state which reduces the photon flux impinging on the detector which in turn decreases the SNR for the considered pixels.So we carried out simulations with the following noise variances ( )  Estimated variances were found to be ( ) ˆ0.1, 0.01, 0.049, and 0.203 Figure 7 shows the location of the estimated polarization states and the true ones over the Poincaré sphere.Again, we observe the excellent behaviour of our model for the case of different noises that reach each class.

True measurement case
Here, the developed model is applied and validated with real measured images acquired by a Stokes imaging polarimeter.The used object to carry out these measurements consists of two dichroic patches (Polaroid) contacted on a diffusing glass slide with a drop of water.The overall mount was backlight illuminated by a vertical polarized beam and the intensity images were sensed by a four probing states rotating quarter-wave-plate analyzer in front of a 12-bits CCD scientific camera.We note that the different orientations were given to transmission axes of the Polaroid shapes to obtain different signatures at the output.Figure 8 shows the intensity images issued from this configuration; Fig. 9 shows the minnorm estimate of Stokes images; and Fig. 10   As expected, the noise that is present in raw intensity data propagates into the Stokes channels.Additional treatments (clustering, validations, etc.) are needed to make the most of these images.Figure 11 illustrates clearly this fact.Indeed, the clutters of points that lie inside the Poincaré sphere do not permit a clear interpretation of polarization properties of the scene.On the other side, our Bayesian approach yields a robust and accurate solution as illustrated by Figs.(10)  ŝ is located near the center (Fig. 10) indicates that the transmission axis of this patch is oriented horizontally which extinguish a large amount of the vertically polarized incident beam.

•
The polarization state corresponding to 2 ŝ lies in the (S 1 -S 2 ) plan indicating a linear state as expected.

•
The incident beam undergoes isotropic depolarization induced by the diffusing glass slide.This is confirmed by observing that the degree of polarization (DOP) of 2 ŝ and 3 ŝ are nearly equal ( ) ( ) ˆDOP DOP ≈ s s .

Conclusion and future extensions
In this work we addressed the problem of estimating Stokes vectors from indirect noisy intensity measurements.Min-norm least squares estimate represents the state of the art method in the field.Because of the influence of measurement noise in the procedure, the estimate may violate the physical admissibility conditions.Moreover, the Stokes vectors vary in a given class, which is in contradiction with the prior knowledge.These shortcomings may lead to severe degradation of the estimated Stokes images, thereby hindering their interpretation and exploitation.
We propose here an alternative markovian Bayesian model accounting for both prior informations (the physical admissibility constraints).The MAP estimate was retained.The proposed optimization algorithm yields robust and accurate estimates on both synthetic and real data considered here.
Work is in progress to incorporate of such priors in the Mueller imaging framework.

Appendix: hyperparameters choice
The hyperparameters are chosen so as to yield weakly informative priors (see for example Ref. [15], p. 735 where equivalent choices are made).We set: ( We remind the reader that ˆLS S is the min-norm estimate of the Stokes image.

Fig. 1 .
Fig.1.The directed acyclic graph (DAG) that represents the dependence relations between the variables.According to usual conventions, square boxes represent fixed or observed quantities whereas circles represent quantities to be estimated.

,
function, taking values 0 or 1 whether the condition in parentheses is false or true.The condition k C is the conjunction of conditions ( ) admissibility of the estimated Stokes vectors.
to the unconstrained problem, e.g., not accounting for C k , violates the conditions C k .The number K of classes is not updated during the procedure but optimizations considering all values of It is well known that the initialization stage is very important for such procedures.The DAG is here initialized as the min-norm estimate ˆin S obtained by Eq. (3).The initialization of I 0 relies on the ˆin S clustering by a variant of the C-means algorithms family where the used distances were redefined in relation with the Stokes images specificities[8].Finally, Φ is initialized by the mean values of ˆin S corresponding to each class.Initializing the other variables follows straightforwardly as given in the Appendix.
pixels in the label map.Fig.2(b)shows these polarization state locations on the Poincaré sphere.

Fig. 2 .
Fig. 2. (a).The label map that was used to generate the simulated Stokes image.The s 1 Stokes vector was assigned to black pixels (1215 pixels), s 2 was assigned to heavy gray pixels (436 pixels), s 3 was assigned to light gray pixels (997 pixels), and S 4 was assigned to white pixels (1448 pixels).(b) Red circles indicate the positions of the four polarization states corresponding to the s i vectors on the Poincaré sphere.A look at Figs 4(a) and 4(b) shows that even at low level noises, the min-norm inversion yields unphysical estimated Stokes vectors (Degree of Polarization (DOP) greater than one).Moreover, if the points remain well grouped into four different classes for a 20 dB SNR [Fig.4(a)], this is no longer the case for a 10 dB SNR.This suggests that classical clustering algorithms may fail to provide a workable label map [see Fig. 4(b)] that allows an accurate physical interpretation of the Stokes image content.We apply now the Bayesian approach introduced in the preceding sections to the same estimation problem.We note that our model provides one Stokes estimate per class.The noise estimated variances are 0.0103 for 0.01 and 0.1015 for 2 0.1

Fig. 3 .Fig. 4 .
Fig. 3. From upper left to bottom right, the four Stokes channels (s 0 (total intensity), s 1 , s 2 , and s 3 ) obtained by the min-norm solution ˆLS S for a noise variance of 0.01.

Fig. 6 .
Fig. 6.Estimated noise variance vs. true noise variance.Circles represent the estimates that correspond to the largest class.Horizontal bars correspond to the variations of the estimates over the four classes.

Fig. 7 .
Fig. 7. Locations of the true polarization states (circles) and the estimated ones (crosses) where different noises reach each class.

Fig. 8 .
Fig. 8. Intensity images corresponding to four polarization probing states used to retrieve the Stokes image.Gray values are scaled for display purposes.

Fig. 9 .
Fig. 9. Stokes channels images issued from intensity images of Fig. 8 by min-norm solution.Gray values are scaled for display purposes.

Fig. 10 .Fig. 11 .
Fig. 10.The label map obtained with the proposed Bayesian approach.Different gray values correspond to different polarization signatures.
and 11(b).The obtained polarization-encoded images are noise free and a neat solution is provided to the polarization states estimation problem.The three locations inside the Poincaré sphere [Fig.11(b)]correspond to the three physical classes of the object, i.e., the two dichroic patches and the diffusing glass slide.The three estimated states are normalized to the s 0 element and are given by: the smallest patch, 2 ŝ to the trapezoidal shape, and 3 ŝ to the background.By looking at Eq. (18) and Figs.10 and 11, one can observe the following:•The polarization state corresponding to 1 January 2007 / Vol. 15, No. 1 / OPTICS EXPRESS 90 estimated Stokes vectors for the different classes.The last two rows lists the RMSE defined as 100 /

Table 1 .
Estimated Stokes vectors of each class.The last two rows lists the RMSE and the ratio of misclassified pixels for each class.