# Masses of Negative Multinomial Distributions: Application to Polarimetric Image Processing

Authors: Philippe Bernardoff; Florent Chatelain; Jean-Yves Tourneret
Journal: Journal of Probability and Statistics
 (2013)

---

## Abstract

This paper derives new closed-form expressions for the masses of negative multinomial distributions. These masses can be maximized to determine the maximum likelihood estimator of its unknown parameters. An application to polarimetric image processing is investigated. We study the maximum likelihood estimators of the polarization degree of polarimetric images using different combinations of images.


---

## Body

## 1. Introduction

The univariate negative binomial distribution is uniquely defined in many statistical textbooks. However, extensions defining multivariate negative multinomial distributions (NMDs) are more controversial. Most definitions are based on the probability generating function (PGF) of these distributions. Doss [1] proposed to define the PGF of an NMD as the inverse λth power of a polynomial linear in each of its variables. This definition can also be found in the famous textbook [2, page 93] and the computation of its modes has been investigated in [3]. A more general class of NMDs introduced in [4] was characterized by PGFs of the form |I-Q|λ|I-QZ|-λ, where λ>0, Q is an n×n matrix, and Z=diag(z1,…,zn). In particular, matrices Q yielding infinitely divisible PGFs were derived. Finally, Bar-Lev et al. [5] introduced NMDs whose PGFs are defined as the inverse λth power of any affine polynomial. Necessary and sufficient conditions on the coefficients of this affine polynomial were derived to obtain the PGF of a multivariate distribution defined on ℕ0n (where ℕ0 is the set of nonnegative integers) [6]. These very general multivariate NMDs were recently used for image processing applications in [7].The family of NMDs introduced in [5] can be defined as follows. Let us denote [n]={1,…,n} the set of the n first nonzero integers. We denote  zT=∏t∈T‍zt as the monomial obtained by multiplying all the entries of the vector z=(z1,…,zn)∈ℝn whose indexes belong to T, where T⊂[n] stands for any subset of the indexes. Let Pn(z)=∑T⊂[n],T≠∅‍pTzT be an affine polynomial with respect to the n variables (z1,…,zn) such that 1-Pn(1)≠0. The NMD distribution defined at pair (n,Pn) is represented by its PFG which is given by
(1)𝒢𝒩ℳ(n,Pn)(z)=(1-Pn(z))-λ(1-Pn(1))λ.
Such laws are denoted as 𝒩ℳ(n,Pn). However, as explained in [6], all couples (n,Pn) do not provide a valid NMD. More specifically, Bernardoff has derived a finite number of conditions over Pn such that (1-Pn(z))-λ(1-Pn(1))λ is the PGF of an NMD for all positive integer n. The corresponding expression of the coefficient of zα in the Taylor expansion of (1-Pn(z))-λ is given by the formula
(2)cα(λ,Pn)=∑k∈Kα‍(λ)|k|pkk!,
where Kα={k:𝒫n→ℕ} and 𝒫n is the set of all subsets of [n]. However, this expression of cα(λ,Pn) does not allow us to explicitly compute the masses of NMDs in the general case.As a first goal of this paper, we propose a way of computing the masses of multivariate NMDs𝒩ℳ(n,Pn) defined above. A specific attention is devoted to bivariate and trivariate cases. In particular, it allows us to retrieve the results of [7] obtained for bivariate NMDs. The second part of the paper is devoted to the application of NMDs to image processing, more specifically to polarimetric image processing [8, 9]. Polarimetric image processing has received a considerable attention in the image processing and optical communities (see for instance [10–12] and references therein). The state of polarization of a polarimetric image is classically characterized by the degree of polarization (DoP) whose estimation is of major importance [13, 14]. The DoP of polarimetric images can be classically estimated by using four images associated to four different polarizations [15]. However, estimating the DoP using less than four images is interesting since it allows one to reduce the acquisition time and the resulting cost of the imaging system. As a consequence, there has been recently an increasing interest in deriving estimators of the DoP based on a reduced number of polarimetric images. Depending on the intensity of the acquired images, polarimetric images are referred to as low flux or high flux images (low flux corresponding to a small intensity and high flux to a larger intensity). DoP estimation based on a single polarimetric image was considered in [16, 17] under high flux and low flux assumptions. DoP estimators derived from two intensity images degraded by fully developed speckle noise were studied in [18, 19]. Finally, imaging systems using three polarimetric images were studied in [20, 21], under high flux and low flux assumptions.This paper studies the maximum likelihood estimators (MLEs) of the square DoP based on two or three polarimetric images. These estimators are computed by maximizing the masses of bivariate or trivariate NMDs derived in the first part of this work.The paper is organized as follows. Section2 recalls important results on NMDs. Section 3 proposes a new way of computing masses of NMDs. A particular attention is devoted to bivariate and trivariate cases. Section 4 addresses the problem of estimating the square DoP of low flux polarimetric images using the maximum likelihood (ML) principle. Different MLEs are constructed depending on the number of available polarimetric images. Simulation results are presented in Section 5.


## 2. Negative Multinomial Distributions

Ann-variate NMD is the distribution of a random vector N=(N1,…,Nn) taking its values in ℕ0n whose PGF is
(3)GN(z)=𝔼(∏k=1nzkNk)=[Pn(z)]-λ,
where 𝔼 denotes the mathematical expectation, z=(z1,…,zn),λ>0, and Pn(z) is an affine polynomial of order n. (A polynomial Pn(z) with respect to z=(z1,…,zn) is affine if the one variable polynomial zj↦Pn(z) can be written as A(-j)zj+B(-j) (for any j=1,…,d), where A(-j) and B(-j) are polynomials with respect to the zi's with i≠j.) These discrete distributions have received much interest in the literature (see for instance [2] and the references therein). Of course, the affine polynomial Pn has to satisfy appropriate conditions to ensure that GN(z) is a PGF. These conditions include the trivial equality Pn(1,…,1)=1. However, determining all pairs (Pn,λ) such that GN(z) is a PGF is still an open problem (see [6], for discussions related to this problem). As explained in [6], the affine polynomial Pn(z) can be rewritten
(4)Pn(z)=An(a1z1,…,anzn)An(a1,…,an),
where a1,…,an are positive numbers and An is an affine polynomial such that An(0,…,0)=1. The Taylor expansions of [An(z)]-λ and [Pn(z)]-λ in the neighborhood of (0,…,0) will be denoted as follows:
(5)[An(z)]-λ=∑α∈ℕ0ncα(λ,An)zα,[Pn(z)]-λ=∑α∈ℕ0ncα(λ,Pn)zα,
where α=(α1,…,αn) and zα=∏i=1n‍ziαi. Equations (3) and (4) clearly show that the masses of multivariate NMDs denoted as cα(λ,Pn) can be expressed as follows:
(6)cα(λ,Pn)=cα(λ,An)∏i=1naiαiAn(a1,…,an)-λ.


## 3. Masses of Negative Multinomial Distributions

In this section, we derive new expressions for the coefficientscα(λ,An) that will be used to compute the masses of NMDs. The particular cases of bivariate and trivariate NMDs will play an important role for the estimation of the DoP on polarimetric images. In order to compute the cα(λ,An), we derive several results summarized in this section whereas all demonstrations are reported in the appendix.Theorem 1.
Denote𝔓n* as the set of nonempty subsets of [n]={1,…,n}. Any affine polynomial An such that An(0)=1 denoted as
(7)An(z)=1-∑T∈𝔓n*aTzT
can be expressed as follows:
(8)An(z)=∏i∈[n](1-aizi)-∑T∈𝔓n*,|T|⩾2dTnzT∏i∈[n]∖T(1-aizi),
where |T| is the cardinal of the set T. Moreover
(9)An(z)=[∏i∈[n]‍(1-aizi)]×(1-∑T∈𝔓n*,|T|⩾2‍dTnzT∏i∈T(1-aizi))(10)=[∏i∈[n]‍(1-aizi)]×(1-Qn(z11-a1z1,…,zn1-anzn)),
where Qn is the polynomial defined by Qn(z)=∑T∈𝔓n*,|T|⩾2‍dTnzT and dTn is related to the 2|T|-1 variables aS,S∈𝔓T* as follows:
(11)dTn=∑T∈𝔓n|T|>1|T|‍aTa[n]∖T+(|T|-1)∏i∈T‍ai.Remark 2.
In the trivariate case defined byn=3, the polynomial A3(z) can be expressed as
(12)A3(z)=(1-a1z1)(1-a2z2)(1-a3z3)×[1-Q3(z11-z1,z21-z2,z31-z3)],
where the coefficients of the polynomial
(13)Q3(z)=b1,2z1z2+b1,3z1z3+b2,3z2z3+b1,2,3z1z2z3
can be determined using the relations
(14)bi,j=ai,j+aiaj,i,j∈{1,2,3},i≠j,b1,2,3=a1,2,3+a1a2,3+a2a1,3+a3a1,2+2a1a2a3.
The next theorem provides a relation between the coefficients of the polynomials An(z) and Qn(z) introduced above.Theorem 3.
LetAn(z)=1-∑T∈𝔓n*‍aTzT,  a=(a1,…,an), and Qn  be the affine polynomial defined in (10) and (11). For any α and γ in ℕn, denote as cγ(λ,An) the coefficient of zγ in the Taylor expansion of [An(z)]-λ and as cα(λ,1-Qn) the coefficient of zα in the Taylor expansion [1-Qn(z)]-λ. The following relation can be obtained:(15)cγ(λ,An)=∑α+β=γ‍cα(λ,1-Qn)(λ1+α)βaββ!(16)=∑0⩽βi⩽γi,i=1,…,n‍cγ-β(λ,1-Qn)×∏i=1n‍(λ+γi-βi)βiaiβiβi!(17)=∑0⩽αi⩽γi,i=1,…,ncα(λ,1-Qn)×∏i=1n(λ+αi)βiaiγi-αi(γi-αi)!.The masses of NMDs can be directly obtained from this theorem. The particular cases of bivariate and trivariate NMDs are considered in the following subsections since the corresponding masses will be useful in the application considered in the second part of this paper.


### 3.1. Bivariate NMDs

Theorem 4.
Consider the affine polynomial of order 2 with variablesz=(z1,z2) defined by
(18)A2(z)=1-∑T∈𝔓2*‍aTzT=1-a1z1-a2z2-a1,2z1z2.
The coefficient of zγ in the Taylor expansion of [A2(z)]-λ, can be computed as follows
(19)cγ(λ,P2)=(λ)max(γ1,γ2)∑ℓ=0min(γ1,γ2)‍(λ+ℓ)min(γ1,γ2)-ℓ(γ1-ℓ)!(γ2-ℓ)!ℓ!×a1γ1-ℓa2γ2-ℓb1,2ℓ.Remark 2.
The result(A.11) was mentioned in [7] without the factorization leading to (19). If a1≠0 and a2≠0, an equivalent formulation of (19) is
(20)cγ(λ,P2)=a1γ1a2γ2γ1!γ2!(λ)max(γ1,γ2)×∑ℓ=0min(γ1,γ2)‍(λ+ℓ)min(γ1,γ2)-ℓ×(γ1ℓ)(γ2ℓ)ℓ!(b1,2a1a2)ℓ.


### 3.2. Trivariate NMDs

Theorem 6.
Consider the affine polynomial with the three variablesz=(z1,z2,z3) defined by
(21)A3(z)=1-∑T∈𝔓3*aTzT.
The coefficient of zγ in the Taylor expansion of [A3(z)]-λ are
(22)cγ(λ,P3)=∑β1=0γ1‍∑β2=0γ2‍∑β3=0γ3‍∑v=∥γ-β∥⌊|γ-β|/2⌋(λ)v‍×b2,3v-γ1+β1b1,3v-γ2+β2b1,2v-γ3+β3∏i=13(v-γi+βi)!b1,2,3|γ-β|-2v(|γ-β|-2v)!×(λ+γ1-β1)β1β1!(λ+γ2-β2)β2β2!×(λ+γ3-β3)β3β3!a1β1a2β2a3β3.
When ai≠0, bi,j≠0, i≠j, i=1,2,3, j=1,2,3, and b1,2,3≠0, an equivalent expression is
(23)cγ(λ,P3)=a1γ1a2γ2a3γ3∑β1=0γ1‍∑β2=0γ2‍∑β3=0γ3‍∑v=∥γ-β∥⌊|γ-β|/2⌋(λ)v‍×(b2,3b1,3b1,2/b1,2,32)v(|γ-β|-2v)!∏i=13(v-γi+βi)!×[∏i=13‍(λ+γi-βi)βiβi!](b1,2,3a1b2,3)γ1-β1×(b1,2,3a2b1,3)γ2-β2(b1,2,3a3b1,2)γ3-β3(24)=(λ)∥γ∥aγγ!∑α1=0γ1‍∑α2=0γ2‍∑α3=0γ3‍∑v=∥α∥⌊|α|/2⌋(λ+∥α∥)v-∥α∥‍(|α|-2v)!∏i=13(v-αi)!×(b2,3b1,3b1,2b1,2,32)vα!(b1,2,3a1b2,3)α1(b1,2,3a2b1,3)α2(b1,2,3a3b1,2)α3×∏i=13‍[(λ+αi)γi-αi(γiαi)].


## 3.1. Bivariate NMDs

Theorem 4.
Consider the affine polynomial of order 2 with variablesz=(z1,z2) defined by
(18)A2(z)=1-∑T∈𝔓2*‍aTzT=1-a1z1-a2z2-a1,2z1z2.
The coefficient of zγ in the Taylor expansion of [A2(z)]-λ, can be computed as follows
(19)cγ(λ,P2)=(λ)max(γ1,γ2)∑ℓ=0min(γ1,γ2)‍(λ+ℓ)min(γ1,γ2)-ℓ(γ1-ℓ)!(γ2-ℓ)!ℓ!×a1γ1-ℓa2γ2-ℓb1,2ℓ.Remark 2.
The result(A.11) was mentioned in [7] without the factorization leading to (19). If a1≠0 and a2≠0, an equivalent formulation of (19) is
(20)cγ(λ,P2)=a1γ1a2γ2γ1!γ2!(λ)max(γ1,γ2)×∑ℓ=0min(γ1,γ2)‍(λ+ℓ)min(γ1,γ2)-ℓ×(γ1ℓ)(γ2ℓ)ℓ!(b1,2a1a2)ℓ.


## 3.2. Trivariate NMDs

Theorem 6.
Consider the affine polynomial with the three variablesz=(z1,z2,z3) defined by
(21)A3(z)=1-∑T∈𝔓3*aTzT.
The coefficient of zγ in the Taylor expansion of [A3(z)]-λ are
(22)cγ(λ,P3)=∑β1=0γ1‍∑β2=0γ2‍∑β3=0γ3‍∑v=∥γ-β∥⌊|γ-β|/2⌋(λ)v‍×b2,3v-γ1+β1b1,3v-γ2+β2b1,2v-γ3+β3∏i=13(v-γi+βi)!b1,2,3|γ-β|-2v(|γ-β|-2v)!×(λ+γ1-β1)β1β1!(λ+γ2-β2)β2β2!×(λ+γ3-β3)β3β3!a1β1a2β2a3β3.
When ai≠0, bi,j≠0, i≠j, i=1,2,3, j=1,2,3, and b1,2,3≠0, an equivalent expression is
(23)cγ(λ,P3)=a1γ1a2γ2a3γ3∑β1=0γ1‍∑β2=0γ2‍∑β3=0γ3‍∑v=∥γ-β∥⌊|γ-β|/2⌋(λ)v‍×(b2,3b1,3b1,2/b1,2,32)v(|γ-β|-2v)!∏i=13(v-γi+βi)!×[∏i=13‍(λ+γi-βi)βiβi!](b1,2,3a1b2,3)γ1-β1×(b1,2,3a2b1,3)γ2-β2(b1,2,3a3b1,2)γ3-β3(24)=(λ)∥γ∥aγγ!∑α1=0γ1‍∑α2=0γ2‍∑α3=0γ3‍∑v=∥α∥⌊|α|/2⌋(λ+∥α∥)v-∥α∥‍(|α|-2v)!∏i=13(v-αi)!×(b2,3b1,3b1,2b1,2,32)vα!(b1,2,3a1b2,3)α1(b1,2,3a2b1,3)α2(b1,2,3a3b1,2)α3×∏i=13‍[(λ+αi)γi-αi(γiαi)].


## 4. Estimating the Polarization Degree of Low Flux Polarimetric Images Using Maximum Likelihood Methods

### 4.1. Low Flux Polarimetric Images

The state of the polarization of the light can be described by the random behavior of a complex vectorA=(AX,AY), called the Jones vector, whose covariance matrix, called the polarization matrix, is
(25)Γ=(𝔼[AXAX*]𝔼[AXAY*]𝔼[AYAX*]𝔼[AYAY*])≜(a1a3+ia4a3-ia4a2),
where* denotes the complex conjugate. The covariance matrix Γ is a nonnegative Hermitian matrix whose diagonal terms are the intensity components in the X and Y directions. The cross terms of Γ are the correlations between the Jones components. If we assume a fully developed speckle, the Jones vector A is distributed according to a complex Gaussian distribution whose probability density function (pdf) is [15]:
(26)p(A)=1π2|Γ|exp(-A†Γ-1A),
where |Γ| is the determinant of the matrix Γ and † denotes the conjugate transpose operator. As a consequence, the statistical properties of A are fully characterized by the covariance matrix Γ. The different components of Γ can be classically estimated by using four intensity images that are related to the components of the Jones vector as follows (see [20], for more details):
(27)I1=|AX|2,I2=|AY|2,I3=12|AX|2+12|AY|2+Re(AXAY*),I4=12|AX|2+12|AY|2+Im(AXAY*).
The state of the polarization of the light is classically characterized by the square DoP defined by [15, pages 134–136]
(28)P2=1-4|Γ|[trace(Γ)]2=1-4[a1a2-(a32+a42)](a1+a2)2,
where trace(Γ) is the trace of the matrix Γ. The light is totally depolarized for P=0, totally polarized for P=1, and partially polarized when P∈]0,1[. As a consequence, estimating the square DoP of a polarimetric image is important in many practical applications. Different estimation methods of P2 using several combinations of intensity images were studied in [20]. Since only one realization of the random vector I=(I1,…,I4)T was available for a given pixel of a polarimetric image, the image was supposed to be locally stationary and ergodic. These assumptions were used to derive square DoP estimators using several neighbor pixels belonging to a so-called estimation window.This paper considers practical applications where the intensity level of the reflected light is very low (low flux assumption), which leads to an additional source of fluctuations on the detected signal. Under the low flux assumption, the quantum nature of the light leads to a Poisson-distributed noise which can become very important relatively to the mean value of the signal at a low photon level. As a consequence, the observed pixels of the low flux polarimetric image are discrete random variables contained in the vectorN=(N1,…,N4) such that the conditional distributions of the random variables Nl∣Il, for l=1,…,4 are independent and distributed according to Poisson distributions with means Il, for l=1,…,4. The resulting joint distribution of N is a multivariate mixed Poisson distribution [22]:
(29)P(N=k)=∫‍⋯(ℝ+)4∫‍∏l=14‍Ilklkl!exp(-Il)f(I)dI,
where k=(k1,…,k4), ki∈ℕ, and f(I) is the joint pdf of the intensity vector. This section studies estimators of the square DoP P2 defined in (28) based on several vectors N1,…,Nn belonging to the estimation window. These estimators are constructed from estimates of the covariance matrix elements ai, i=1,…,4. As explained in the introduction, several studies have been recently devoted to the estimation of the square DoP using less than four polarimetric images. This paper goes into this direction by deriving estimators based on the observation of 2 or 3 polarimetric images.The joint distribution of the intensity vectorI is known to be a multivariate gamma distribution whose Laplace transform is [20]
(30)E[exp(∑k=14‍zkIk)]=1P4(z),
where the affine polynomial P4 is
(31)P4(z)=1+zμ+ka[2z1z2+z3z4+(z1+z2)(z3+z4)]
with z=(z1,…,z4) and
(32)ka=12(a1a2-a32-a42),μ=(a1,a2,a3+a1+a22,a4+a1+a22)T.
As a consequence, the distribution of N is an NMD whose PGF can be written as (the interested reader is invited to consult [22] for more details)(33)GN(z)=1P4(z1-1,z2-1,z3-1,z4-1).
The results of Section 2 allow us to compute the masses of N that will be useful for studying the maximum likelihood estimator (MLE) of the square DoP.


### 4.2. MLE Using Three Polarimetric Images

The PGF ofN~=(N1,N2,N3) can be computed from (33) by setting z4=1. The following result can be obtained:
(34)GN~(z)=1P3(z)
with
(35)P3(z)=P3(0)+z1(μ1-3ka)+z2(μ2-3ka)+z3(μ3-2ka)+ka(2z1z2+z1z3+z2z3),z=(z1,z2,z3), and P3(0)=1-∑i=13‍μi+4ka. The results of Section 3.2 can then be used to express the masses of N~ as a function of θ=(a1,a2,a3,a42)T.The ML estimator ofθ based on several vectors N~k belonging to the estimation window (where k=1,…,K and K is the number of pixels of the observation window) is obtained by maximizing the log-likelihood
(36)l3(N~1,…,N~K∣θ)=∑k=1K‍log[P(N~k)]
with respect to θ (note again that P(N~k) is the mass of a trivariate NMD that has been computed in Section 3.2). The practical determination of the ML estimator of θ is achieved by using a Newton-Raphson procedure. The ML estimators of the vector θ, denoted as θ~=(a~1,a~2,a~3,a42~)T, are then plugged into (28) to provide the ML estimator of the square DoP based on three polarimetric images: (37)P~2=1-4[a~1a~2-(a~32+a42~)](a~1+a~2)2.


### 4.3. MLE Using Two Polarimetric Images

The PGF ofN_=(N1,N2) can be computed from (34) by setting z3=1. The following result can be obtained:
(38)GN_(z)=1P2(z)
with
(39)P2(z)=P2(0)+z1(μ1-2ka)+z2(μ2-2ka)+2kaz1z2,z=(z1,z2) and P2(0)=1-∑i=12‍μi+2ka. The results of Section 3.1 can then be used to express the masses of N_ as a function of θ=(a1,a2,ka)T.The MLE ofθ based on several vectors N_k belonging to the estimation window is obtained by maximizing the log-likelihood
(40)l2(N_1,…,N_K∣θ)=∑k=1K‍log[P(N_k)]
with respect to θ (note that P(N_k) is the mass of a bivariate NMD that has been computed in Section 3.1). The practical determination of the ML estimator of θ is achieved by using a Newton-Raphson procedure. The ML estimators of the vector θ elements, denoted as θ_=(a_1,a_2,ka_)T, are then plugged into (28) to provide the MLE of the square DoP based on two polarimetric images:
(41)P_2=1-8ka_(a_1+a_2)2.


## 4.1. Low Flux Polarimetric Images

The state of the polarization of the light can be described by the random behavior of a complex vectorA=(AX,AY), called the Jones vector, whose covariance matrix, called the polarization matrix, is
(25)Γ=(𝔼[AXAX*]𝔼[AXAY*]𝔼[AYAX*]𝔼[AYAY*])≜(a1a3+ia4a3-ia4a2),
where* denotes the complex conjugate. The covariance matrix Γ is a nonnegative Hermitian matrix whose diagonal terms are the intensity components in the X and Y directions. The cross terms of Γ are the correlations between the Jones components. If we assume a fully developed speckle, the Jones vector A is distributed according to a complex Gaussian distribution whose probability density function (pdf) is [15]:
(26)p(A)=1π2|Γ|exp(-A†Γ-1A),
where |Γ| is the determinant of the matrix Γ and † denotes the conjugate transpose operator. As a consequence, the statistical properties of A are fully characterized by the covariance matrix Γ. The different components of Γ can be classically estimated by using four intensity images that are related to the components of the Jones vector as follows (see [20], for more details):
(27)I1=|AX|2,I2=|AY|2,I3=12|AX|2+12|AY|2+Re(AXAY*),I4=12|AX|2+12|AY|2+Im(AXAY*).
The state of the polarization of the light is classically characterized by the square DoP defined by [15, pages 134–136]
(28)P2=1-4|Γ|[trace(Γ)]2=1-4[a1a2-(a32+a42)](a1+a2)2,
where trace(Γ) is the trace of the matrix Γ. The light is totally depolarized for P=0, totally polarized for P=1, and partially polarized when P∈]0,1[. As a consequence, estimating the square DoP of a polarimetric image is important in many practical applications. Different estimation methods of P2 using several combinations of intensity images were studied in [20]. Since only one realization of the random vector I=(I1,…,I4)T was available for a given pixel of a polarimetric image, the image was supposed to be locally stationary and ergodic. These assumptions were used to derive square DoP estimators using several neighbor pixels belonging to a so-called estimation window.This paper considers practical applications where the intensity level of the reflected light is very low (low flux assumption), which leads to an additional source of fluctuations on the detected signal. Under the low flux assumption, the quantum nature of the light leads to a Poisson-distributed noise which can become very important relatively to the mean value of the signal at a low photon level. As a consequence, the observed pixels of the low flux polarimetric image are discrete random variables contained in the vectorN=(N1,…,N4) such that the conditional distributions of the random variables Nl∣Il, for l=1,…,4 are independent and distributed according to Poisson distributions with means Il, for l=1,…,4. The resulting joint distribution of N is a multivariate mixed Poisson distribution [22]:
(29)P(N=k)=∫‍⋯(ℝ+)4∫‍∏l=14‍Ilklkl!exp(-Il)f(I)dI,
where k=(k1,…,k4), ki∈ℕ, and f(I) is the joint pdf of the intensity vector. This section studies estimators of the square DoP P2 defined in (28) based on several vectors N1,…,Nn belonging to the estimation window. These estimators are constructed from estimates of the covariance matrix elements ai, i=1,…,4. As explained in the introduction, several studies have been recently devoted to the estimation of the square DoP using less than four polarimetric images. This paper goes into this direction by deriving estimators based on the observation of 2 or 3 polarimetric images.The joint distribution of the intensity vectorI is known to be a multivariate gamma distribution whose Laplace transform is [20]
(30)E[exp(∑k=14‍zkIk)]=1P4(z),
where the affine polynomial P4 is
(31)P4(z)=1+zμ+ka[2z1z2+z3z4+(z1+z2)(z3+z4)]
with z=(z1,…,z4) and
(32)ka=12(a1a2-a32-a42),μ=(a1,a2,a3+a1+a22,a4+a1+a22)T.
As a consequence, the distribution of N is an NMD whose PGF can be written as (the interested reader is invited to consult [22] for more details)(33)GN(z)=1P4(z1-1,z2-1,z3-1,z4-1).
The results of Section 2 allow us to compute the masses of N that will be useful for studying the maximum likelihood estimator (MLE) of the square DoP.


## 4.2. MLE Using Three Polarimetric Images

The PGF ofN~=(N1,N2,N3) can be computed from (33) by setting z4=1. The following result can be obtained:
(34)GN~(z)=1P3(z)
with
(35)P3(z)=P3(0)+z1(μ1-3ka)+z2(μ2-3ka)+z3(μ3-2ka)+ka(2z1z2+z1z3+z2z3),z=(z1,z2,z3), and P3(0)=1-∑i=13‍μi+4ka. The results of Section 3.2 can then be used to express the masses of N~ as a function of θ=(a1,a2,a3,a42)T.The ML estimator ofθ based on several vectors N~k belonging to the estimation window (where k=1,…,K and K is the number of pixels of the observation window) is obtained by maximizing the log-likelihood
(36)l3(N~1,…,N~K∣θ)=∑k=1K‍log[P(N~k)]
with respect to θ (note again that P(N~k) is the mass of a trivariate NMD that has been computed in Section 3.2). The practical determination of the ML estimator of θ is achieved by using a Newton-Raphson procedure. The ML estimators of the vector θ, denoted as θ~=(a~1,a~2,a~3,a42~)T, are then plugged into (28) to provide the ML estimator of the square DoP based on three polarimetric images: (37)P~2=1-4[a~1a~2-(a~32+a42~)](a~1+a~2)2.


## 4.3. MLE Using Two Polarimetric Images

The PGF ofN_=(N1,N2) can be computed from (34) by setting z3=1. The following result can be obtained:
(38)GN_(z)=1P2(z)
with
(39)P2(z)=P2(0)+z1(μ1-2ka)+z2(μ2-2ka)+2kaz1z2,z=(z1,z2) and P2(0)=1-∑i=12‍μi+2ka. The results of Section 3.1 can then be used to express the masses of N_ as a function of θ=(a1,a2,ka)T.The MLE ofθ based on several vectors N_k belonging to the estimation window is obtained by maximizing the log-likelihood
(40)l2(N_1,…,N_K∣θ)=∑k=1K‍log[P(N_k)]
with respect to θ (note that P(N_k) is the mass of a bivariate NMD that has been computed in Section 3.1). The practical determination of the ML estimator of θ is achieved by using a Newton-Raphson procedure. The ML estimators of the vector θ elements, denoted as θ_=(a_1,a_2,ka_)T, are then plugged into (28) to provide the MLE of the square DoP based on two polarimetric images:
(41)P_2=1-8ka_(a_1+a_2)2.


## 5. Simulation Results

### 5.1. Estimation Performance

The performance of the ML estimators of the square DoP based on two or three polarimetric images has been evaluated via several experiments. The first simulations compare the log Mean Square Errors (MSEs) of the square DoP estimators constructed from two or three images. Eleven different covariance matrices of the Jones vector have been considered in order to define typical values of the DoP. The values ofai, for i=1,…,4, (defining the covariance matrix elements of the Jones vector) and the corresponding values of the square DoPs are reported in Table 1. Note that all the covariance matrices Γj, j=0,…,10, have been normalized so that the mean of the total intensity E[I1+I2]=E[N1+N2]=a1+a2 is equal to 4. Thus, the average number of photons collected on each pixel equals 4 for each matrix of the considered set. This point is in agreement with the low flux assumption.Table 1
Covariance matrix elements and square DoP values for the Jones vector.
Γ
0
Γ
1
Γ
2
Γ
3
Γ
4
Γ
5
Γ
6
Γ
7
Γ
8
Γ
9
Γ
10
a
1
2.00
1.71
2.40
1.82
2.00
2.67
2.40
2.24
2.74
2.00
2.00
a
2
2.00
2.29
1.60
2.18
2.00
1.33
1.60
1.76
1.26
2.00
2.00
a
3
0.00
0.40
0.48
0.91
0.89
0.96
1.07
1.05
1.46
0.60
1.41
a
4
0.00
0.40
0.64
0.58
0.89
0.80
1.05
1.28
0.73
1.80
1.41
P
2
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0Figure1(a) display the empirical log MSEs of the square DoP estimators for the set of covariance matrices defined in Table 1 as a function of the true square DoP value. The red plus markers + correspond to the estimators obtained for two polarimetric images (2D MLE), whereas the blue cross markers × correspond to the MLE obtained for three polarimetric images (3D MLE). Note that these empirical MSEs have been computed for a square observation window of size n=51×51 pixels, based on 1000 Monte-Carlo runs. The theoretical asymptotic log MSEs of the MLE are also displayed in Figure 1(a) with continuous lines. These asymptotic values correspond to the Cramer-Rao Lower Bound (CRLB) for the parameter P2. The MLE is known to be asymptotically unbiased and efficient under mild regularity conditions (that are satisfied for P2). Thus, the MSE of the estimates can be approximated a for large sample by the CRLB. More details about the way of computing the square DoP CLRBs can be found in [23]. Figure 1(a) indicates that the empirical MSEs are in good agreement with the corresponding CRLBs, except for the matrices Γ0 and Γ10. Indeed, the CRLBs for these two matrices cannot be computed since the true value of the parameters belongs to the boundary of its definition domain. The empirical bias, standard deviations (“std”), MSEs, and asymptotic variances (“avar”) of the estimators of P2 are also reported in Table 2. It is interesting to note that the MLE obtained using 3 images is slightly more biased than the one obtained using 2 images. However, the MLE based on 3 images provides lower MSEs than the estimator based on 2 images, as expected.Table 2
Simulation results for the estimation ofP2 using 2 or 3 images, obtained from 1000 Monte-Carlo runs (n=51×51).
2 images MLE
3 images MLE
Bias
std
MSE
avar
Bias
std
MSE
avar
Γ
0
3.57
e
-
03
1.31
e
-
02
1.85
e
-
04
—
4.12
e
-
03
9.08
e
-
03
9.93
e
-
05
—
Γ
1
4.28
e
-
04
3.14
e
-
02
9.83
e
-
04
9.78
e
-
04
1.71
e
-
03
2.57
e
-
02
6.62
e
-
04
7.25
e
-
04
Γ
2
9.49
e
-
04
3.23
e
-
02
1.04
e
-
03
9.82
e
-
04
1.20
e
-
03
2.86
e
-
02
8.18
e
-
04
7.83
e
-
04
Γ
3
-
2.15
e
-
04
3.13
e
-
02
9.81
e
-
04
9.67
e
-
04
4.84
e
-
04
2.88
e
-
02
8.28
e
-
04
8.07
e
-
04
Γ
4
2.49
e
-
03
3.02
e
-
02
9.19
e
-
04
8.91
e
-
04
1.71
e
-
03
2.81
e
-
02
7.91
e
-
04
7.67
e
-
04
Γ
5
2.86
e
-
04
2.85
e
-
02
8.10
e
-
04
7.54
e
-
04
4.80
e
-
04
2.65
e
-
02
7.03
e
-
04
6.69
e
-
04
Γ
6
5.25
e
-
04
2.59
e
-
02
6.72
e
-
04
6.58
e
-
04
6.32
e
-
04
2.47
e
-
02
6.11
e
-
04
5.95
e
-
04
Γ
7
1.05
e
-
03
2.37
e
-
02
5.60
e
-
04
5.43
e
-
04
1.04
e
-
03
2.26
e
-
02
5.11
e
-
04
4.95
e
-
04
Γ
8
1.59
e
-
04
2.01
e
-
02
4.03
e
-
04
4.12
e
-
04
3.17
e
-
04
1.90
e
-
02
3.60
e
-
04
3.73
e
-
04
Γ
9
-
2.88
e
-
04
1.76
e
-
02
3.10
e
-
04
3.20
e
-
04
-
4.46
e
-
04
1.65
e
-
02
2.73
e
-
04
2.86
e
-
04
Γ
10
-
5.98
e
-
03
8.98
e
-
03
1.16
e
-
04
—
-
5.73
e
-
03
8.58
e
-
03
1.06
e
-
04
—log MSEs of the square DoP estimates using 2 and 3 images versus P2 for the set of polarization matrices defined in Table 1 under (a) low flux and (b) high flux assumptions (n=51×51, ML: maximum likelihood estimators, and Asympt.: theoretical asymptotic value of the log MSE for a given estimator).
(a)
Low flux
(b)
High fluxIn order to appreciate the influence of the Poisson noise due to the low flux assumption, experiments have been conducted using the high flux assumption. In this case, the intensity vectorI is assumed to be known. Thus, the high flux MLEs using two and three images can be derived from the intensity vectors I_k=(I1k,I2k)T and I~k=(I1k,I2k,I3k)T [20]. The results are depicted in Figure 1(b). A comparison between Figures 1(a) and 1(b) allows one to appreciate a similar global behaviour for all the estimators, with a maximum MSE near P2=1/3 and decreasing MSEs as P2 goes to 0 or 1. The degradation of the estimation performance due to the presence of Poisson noise (due to the low flux assumption) can also be clearly noticed.The next set of simulations studies the performance of the different estimators as a function of the sample sizen. Figures 2(a) and 2(b) show the log MSEs of the square DoP estimates obtained for 2 and 3 images (for the two particular matrices Γ2 and Γ8). The empirical bias, standard deviations (“std”), MSEs, and asymptotic variances (“avar”) are also reported in Tables 3 and 4. One can see that the empirical MSEs are in good agreement with their theoretical asymptotic values for a large enough sample size, that is, n>25×25. Moreover, the gain of a performance using 3 images instead of 2 is more significant for small values of P2 (indeed, the difference between the different curves is more pronounced in the left figure corresponding to P2=0.2 than in the right figure corresponding to P2=0.8.)Table 3
Simulation results for the estimation ofP2 using 2 or 3 images, obtained from 1000 Monte-Carlo runs (matrix Γ2).
n
2 images MLE
3 images MLE
Bias
std
MSE
avar
Bias
std
MSE
avar
5
2
9.03
e
-
02
2.49
e
-
01
7.01
e
-
02
1.02
e
-
01
1.20
e
-
01
2.20
e
-
01
6.27
e
-
02
8.15
e
-
02
7
2
5.10
e
-
02
1.97
e
-
01
4.13
e
-
02
5.21
e
-
02
7.04
e
-
02
1.69
e
-
01
3.36
e
-
02
4.16
e
-
02
1
1
2
2.60
e
-
02
1.39
e
-
01
2.01
e
-
02
2.11
e
-
02
3.22
e
-
02
1.14
e
-
01
1.39
e
-
02
1.68
e
-
02
1
5
2
7.68
e
-
03
1.02
e
-
01
1.04
e
-
02
1.14
e
-
02
1.59
e
-
02
8.48
e
-
02
7.43
e
-
03
9.05
e
-
03
2
5
2
-
7.42
e
-
04
6.25
e
-
02
3.90
e
-
03
4.09
e
-
03
2.16
e
-
03
5.36
e
-
02
2.87
e
-
03
3.26
e
-
03
3
7
2
-
5.19
e
-
04
4.32
e
-
02
1.87
e
-
03
1.87
e
-
03
6.87
e
-
05
3.79
e
-
02
1.44
e
-
03
1.49
e
-
03
5
1
2
9.49
e
-
04
3.23
e
-
02
1.04
e
-
03
9.82
e
-
04
1.20
e
-
03
2.86
e
-
02
8.18
e
-
04
7.83
e
-
04Table 4
Simulation results for the estimation ofP2 using 2 or 3 images, obtained from 1000 Monte-Carlo runs (matrix Γ8).
n
2 images MLE
3 images MLE
Bias
std
MSE
avar
Bias
std
MSE
avar
5
2
-
3.91
e
-
03
2.08
e
-
01
4.31
e
-
02
4.29
e
-
02
2.80
e
-
02
1.61
e
-
01
2.67
e
-
02
3.88
e
-
02
7
2
4.04
e
-
03
1.44
e
-
01
2.08
e
-
02
2.19
e
-
02
1.89
e
-
02
1.24
e
-
01
1.56
e
-
02
1.98
e
-
02
1
1
2
2.14
e
-
03
9.38
e
-
02
8.80
e
-
03
8.86
e
-
03
8.70
e
-
03
8.47
e
-
02
7.24
e
-
03
8.01
e
-
03
1
5
2
-
1.27
e
-
03
7.14
e
-
02
5.10
e
-
03
4.76
e
-
03
2.70
e
-
03
6.63
e
-
02
4.39
e
-
03
4.31
e
-
03
2
5
2
5.23
e
-
04
3.98
e
-
02
1.58
e
-
03
1.71
e
-
03
6.87
e
-
04
3.75
e
-
02
1.40
e
-
03
1.55
e
-
03
3
7
2
2.24
e
-
03
2.89
e
-
02
8.42
e
-
04
7.83
e
-
04
2.54
e
-
03
2.75
e
-
02
7.61
e
-
04
7.08
e
-
04
5
1
2
1.59
e
-
04
2.01
e
-
02
4.03
e
-
04
4.12
e
-
04
3.17
e
-
04
1.90
e
-
02
3.60
e
-
04
3.73
e
-
04log MSE of the estimated square DoP P2 using 2 or 3 intensity images versus the logarithm of the sample size for the matrices (a) Γ2 and (b) Γ8 (ML: maximum likelihood estimators, and Asympt. theoretical asymptotic value of the log MSE for a given estimator).
(a)
MatrixΓ2
(b)
MatrixΓ8Figures3(a) and 3(b) display the MSEs of the MLE under the high flux assumption. By comparing Figures 2 and 3, one can observe that the gain of performance using 3 images instead of 2 is more important in the high flux scenario than under a low flux assumption.log MSE of the estimated square DoP P2 using 2 or 3 intensity images versus the logarithm of the sample size for the matrices (a) Γ2 and (b) Γ8 under high flux assumption (ML: maximum likelihood estimators, Asympt.: theoretical asymptotic value of the log MSE for a given estimator).
(a)
MatrixΓ2
(b)
MatrixΓ8


### 5.2. Application to Synthetic Polarimetric Images

In order to appreciate the estimation performance on polarimetric images, we consider a synthetic polarimetric image of size512×512 composed of three distinct objects located on a homogeneous background depicted in Figure 4 (see also [20]). The polarimetric properties of these objects and background (i.e., the covariance matrix of the Jones vector and the square DoPs) are reported in Table 5. The polarimetric low flux images generated according to this model are also represented in Figure 5 in negative colors (bright pixels correspond to a small number of photons, dark ones correspond to a large number of photons). Note that these images represent the values of the vector N=(N1,N2,N3,N4)T for each pixel. The square DoP of each pixel x(i,j) (for i,j=1,…,512) has been estimated from vectors belonging to windows of size n=15×15 centered around the pixel of coordinates (i,j) in the analyzed image. The estimated square DoPs are depicted in Figures 6(a) and 6(b) for the MLEs using 2 and 3 images, respectively. One can see that the ML method for 3 images gives more homogeneous results on each object than the ML method derived for 2 images. This result confirms that the MLE using 3 images performs better than the MLE using only 2 images. However, the polarimetric properties of a scene can be clearly recovered with 2 or 3 low flux images.Table 5
Polarimetric properties of elements that compose the scene displayed in Figure4.
Object
Polarization matrixΓ
P
2
Remarks
Background
(
0.79
0
0
0.98
)
0.0115
Very depolarizing and dark background
1
(
3.6
0
0
0.22
)
0.783
Very bright and weakly depolarizing object (typically steel)
2
(
3
0.1
0.1
0.6
)
0.447
Bright object quite depolarizing
3
(
0.7
0.5
+
0.2
i
0.5
-
0.2
i
1.07
)
0.414
Dark object whose mean total intensity is the same as the backgroundComposition of the scene used to generate synthetic polarimetric low flux images and associated theoretical squared DoP.
(a)
Scene
(b)
Theoretical squared DoPSynthetic intensity images (negative colors) for the scene depicted in Figure4 and described in Table 5.
(a)
Low flux intensityN1
(b)
Low flux intensityN2
(c)
Low flux intensityN3
(d)
Low flux intensityN4Estimates ofP2 using 2 or 3 low flux intensity images for the synthetic polarimetric images for an estimation window of size n=15×15. MLE: maximum likelihood estimator.
(a)
MLE2 images
(b)
MLE3 images


### 5.3. Estimation Results on Real-Word Polarimetric Images

Finally, the ML estimator based on three images is applied on real polarimetric data. These images are acquired by using a laser as a coherent illumination source. The scene consists of two disks. The first one, intended to provide low DoP, is a grey diffuse material (left object in Figure7). The second one is made of sand blasted aluminium providing high DoP (right object in Figure 7). Due to the experimental conditions, the measured intensities are quite low. As a consequence, these intensities are assumed to be distributed according to an NMD. The intensity images corresponding to N1, N2, N3, and N4 are depicted in Figure 7. The interested reader is invited to read [20] for more details on these data. It is important to note that the two disks exhibit the similar level of total reflectivity N1+N2 and can hardly be distinguished without a polarimetric processing.Real-world polarimetric intensity images of a scene composed of a plastic disk (left) and a steel disk (right).
(a)
Low flux intensityN1
(b)
Low flux intensityN2
(c)
Low flux intensityN3
(d)
Low flux intensityN4Figure8 shows the ML estimates of the square DoPs P2 for 3 images and an estimation window of size n=9×9 pixels. As expected, the values of the estimates are quite different on each disk: quite higher on the metal than on the plastic disk. This result is in good agreement with the theoretical properties of the considered material. This emphasizes the interest in considering efficient estimators based on NMDs for polarimetric images.Figure 8
Estimates ofP2 using the 3 low flux intensity images N1, N2, and N3 for the real polarimetric images for an estimation window of size n=9×9. “MLE”: maximum likelihood estimator.


## 5.1. Estimation Performance

The performance of the ML estimators of the square DoP based on two or three polarimetric images has been evaluated via several experiments. The first simulations compare the log Mean Square Errors (MSEs) of the square DoP estimators constructed from two or three images. Eleven different covariance matrices of the Jones vector have been considered in order to define typical values of the DoP. The values ofai, for i=1,…,4, (defining the covariance matrix elements of the Jones vector) and the corresponding values of the square DoPs are reported in Table 1. Note that all the covariance matrices Γj, j=0,…,10, have been normalized so that the mean of the total intensity E[I1+I2]=E[N1+N2]=a1+a2 is equal to 4. Thus, the average number of photons collected on each pixel equals 4 for each matrix of the considered set. This point is in agreement with the low flux assumption.Table 1
Covariance matrix elements and square DoP values for the Jones vector.
Γ
0
Γ
1
Γ
2
Γ
3
Γ
4
Γ
5
Γ
6
Γ
7
Γ
8
Γ
9
Γ
10
a
1
2.00
1.71
2.40
1.82
2.00
2.67
2.40
2.24
2.74
2.00
2.00
a
2
2.00
2.29
1.60
2.18
2.00
1.33
1.60
1.76
1.26
2.00
2.00
a
3
0.00
0.40
0.48
0.91
0.89
0.96
1.07
1.05
1.46
0.60
1.41
a
4
0.00
0.40
0.64
0.58
0.89
0.80
1.05
1.28
0.73
1.80
1.41
P
2
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0Figure1(a) display the empirical log MSEs of the square DoP estimators for the set of covariance matrices defined in Table 1 as a function of the true square DoP value. The red plus markers + correspond to the estimators obtained for two polarimetric images (2D MLE), whereas the blue cross markers × correspond to the MLE obtained for three polarimetric images (3D MLE). Note that these empirical MSEs have been computed for a square observation window of size n=51×51 pixels, based on 1000 Monte-Carlo runs. The theoretical asymptotic log MSEs of the MLE are also displayed in Figure 1(a) with continuous lines. These asymptotic values correspond to the Cramer-Rao Lower Bound (CRLB) for the parameter P2. The MLE is known to be asymptotically unbiased and efficient under mild regularity conditions (that are satisfied for P2). Thus, the MSE of the estimates can be approximated a for large sample by the CRLB. More details about the way of computing the square DoP CLRBs can be found in [23]. Figure 1(a) indicates that the empirical MSEs are in good agreement with the corresponding CRLBs, except for the matrices Γ0 and Γ10. Indeed, the CRLBs for these two matrices cannot be computed since the true value of the parameters belongs to the boundary of its definition domain. The empirical bias, standard deviations (“std”), MSEs, and asymptotic variances (“avar”) of the estimators of P2 are also reported in Table 2. It is interesting to note that the MLE obtained using 3 images is slightly more biased than the one obtained using 2 images. However, the MLE based on 3 images provides lower MSEs than the estimator based on 2 images, as expected.Table 2
Simulation results for the estimation ofP2 using 2 or 3 images, obtained from 1000 Monte-Carlo runs (n=51×51).
2 images MLE
3 images MLE
Bias
std
MSE
avar
Bias
std
MSE
avar
Γ
0
3.57
e
-
03
1.31
e
-
02
1.85
e
-
04
—
4.12
e
-
03
9.08
e
-
03
9.93
e
-
05
—
Γ
1
4.28
e
-
04
3.14
e
-
02
9.83
e
-
04
9.78
e
-
04
1.71
e
-
03
2.57
e
-
02
6.62
e
-
04
7.25
e
-
04
Γ
2
9.49
e
-
04
3.23
e
-
02
1.04
e
-
03
9.82
e
-
04
1.20
e
-
03
2.86
e
-
02
8.18
e
-
04
7.83
e
-
04
Γ
3
-
2.15
e
-
04
3.13
e
-
02
9.81
e
-
04
9.67
e
-
04
4.84
e
-
04
2.88
e
-
02
8.28
e
-
04
8.07
e
-
04
Γ
4
2.49
e
-
03
3.02
e
-
02
9.19
e
-
04
8.91
e
-
04
1.71
e
-
03
2.81
e
-
02
7.91
e
-
04
7.67
e
-
04
Γ
5
2.86
e
-
04
2.85
e
-
02
8.10
e
-
04
7.54
e
-
04
4.80
e
-
04
2.65
e
-
02
7.03
e
-
04
6.69
e
-
04
Γ
6
5.25
e
-
04
2.59
e
-
02
6.72
e
-
04
6.58
e
-
04
6.32
e
-
04
2.47
e
-
02
6.11
e
-
04
5.95
e
-
04
Γ
7
1.05
e
-
03
2.37
e
-
02
5.60
e
-
04
5.43
e
-
04
1.04
e
-
03
2.26
e
-
02
5.11
e
-
04
4.95
e
-
04
Γ
8
1.59
e
-
04
2.01
e
-
02
4.03
e
-
04
4.12
e
-
04
3.17
e
-
04
1.90
e
-
02
3.60
e
-
04
3.73
e
-
04
Γ
9
-
2.88
e
-
04
1.76
e
-
02
3.10
e
-
04
3.20
e
-
04
-
4.46
e
-
04
1.65
e
-
02
2.73
e
-
04
2.86
e
-
04
Γ
10
-
5.98
e
-
03
8.98
e
-
03
1.16
e
-
04
—
-
5.73
e
-
03
8.58
e
-
03
1.06
e
-
04
—log MSEs of the square DoP estimates using 2 and 3 images versus P2 for the set of polarization matrices defined in Table 1 under (a) low flux and (b) high flux assumptions (n=51×51, ML: maximum likelihood estimators, and Asympt.: theoretical asymptotic value of the log MSE for a given estimator).
(a)
Low flux
(b)
High fluxIn order to appreciate the influence of the Poisson noise due to the low flux assumption, experiments have been conducted using the high flux assumption. In this case, the intensity vectorI is assumed to be known. Thus, the high flux MLEs using two and three images can be derived from the intensity vectors I_k=(I1k,I2k)T and I~k=(I1k,I2k,I3k)T [20]. The results are depicted in Figure 1(b). A comparison between Figures 1(a) and 1(b) allows one to appreciate a similar global behaviour for all the estimators, with a maximum MSE near P2=1/3 and decreasing MSEs as P2 goes to 0 or 1. The degradation of the estimation performance due to the presence of Poisson noise (due to the low flux assumption) can also be clearly noticed.The next set of simulations studies the performance of the different estimators as a function of the sample sizen. Figures 2(a) and 2(b) show the log MSEs of the square DoP estimates obtained for 2 and 3 images (for the two particular matrices Γ2 and Γ8). The empirical bias, standard deviations (“std”), MSEs, and asymptotic variances (“avar”) are also reported in Tables 3 and 4. One can see that the empirical MSEs are in good agreement with their theoretical asymptotic values for a large enough sample size, that is, n>25×25. Moreover, the gain of a performance using 3 images instead of 2 is more significant for small values of P2 (indeed, the difference between the different curves is more pronounced in the left figure corresponding to P2=0.2 than in the right figure corresponding to P2=0.8.)Table 3
Simulation results for the estimation ofP2 using 2 or 3 images, obtained from 1000 Monte-Carlo runs (matrix Γ2).
n
2 images MLE
3 images MLE
Bias
std
MSE
avar
Bias
std
MSE
avar
5
2
9.03
e
-
02
2.49
e
-
01
7.01
e
-
02
1.02
e
-
01
1.20
e
-
01
2.20
e
-
01
6.27
e
-
02
8.15
e
-
02
7
2
5.10
e
-
02
1.97
e
-
01
4.13
e
-
02
5.21
e
-
02
7.04
e
-
02
1.69
e
-
01
3.36
e
-
02
4.16
e
-
02
1
1
2
2.60
e
-
02
1.39
e
-
01
2.01
e
-
02
2.11
e
-
02
3.22
e
-
02
1.14
e
-
01
1.39
e
-
02
1.68
e
-
02
1
5
2
7.68
e
-
03
1.02
e
-
01
1.04
e
-
02
1.14
e
-
02
1.59
e
-
02
8.48
e
-
02
7.43
e
-
03
9.05
e
-
03
2
5
2
-
7.42
e
-
04
6.25
e
-
02
3.90
e
-
03
4.09
e
-
03
2.16
e
-
03
5.36
e
-
02
2.87
e
-
03
3.26
e
-
03
3
7
2
-
5.19
e
-
04
4.32
e
-
02
1.87
e
-
03
1.87
e
-
03
6.87
e
-
05
3.79
e
-
02
1.44
e
-
03
1.49
e
-
03
5
1
2
9.49
e
-
04
3.23
e
-
02
1.04
e
-
03
9.82
e
-
04
1.20
e
-
03
2.86
e
-
02
8.18
e
-
04
7.83
e
-
04Table 4
Simulation results for the estimation ofP2 using 2 or 3 images, obtained from 1000 Monte-Carlo runs (matrix Γ8).
n
2 images MLE
3 images MLE
Bias
std
MSE
avar
Bias
std
MSE
avar
5
2
-
3.91
e
-
03
2.08
e
-
01
4.31
e
-
02
4.29
e
-
02
2.80
e
-
02
1.61
e
-
01
2.67
e
-
02
3.88
e
-
02
7
2
4.04
e
-
03
1.44
e
-
01
2.08
e
-
02
2.19
e
-
02
1.89
e
-
02
1.24
e
-
01
1.56
e
-
02
1.98
e
-
02
1
1
2
2.14
e
-
03
9.38
e
-
02
8.80
e
-
03
8.86
e
-
03
8.70
e
-
03
8.47
e
-
02
7.24
e
-
03
8.01
e
-
03
1
5
2
-
1.27
e
-
03
7.14
e
-
02
5.10
e
-
03
4.76
e
-
03
2.70
e
-
03
6.63
e
-
02
4.39
e
-
03
4.31
e
-
03
2
5
2
5.23
e
-
04
3.98
e
-
02
1.58
e
-
03
1.71
e
-
03
6.87
e
-
04
3.75
e
-
02
1.40
e
-
03
1.55
e
-
03
3
7
2
2.24
e
-
03
2.89
e
-
02
8.42
e
-
04
7.83
e
-
04
2.54
e
-
03
2.75
e
-
02
7.61
e
-
04
7.08
e
-
04
5
1
2
1.59
e
-
04
2.01
e
-
02
4.03
e
-
04
4.12
e
-
04
3.17
e
-
04
1.90
e
-
02
3.60
e
-
04
3.73
e
-
04log MSE of the estimated square DoP P2 using 2 or 3 intensity images versus the logarithm of the sample size for the matrices (a) Γ2 and (b) Γ8 (ML: maximum likelihood estimators, and Asympt. theoretical asymptotic value of the log MSE for a given estimator).
(a)
MatrixΓ2
(b)
MatrixΓ8Figures3(a) and 3(b) display the MSEs of the MLE under the high flux assumption. By comparing Figures 2 and 3, one can observe that the gain of performance using 3 images instead of 2 is more important in the high flux scenario than under a low flux assumption.log MSE of the estimated square DoP P2 using 2 or 3 intensity images versus the logarithm of the sample size for the matrices (a) Γ2 and (b) Γ8 under high flux assumption (ML: maximum likelihood estimators, Asympt.: theoretical asymptotic value of the log MSE for a given estimator).
(a)
MatrixΓ2
(b)
MatrixΓ8


## 5.2. Application to Synthetic Polarimetric Images

In order to appreciate the estimation performance on polarimetric images, we consider a synthetic polarimetric image of size512×512 composed of three distinct objects located on a homogeneous background depicted in Figure 4 (see also [20]). The polarimetric properties of these objects and background (i.e., the covariance matrix of the Jones vector and the square DoPs) are reported in Table 5. The polarimetric low flux images generated according to this model are also represented in Figure 5 in negative colors (bright pixels correspond to a small number of photons, dark ones correspond to a large number of photons). Note that these images represent the values of the vector N=(N1,N2,N3,N4)T for each pixel. The square DoP of each pixel x(i,j) (for i,j=1,…,512) has been estimated from vectors belonging to windows of size n=15×15 centered around the pixel of coordinates (i,j) in the analyzed image. The estimated square DoPs are depicted in Figures 6(a) and 6(b) for the MLEs using 2 and 3 images, respectively. One can see that the ML method for 3 images gives more homogeneous results on each object than the ML method derived for 2 images. This result confirms that the MLE using 3 images performs better than the MLE using only 2 images. However, the polarimetric properties of a scene can be clearly recovered with 2 or 3 low flux images.Table 5
Polarimetric properties of elements that compose the scene displayed in Figure4.
Object
Polarization matrixΓ
P
2
Remarks
Background
(
0.79
0
0
0.98
)
0.0115
Very depolarizing and dark background
1
(
3.6
0
0
0.22
)
0.783
Very bright and weakly depolarizing object (typically steel)
2
(
3
0.1
0.1
0.6
)
0.447
Bright object quite depolarizing
3
(
0.7
0.5
+
0.2
i
0.5
-
0.2
i
1.07
)
0.414
Dark object whose mean total intensity is the same as the backgroundComposition of the scene used to generate synthetic polarimetric low flux images and associated theoretical squared DoP.
(a)
Scene
(b)
Theoretical squared DoPSynthetic intensity images (negative colors) for the scene depicted in Figure4 and described in Table 5.
(a)
Low flux intensityN1
(b)
Low flux intensityN2
(c)
Low flux intensityN3
(d)
Low flux intensityN4Estimates ofP2 using 2 or 3 low flux intensity images for the synthetic polarimetric images for an estimation window of size n=15×15. MLE: maximum likelihood estimator.
(a)
MLE2 images
(b)
MLE3 images


## 5.3. Estimation Results on Real-Word Polarimetric Images

Finally, the ML estimator based on three images is applied on real polarimetric data. These images are acquired by using a laser as a coherent illumination source. The scene consists of two disks. The first one, intended to provide low DoP, is a grey diffuse material (left object in Figure7). The second one is made of sand blasted aluminium providing high DoP (right object in Figure 7). Due to the experimental conditions, the measured intensities are quite low. As a consequence, these intensities are assumed to be distributed according to an NMD. The intensity images corresponding to N1, N2, N3, and N4 are depicted in Figure 7. The interested reader is invited to read [20] for more details on these data. It is important to note that the two disks exhibit the similar level of total reflectivity N1+N2 and can hardly be distinguished without a polarimetric processing.Real-world polarimetric intensity images of a scene composed of a plastic disk (left) and a steel disk (right).
(a)
Low flux intensityN1
(b)
Low flux intensityN2
(c)
Low flux intensityN3
(d)
Low flux intensityN4Figure8 shows the ML estimates of the square DoPs P2 for 3 images and an estimation window of size n=9×9 pixels. As expected, the values of the estimates are quite different on each disk: quite higher on the metal than on the plastic disk. This result is in good agreement with the theoretical properties of the considered material. This emphasizes the interest in considering efficient estimators based on NMDs for polarimetric images.Figure 8
Estimates ofP2 using the 3 low flux intensity images N1, N2, and N3 for the real polarimetric images for an estimation window of size n=9×9. “MLE”: maximum likelihood estimator.



---

*Source: 170967-2013-02-24.xml*