Clustering of Mueller matrix images for skeletonized structure detection

This paper extends and refines previous work on clustering of polarization-encoded images. The polarization-encoded images used in this work are considered as multidimensional parametric images where a clustering scheme based on Markovian Bayesian inference is applied. Hidden Markov Chains Model (HMCM) and Hidden Hierarchical Markovian Model (HHMM) show to handle effectively Mueller images and give very good results for biological tissues (vegetal leaves). Pretreatments attempting to reduce the image dimensionality based on the Principal Component Analysis (PCA) turns out to be useless for Mueller matrix images. © 2004 Optical Society of America OCIS codes: (120.5410) Polarimetry; (110.2960) Image analysis; (100.5010) Pattern recognition and feature extraction. ___________________________________________________________________________ References and links 1. J. Zallat, C. Collet, and Y. Takakura, “Clustering of polarization-encoded images,” Appl. Opt. 43, 283-292 (2004). 2. J.M. Bueno, and P. Artal, “Double-pass imaging polarimetry in the human eye,” Opt. Lett. 24, 64-66 (1999). 3. P.Y. Gerligand, M.H. Smith, and R.A. Chipman, “Polarimetric images of a cone,” Opt. Express 4, 420-430 (1999). http://www.opticsexpress.org/abstract.cfm?URI=OPEX-4-10-420 4. G.D. Lewis, D.L. Jordan, and P.J. Roberts, “Backscattering target detection in a turbid medium by polarization discrimination,” Appl. Opt. 38, 3937-3944 (1999). 5. J.-N. Provost, C. Collet, P. Rostaing, P. Pérez, and P. Bouthemy, “Hierarchical Markovian Segmentation of Multispectral Images for the Reconstruction of Water Depth Maps,” Computer Vision and Image Understanding 93, 155-174 (2004). 6. Collet, C., M. Louys, C. Bot, and A. Oberto, “Markov Model for Multispectral Image analysis : application to Small Magellanic Cloud segmentation,” in International Conference on Image Processing ICIP'03. 2003. Barcelona, Spain. 7. P.A. Devijver , “Baum's forward-backward algorithm revisited,” Pattern Recognition Lett. 39, 369-373 (1985). 8. N. Giordana, and W. Pieczynski, “Estimation of generalized multisensor hidden Markov chains and unsupervised image segmentation,” IEEE Transactions on Pattern Analysis and Machine Intelligence 19, 465-475 (1997). 9. J. M. Laferté, P. Pérez, and F. Heitz, “Discrete Markov Image Modeling and Inference on the Quad-tree,” IEEE Transactions on Image Processing 9, 390-404 (2000). 10. C. H. Chen, L.F. Pau, and P.S.P. Wang, Handbook of Pattern Recognition and Computer Vision. WORLD SCIENTIFIC. 1044 (2000.) 11. A. P. Dempster, A.P., N.M. Laird, and D.B. Rubin, “Maximum likelihood from incomplete data via the EM algorithm,” Royal Statistical Society, 1-38 1976. 12. R. O. Duda, P.E. Hart, and D.G. Stork, Pattern Classification, 2nd Edition Ed. Cloth. 680 (2000). 13. Movie of Mueller images after PCA transform. 2004. ftp://picabia.u-strasbg.fr/pub/www/collet/Publis/Animations/acp.gif 14. S. D. Baker, Unsupervised pattern recognition, PhD Thesis Antwerpen University (2002). http://www.ruca.ua.ac.be/visielab/theses/debacker/SteveThesis.pdf 15. A. Hyvärinen, J. Karhunen, and E. Oja, Independent Component Analysis, John Wiley & Sons (2001). 16. Movie of HMCM segmentation. 2004. ftp://picabia.u-strasbg.fr/pub/www/collet/Publis/Animations/Iteration_map.gif #3811 $15.00 US Received 11 February 2004; revised 12 March 2004; accepted 23 March 2004 (C) 2004 OSA 5 April 2004 / Vol. 12, No. 7 / OPTICS EXPRESS 1271 17. Movie of HMCM segmentation map with different number of classes. 2004. ftp://picabia.u-strasbg.fr/pub/www/collet/Publis/Animations/3-6classes_resample.gif 18. Movie of HHMM segmentation. 2004. ftp://picabia.u-strasbg.fr/pub/www/collet/Publis/Animations/HHMM_animation.gif 19. S. R. Cloude and E. Pottier, “Concept of polarization entropy in optical scattering,” Opt. Eng. 6 1599-610 (1995). 20. S. Y. Lu and R.A. Chipman, “Interpretation of Mueller matrices based on polar decomposition,” J. Opt. Soc. Am. A 13 1-8 (1996).


Introduction and objectives
We have previously reported on Mueller images clustering algorithms [1] after noticing that there was no available literature that take into account the image structure of the measurements.Indeed, all available references use only a physical pixel-based, processing approach.In that study, the aim was to consistently consider the two-dimensional structure of polarization-encoded images.The clustering procedure based on Markovian algorithm was able to segment properly the Mueller matrix image under poor illumination conditions.The task of segmenting skeletonized structures remained pending, due to possible block effects induced by HHMM algorithms based on a quad-tree topology.In this paper, we present our approach to deal with this peculiar case.
In the following, the Mueller matrix image is arranged in a three dimensional structure ( N M p × × ) where N M × defines the image size in pixels and 1,16 p = indexes the Mueller matrix elements.We will use indifferently the 'Mueller matrix image' or the 'Mueller image cube' terms.Due to computing time considerations, we use Hidden Markov Chains Model (HMCM) or Hidden Hierarchical Markovian Model (HHMM) instead of Markov Field Models for their fast convergence features.Both methods (based on chain or quad-tree model) are unsupervised (model parameters estimation stand alone) and robust: segmentation maps obtained are quite similar.These approaches are validated on raw Mueller images of leaves where one can easily observe the great robustness and the capacity of these Markov models to extract finest details, with low computing time.
In this work, we address the problem of analyzing Mueller images of physical objects and explore the potential of this technique for classification issues.We describe shortly our image acquisition system and, after recalling some definitions and a survey of the actual state of art concerning the interpretation of Mueller matrices, we propose an analyzing procedure of fully polarimetric images based on Markovian assumption in a Bayesian inference framework.

Experimental setup
A fully polarimetric imaging system that records spatially dependent intensity patterns of polarized light that is diffusely scattered from a target allows to acquire the whole Mueller matrix image [1].We define the Mueller matrix image as the two-dimensional measurement of the Mueller matrix attached to each pixel.Many such imaging systems have been designed and built for a wide variety of applications, including, medical [2], metrologic [3] and remote sensing [4] among others.This study was made possible by a fully polarimetric, high precision, and well calibrated imaging system that we developed.The experimental setup is a classical dual rotating quarter wave plates.It uses an incoherent source coupled with an interferential filter centered at λ = 632.8nm with a 10 nm bandwidth, positioned in front of the camera.The incident beam is polarized through a polarizing optic before impinging on the target.The object is imaged in backlight illumination (see Fig. 1) through the polarization analysis optic, and onto a digital CCD camera (12 bits image definition).α , this intensity can be written as where ij M is the pixel's ( ) , i j 4×4 Mueller matrix.

To fully determine
where Appropriate values of 1 k α and 2 l α ( , 1,4) k l = that maximize the determinants of A and P are obtained by minimization to avoid matrix singularities.
Calibration procedures ensure a maximum relative error over the whole image of less than 0.1% for the free-space Mueller matrix, and less than 1% for the Fresnel reflection matrix of a BK7 prism's surface.
To enhance the signal level, every acquired intensity image is the result of six frame integrations.

Physical analysis of Mueller matrices
Mathematically, Mueller matrices act as linear operators on the incoming Stokes vectors yielding the outgoing Stokes vectors.This provides a mathematical formalism known as "Stokes-Mueller formalism" to treat polarized interaction in a coherent way.The main interest of this approach lies in the addition theorem of Stokes vectors and the optical equivalence principle, allowing decomposition of a general matrix as the incoherent sum of matrices representing "simpler" processes.Furthermore, Mueller matrices are able to represent depolarizing systems in contrast with Jones vector-matrix formalism.Although it is always possible to associate any (non-depolarizing) system's Jones matrix to its Mueller matrix, the reverse is not always true.This is due to the fact that Mueller matrices have sixteen degrees of freedom while Jones matrices admit only eight.The extra degrees of freedom come from the existent correlations between the elements of the Jones matrices.Analyzing these correlations may be of great interest to understand the underlying physical interactions.This task is simplified by using the coherency approach introduced by Cloude [5].
Cloude shows that there exists a one to one correspondence between any Mueller matrix ( ) M and a 4 4 × complex hermitian semi-definite positive matrix ( ) T called "system coherency matrix".In this paper, we shall note ( ) O the linear operator that allows to pass from a Muller matrix to the system coherency matrix and ( ) O − its inverse.Cloude showed that a necessary and sufficient condition for ( ) M to be physically realizable is that all eigenvalues ( ) ; 0 3 i i λ = L of its associated coherency matrix ( ) T are real positive numbers.

Moreover, ( )
M is non-depolarizing, if and only if the rank of ( ) T equals one.One can combine the coherency matrix eigenvalues in a scalar quantity H called "entropy" that characterizes the polarimetric disorder in the physical system.H is defined as The system is non-depolarizing if its entropy is zero and a perfect depolarizer has its entropy equal to one.Practically, systems showing low entropy can be approximated by a nondepolarizing Mueller matrix associated to the largest eigenvalue and its corresponding eigenvector.
Another interesting approach to physically analyze experimentally-obtained Mueller matrices is by using the polar decomposition developed by Lu and Chipman [6] From these matrices the diattenuation, retardance, and depolarization properties are readily determined.From these matrices the diattenuation, retardance, and depolarization properties are readily determined.
In order to extract high-level information from the Mueller matrix image, one may need to segment the image properly, prior to any post-processing algorithms.Image segmentation is one of the first fundamental computer vision's problems, it consists in the partitioning process of an image into non-overlapping regions, where each one exhibits homogeneity features (intensity distribution, spatial clustering, spectral or in-scale homogeneity) that depend on the application needs.We show here that the recent Markovian developments in image segmentation, applied to Mueller matrix images, prove to be robust and efficient for skeletonized structure detection.

Markovian Model and Bayesian Inference
Our motivations for using Hidden Markov Model are to provide fast computations and efficient structures to analyze parameter data cubes.Indeed, different communities needs tools to interpret large data corresponding to multispectral observations (e.g., remote sensing [7], astronomical images [8]), parameter observations (e.g., polarimetric imaging [1]) or volume information (e.g., multimodal slices of RMI, SPECT, X cubes for medical applications), etc.
The segmentation problem is to estimate the unobserved realization X x = from the observed realization Y y = , where , , ,..., K ω ω ω Ω = and n , respectively.One assumes that each observation s y is a filtered and noisy measurement of the observed scene : • the segmentation process tries to determine which label (or class { } . The segmentation step is thus of great importance to cluster the pixels into sets of similar spectral, statistical or geometric features, in order to be able to extract their characteristics with further discriminant analysis (classification step).In the context of unsupervised Bayesian segmentation, we need to estimate all the parameters defining the joint distribution ( , ) p X Y between the label field X that have to be estimated and the noisy observed field Y.In other words, X p is the marginal distribution of X, and ( / ) p Y X is the distribution of Y conditional to X. Denoting by C the set of cliques (a clique being either a singleton or a set of sites mutually neighbors), the X distribution is where Z is a normalization coefficient and c Ψ stand for a potential function [10], [11].Moreover, a trade-off with the data-driven term is necessary to maximize the global joint distribution ( , ) P X Y :  [12]) and the goal consists in minimizing the average risk of misclassification [12].This trade-off converges slowly in the case of Markov field because optimization algorithm requires Monte Carlo samplings and a large number of iterations before convergence is reached (simulated annealing algorithm).More recent studies in image analysis have suggested replacing the purely spatial priors used in Markov field with hierarchical priors, in which interaction between variables are not supported by the grid over the image pixels but are defined from scale to scale.Thus, the HHMM modeling scheme captures, over a label pyramidal lattice (quad-tree, Fig. 2(a)), significant interscale statistical dependencies whereas the intrascale dependencies with the observation vector are statistically defined at the finest resolution by a density probability function linking label and observation.This approach, described in [1], [7], may induce block effects due to quad-tree structure.
The usefulness of HMCM stems from their ability to learn Hidden Markov Model parameters through the Baum-Welch re-estimation procedure [9].Firstly, the images are scanned by a Hilbert Peano path [10] (cf.Fig. 2(a)) in order to transform the cube into a chain of Mueller vector, then a Markov chain, defined by its transition matrix, is estimated.The Baum-Welch algorithm is an iterative update algorithm which re-estimates parameters of a given hidden Markov model to produce a new model which has a higher probability of generating the given observation sequence.This re-estimation procedure is continued until no more significant improvement in probability can be obtained: the local minimum is thus found.
The foremost importance of HHMM and HMCM models consists in performing exactly, in a non-iterative way, the Maximum a Posteriori Mode (MPM) inference [11].Such noniterative approaches can also provide substantial gains in speed and result quality.
These models are justified by the fact that for many natural scenes, neighboring pixels are more likely to belong to the same class than pixels that are farther away from each other : this property is translated on a Markov 1D-chain (HMCM, cf.Fig. 2(a)) [10] or in scale on a Markovian quad-tree (HHMM, cf.Fig. 2(b)) [7].One important contribution in this paper consists in showing that these two approaches give segmentation maps of similar quality on real data cubes.From the observed cube Y to the segmented image X, the HHMM and HMCM algorithms can be decomposed into three main phases: Initialization: this phase provides a first estimation of data-driven parameters.K-means algorithm can be used [12].
Segmentation: the label map is then achieved using the Maximum a Posteriori Mode (MPM) segmentation rule [14].Return to step "parameter iteration" until convergence (weak label variations on map X between two iterations).
The number of desired classes for the segmentation map is the single parameter given by the user.The software of complete segmentation chain is available on line on http://voltairemiv.u-strasbg.fr/.m Mueller element image.It corresponds to conventional intensity image.We observe that this image is quite hard to segment since the finest details do not appear clearly (finest veins are not observed).This observation was pointed out in [1].

Experimental results and discussions
Figure 4 While the conventional image (Fig. 3(a)) seems to be unsuitable for segmentation, the animation (normalized to the m 00 element image), Fig. 4(a) shows clearly that the polarimetric response of the object provides the appropriate candidate to the segmentation task.The goal is to segment this data cube into different classes corresponding to physical properties of the leaf.We observe further that the images of m 11 /m 00 , m 22 /m 00 , and m 33 /m 00 , Mueller elements provide a better contrast between the leaves constituents (nervures and tissue).This indicates a smoothing effect in the transmittance heterogeneity over the leaf surface.The segmentation map obtained with HMCM analysis (Fig. 3(b)) reveal the finest details of the object.This can not be reached with classical intensity based imaging.Some authors suggest to reduce the data cube by means of a Principal Component Analysis (PCA) transform [14].PCA consists in projecting the 16 Mueller images on different axes (inner vectors) corresponding to the largest energy.The goal is to project the image on axes maximizing the variance, in order to decorrelate the different parameter bands.Movie of Mueller images after PCA transform is available on-line [15].This transform, used generally for data reduction (i.e., multispectral or hyperspectral imagery [16]) is useless in our case.Indeed, the segmentation map obtained by feeding the Markovian algorithms with such pictures does not improve the result (Fig. 5(a)).One may assume that an energy-based projection is not adapted to Mueller parameters.Another way for the data reduction we explore is based on Independent Component Analysis (ICA) [17].Finally, some results are displayed on-line : HMCM [18] map (Fig. 4(b)) with a different number of classes [19] and HHMM results [20].Similar results for both Markovian algorithms (weak block effect Fig. Once the label maps are obtained from the Markovian procedure, different classes can be merged properly to provide appropriate masks corresponding to nervures and tissue.These masks were used to extract the mean Mueller matrices of nervures N M and tissue T M from the Mueller image cube.We found 1.0000 0.0227 0.0031 0.0028 0.0077 0.2066 0.0038 0.0096 0.0009 0.0121 0.2225 0.0024 0.0035 0.0118 0.0082 0.1306 1.0000 0.0269 0.0021 0.0018 0.0101 0.3236 0.0087 0.0023 0.0008 0.0024 0.3276 0.0009 0.0026 0.0023 0.0029 0.2754  M and T M are found to be nearly equal 0.96 and 0.91, suggesting that the two classes reach a close-to-ideal depolarizer.Figure 6 shows the eigenvalues distribution.We observe a dominant eigenvalue while the three remaining eigenvalues have almost the same order of magnitude.This makes it reasonable to approximate the Mueller matrices as the sum of a main interaction mechanism and an isotropic depolarizer (i.e.0 iso ≈ + M M M ).The Mueller matrix of an isotropic depolarizer has all their elements bring to naught except m 00 ,which is equal to d (depolarizer strength).Its associated coherency matrix is given by 1 2 4 iso d = T I , where I 4 is the 4×4 identity matrix.These two matrices can be constructed using the following equations: where v 0 is the eigenvector corresponding to λ 0 .Table 1 summarizes the results obtained by using Eq. ( 11).These results indicate the following: i) the two classes composing our object depolarize a large amount of incoming polarization states, ii) They act equally on all polarization states in magnitude.However, the handedness and the tilt of the incoming wave are inverted.We point out the performance of the clustering based on Markovian Bayesian inference in view of the very low signal level (10 nm wideband for the interferential filter) and the 20% (nervures) and 30% (tissue) remaining unpolarized fraction of light that escapes the vegetal leaf.

3.3
The use of the polar decomposition for the two Mueller matrices under study shows that the pure depolarizer and pure attenuator matrices are nearly equal to identity matrix.The two classes act primarily as general pure depolarizers which confirm the conclusion of the analysis based on Cloude's theory and validate the Markovian segmentation approach : the obtained classes given by this process have physical significance.

Conclusion
A new analysis framework for Mueller matrices images is introduced which allows an original study of the polarimetric properties of skeletonized shapes.The major interest of this study comes from the association of up-to-date image processing algorithms with physical interpretation of the results.Markovian based algorithms turn out to be of great benefit for such inference tasks: the only image dependent parameter that must be provided by the user is the number of classes, all other parameters are computed automatically.Moreover, this approach is quite general and can be extended to a large variety of samples, e.g.Mueller images of biological tissues.Work in progress pursues two lines of inquiries: i) Bayesian tree-structure image using wavelet domain hidden Markov models, ii) noise features variations and its impact on the clustering algorithms performances.

Fig. 1 .
Fig. 1.Schematic of the Mueller matrix imaging polarimeter used in this study.
which yields three Mueller matrices for a pure diattenuator, D M , a pure retarder, R M , and a pure depolarizer, dep M , related to M by matrix multiplication s y y s S = ∈ and S is the set of voxel positions in the data cube to be segmented.The variables s X and s Y take their values in { } 1 2

Fig. 2 .
Fig. 2. (a) Markov chain model and Hilbert Peano paths for different image sizes (2 2 , 4 2 , 8 2 and 16 2 pixels).The image becomes a vector after Hilbert-Peano scan.Markovian prior model describes the transitions within the chain between labels X n and X n+1 whereas the data-driven parameter links X n+1 and Y n+1 .(b) Markovian quadtree with inter-scale transition probabilities a ij (s -stands for the father of site s in the tree).Data-driven probability density function between observations (white circle) and labels (black circle) equals ( ) ( ) / i s s i f l p y l x ω = = =

Figure 3 (
Figure 3(a) represents the 00m Mueller element image.It corresponds to conventional intensity image.We observe that this image is quite hard to segment since the finest details do not appear clearly (finest veins are not observed).This observation was pointed out in[1].Figure4(a) contains a movie of raw Mueller images (16 parameters on each pixel: ; , 0..3

Fig. 3 .
Figure 3(a) represents the 00m Mueller element image.It corresponds to conventional intensity image.We observe that this image is quite hard to segment since the finest details do not appear clearly (finest veins are not observed).This observation was pointed out in[1].Figure4(a) contains a movie of raw Mueller images (16 parameters on each pixel: ; , 0..3 ij m i j = ).These parametric images form a data cube of size 256 by 256 by 16 of floating values ( ; , 0..3; , 0..255 kl ij m i j k l = =).

Fig. 4 .
Fig. 4.. (a) (1.54 Mb) movie of raw Mueller images (16 parameters for each pixel: ( ; , 0 3 ij m i j = L ).(b) (1.04 Mb) Sequence of HMCM maps from initialization (based on Kmean algorithm) up to label map: convergence is performed in 6 iterations, taking less than 1 mn on a PC Pentium IV, 2 GHz, 4 Gb RAM..Each class is displayed with a random gray level.

Fig. 5 .
Fig. 5. (a) Segmentation map PCA transform (b) HHMM map : convergence is performed in 6 iterations, taking less than 1 mn on a PC Pentium IV, 2 GHz, 4 Gb RAM.One observes weak blocky effects.This label map is similar to HMCM map obtained on Fig. 4(b).