1. 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 p=1,16 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.

2. Experimental and theoretical considerations

2.1 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.8 nm 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).

 

Fig. 1. Schematic of the Mueller matrix imaging polarimeter used in this study.

Download Full Size | PDF

It can be shown that the measured intensity at the location of pixel I_ij (i, j) is a function of the angular positions of quarter-wave plates α ₁ and α ₂. For fixed angle positions of the two quarter-wave plates, ${\alpha }_1^k$ and ${\alpha }_2^l$, this intensity can be written as

where M _ij is the pixel’s (i, j) 4×4 Mueller matrix.

To fully determine M _ij requires 16 intensity measurements that correspond to four positions of both angles ${\alpha }_1^k$ and ${\alpha }_2^l$(k,l=1,4), yielding

where

Appropriate values of ${\alpha }_1^k$ and ${\alpha }_2^l$ (k,l=1,4) 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.

2.2 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 ^-1) its inverse. Cloude showed that a necessary and sufficient condition for (M) to be physically realizable is that all eigenvalues (λ_i ;i=0…3) 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 non-depolarizing 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] which yields three Mueller matrices for a pure diattenuator, M _D , a pure retarder, M _R , and a pure depolarizer, M _dep , related to M by matrix multiplication

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.

2.3 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 y=y_s , s∈S and S is the set of voxel positions in the data cube to be segmented. The variables X_s and Y_s take their values in Ω={ω ₁,ω ₂,…,ω_K } and ▯ ⁿ , respectively. One assumes that each observation y_s is a filtered and noisy measurement of the observed scene :

• the segmentation process tries to determine which label (or class {ω ₁, ω ₂,…,ω_K }) corresponds to each observation Y_s ;

• the classification process tries to determine which label (or class {ω ₁, ω ₂,…,ω_K }) corresponds to the real nature of observation (this may be, nomenclature in remote sensing applications, i.e., forests {ω ₁}, fields {ω ₂}, rivers {ω ₃},… lakes {ω_K }).

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, p_X is the marginal distribution of X, and p(Y/X) is the distribution of Y conditional to X. p(Y/X) is called the data driven term and $p(Y⁄X)={\displaystyle \prod _{s\in S}}p(Y_s⁄X_s)$ where p(Y_s/X_s ) follows a Gaussian distribution. Under Markovian assumption, the X distribution takes into account neighborhood labels (spatially or in scale) on each site in order to provide a form of context handling in the labeling process. This a priori favors homogeneity in the labeling process. Let us assume that the X distribution is markovian, i.e., p(x_s /x_t ,t≠s)=p (x_s /x_t , t∈Λ _s ) then X follows a Gibbs distribution with respect to a neighboring system (Λ _s ) _s∝S , for example the set of the four nearest neighbors of s∈S. 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):

where Z_Y is a normalization coefficient depending on Y. In the Bayesian framework, a cost function defines the segmentation criterion (Maximum A Posteriori, Modes of Posterior Marginals, Mean Field,… [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 non-iterative 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.

 

Fig. 2. (a) Markov chain model and Hilbert Peano paths for different image sizes (2², 4², 8² and 16² 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 f_i (l)=p(y_s =l/x_s =ω_i )

Download Full Size | PDF

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].

Parameter estimation: Both approaches (HHMM or HMCM) use, roughly speaking, the same parameter estimation, which is based on Expectation Method to perform iteratively the maximum likelihood estimation with incompletely observed data[13],[10],[7],[11].

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/.

3. Experimental results and discussions

Figure 3(a) represents the 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].

Figure 4(a) contains a movie of raw Mueller images (16 parameters on each pixel: m_ij ; i,j=0..3). These parametric images form a data cube of size 256 by 256 by 16 of floating values ($m_{\mathit{\text{ij}}}^{\mathit{\text{kl}}}$; i, j=0..3; k,l=0..255).

 

Fig. 3. (a) m ₀₀ Mueller element image. This element image corresponds to a conventional intensity image with which results must be compared. (b) HMCM analysis with 6 classes

Download Full Size | PDF

While the conventional image (Fig. 3(a)) seems to be unsuitable for segmentation, the animation (normalized to the m ₀₀ 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 ₁₁/m ₀₀, m ₂₂/m ₀₀, and m ₃₃/m ₀₀, 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. 5(b)) are observed.

 

Fig. 4. . (a) (1.54 Mb) movie of raw Mueller images (16 parameters for each pixel: (m_ij; i,j=0…3). (b) (1.04 Mb) Sequence of HMCM maps from initialization (based on K-mean 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.

Download Full Size | PDF

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 M _N and tissue M _T from the Mueller image cube. We found

 

Fig. 5. (a) Segmentation map after 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).

Download Full Size | PDF

 

Fig. 6. Eigenvalues of O(M _N ) and O(M _T )

Download Full Size | PDF

Entropy values of, M _N and M _T 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. M≈M ₀+M _iso ). The Mueller matrix of an isotropic depolarizer has all their elements bring to naught except m₀₀ ,which is equal to d (depolarizer strength). Its associated coherency matrix is given by T _iso =½d I ₄, where I ₄ is the 4×4 identity matrix. These two matrices can be constructed using the following equations:

where v₀ is the eigenvector corresponding to λ₀. 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.

Table 1. Main interaction mechanism + isotropic depolarizer results

View Table

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.

4. 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.