Abstract

The recently proposed restoration-segmentation algorithms dedicated to polarization encoded images suffer two important limitations: the number of classes into which the image is segmented is not obtained automatically, and more importantly the quality of the segmentation is affected by the nonuniformity of the illumination of the scene. We propose here a new method addressing these issues. It is based on a global estimation-segmentation approach, explicitly modeling the nonuniform illumination. The physical admissibility of the retrieved Mueller matrices is ensured. Results stemming from synthetic and real data are provided and support the proposed approach.

© 2015 Optical Society of America

1. Introduction

Polarization imaging captures polarization properties of light and/or scenes in spatially distributed measurements. It gives rise to a multicomponent information attached to each pixel. In this context, Mueller imaging is the estimation of the Mueller matrix attached to each pixel in the measured images.

Polarisation-sensitive imaging is a form of optical inspection that can reveal features in a sample that appear invisible to intensity and/or phase detection systems. Many polarization-sensitive imaging systems have been designed for several fields including metrology and medical [1–6 ] applications. This is mainly based on the ability to build an effective Polarization State Analyzer (PSA) that acquires polarized radiances yielding polarization parameters corresponding to each pixel in the image. This is done by using birefringent optical elements such as modulators in the PSA architecture. Furthermore, if one can control the polarization state of the incoming light through a Polarization State Generator (PSG), the system allows acquisition of images from which a Mueller matrix can be estimated at each pixel.

A key issue that arises at early stages in imaging polarimeters design and operation is proper calibration [7–9 ] that allows to define precisely the wavelength-dependent polarization measurement matrix that relates the measured polarized spectral radiances to the sought polarization parameters at different wavelengths and permits to get rid of all types of systematic errors attached to the imaging instrument. Systematic errors will not be considered in the paper since an efficient calibration procedure of polarimetric imaging systems is employed in order to extract the desired polarization-encoded images effectively while minimizing/eliminating the impact of systematic errors.

In this paper we deal with random errors, i.e. noise and with non-uniform illumination. First, the term noise includes different contributions that arise from the use of a pixelated photosensitive device as a detector. To be more specific, this includes readout noise, thermal noise, dark current noise, and photon noise, among others. These quantities are device dependent and could be reduced using a high quality camera (High quantum efficiency, high-performance electronics, etc.) and a proper choice of acquisition parameters (exposure (integration) time, multi-frame averaging, etc.) but can never be totally eliminated. Thus, one needs to account for noise in polarization imaging. Moreover, practical working conditions (in vivo imaging) constrain to consider short time acquisition or narrow band imaging (to have a good conditioning of the polarimetric system) such that the Signal-to-noise ratio (SNR) cannot be as high as desired. As a consequence, the measured polarized radiances are contaminated by noise and may prevent the admissibility of the calculated Mueller matrices ([10] or [11] for the case of Stokes vectors). A common approximation for modeling noise is to consider it as Gaussian, additive, independent at each pixel, and independent of the signal intensity. Even if it is physically unrealistic [12], we assume this model in the following since it has shown to work well with our real data.

Then, we also deal with non-uniform illumination. Illumination uniformity can only be approximated in laboratory measurement conditions since a strictly parallel uniform beam cannot be obtained. Working with a “strictly” parallel uniform beam imposes to use complicated illumination systems and to image only planar samples at a fixed focal distance. Every departure from this ideal situation would necessitate a whole recalibration. A possibility to deal with nonuniform illumination is to normalize the estimated Mueller matrices by the first Mueller element (mean transmittance image) or equivalently to normalize the measured data by the same element. Such an approach has two major defaults. First, it is clear that the spectral polarimetric signature will be hindered by the normalization operation since it eliminates the energy relationship at different wavelengths. Then, under Gaussian assumption of the noise, the elements of the Mueller matrix follow a multivariate normal distribution when using the pseudo-inverse for the estimation. By normalizing each element of the Mueller matrix by the first one, each element of the normalized Mueller matrix follows then the distribution of the ratio of jointly normal variables. The same conclusion holds for the normalized measured data. As explained in [13], such a distribution is a product of two factors, the first being a Cauchy density. This distribution is heavy tailed and can be bimodal. These properties yield inappropriate segmentation as we will show in Sec. 5. To solve this problem, we deal with non-uniform illumination by defining a model that explicitly accounts for illumination variation.

The data acquired using imaging polarimeters are noisy and indirect versions of the parameters of interest. The central issues of any inversion procedure are matching the physical admissibility conditions of the estimated parameters and reducing the effect of noise. The pseudo-inverse solution is particularly ill-suited in this respect: noise is propagated and each estimate is not necessarily physically admissible. One should take advantage of the image structure so as to efficiently reduce the effect of noise. In this framework and addressing both aforementioned issues, we have proposed a specific algorithm that segments the Mueller images in a given number of classes of constant signature [14], and we have also proposed an approach attenuating the effect of noise for spatially varying Mueller parameters, resorting to non-local means (NLM) [10]. The former algorithm assumes that the illumination of the scene is uniform and cannot adequately estimate the number of classes into which the images are segmented, while the latter algorithm does not provide a segmentation map of the Mueller image though this map could be of great interest to analyze the image content. We propose here an appropriate method for the estimation of piecewise constant polarization signatures in Mueller imagery, where nonuniformity in the illumination is explicitly modeled. Non-uniform illumination is a major issue, since it occurs frequently in practice and largely degrades the quality of the segmentation map when ignored in the model.

The article is organized as follows. The intensity images obtained with a Mueller polarimeter in the presence of illumination inhomogeneities are modeled in section 2. Section 3 addresses the estimation of the parameters of the model, in the maximum likelihood framework. In section 4, estimation of the segmentation map resorting to Potts’ model is dealt with. Results are presented in section 5. Conclusions are drawn and perspectives are given in section 6.

2. Modeling of Mueller intensity images under non–uniform illumination

At any pixel location s of the image, the N intensity measurements I(s) (N ≥ 16) are related to the Mueller matrix M(s) by:

I(s)=PM(s)_(+noise),
where the underline symbol is the operator which maps the 4 × 4 matrix M(s) into a 16 × 1 vector formed by stacking up its columns. The known N × 16 matrix P is the Polarization Measurement Matrix (PMM).

As in [14], the objects of the scene are supposed to have a constant Mueller matrix signature, which corresponds to many situations of practical interest. If we account for lighting nonuniformity and if we suppose that pixel s is associated to class k, we may write

I(s)=κ(s)PMk_+ε(s)=κ(s)Ψk+ε(s),
where M k is the Mueller matrix associated to class k and Ψ k is the related ideal intensity observation ( Ψk=PMk_). The noise ε(s) is supposed to follow a Gaussian distribution ε(s) ∼ 𝒩 (0, Σ k) where Σ k is the covariance matrix of class k. The covariance matrix is supposed to be diagonal in the sequel. Note that the model of image noise that is used here corresponds to a common approximation since it is Gaussian, additive, independent at each pixel, and independent of the signal intensity. The term κ(s) models the lighting distribution and is supposed to take real positive values (κ(s) has to be positive for M k to be a Mueller matrix). We consider κ(s) ∼ 𝒰 (0, +∞), where 𝒰 is the uniform distribution. This prior is not a proper distribution but the related posterior distribution will be a proper one.

With this model, M k can only be determined up to a multiplicative constant. To remove this ambiguity, M k is normalized such that its upper-left element m 00 is set to 1. The term κ(s) accounts also for the mean intensity transmittance.

The model is parametrized as follows: K denotes the number of classes and we define Θ = {Θ k}k=1...K = {πk, M k, Ψ k, Σ k}k=1...K. The term πk is the prior of class k, with πk = P(y(s) = k|Θ), where y(s) = k means that pixel s is associated to class k. Since M k and Ψ k are linked by the relation Ψk=PMk_, one of the two variables suffices but both are introduced for the sake of clarity (equations may be written considering either M k or Ψ k).

By using the definition of conditional probability and the assumption that y(s) and κ(s) are conditionally independent given Θ, we have :

P(I(s),y(s)=k,κ(s)|Θ)=P(I(s)|y(s)=k,κ(s),Θ)P(y(s)=k|Θ)P(κ(s)|Θ).

As a reminder, the second term on the right hand side of Eq. (3) is πk and the third one is modeled with κ(s) ∼ 𝒰 (0, +∞). Thanks to Eq. (2) and the related assumptions (ε(s) ∼ 𝒩 (0, Σ k)), the first one can be computed using:

I(s)|y(s)=k,κ(s),Θ~𝒩(κ(s)Ψk,Σk).

The likelihood at pixel s can be obtained by marginalizing Eq. (3) over κ(s) and y(s):

P(I(s)|Θ)=k=1Kκ(s)P(I(s),y(s)=k,κ(s)|Θ)dκ(s),
from which we derive (thanks to Eq. (3))
P(I(s)|Θ)=k=1Kκ(s)πkP(κ(s)|Θ)P(I(s)|y(s)=k,κ(s),Θ)dκ(s).
Finally, under pixel-independence assumptions considered in Eq. (8), the likelihood L(Θ) writes:
L(Θ)=P(I|Θ)=sP(I(s)|Θ).
The expression of the likelihood as a product of pixelwise factors stems from:
P(I,y,κ|Θ)=P(I,y,κ|Θ)P(y,κ|Θ)=[sP(I(s)|y(s),κ(s)|Θ)]P(y|Θ)P(κ|Θ)=[sP(I(s)|y(s),κ(s),Θ)][sP(y(s)|Θ)][sP(κ(s)|Θ)]=sP(I(s),y(s),κ(s)|Θ),
with y = {y(s)} and κ = {κ(s)}, and where the justification is completed by a marginalization over y and κ.

3. Estimation of the model parameters and of the number of classes

The estimation of the model parameters and the number of classes is achieved through an iterative process as depicted in Fig. 1. We start with one class (K = 1). Variants of the Expectation Maximization (EM) algorithm are used to estimate the parameters of the model (subsection 3.1). The initialization of the EM algorithm is addressed in appendix A.1. The number of classes is increased by splitting a class in two if at least one class among the K classes does not meet a criterion based on statistical tests described in subsection 3.2. If the number of classes is to be increased, the class with the worst criterion is split. Appendix A.2 describes how a model with K + 1 classes is initialized from the model with K classes. The parameters are then estimated as previously with the EM algorithm. The procedure is repeated until all classes meet the criterion.

 figure: Fig. 1

Fig. 1 Steps related to the estimation of the parameters Θ and of the number K of classes.

Download Full Size | PPT Slide | PDF

3.1. Maximum likelihood parameter estimation

The model parameters are estimated using maximum likelihood. Optimization of the likelihood (Eq. (7)) is cumbersome, and is achieved using the EM algorithm, which is a widely used procedure for estimating Gaussian mixture densities [15, 16]. It introduces an auxiliary function Q(Θ, Θ′), such that if Q(Θ, Θ′) ≥ Q(Θ′, Θ′), then L(Θ) ≥ L(Θ′). Based on this property, the EM algorithm proceeds by defining recursively a series of models {Θ 1, Θ 2,...} verifying Q(Θ i+1, Θ i) ≥ Q(Θ i, Θ i) and thereby L(Θ i+1) ≥ L(Θ i).

In the present case, the auxiliary function writes (see Appendix B.1)

Q(Θ,Θ)=sk=1Kwk(s)κ(s)=0+log[P(I(s),y(s)=k,κ(s)|Θk)]×P(κ(s)|y(s)=k,I(s),Θ)dκ(s),
with w′k(s) = P(y(s) = k|I(s), Θ′).

To the best of our knowledge, the M step (ie, the maximization of Q(Θ, Θ′) with respect to Θ) of the EM algorithm is not tractable because of the integration in κ(s) and the subsequent optimization of Q. To tackle this problem, we resort to a stochastic EM (SEM) algorithm [17]. In our case, the SEM works as follows: for each pixel s and for each class k, we generate the missing data κk(s) according to its posterior distribution [18]:

κk(s)|y(s)=k,I(s),Θ~𝒩+(αkβk,αk),
where 𝒩 + is a truncated normal distribution (the support of 𝒩 + is ℝ+) and where α′k and β′k are defined in Appendix B.1. We get:
Θ^=argmaxΘsk=1Kwk(s)log[P(I(s),y(s)=k,kκ(s)|Θk)]=argmaxΘsk=1Kwk(s)log[πkP(I(s)|Θk,y(s)=k,kκ(s))]
As in the case of the Gaussian mixture density case, we get:
π^k=1#pixelsswk(s),
and have to solve
argmin{Ψ,Σ}swk(s)[log(|Σ|)+(I(s)κk(s)Ψ)TΣ1(I(s)κk(s)Ψ)].

As an under-optimal solution to this optimization problem, we finally define (see the rationale and the computation in Appendix B.2)

I¯k=1swk(s)κk2(s)swk(s)κk(s)I(s).
We retain:
M^k=argminMΣk1/2(I¯kPM)2,
and
Σ^k=(1swk(s)swk(s)(I(s)κk(s)Ψ^k)(I(s)κk(s)Ψ^k)T)IN×N,
where ⊙ stands for the Hadamard product and where I N×N is the identity matrix of size N × N. The cost function in Eq. (15) is optimized with the approach developed in [10] under physical admissibility constraints (a Mueller matrix is admissible if and only if the eigenvalues of the related coherency matrix are positive or null [19]). Finally, at the end of each step of the EM algorithm, the matrices k are normalized so that their mean intensity transmittance are normalized (m 00 = 1).

3.2. Criterion for potentially increasing the number of classes

As shown in Fig. 1, the number of classes is increased by one if at least one class among the K classes does not meet a certain criterion. For a given class k, the residual k(s) at pixel s writes:

r^k(s)=I(s)κ^k(s)Ψ^kN
where κ̂k(s) maximizes the posterior distribution defined Eq. (10). For a given class k, we have N distributions of scalar residuals, and N values of the criterion are computed. If the number of classes is too small to explain the data, we expect at least one distribution out of N to be non-Gaussian or non-unimodal (for at least a class).

Let pk,n denote the p-value associated to the test of the n-th residual component of class k. Two different tests are used: we test either the normality using the Jarque–Bera test [20], or the unimodality using the dip test [21]. From class k point of view, all pixels do not share the same weights wk(s) = P(y(s) = k, |I(s), Θ). To account for this, we use the strategy proposed in [22]: each pixel is associated to a class, the probability of choosing class k for pixel s being wk(s). The class being drawn randomly, the residual is computed for each pixel. For both tests, the p–value is uniformly distributed under the null hypothesis 0 (all residual components are Gaussian or all residual components are unimodal, depending on the test considered) and is expected to be small if the null hypothesis is false. The choice to increase (or not) the number of classes is made using the mimimum p–value pmin:

pmin=min{pk,n}k=1K,n=1N.

The problem of multiple comparisons is tackled here with the well–known Bonferroni correction. The null hypothesis is rejected with the α–level if

pmin1(1α)1NK.
Under the hypothesis that all statistical tests are independent, the probability to reject wrongly the null hypothesis is α. If the null hypothesis is rejected, the number K of classes is increased by splitting a class in two. Appendix A.2 describes how a model composed of K + 1 classes is initialized from a model of size K.

4. Mueller image segmentation

Spatial regularization is performed, as in [14], by modeling the segmentation map as a Markov random field with a Potts model of parameter β. The prior probability of the segmentation map I 0 writes:

P(I0|β)=1Z(β)exp(βU(I0)),
where
U(I0)=[t,s]𝕀(I0(t)=I0(s))
is the energy associated to the field and Z is the partition function. Variables t and s denote pixel sites of the image, and [t, s] denotes summation over neighboring sites (we used eight neighbors in the present work). 𝕀(.) is the indicator function taking values 1 or 0 whether the condition between parentheses is true or false.

The segmentation of I can be achieved as follows:

{I^0,β^}=argmax{I0,β}logP(I0,β|I,Θ)=argmax{I0,β}logP(I0|β)+slogP(I(s)|I0(s),Θ),
where a float prior is taken for β. To address this optimization problem, we resort to an alternating optimization approach. It consists in building a series of segmentation maps {I0(1),I0(2),} with the associated set {β (1), β (2),...} of β’s until convergence. The segmentation map I0(n) is determined from Eq. (22) by setting β equal to β (n), whereas β (n) is determined from Eq. (22) by setting I 0 equal to I0(n1). For the initialization, β (1) is set to 0 (no regularization). We explain hereafter how Eq. (22) can be optimized, when either β or I 0 is fixed.

Given β, the criterion of Eq. (22) can be optimized using graph–cut algorithms [23, 24, 25]. The last term of Eq. (22), namely P(I(s)|I 0(s), Θ), is computed using:

{P(I(s),κ(s)|I0(s),Θ)=P(I(s)|I0(s),κ(s),Θ)P(κ(s)),I(s)|y(s)=k,κ(s),Θ~𝒩(κ(s)Ψk,Σk),P(I(s)|I0(s),Θ)=P(I(s),κ(s)|I0(s),Θ)dκ(s),
and the integral in the previous equation can be computed exactly (see Eq. (32)).

Given I 0, the estimation of β is achieved in a manner similar to the one in [14].

5. Algorithm applications

5.1. Results on synthetic data

5.1.1. Evaluation protocol

We synthesized a hard 512 × 512 pixel Mueller image test case according to Eq. (2). The number of classes has been set to 4, each class being assigned to a region of the image (Fig. 2(a)).

 figure: Fig. 2

Fig. 2 (a) Simulated segmentation map, (b) simulated lighting map, and (c) example of a mean transmittance image obtained for a 10 dB SNR. Representative segmentation maps obtained (SNR=10 dB) with the proposed approach (d), and with the approach of [14] (without normalization (e), by using the pseudo-inverse for normalization (thresholding (f) and polynomial modeling (g)), and by using the ground truth for normalization (h).

Download Full Size | PPT Slide | PDF

The Mueller matrix of class 1, M 1, has been generated using the λ-parametrization [14] which assures a valid Mueller matrix parametrized with 16 real parameters {λi}i=1...16. Each λi related to matrix M 1 has been drawn according to a normal distribution with standard deviation of 1. The parametrizations of M 2, M 3 and M 4 have been deduced from the one of M 1 as {λ 1, λ 2, λ 3, λ42, λ 5 ...λ 16}, {λ 1, λ 2, λ 3, 0, λ 5 ...λ 16}, {λ 1, λ 2, λ32, 0, λ 5 ...λ 16}, respectively. Such a procedure provides Mueller matrices with similar characteristics. The corresponding Mueller matrices are displayed in Table 1.

Tables Icon

Table 1. Mueller matrices used for the creation of the synthetic data.

The PMM P of Eq. (2) has been defined as the PMM of an optimal rotating–retarder Mueller polarimeter with the optimal parameter vector ηopt as defined in [9]. The non–uniform lighting map κ(s) has been defined as an isotropic truncated Gaussian distribution, with values varying from 1 to 5 (see Fig. 2(b)).

It can be noticed that the term κ(s) being estimated with the proposed approach cannot be compared with the one used for generating the data since the element m 00 of each Mueller matrix is set to 1 during the estimation, so that the estimated κ(s) accounts for the mean intensity transmittance as well as for the lighting.

Finally, a zero mean white Gaussian noise has been added to the signal of interest with a diagonal covariance matrix {1σ2Σk}k=14 depending on the class k. Each diagonal element of {Σk}k=1...4 has been drawn according to a uniform distribution on the unit interval. The standard deviation σ is estimated so as to correspond to a signal to noise ratio (SNR) of fixed value computed on the mean transmittance image (first Mueller channel). Figure 3 represents the 16 intensity images obtained for a 10 dB SNR.

 figure: Fig. 3

Fig. 3 Exemple of synthetic intensity images obtained with a SNR of 10 dB.

Download Full Size | PPT Slide | PDF

The proposed algorithm has been tested on data obtained with a 10 dB and a 3 dB SNR. For each SNR, 100 different datasets have been generated with different noise realizations. Each dataset has been processed several times, using three different α-levels for the statistical tests (5%, 1%, and 0.1%), and two different tests, the Jarque–Bera test [20] and the dip test [21] (Sec. 3.2). We propose to evaluate (i) the ability of the method to detect the correct number of classes, and (ii) the accuracy of the estimation of the Mueller matrices and the quality of the segmentation maps. It may be noted that segmentation of these images is a challenging task because the different classes have similar Mueller matrices, the SNR is low, and the lighting is not uniform.

5.1.2. Estimation of the number of classes

The number K of classes estimated over 100 realizations with the Jarque–Bera test is given in Tab. 2 for three different α’s and two different SNR values.

Tables Icon

Table 2. Number of occurrences of the estimated number K of classes for 100 experiments, for different values of α and of the SNR, with the Jarque-Bera test.

For iterations where K is less than the expected number of classes, we note the ability of the Jarque–Bera test to detect the non–Gaussianity of the residuals. Indeed, in all conducted experiments, the estimated number of classes is never less than the true value.

The number of type I errors ( 0 is erroneously rejected leading to an estimation of K higher than 4) should be approximately 5, 1, and 0.1 for α–levels respectively equal to 5%, 1%, and 0.1%. This means that the effective rate of the type I errors is slightly higher than expected. This comes from the fact that the estimated residuals do not strictly adhere to a Gaussian distribution, because of slight error estimation in κ(s) and M k and subsequent error propagation.

A similar study has been conducted with the dip test. This test evaluates the modality of a distribution (unimodal vs. multimodal) and is used for splitting classes. In this case, the algorithm returns a unique class for every experiment. The number of classes is underestimated due to the inability of the test to reject 0. This result is not surprising because the generated Mueller matrices are very similar and the addition of noise leads to a unimodal histogram (adding a Gaussian noise comes down to convolve the noise-free histogram with a Gaussian kernel) for I, and also for the histogram of the estimated residuals. When using a higher SNR or when dealing with less similar Mueller matrices, the dip test provides a better rejection of 0, leading to a better estimation of the number of classes. Finally, as could have been expected, since a Gaussian distribution is unimodal while a unimodal distribution may not be Gaussian, the Jarque–Bera test provides a better rejection of 0 than the dip test, with the risk of increasing the type I error since there is no guarantee that the estimated residuals are perfectly Gaussian.

5.1.3. Accuracy of the parameter estimation and quality of the segmentation map

The α–level has no influence on the parameter estimation but merely on the choice of the number of classes. In the following, we consider only results that have been obtained with an α–level of 0.1% and that have a proper estimation of the number of classes.

For each class k, the relative estimation errors REMk and RE Σ k, are computed as :

REMk=M^kMkMk,REΣk=Σ^kΣkΣk,
where ‖.‖ is the Frobenius norm. In Eq. (24) M k is normalized so that its left upper element is 1. Both quantities in percent are shown in Table 3 for the proposed method.

Tables Icon

Table 3. Mean values of the relative estimation error RE (in percent) for the four Mueller matrices and the corresponding covariance matrices. Averaging has been done over 100 realizations. Numbers between parentheses represent standard-deviation of related quantities.

We can first observe that the model parameters are well–estimated. This validates the estimation procedure, and in particular the EM algorithm. Note that the relative estimation errors RE Σk do not vary much with the SNR. This is not surprising since this is a relative estimation error and reducing the SNR leads to an increase in the norm of Σ k.

Finally, the quality of the segmentation map can be assessed with the number of misclassified pixels. The average number of misclassified pixels for the proposed approach is about 27 pixels (for 3 dB) and 18 pixels (for 10 dB), among 512 × 512 = 262 144 pixels (less than 10−2 %). An example of segmentation map is shown in Fig. 2(d) for a 10 dB SNR. The quality of the results stems from the accuracy of the stochastic model proposed and the related parameter estimation scheme.

We propose to compare this approach with the one of [14]. The segmentation map obtained with [14] is shown in Fig. 2(e) (in [14], the number of classes is user-defined). We observe that the first class has been split in three, because of illumination variation over this class. The segmentation proposed in [14] works well if the lighting map is constant. One can think to normalize the data in order to make the lighting map constant, and then to use the algorithm of [14] on the normalized data. Normalization could be achieved by estimating the mean transmittance image m 00 with the pseudo-inverse approach. The normalized data are computed as Inorm(s)=I(s)m00(s), and the observation model (Eq. (2)) then writes:

Inorm(s)=κ(s)m00(s)PMk_+ε(s)m00(s),PMk_+ε(s)m00(s).
In order to prevent the division by values close to 0, several thresholding procedures have been tested for m 00. In some tests, all values of m 00 lower than a threshold have been set to this threshold. In other tests, all positive values of m 00 lower than a positive threshold have been set to the threshold, whereas all negative values of m 00 larger than a negative threshold have been set to this threshold. In all cases, results obtained are not satisfactory, with relative estimation errors varying from 17% to 85%, and with the number of misclassified pixels varying from 37% to 70% (see for example Fig. 2(f)). The poor quality is due to the fact that the estimation of m 00(s) is corrupted by a Gaussian noise when using the pseudo–inverse solution. As an example, we can show that the standard deviation related to the estimation of m 00(s) (for each s) is equal to 1.97 for 3 dB and 0.94 for 10 dB, whereas the ground truth (the generated m 00 channel) varies from 1 to 5 approximately. Consequently, the normalization of the model, as written in Eq. (25), leads to a new residual ε(s)m00(s), which is no longer distributed according to a Gaussian distribution, but to a Cauchy-like distribution (please see Introduction), and the model in [14] is no longer appropriate.

Finally, spatial information could be used to obtain a better estimation of m 00(s) (or κ(s)), an so, to achieve a better normalization of the data. We tried to approximate m 00 using polynomials and to normalize the data with the ground truth κ(s). The segmentation maps obtained with polynomial approximation and with ground–truth normalization are presented in Fig. 2(g) and 2(h), respectively. In both cases, the background was split into three classes, and the three objects (triangle, rectangle and circle) are associated to a unique class. The three classes associated with the background have similar Mueller matrices, but with three different covariance matrices that account for the non–stationarity of the noise. Indeed, in the case of ground-truth normalization, Eq. (25) writes:

Inorm(s)=PMk_+ε(s)κ(s),
with varying κ(s). This clearly demonstrates that segmentation of Mueller matrices in the context of non-uniform illumination cannot be achieved by a two step approach relying on a normalization step and a segmentation step. This validates also the proposed approach that accounts for lighting non–uniformity directly in the core of the segmentation process.

5.2. Results on real data

The proposed model is applied to real images acquired by a dual rotating quarter–wave plate Mueller imaging polarimeter in a transmission configuration. A sketch of the object that has been used to carry out these measurements is represented in Fig. 4(a). The object is made with a transparent thin film (cellophane) on which we placed three dichroic patches (Polaroid) and a piece of triangular paper. The first dichroic (dichroic 1 in Fig. 4(a)) is approximately orthogonally oriented to the other two (dichroic 2 and 3). The acquired images are shown in Fig. 5: we can see the different elements previously presented as well as the presence of a halo-like background due to non–uniform lighting of the object.

 figure: Fig. 4

Fig. 4 Sketch of the scene composition (a), automatic segmentation with the proposed approach using the Jarque–Bera test (7 classes) (b), and using the dip test (c).

Download Full Size | PPT Slide | PDF

 figure: Fig. 5

Fig. 5 The 16 intensity images of the object acquired by a rotating quarter–wave plate Mueller imaging polarimeter in transmission configuration.

Download Full Size | PPT Slide | PDF

We have applied the proposed algorithm with the two statistical tests. When using the Jarque–Bera test, the proposed approach fails to converge to the correct number of classes. The estimated number of classes reaches the maximal authorized number of 20 and the algorithm is stopped. In this case, the main elements are found, but additional classes are created at the edges of the polaroids and on the background. The Gaussianity criterion is not satisfactory for the presented images: first, as written in Section 5.1.2, the distribution of lighting can yield non-Gaussian estimated residuals. In the case of real imaging, other reasons may lead to the non-Gaussianity of the residuals such as the intrinsic variability of a class or systematic errors that affect the PMM.

With the proposed algorithm, we can easily obtain a segmentation map for each K. Figure 4(b) represents the segmentation map with 7 classes. The main elements (except the cellophane) are relatively well-segmented.

When using the dip test, the proposed approach converges into the correct number of classes (K = 5). The segmentation map is presented in Fig. 4(c). Results obtained are in good accordance with the known physical composition of the scene (Fig. 4(a)). One can observe misclassification for very few pixels at the border of the patches due maybe to physical deformations (polaroids were cut with scissors, so that they have been distorted locally). Moreover, the very low level of illumination at the right border of the image makes the segmentation very difficult so that a part of the cellophane is considered as paper (opaque screen).

In the following, we compare the estimated Mueller matrices of the dichroic patches (obtained by the proposed approach with the dip test) with the matrices of a linear polarizer. Since the patches are placed on the cellophane, it is necessary to remove the contribution of cellophane by right–multiplying the estimated matrices by the inverse of the Mueller matrix associated to the cellophane. Table 4 presents for each dichroic patch the estimated Mueller matrix (after removal of the contribution of the cellophane), and the Mueller matrix of the related linear polarizer (its orientation angle ζ has been estimated using a least square approach).

Tables Icon

Table 4. Mueller matrices associated to dichroic patches and to their related linear polarizer.

As expected, the estimated Mueller matrices are very similar to those of a linear polarizer. The relative error (with the Frobenius norm) between the estimated matrix and the one of a linear polarizer is about 3% for the first class, 4% for the second and 2% for the last one. The dichroics 2 and 3 are about in the same orientation (118.7° and 123.3°), whereas the dichroic 1 (28.7°) is approximately orthogonal to the other two (118.7° − 28.7° = 90° and 123.3° − 28.7° = 94.6°).

6. Conclusion

We have presented a segmentation algorithm dedicated to Mueller images in the presence of severe illumination inhomogeneities. These images are modeled with a Gaussian mixture distribution. Lighting non-uniformity is incorporated in the model. The parameters are estimated in the maximum likelihood sense with an EM algorithm. The EM algorithm is used in an iterative framework based on three steps: an estimation step, a statistical test, and if needed, the splitting of a class in two. This framework allows the estimation of the number of classes. Results obtained with synthetic and real data validate the proposed approach. They clearly demonstrate that segmentation of Mueller matrices in the context of lighting non–uniformity cannot be achieved by a two step approach that relies on a normalization step followed by a segmentation step.

A. Initialization of the EM algorithm

This appendix addresses the initialization of Θ, for K = 1 and for K ≠ 1.

A.1. Initialization for K = 1

For K = 1 and with Θ 1 = {π 1 = 1, M 1, Ψ 1, Σ 1}, κ̂(s) is first set equal to the upper-left element of the Mueller matrix estimated at pixel s using the pseudo–inverse estimate. If κ̂(s) is smaller than a threshold (2.10−16), κ̂(s) is set to this threshold. Ψ 1 and M 1 are then computed as follows:

Ψ1=sI(s)sκ^(s)andM1=(PtP)1PtΨ1.
Matrix M 1 may not be a Mueller matrix and Ψ 1 may not be equal to PM 1 (if N > 16). These constraints will be taken into account in the EM algorithm. Matrix Σ 1 is then estimated from the residuals I(s) − κ̂(s)Ψ 1.

A.2. Initialization from a previous model

We explain here how a model with K + 1 classes and Θ = { Θ k}k=1...K+1 = {πk, M k, Ψ k, Σ k}k=1...K+1 can be derived from a model Θ(p)={Θk(p)}k=1K={πk(p),Mk(p),Ψk(p),Σk(p)}k=1K with K classes.

Let k and n be such that pmin = p k,n (see Eq. (18)). The initialization of the new model is achieved by splitting the class k of Θ (p) in two so that Θk=Θk(p) for k = 1...K and kk . The parameters Θ K+1 and Θ k are computed in two steps.

First, the entire set of pixels is split in two classes by thresholding the n -th component of the residuals k (s). Several methods can be used for splitting. We used the one proposed in [22].

The threshold t is estimated such that the evaluation at t of the difference between the empirical cumulative distribution function (cdf) of the n -th component of the residual r^k(s) and the associated Gaussian cdf (its parameters are estimated by the sample mean and variance) is as large as possible. Since all points do not contribute equally to the estimation of the cdfs, both estimations have to account for the weights wk(p)(s)=P(y(s)=k|I(s),Θ(p)). Comparing the n -th component of the residual k (s) with the threshold t enables finally to create two clusters C K+1 and C k of pixels, one with residuals larger than t, one with residuals lower than or equal to t.

Finally, the parameters Θ K+1 and Θ k are determined in a manner similar to the case k = 1, with three main differences. First, the whole set of pixels is not used for the estimation: only pixels of C K+1 and C k are used for the estimation of Θ K+1 and Θ k, respectively. Then, weighted means and variances are considered for the estimation using wk(p)(s) as the weights. This enables to reduce the contribution of pixels that belong to other classes. Finally, the priors are set as follows:

πK+1=πk(p)=sCK+1wk(p)(s)sCK+1Ckwk(p)(s),andπk=πk(p)sCkwk(p)(s)sCK+1Ckwk(p)(s)

B. Details of the computations at stake in the EM algorithm

In this section, we detail the computations yielding the different expressions used in the EM algorithm.

B.1. Derivation of the auxiliary function

This subsection addresses the computation of the auxiliary function Q(Θ, Θ′). For Gaussian mixture densities estimation, the auxiliary function is based on unobserved data y(s) associating each data to a class. Equality y(s) = k means that the sample at pixel s corresponds to the k–th class. In our case, we also consider the lighting non-uniformity κ(s) at pixel s as an unobserved data.

The complete log–likelihood (ie, the likelihood of the observed and unobserved data) writes:

logP(I,y,κ|Θ)=slog[P(I(s),y(s),κ(s)|Θ)].
The auxiliary function is expressed as:
Q(Θ,Θ)=yκlog[P(I,y,κ|Θ)]P(y,κ|I,Θ)dκ.

In Eq. (30), y stands for a segmentation map associating each pixel to its class. The sum over y corresponds to a sum over KM segmentation maps, where M is the number of pixels. In the same way, κ(s) represents a continuous map associating a value κ to each pixel. The symbol κ is consequently a multiple integration over M variables. Given Θ′, we aim at determining the optimizer of Q(., Θ′). To achieve this goal, Eq. (30) has to be simplified. According to computations detailed in [16], we get:

Q(Θ,Θ)=sk=1Kκ(s)=0+log[P(I(s),y(s)=k,κ(s)|Θk)]P(y(s)=k,κ(s)|I(s),Θ)dκ(s).
We then decompose
P(y(s)=k,κ(s)|I(s),Θ)=P(κ(s)|y(s)=k,I(s),Θ)P(y(s)=k|I(s),Θ),
and extract the rightmost factor of the right-hand side out of the integral of Eq. (31). This factor may be computed as:
wk(s)=P(y(s)=k|I(s),Θ)P(y(s)=k,I(s)|Θ)P(I(s)|y(s)=k,Θ)P(y(s)=k|Θ)[P(I(s),κ(s)|y(s)=k,Θ)dκ(s)]πkπkκ(s)=0+P(I(s)|y(s)=k,κ(s),Θk)dκ(s)πkαk|Σk|exp(αkβk2γk2)(1+erf(βkαk2)),
with
{αk=(ΨkTΣk1Ψk)1βk=ΨkTΣk1I(s)γk=I(s)Σk1I(s).
The proportionality coefficient in Eq. (32) can be determined with the relation: k=1Kwk(s)=1. Finally, the auxiliary function writes:
Q(Θ,Θ)=sk=1Kwk(s)κ(s)=0+log[P(I(s),y(s)=k,κ(s)|Θk)]×P(κ(s)|y(s)=k,I(s),Θ)dκ(s).

B.2. Suboptimal solution for the SEM cost function

By setting the gradient of Eq. (13) with respect to Ψ equal to the null vector, we get the optimal value of Ψ k denoted Ī k:

I¯k=1swk(s)κk2(s)swk(s)κk(s)I(s).
Ī k enables to reduce the cost function of Eq. (13) as much as possible, but there is no guarantee that Ī k is an acceptable solution. To be an acceptable solution, there must exist a Mueller matrix M k verifying Ψk=PMk_. To overcome this problem, we show in Appendix B.3 that the criterion of Eq. (13) can be written:
{M^k,Σ^k}=argminΘ={M,Σ}(swk(s))log(|Σ|)+(swk(s)κk2(s))Σ1/2(I¯kPM)2.
To our knowledge, there is no closed–form expression of the minimizer of Eq. (36). However, there is a variant of the EM algorithm (the Generalized EM) that consists in determining Θ so as to increase Q(Θ, Θ′) (but Θ is not necessarily the optimizer of Q(., Θ′)) [16]. Consequently, we propose to determine k from Eq. (36) using Σ′ k as a surrogate value to Σ k, and then to determine Σ̂ k from k. Using Σ′ k as a surrogate value to Σ k, the criterion of Eq. (36) writes:
M^k=argminMΣk1/2(I¯kPM)2,
whose solution (pseudo–inverse, PI) is (PtΣk1/2P)1PtΣk1I¯k. If the PI solution is a Mueller matrix (the eigenvalues of the related coherency matrix are positive or null [19]), k is set to the PI solution. Otherwise, criterion of Eq. (37) is optimized using a constrained optimization procedure described in [10].

Finally, since Σ k is supposed to be diagonal, Σ̂ k can be determined as follows from the value of k (and consequently Ψ̂ k):

Σ^k=(1swk(s)swk(s)(I(s)κk(s)Ψ^k)(I(s)κk(s)Ψ^k)T)IN×N,
where ○ stands for the Hadamard product and where I N×N is the identity matrix of size N × N. Finally, at the end of each step of the EM algorithm, the matrices M k are normalized so that their mean intensity transmittance are normalized.

B.3. Derivation of the criterion of Eq. (36)

We detail here the derivation of Eq. (36) (see Appendix B.2). We note:

ck(Ψ,Σ)=swk(s)(I(s)κk(s)Ψ)TΣ1(I(s)κk(s)Ψ).

Since Σ is a diagonal matrix, we have:

ck(Ψ,Σ)=swk(s)κk2(s)(I(1)(s)Ψ(1))T(I(1)(s)Ψ(1))),
with I(1)(s)=Σ1/2I(s)κk(s) and Ψ (1) = Σ −1/2 Ψ.

If

D(s)=wk(s)κk2(s)swk(s)κk2(s)IN×N,
where I N×N is the identity matrix of size N × N, we can write:
ck(Ψ,Σ)=(swk(s)κk2(s))s(I(1)(s)Ψ(1))TD(s)(I(1)(s)Ψ(1)).

By denoting P (1) = Σ −1/2 P, we get Ψ (1) = Σ −1/2 Ψ = Σ −1/2 PM = P (1) M and:

s(I(1)(s)P(1)M)TD(s)(I(1)(s)P(1)M)=I¯P(1)M2(+constantterm),
with I¯=sD(s)I(1)(s). This is due to the fact that the two terms of Eq. (43) have the same gradient according to M (because sD(s) is the identity matrix and the matrices D(s) are diagonal). Finally, based on the definition of Ī k (see Eq. (35)), we have finally:
ck(Ψ,Σ)=(swk(s)κk2(s))Σ1/2(I¯kPM)2(+constantterm)
The criterion of Eq. (36) is then easily derived from Eq. (13) using Eq. (44).

Acknowledgments

These developments were achieved in the framework of the “Imagerie polarimétrique pour le diagnostic et le suivi du traittement du cancer du col utérin” project supported by the Institut National du Cancer.

References and links

1. D. Miyasaki, 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). [CrossRef]  

2. J. Chung, W. Jung, M. J. Hammer-Wilson, P. Wilder-Smith, and Z. Chen, “Use of polar decomposition for the diagnosis of oral precancer,” Appl. Opt. 46, 3038–3045 (2007). [CrossRef]   [PubMed]  

3. M. F. G. Wood, N. Ghosh, X. Guo, and I. A. Vitkin, “Towards noninvasive glucose sensing using polarization analysis of multiply scattered light,” Chap. 17, Handbook of Optical Sensing of Glucose in Biological Fluids and Tissues, V. V. Tuchin, ed., Series in Medical Physics and Biomedical Engineering, Vol. 12 (Taylor and Francis Publishing, 2008).

4. M. F. G. Wood, N. Ghosh, M. A. Wallenburg, S. H. Li, R. D. Weisel, B. C. Wilson, R. K. Li, and I. A. Vitkin, “Polarization birefringence Measurements for characterizing the myocardium, including healthy, infracted, and stem cell treated regenerating cardiac tissues,” J. Biomed Opt. 15, 047009 (2010). [CrossRef]  

5. P. Shukla and A. Pradhan, “Mueller decomposition images for cervical tissue: Potential for discriminating normal and dysplastic states,” Opt. Express 17, 1600–1609 (2009). [CrossRef]   [PubMed]  

6. 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, 1582–1593 (2011). [CrossRef]   [PubMed]  

7. J. S. Tyo, “Design of optimal polarimeters: maximization of signal-to-noise ratio and minimization of systematic errors,” Appl. Opt. 41619–630 (2002). [CrossRef]   [PubMed]  

8. S. N. Savenkov, “Optimization and structuring of the instrument matrix for polarimetric measurements,” Opt. Eng. 41, 965–972 (2002). [CrossRef]  

9. J. Zallat, S. Anouz, and M. P. Stoll, “Optimal configurations for imaging polarimeters: impact of image noise and systematic errors,” J. Opt. A-Pure Appl. Op. 8, 807–814 (2006). [CrossRef]  

10. S. Faisan, C. Heinrich, G. Sfikas, and J. Zallat, “Estimation of Mueller matrices using non-local means filtering,” Opt. Express 21, 4424–4438 (2013). [CrossRef]   [PubMed]  

11. J. Zallat and C. Heinrich, “Polarimetric data reduction: a Bayesian approach,” Opt. Express 15, 83–96 (2007). [CrossRef]   [PubMed]  

12. S. W. Hasinoff, “ Photon, Poisson Noise,” in Computer Vision, Katsushi Ikeuchi, ed. (SpringerUS, 2014) [CrossRef]  

13. A. Cedilnik, K. Košmelj, and A. Blejec, “The distribution of the ratio of jointly normal variables,” Metodološki Zvezki – Advances in Methodology and Statistics 1, 99–108 (2004).

14. J. Zallat, C. Heinrich, and M. Petremand, “A Bayesian approach for polarimetric data reduction: the Mueller imaging case,” Opt. Express 16, 7119–7133 (2008). [CrossRef]   [PubMed]  

15. G. McLachlan and T. Krischnan, The EM Algorithm and Extensions, 2nd edition (Wiley, 2008). [CrossRef]  

16. J. Bilmes, “A gentle tutorial of the EM algorithm and its application to parameter estimation for gaussian mixture and hidden markov models,” Tech. rep., International Computer Science Institute (1998).

17. G. Celeux and J. Diebolt, “The SEM algorithm: a probabilistic teacher algorithm derived from the EM algorithm for the mixture problem,” Comp. Statis. Quaterly 2, 73–82 (1985).

18. W. H. Greene, Econometric Analysis, 7th edition (Prentice Hall, 2012), Chap. 9.

19. S. Cloude and E. Pottier, “Concept of polarization entropy in optical scattering: polarization analysis and measurement,” Opt. Eng. 34, 1599–1610 (1995). [CrossRef]  

20. C. M. Jarque and A. K. Bera, “A test for normality of observations and regression residuals,” Int. Stat. Rev 55, 163–172 (1987). [CrossRef]  

21. J. A. Hartigan and P. M. Hartigan, “The dip test of unimodality,” Ann. Stat. 13, 70–84 (1985). [CrossRef]  

22. D. Ververidis and C. Kotropoulos, “Gaussian mixture modeling by exploiting the Mahalanobis distance,” IEEE T. Signal Process 56, 2797–2811 (2008). [CrossRef]  

23. Y. Boykov and V. Kolmogorov, “An experimental comparison of min-cut/max-flow algorithms for energy minimization in vision,” IEEE T. Pattern Anal. 26, 1124–1137 (2004). [CrossRef]  

24. Y. Boykov, O. Veksler, and R. Zabih, “Fast approximate energy minimization via graph cuts,” IEEE T. Pattern Anal. 23, 1222–1239 (2001). [CrossRef]  

25. V. Kolmogorov and R. Zabih, “What energy functions can be minimized via graph cuts?,,” IEEE T. Pattern Anal. 26, 147–159 (2004). [CrossRef]  

References

  • View by:
  • |
  • |
  • |

  1. D. Miyasaki, 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).
    [Crossref]
  2. J. Chung, W. Jung, M. J. Hammer-Wilson, P. Wilder-Smith, and Z. Chen, “Use of polar decomposition for the diagnosis of oral precancer,” Appl. Opt. 46, 3038–3045 (2007).
    [Crossref] [PubMed]
  3. M. F. G. Wood, N. Ghosh, X. Guo, and I. A. Vitkin, “Towards noninvasive glucose sensing using polarization analysis of multiply scattered light,” Chap. 17, Handbook of Optical Sensing of Glucose in Biological Fluids and Tissues, V. V. Tuchin, ed., Series in Medical Physics and Biomedical Engineering, Vol. 12 (Taylor and Francis Publishing, 2008).
  4. M. F. G. Wood, N. Ghosh, M. A. Wallenburg, S. H. Li, R. D. Weisel, B. C. Wilson, R. K. Li, and I. A. Vitkin, “Polarization birefringence Measurements for characterizing the myocardium, including healthy, infracted, and stem cell treated regenerating cardiac tissues,” J. Biomed Opt. 15, 047009 (2010).
    [Crossref]
  5. P. Shukla and A. Pradhan, “Mueller decomposition images for cervical tissue: Potential for discriminating normal and dysplastic states,” Opt. Express 17, 1600–1609 (2009).
    [Crossref] [PubMed]
  6. 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, 1582–1593 (2011).
    [Crossref] [PubMed]
  7. J. S. Tyo, “Design of optimal polarimeters: maximization of signal-to-noise ratio and minimization of systematic errors,” Appl. Opt. 41619–630 (2002).
    [Crossref] [PubMed]
  8. S. N. Savenkov, “Optimization and structuring of the instrument matrix for polarimetric measurements,” Opt. Eng. 41, 965–972 (2002).
    [Crossref]
  9. J. Zallat, S. Anouz, and M. P. Stoll, “Optimal configurations for imaging polarimeters: impact of image noise and systematic errors,” J. Opt. A-Pure Appl. Op. 8, 807–814 (2006).
    [Crossref]
  10. S. Faisan, C. Heinrich, G. Sfikas, and J. Zallat, “Estimation of Mueller matrices using non-local means filtering,” Opt. Express 21, 4424–4438 (2013).
    [Crossref] [PubMed]
  11. J. Zallat and C. Heinrich, “Polarimetric data reduction: a Bayesian approach,” Opt. Express 15, 83–96 (2007).
    [Crossref] [PubMed]
  12. S. W. Hasinoff, “ Photon, Poisson Noise,” in Computer Vision, Katsushi Ikeuchi, ed. (SpringerUS, 2014)
    [Crossref]
  13. A. Cedilnik, K. Košmelj, and A. Blejec, “The distribution of the ratio of jointly normal variables,” Metodološki Zvezki – Advances in Methodology and Statistics 1, 99–108 (2004).
  14. J. Zallat, C. Heinrich, and M. Petremand, “A Bayesian approach for polarimetric data reduction: the Mueller imaging case,” Opt. Express 16, 7119–7133 (2008).
    [Crossref] [PubMed]
  15. G. McLachlan and T. Krischnan, The EM Algorithm and Extensions, 2nd edition (Wiley, 2008).
    [Crossref]
  16. J. Bilmes, “A gentle tutorial of the EM algorithm and its application to parameter estimation for gaussian mixture and hidden markov models,” Tech. rep., International Computer Science Institute (1998).
  17. G. Celeux and J. Diebolt, “The SEM algorithm: a probabilistic teacher algorithm derived from the EM algorithm for the mixture problem,” Comp. Statis. Quaterly 2, 73–82 (1985).
  18. W. H. Greene, Econometric Analysis, 7th edition (Prentice Hall, 2012), Chap. 9.
  19. S. Cloude and E. Pottier, “Concept of polarization entropy in optical scattering: polarization analysis and measurement,” Opt. Eng. 34, 1599–1610 (1995).
    [Crossref]
  20. C. M. Jarque and A. K. Bera, “A test for normality of observations and regression residuals,” Int. Stat. Rev 55, 163–172 (1987).
    [Crossref]
  21. J. A. Hartigan and P. M. Hartigan, “The dip test of unimodality,” Ann. Stat. 13, 70–84 (1985).
    [Crossref]
  22. D. Ververidis and C. Kotropoulos, “Gaussian mixture modeling by exploiting the Mahalanobis distance,” IEEE T. Signal Process 56, 2797–2811 (2008).
    [Crossref]
  23. Y. Boykov and V. Kolmogorov, “An experimental comparison of min-cut/max-flow algorithms for energy minimization in vision,” IEEE T. Pattern Anal. 26, 1124–1137 (2004).
    [Crossref]
  24. Y. Boykov, O. Veksler, and R. Zabih, “Fast approximate energy minimization via graph cuts,” IEEE T. Pattern Anal. 23, 1222–1239 (2001).
    [Crossref]
  25. V. Kolmogorov and R. Zabih, “What energy functions can be minimized via graph cuts?,,” IEEE T. Pattern Anal. 26, 147–159 (2004).
    [Crossref]

2013 (1)

2011 (1)

2010 (1)

M. F. G. Wood, N. Ghosh, M. A. Wallenburg, S. H. Li, R. D. Weisel, B. C. Wilson, R. K. Li, and I. A. Vitkin, “Polarization birefringence Measurements for characterizing the myocardium, including healthy, infracted, and stem cell treated regenerating cardiac tissues,” J. Biomed Opt. 15, 047009 (2010).
[Crossref]

2009 (1)

2008 (2)

J. Zallat, C. Heinrich, and M. Petremand, “A Bayesian approach for polarimetric data reduction: the Mueller imaging case,” Opt. Express 16, 7119–7133 (2008).
[Crossref] [PubMed]

D. Ververidis and C. Kotropoulos, “Gaussian mixture modeling by exploiting the Mahalanobis distance,” IEEE T. Signal Process 56, 2797–2811 (2008).
[Crossref]

2007 (2)

2006 (1)

J. Zallat, S. Anouz, and M. P. Stoll, “Optimal configurations for imaging polarimeters: impact of image noise and systematic errors,” J. Opt. A-Pure Appl. Op. 8, 807–814 (2006).
[Crossref]

2004 (3)

V. Kolmogorov and R. Zabih, “What energy functions can be minimized via graph cuts?,,” IEEE T. Pattern Anal. 26, 147–159 (2004).
[Crossref]

Y. Boykov and V. Kolmogorov, “An experimental comparison of min-cut/max-flow algorithms for energy minimization in vision,” IEEE T. Pattern Anal. 26, 1124–1137 (2004).
[Crossref]

A. Cedilnik, K. Košmelj, and A. Blejec, “The distribution of the ratio of jointly normal variables,” Metodološki Zvezki – Advances in Methodology and Statistics 1, 99–108 (2004).

2002 (3)

2001 (1)

Y. Boykov, O. Veksler, and R. Zabih, “Fast approximate energy minimization via graph cuts,” IEEE T. Pattern Anal. 23, 1222–1239 (2001).
[Crossref]

1995 (1)

S. Cloude and E. Pottier, “Concept of polarization entropy in optical scattering: polarization analysis and measurement,” Opt. Eng. 34, 1599–1610 (1995).
[Crossref]

1987 (1)

C. M. Jarque and A. K. Bera, “A test for normality of observations and regression residuals,” Int. Stat. Rev 55, 163–172 (1987).
[Crossref]

1985 (2)

J. A. Hartigan and P. M. Hartigan, “The dip test of unimodality,” Ann. Stat. 13, 70–84 (1985).
[Crossref]

G. Celeux and J. Diebolt, “The SEM algorithm: a probabilistic teacher algorithm derived from the EM algorithm for the mixture problem,” Comp. Statis. Quaterly 2, 73–82 (1985).

Anouz, S.

J. Zallat, S. Anouz, and M. P. Stoll, “Optimal configurations for imaging polarimeters: impact of image noise and systematic errors,” J. Opt. A-Pure Appl. Op. 8, 807–814 (2006).
[Crossref]

Antonelli, M. R.

Benali, A.

Bera, A. K.

C. M. Jarque and A. K. Bera, “A test for normality of observations and regression residuals,” Int. Stat. Rev 55, 163–172 (1987).
[Crossref]

Bilmes, J.

J. Bilmes, “A gentle tutorial of the EM algorithm and its application to parameter estimation for gaussian mixture and hidden markov models,” Tech. rep., International Computer Science Institute (1998).

Blejec, A.

A. Cedilnik, K. Košmelj, and A. Blejec, “The distribution of the ratio of jointly normal variables,” Metodološki Zvezki – Advances in Methodology and Statistics 1, 99–108 (2004).

Boykov, Y.

Y. Boykov and V. Kolmogorov, “An experimental comparison of min-cut/max-flow algorithms for energy minimization in vision,” IEEE T. Pattern Anal. 26, 1124–1137 (2004).
[Crossref]

Y. Boykov, O. Veksler, and R. Zabih, “Fast approximate energy minimization via graph cuts,” IEEE T. Pattern Anal. 23, 1222–1239 (2001).
[Crossref]

Cedilnik, A.

A. Cedilnik, K. Košmelj, and A. Blejec, “The distribution of the ratio of jointly normal variables,” Metodološki Zvezki – Advances in Methodology and Statistics 1, 99–108 (2004).

Celeux, G.

G. Celeux and J. Diebolt, “The SEM algorithm: a probabilistic teacher algorithm derived from the EM algorithm for the mixture problem,” Comp. Statis. Quaterly 2, 73–82 (1985).

Chen, Z.

Chung, J.

Cloude, S.

S. Cloude and E. Pottier, “Concept of polarization entropy in optical scattering: polarization analysis and measurement,” Opt. Eng. 34, 1599–1610 (1995).
[Crossref]

De Martino, A.

Diebolt, J.

G. Celeux and J. Diebolt, “The SEM algorithm: a probabilistic teacher algorithm derived from the EM algorithm for the mixture problem,” Comp. Statis. Quaterly 2, 73–82 (1985).

Faisan, S.

Gayet, B.

Ghosh, N.

M. F. G. Wood, N. Ghosh, M. A. Wallenburg, S. H. Li, R. D. Weisel, B. C. Wilson, R. K. Li, and I. A. Vitkin, “Polarization birefringence Measurements for characterizing the myocardium, including healthy, infracted, and stem cell treated regenerating cardiac tissues,” J. Biomed Opt. 15, 047009 (2010).
[Crossref]

M. F. G. Wood, N. Ghosh, X. Guo, and I. A. Vitkin, “Towards noninvasive glucose sensing using polarization analysis of multiply scattered light,” Chap. 17, Handbook of Optical Sensing of Glucose in Biological Fluids and Tissues, V. V. Tuchin, ed., Series in Medical Physics and Biomedical Engineering, Vol. 12 (Taylor and Francis Publishing, 2008).

Greene, W. H.

W. H. Greene, Econometric Analysis, 7th edition (Prentice Hall, 2012), Chap. 9.

Guo, X.

M. F. G. Wood, N. Ghosh, X. Guo, and I. A. Vitkin, “Towards noninvasive glucose sensing using polarization analysis of multiply scattered light,” Chap. 17, Handbook of Optical Sensing of Glucose in Biological Fluids and Tissues, V. V. Tuchin, ed., Series in Medical Physics and Biomedical Engineering, Vol. 12 (Taylor and Francis Publishing, 2008).

Hammer-Wilson, M. J.

Hartigan, J. A.

J. A. Hartigan and P. M. Hartigan, “The dip test of unimodality,” Ann. Stat. 13, 70–84 (1985).
[Crossref]

Hartigan, P. M.

J. A. Hartigan and P. M. Hartigan, “The dip test of unimodality,” Ann. Stat. 13, 70–84 (1985).
[Crossref]

Hasinoff, S. W.

S. W. Hasinoff, “ Photon, Poisson Noise,” in Computer Vision, Katsushi Ikeuchi, ed. (SpringerUS, 2014)
[Crossref]

Heinrich, C.

Ikeuchi, K.

Jarque, C. M.

C. M. Jarque and A. K. Bera, “A test for normality of observations and regression residuals,” Int. Stat. Rev 55, 163–172 (1987).
[Crossref]

Jung, W.

Kolmogorov, V.

Y. Boykov and V. Kolmogorov, “An experimental comparison of min-cut/max-flow algorithms for energy minimization in vision,” IEEE T. Pattern Anal. 26, 1124–1137 (2004).
[Crossref]

V. Kolmogorov and R. Zabih, “What energy functions can be minimized via graph cuts?,,” IEEE T. Pattern Anal. 26, 147–159 (2004).
[Crossref]

Košmelj, K.

A. Cedilnik, K. Košmelj, and A. Blejec, “The distribution of the ratio of jointly normal variables,” Metodološki Zvezki – Advances in Methodology and Statistics 1, 99–108 (2004).

Kotropoulos, C.

D. Ververidis and C. Kotropoulos, “Gaussian mixture modeling by exploiting the Mahalanobis distance,” IEEE T. Signal Process 56, 2797–2811 (2008).
[Crossref]

Krischnan, T.

G. McLachlan and T. Krischnan, The EM Algorithm and Extensions, 2nd edition (Wiley, 2008).
[Crossref]

Li, R. K.

M. F. G. Wood, N. Ghosh, M. A. Wallenburg, S. H. Li, R. D. Weisel, B. C. Wilson, R. K. Li, and I. A. Vitkin, “Polarization birefringence Measurements for characterizing the myocardium, including healthy, infracted, and stem cell treated regenerating cardiac tissues,” J. Biomed Opt. 15, 047009 (2010).
[Crossref]

Li, S. H.

M. F. G. Wood, N. Ghosh, M. A. Wallenburg, S. H. Li, R. D. Weisel, B. C. Wilson, R. K. Li, and I. A. Vitkin, “Polarization birefringence Measurements for characterizing the myocardium, including healthy, infracted, and stem cell treated regenerating cardiac tissues,” J. Biomed Opt. 15, 047009 (2010).
[Crossref]

McLachlan, G.

G. McLachlan and T. Krischnan, The EM Algorithm and Extensions, 2nd edition (Wiley, 2008).
[Crossref]

Miyasaki, D.

Novikova, T.

Petremand, M.

Pierangelo, A.

Pottier, E.

S. Cloude and E. Pottier, “Concept of polarization entropy in optical scattering: polarization analysis and measurement,” Opt. Eng. 34, 1599–1610 (1995).
[Crossref]

Pradhan, A.

Saito, M.

Sato, Y.

Savenkov, S. N.

S. N. Savenkov, “Optimization and structuring of the instrument matrix for polarimetric measurements,” Opt. Eng. 41, 965–972 (2002).
[Crossref]

Sfikas, G.

Shukla, P.

Stoll, M. P.

J. Zallat, S. Anouz, and M. P. Stoll, “Optimal configurations for imaging polarimeters: impact of image noise and systematic errors,” J. Opt. A-Pure Appl. Op. 8, 807–814 (2006).
[Crossref]

Tyo, J. S.

Validire, P.

Veksler, O.

Y. Boykov, O. Veksler, and R. Zabih, “Fast approximate energy minimization via graph cuts,” IEEE T. Pattern Anal. 23, 1222–1239 (2001).
[Crossref]

Ververidis, D.

D. Ververidis and C. Kotropoulos, “Gaussian mixture modeling by exploiting the Mahalanobis distance,” IEEE T. Signal Process 56, 2797–2811 (2008).
[Crossref]

Vitkin, I. A.

M. F. G. Wood, N. Ghosh, M. A. Wallenburg, S. H. Li, R. D. Weisel, B. C. Wilson, R. K. Li, and I. A. Vitkin, “Polarization birefringence Measurements for characterizing the myocardium, including healthy, infracted, and stem cell treated regenerating cardiac tissues,” J. Biomed Opt. 15, 047009 (2010).
[Crossref]

M. F. G. Wood, N. Ghosh, X. Guo, and I. A. Vitkin, “Towards noninvasive glucose sensing using polarization analysis of multiply scattered light,” Chap. 17, Handbook of Optical Sensing of Glucose in Biological Fluids and Tissues, V. V. Tuchin, ed., Series in Medical Physics and Biomedical Engineering, Vol. 12 (Taylor and Francis Publishing, 2008).

Wallenburg, M. A.

M. F. G. Wood, N. Ghosh, M. A. Wallenburg, S. H. Li, R. D. Weisel, B. C. Wilson, R. K. Li, and I. A. Vitkin, “Polarization birefringence Measurements for characterizing the myocardium, including healthy, infracted, and stem cell treated regenerating cardiac tissues,” J. Biomed Opt. 15, 047009 (2010).
[Crossref]

Weisel, R. D.

M. F. G. Wood, N. Ghosh, M. A. Wallenburg, S. H. Li, R. D. Weisel, B. C. Wilson, R. K. Li, and I. A. Vitkin, “Polarization birefringence Measurements for characterizing the myocardium, including healthy, infracted, and stem cell treated regenerating cardiac tissues,” J. Biomed Opt. 15, 047009 (2010).
[Crossref]

Wilder-Smith, P.

Wilson, B. C.

M. F. G. Wood, N. Ghosh, M. A. Wallenburg, S. H. Li, R. D. Weisel, B. C. Wilson, R. K. Li, and I. A. Vitkin, “Polarization birefringence Measurements for characterizing the myocardium, including healthy, infracted, and stem cell treated regenerating cardiac tissues,” J. Biomed Opt. 15, 047009 (2010).
[Crossref]

Wood, M. F. G.

M. F. G. Wood, N. Ghosh, M. A. Wallenburg, S. H. Li, R. D. Weisel, B. C. Wilson, R. K. Li, and I. A. Vitkin, “Polarization birefringence Measurements for characterizing the myocardium, including healthy, infracted, and stem cell treated regenerating cardiac tissues,” J. Biomed Opt. 15, 047009 (2010).
[Crossref]

M. F. G. Wood, N. Ghosh, X. Guo, and I. A. Vitkin, “Towards noninvasive glucose sensing using polarization analysis of multiply scattered light,” Chap. 17, Handbook of Optical Sensing of Glucose in Biological Fluids and Tissues, V. V. Tuchin, ed., Series in Medical Physics and Biomedical Engineering, Vol. 12 (Taylor and Francis Publishing, 2008).

Zabih, R.

V. Kolmogorov and R. Zabih, “What energy functions can be minimized via graph cuts?,,” IEEE T. Pattern Anal. 26, 147–159 (2004).
[Crossref]

Y. Boykov, O. Veksler, and R. Zabih, “Fast approximate energy minimization via graph cuts,” IEEE T. Pattern Anal. 23, 1222–1239 (2001).
[Crossref]

Zallat, J.

Ann. Stat. (1)

J. A. Hartigan and P. M. Hartigan, “The dip test of unimodality,” Ann. Stat. 13, 70–84 (1985).
[Crossref]

Appl. Opt. (2)

Comp. Statis. Quaterly (1)

G. Celeux and J. Diebolt, “The SEM algorithm: a probabilistic teacher algorithm derived from the EM algorithm for the mixture problem,” Comp. Statis. Quaterly 2, 73–82 (1985).

IEEE T. Pattern Anal. (3)

Y. Boykov and V. Kolmogorov, “An experimental comparison of min-cut/max-flow algorithms for energy minimization in vision,” IEEE T. Pattern Anal. 26, 1124–1137 (2004).
[Crossref]

Y. Boykov, O. Veksler, and R. Zabih, “Fast approximate energy minimization via graph cuts,” IEEE T. Pattern Anal. 23, 1222–1239 (2001).
[Crossref]

V. Kolmogorov and R. Zabih, “What energy functions can be minimized via graph cuts?,,” IEEE T. Pattern Anal. 26, 147–159 (2004).
[Crossref]

IEEE T. Signal Process (1)

D. Ververidis and C. Kotropoulos, “Gaussian mixture modeling by exploiting the Mahalanobis distance,” IEEE T. Signal Process 56, 2797–2811 (2008).
[Crossref]

Int. Stat. Rev (1)

C. M. Jarque and A. K. Bera, “A test for normality of observations and regression residuals,” Int. Stat. Rev 55, 163–172 (1987).
[Crossref]

J. Biomed Opt. (1)

M. F. G. Wood, N. Ghosh, M. A. Wallenburg, S. H. Li, R. D. Weisel, B. C. Wilson, R. K. Li, and I. A. Vitkin, “Polarization birefringence Measurements for characterizing the myocardium, including healthy, infracted, and stem cell treated regenerating cardiac tissues,” J. Biomed Opt. 15, 047009 (2010).
[Crossref]

J. Opt. A-Pure Appl. Op. (1)

J. Zallat, S. Anouz, and M. P. Stoll, “Optimal configurations for imaging polarimeters: impact of image noise and systematic errors,” J. Opt. A-Pure Appl. Op. 8, 807–814 (2006).
[Crossref]

J. Opt. Soc. Am. A (1)

Metodološki Zvezki – Advances in Methodology and Statistics (1)

A. Cedilnik, K. Košmelj, and A. Blejec, “The distribution of the ratio of jointly normal variables,” Metodološki Zvezki – Advances in Methodology and Statistics 1, 99–108 (2004).

Opt. Eng. (2)

S. N. Savenkov, “Optimization and structuring of the instrument matrix for polarimetric measurements,” Opt. Eng. 41, 965–972 (2002).
[Crossref]

S. Cloude and E. Pottier, “Concept of polarization entropy in optical scattering: polarization analysis and measurement,” Opt. Eng. 34, 1599–1610 (1995).
[Crossref]

Opt. Express (5)

Other (5)

M. F. G. Wood, N. Ghosh, X. Guo, and I. A. Vitkin, “Towards noninvasive glucose sensing using polarization analysis of multiply scattered light,” Chap. 17, Handbook of Optical Sensing of Glucose in Biological Fluids and Tissues, V. V. Tuchin, ed., Series in Medical Physics and Biomedical Engineering, Vol. 12 (Taylor and Francis Publishing, 2008).

S. W. Hasinoff, “ Photon, Poisson Noise,” in Computer Vision, Katsushi Ikeuchi, ed. (SpringerUS, 2014)
[Crossref]

G. McLachlan and T. Krischnan, The EM Algorithm and Extensions, 2nd edition (Wiley, 2008).
[Crossref]

J. Bilmes, “A gentle tutorial of the EM algorithm and its application to parameter estimation for gaussian mixture and hidden markov models,” Tech. rep., International Computer Science Institute (1998).

W. H. Greene, Econometric Analysis, 7th edition (Prentice Hall, 2012), Chap. 9.

Cited By

OSA participates in Crossref's Cited-By Linking service. Citing articles from OSA journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (5)

Fig. 1
Fig. 1 Steps related to the estimation of the parameters Θ and of the number K of classes.
Fig. 2
Fig. 2 (a) Simulated segmentation map, (b) simulated lighting map, and (c) example of a mean transmittance image obtained for a 10 dB SNR. Representative segmentation maps obtained (SNR=10 dB) with the proposed approach (d), and with the approach of [14] (without normalization (e), by using the pseudo-inverse for normalization (thresholding (f) and polynomial modeling (g)), and by using the ground truth for normalization (h).
Fig. 3
Fig. 3 Exemple of synthetic intensity images obtained with a SNR of 10 dB.
Fig. 4
Fig. 4 Sketch of the scene composition (a), automatic segmentation with the proposed approach using the Jarque–Bera test (7 classes) (b), and using the dip test (c).
Fig. 5
Fig. 5 The 16 intensity images of the object acquired by a rotating quarter–wave plate Mueller imaging polarimeter in transmission configuration.

Tables (4)

Tables Icon

Table 1 Mueller matrices used for the creation of the synthetic data.

Tables Icon

Table 2 Number of occurrences of the estimated number K of classes for 100 experiments, for different values of α and of the SNR, with the Jarque-Bera test.

Tables Icon

Table 3 Mean values of the relative estimation error RE (in percent) for the four Mueller matrices and the corresponding covariance matrices. Averaging has been done over 100 realizations. Numbers between parentheses represent standard-deviation of related quantities.

Tables Icon

Table 4 Mueller matrices associated to dichroic patches and to their related linear polarizer.

Equations (45)

Equations on this page are rendered with MathJax. Learn more.

I ( s ) = P M ( s ) _ ( + noise ) ,
I ( s ) = κ ( s ) P M k _ + ε ( s ) = κ ( s ) Ψ k + ε ( s ) ,
P ( I ( s ) , y ( s ) = k , κ ( s ) | Θ ) = P ( I ( s ) | y ( s ) = k , κ ( s ) , Θ ) P ( y ( s ) = k | Θ ) P ( κ ( s ) | Θ ) .
I ( s ) | y ( s ) = k , κ ( s ) , Θ ~ 𝒩 ( κ ( s ) Ψ k , Σ k ) .
P ( I ( s ) | Θ ) = k = 1 K κ ( s ) P ( I ( s ) , y ( s ) = k , κ ( s ) | Θ ) d κ ( s ) ,
P ( I ( s ) | Θ ) = k = 1 K κ ( s ) π k P ( κ ( s ) | Θ ) P ( I ( s ) | y ( s ) = k , κ ( s ) , Θ ) d κ ( s ) .
L ( Θ ) = P ( I | Θ ) = s P ( I ( s ) | Θ ) .
P ( I , y , κ | Θ ) = P ( I , y , κ | Θ ) P ( y , κ | Θ ) = [ s P ( I ( s ) | y ( s ) , κ ( s ) | Θ ) ] P ( y | Θ ) P ( κ | Θ ) = [ s P ( I ( s ) | y ( s ) , κ ( s ) , Θ ) ] [ s P ( y ( s ) | Θ ) ] [ s P ( κ ( s ) | Θ ) ] = s P ( I ( s ) , y ( s ) , κ ( s ) | Θ ) ,
Q ( Θ , Θ ) = s k = 1 K w k ( s ) κ ( s ) = 0 + log [ P ( I ( s ) , y ( s ) = k , κ ( s ) | Θ k ) ] × P ( κ ( s ) | y ( s ) = k , I ( s ) , Θ ) d κ ( s ) ,
κ k ( s ) | y ( s ) = k , I ( s ) , Θ ~ 𝒩 + ( α k β k , α k ) ,
Θ ^ = arg max Θ s k = 1 K w k ( s ) log [ P ( I ( s ) , y ( s ) = k , k κ ( s ) | Θ k ) ] = arg max Θ s k = 1 K w k ( s ) log [ π k P ( I ( s ) | Θ k , y ( s ) = k , k κ ( s ) ) ]
π ^ k = 1 # pixels s w k ( s ) ,
arg min { Ψ , Σ } s w k ( s ) [ log ( | Σ | ) + ( I ( s ) κ k ( s ) Ψ ) T Σ 1 ( I ( s ) κ k ( s ) Ψ ) ] .
I ¯ k = 1 s w k ( s ) κ k 2 ( s ) s w k ( s ) κ k ( s ) I ( s ) .
M ^ k = arg min M Σ k 1 / 2 ( I ¯ k P M ) 2 ,
Σ ^ k = ( 1 s w k ( s ) s w k ( s ) ( I ( s ) κ k ( s ) Ψ ^ k ) ( I ( s ) κ k ( s ) Ψ ^ k ) T ) I N × N ,
r ^ k ( s ) = I ( s ) κ ^ k ( s ) Ψ ^ k N
p min = min { p k , n } k = 1 K , n = 1 N .
p min 1 ( 1 α ) 1 N K .
P ( I 0 | β ) = 1 Z ( β ) exp ( β U ( I 0 ) ) ,
U ( I 0 ) = [ t , s ] 𝕀 ( I 0 ( t ) = I 0 ( s ) )
{ I ^ 0 , β ^ } = arg max { I 0 , β } log P ( I 0 , β | I , Θ ) = arg max { I 0 , β } log P ( I 0 | β ) + s log P ( I ( s ) | I 0 ( s ) , Θ ) ,
{ P ( I ( s ) , κ ( s ) | I 0 ( s ) , Θ ) = P ( I ( s ) | I 0 ( s ) , κ ( s ) , Θ ) P ( κ ( s ) ) , I ( s ) | y ( s ) = k , κ ( s ) , Θ ~ 𝒩 ( κ ( s ) Ψ k , Σ k ) , P ( I ( s ) | I 0 ( s ) , Θ ) = P ( I ( s ) , κ ( s ) | I 0 ( s ) , Θ ) d κ ( s ) ,
RE Mk = M ^ k M k M k , RE Σ k = Σ ^ k Σ k Σ k ,
I norm ( s ) = κ ( s ) m 00 ( s ) P M k _ + ε ( s ) m 00 ( s ) , P M k _ + ε ( s ) m 00 ( s ) .
I norm ( s ) = P M k _ + ε ( s ) κ ( s ) ,
Ψ 1 = s I ( s ) s κ ^ ( s ) and M 1 = ( P t P ) 1 P t Ψ 1 .
π K + 1 = π k ( p ) = s C K + 1 w k ( p ) ( s ) s C K + 1 C k w k ( p ) ( s ) , and π k = π k ( p ) s C k w k ( p ) ( s ) s C K + 1 C k w k ( p ) ( s )
log P ( I , y , κ | Θ ) = s log [ P ( I ( s ) , y ( s ) , κ ( s ) | Θ ) ] .
Q ( Θ , Θ ) = y κ log [ P ( I , y , κ | Θ ) ] P ( y , κ | I , Θ ) d κ .
Q ( Θ , Θ ) = s k = 1 K κ ( s ) = 0 + log [ P ( I ( s ) , y ( s ) = k , κ ( s ) | Θ k ) ] P ( y ( s ) = k , κ ( s ) | I ( s ) , Θ ) d κ ( s ) .
P ( y ( s ) = k , κ ( s ) | I ( s ) , Θ ) = P ( κ ( s ) | y ( s ) = k , I ( s ) , Θ ) P ( y ( s ) = k | I ( s ) , Θ ) ,
w k ( s ) = P ( y ( s ) = k | I ( s ) , Θ ) P ( y ( s ) = k , I ( s ) | Θ ) P ( I ( s ) | y ( s ) = k , Θ ) P ( y ( s ) = k | Θ ) [ P ( I ( s ) , κ ( s ) | y ( s ) = k , Θ ) d κ ( s ) ] π k π k κ ( s ) = 0 + P ( I ( s ) | y ( s ) = k , κ ( s ) , Θ k ) d κ ( s ) π k α k | Σ k | exp ( α k β k 2 γ k 2 ) ( 1 + erf ( β k α k 2 ) ) ,
{ α k = ( Ψ k T Σ k 1 Ψ k ) 1 β k = Ψ k T Σ k 1 I ( s ) γ k = I ( s ) Σ k 1 I ( s ) .
Q ( Θ , Θ ) = s k = 1 K w k ( s ) κ ( s ) = 0 + log [ P ( I ( s ) , y ( s ) = k , κ ( s ) | Θ k ) ] × P ( κ ( s ) | y ( s ) = k , I ( s ) , Θ ) d κ ( s ) .
I ¯ k = 1 s w k ( s ) κ k 2 ( s ) s w k ( s ) κ k ( s ) I ( s ) .
{ M ^ k , Σ ^ k } = arg min Θ = { M , Σ } ( s w k ( s ) ) log ( | Σ | ) + ( s w k ( s ) κ k 2 ( s ) ) Σ 1 / 2 ( I ¯ k P M ) 2 .
M ^ k = arg min M Σ k 1 / 2 ( I ¯ k PM ) 2 ,
Σ ^ k = ( 1 s w k ( s ) s w k ( s ) ( I ( s ) κ k ( s ) Ψ ^ k ) ( I ( s ) κ k ( s ) Ψ ^ k ) T ) I N × N ,
c k ( Ψ , Σ ) = s w k ( s ) ( I ( s ) κ k ( s ) Ψ ) T Σ 1 ( I ( s ) κ k ( s ) Ψ ) .
c k ( Ψ , Σ ) = s w k ( s ) κ k 2 ( s ) ( I ( 1 ) ( s ) Ψ ( 1 ) ) T ( I ( 1 ) ( s ) Ψ ( 1 ) ) ) ,
D ( s ) = w k ( s ) κ k 2 ( s ) s w k ( s ) κ k 2 ( s ) I N × N ,
c k ( Ψ , Σ ) = ( s w k ( s ) κ k 2 ( s ) ) s ( I ( 1 ) ( s ) Ψ ( 1 ) ) T D ( s ) ( I ( 1 ) ( s ) Ψ ( 1 ) ) .
s ( I ( 1 ) ( s ) P ( 1 ) M ) T D ( s ) ( I ( 1 ) ( s ) P ( 1 ) M ) = I ¯ P ( 1 ) M 2 ( + constant term ) ,
c k ( Ψ , Σ ) = ( s w k ( s ) κ k 2 ( s ) ) Σ 1 / 2 ( I ¯ k P M ) 2 ( + constant term )

Metrics