Masses of Negative Multinomial Distributions: Application to Polarimetric Image Processing

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.

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 |𝐈−𝐐|^𝜆|𝐈−𝐐𝐙|^−𝜆, where 𝜆>0, 𝐐 is an 𝑛×𝑛 matrix, and 𝐙=diag(𝑧₁,…,𝑧_𝑛). In particular, matrices 𝐐 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 ℕ^𝑛₀ (where ℕ₀ 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 [𝑛]={1,…,𝑛} the set of the 𝑛 first nonzero integers. We denote 𝐳^𝑇=∏_𝑡∈𝑇𝑧_𝑡 as the monomial obtained by multiplying all the entries of the vector 𝐳=(𝑧₁,…,𝑧_𝑛)∈ℝ^𝑛 whose indexes belong to 𝑇, where 𝑇⊂[𝑛] stands for any subset of the indexes. Let 𝑃_𝑛(𝐳)=∑_{𝑇⊂[𝑛],𝑇≠∅}𝑝_𝑇𝐳^𝑇 be an affine polynomial with respect to the 𝑛 variables (𝑧₁,…,𝑧_𝑛) such that 1−𝑃_𝑛(𝟏)≠0. The NMD distribution defined at pair (𝑛,𝑃_𝑛) is represented by its PFG which is given by 𝒢_{𝒩ℳ(𝑛,𝑃𝑛)}(𝐳)=1−𝑃_𝑛(𝐳)^−𝜆1−𝑃_𝑛(𝟏)^𝜆. (1) Such laws are denoted as 𝒩ℳ(𝑛,𝑃_𝑛). However, as explained in [6], all couples (𝑛,𝑃_𝑛) do not provide a valid NMD. More specifically, Bernardoff has derived a finite number of conditions over 𝑃_𝑛 such that (1−𝑃_𝑛(𝐳))^−𝜆(1−𝑃_𝑛(𝟏))^𝜆 is the PGF of an NMD for all positive integer 𝑛. The corresponding expression of the coefficient of 𝐳^𝜶 in the Taylor expansion of (1−𝑃_𝑛(𝐳))^−𝜆 is given by the formula 𝑐_𝛼𝜆,𝑃_𝑛=_{𝑘∈𝐾𝛼}(𝜆)_|𝑘|𝐩^𝑘,𝑘! (2) where 𝐾_𝛼={𝑘∶𝒫_𝑛→ℕ} and 𝒫_𝑛 is the set of all subsets of [𝑛]. However, this expression of 𝑐_𝛼(𝜆,𝑃_𝑛) 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 𝒩ℳ(𝑛,𝑃_𝑛) 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. Section 2 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

An 𝑛-variate NMD is the distribution of a random vector 𝐍=(𝑁₁,…,𝑁_𝑛) taking its values in ℕ^𝑛₀ whose PGF is 𝐺_𝐍(𝐳)=𝔼^𝑛_𝑘=1𝑧^𝑁𝑘_𝑘=𝑃_𝑛(𝐳)^−𝜆, (3) where 𝔼 denotes the mathematical expectation, 𝐳=(𝑧₁,…,𝑧_𝑛),𝜆>0, and 𝑃_𝑛(𝐳) is an affine polynomial of order 𝑛. (A polynomial 𝑃_𝑛(𝐳) with respect to 𝐳=(𝑧₁,…,𝑧_𝑛) is affine if the one variable polynomial 𝑧_𝑗↦𝑃_𝑛(𝐳) can be written as 𝐴^(−𝑗)𝑧_𝑗+𝐵^(−𝑗) (for any 𝑗=1,…,𝑑), where 𝐴^(−𝑗) and 𝐵^(−𝑗) are polynomials with respect to the 𝑧_𝑖's with 𝑖≠𝑗.) These discrete distributions have received much interest in the literature (see for instance [2] and the references therein). Of course, the affine polynomial 𝑃_𝑛 has to satisfy appropriate conditions to ensure that 𝐺_𝐍(𝐳) is a PGF. These conditions include the trivial equality 𝑃_𝑛(1,…,1)=1. However, determining all pairs (𝑃_𝑛,𝜆) such that 𝐺_𝐍(𝐳) is a PGF is still an open problem (see [6], for discussions related to this problem). As explained in [6], the affine polynomial 𝑃_𝑛(𝐳) can be rewritten 𝑃_𝑛(𝐳)=𝐴_𝑛𝑎₁𝑧₁,…,𝑎_𝑛𝑧_𝑛𝐴_𝑛𝑎₁,…,𝑎_𝑛, (4) where 𝑎₁,…,𝑎_𝑛 are positive numbers and 𝐴_𝑛 is an affine polynomial such that 𝐴_𝑛(0,…,0)=1. The Taylor expansions of [𝐴_𝑛(𝐳)]^−𝜆 and [𝑃_𝑛(𝐳)]^−𝜆 in the neighborhood of (0,…,0) will be denoted as follows: 𝐴_𝑛(𝐳)^−𝜆=_{𝜶∈ℕ𝑛0}𝑐_𝜶𝜆,𝐴_𝑛𝐳^𝜶,𝑃_𝑛(𝐳)^−𝜆=_{𝜶∈ℕ𝑛0}𝑐_𝜶𝜆,𝑃_𝑛𝐳^𝜶, (5) where 𝜶=(𝛼₁,…,𝛼_𝑛) and 𝐳^𝜶=∏^𝑛_𝑖=1𝑧^𝛼𝑖_𝑖. Equations (3) and (4) clearly show that the masses of multivariate NMDs denoted as 𝑐_𝜶(𝜆,𝑃_𝑛) can be expressed as follows: 𝑐_𝜶𝜆,𝑃_𝑛=𝑐_𝜶𝜆,𝐴_𝑛∏^𝑛_𝑖=1𝑎^𝛼𝑖_𝑖𝐴_𝑛𝑎₁,…,𝑎_𝑛^−𝜆. (6)

3. Masses of Negative Multinomial Distributions

In this section, we derive new expressions for the coefficients 𝑐_𝜶(𝜆,𝐴_𝑛) 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 𝑐_𝜶(𝜆,𝐴_𝑛), we derive several results summarized in this section whereas all demonstrations are reported in the appendix.

Theorem 1. Denote 𝔓^∗_𝑛 as the set of nonempty subsets of [𝑛]={1,…,𝑛}. Any affine polynomial 𝐴_𝑛 such that 𝐴_𝑛(𝟎)=1 denoted as 𝐴_𝑛(𝐳)=1−_{𝑇∈𝔓∗𝑛}𝑎_𝑇𝐳^𝑇 (7) can be expressed as follows: 𝐴_𝑛(𝐳)=_{[𝑛]𝑖∈}1−𝑎_𝑖𝑧_𝑖−_{𝑇∈𝔓∗𝑛,|𝑇|⩾2}𝑑^𝑛_𝑇𝐳^𝑇_{[𝑛]𝑖∈⧵𝑇}1−𝑎_𝑖𝑧_𝑖, (8) where |𝑇| is the cardinal of the set 𝑇. Moreover 𝐴_𝑛(𝐳)=_{[𝑛]𝑖∈}1−𝑎_𝑖𝑧_𝑖×1−_{𝑇∈𝔓∗𝑛,|𝑇|⩾2}𝑑^𝑛_𝑇𝐳^𝑇∏_𝑖∈𝑇1−𝑎_𝑖𝑧_𝑖=_{[𝑛]𝑖∈}1−𝑎_𝑖𝑧_𝑖×1−𝑄_𝑛𝑧₁1−𝑎₁𝑧₁𝑧,…,^𝑛1−𝑎_𝑛𝑧_𝑛, (9)(10) where 𝑄_𝑛 is the polynomial defined by 𝑄_𝑛(𝐳)=∑_𝑇∈𝔓^∗_{𝑛,|𝑇|⩾2}𝑑^𝑛_𝑇𝐳^𝑇 and 𝑑^𝑛_𝑇 is related to the 2^|𝑇|−1 variables 𝑎_𝑆,𝑆∈𝔓^∗_𝑇 as follows: 𝑑^𝑛_𝑇=^|𝑇|_{𝑇∈𝔓𝑛|𝑇|>1}𝑎_𝑇𝑎^{[𝑛]⧵𝑇}+(|𝑇|−1)_𝑖∈𝑇𝑎_𝑖. (11)

Remark 2. In the trivariate case defined by 𝑛=3, the polynomial 𝐴₃(𝐳) can be expressed as 𝐴₃(𝐳)=1−𝑎₁𝑧₁1−𝑎₂𝑧₂1−𝑎₃𝑧₃×1−𝑄₃𝑧₁1−𝑧₁,𝑧₂1−𝑧₂,𝑧₃1−𝑧₃, (12) where the coefficients of the polynomial 𝑄₃(𝐳)=𝑏_1,2𝑧₁𝑧₂+𝑏_1,3𝑧₁𝑧₃+𝑏_2,3𝑧₂𝑧₃+𝑏_1,2,3𝑧₁𝑧₂𝑧₃ (13) can be determined using the relations 𝑏_𝑖,𝑗=𝑎_𝑖,𝑗+𝑎_𝑖𝑎_𝑗{}𝑏,𝑖,𝑗∈1,2,3,𝑖≠𝑗,_1,2,3=𝑎_1,2,3+𝑎₁𝑎_2,3+𝑎₂𝑎_1,3+𝑎₃𝑎_1,2+2𝑎₁𝑎₂𝑎₃. (14) The next theorem provides a relation between the coefficients of the polynomials 𝐴_𝑛(𝐳) and 𝑄_𝑛(𝐳) introduced above.

Theorem 3. Let 𝐴_𝑛(𝐳)=1−∑_𝑇∈𝔓^∗_𝑛𝑎_𝑇𝑧^𝑇,  𝐚=(𝑎₁,…,𝑎_𝑛), and 𝑄_𝑛  be the affine polynomial defined in (10) and (11). For any 𝛼 and 𝛾 in ℕ^𝑛, denote as 𝑐_𝛾(𝜆,𝐴_𝑛) the coefficient of 𝐳^𝛾 in the Taylor expansion of [𝐴_𝑛(𝐳)]^−𝜆 and as 𝑐_𝛼(𝜆,1−𝑄_𝑛) the coefficient of 𝐳^𝛼 in the Taylor expansion [1−𝑄_𝑛(𝐳)]^−𝜆. The following relation can be obtained:𝑐_𝛾𝜆,𝐴_𝑛=_{𝛼+𝛽=𝛾}𝑐_𝛼𝜆,1−𝑄_𝑛(𝜆𝟏+𝛼)_𝛽𝐚^𝛽𝛽!=_{0⩽𝛽𝑖⩽𝛾𝑖,𝑖=1,…,𝑛}𝑐_𝛾−𝛽𝜆,1−𝑄_𝑛×^𝑛_𝑖=1𝜆+𝛾_𝑖−𝛽_𝑖𝛽𝑖𝑎^𝛽𝑖_𝑖𝛽_𝑖!=_{0⩽𝛼𝑖⩽𝛾𝑖,𝑖=1,…,𝑛}𝑐_𝛼𝜆,1−𝑄_𝑛×^𝑛_𝑖=1𝜆+𝛼_𝑖𝛽𝑖𝑎^{𝛾𝑖−𝛼𝑖}_𝑖𝛾_𝑖−𝛼_𝑖!. (15)(16)(17)

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 variables 𝐳=(𝑧₁,𝑧₂) defined by 𝐴₂(𝐳)=1−_{𝑇∈𝔓∗2}𝑎_𝑇𝑧^𝑇=1−𝑎₁𝑧₁−𝑎₂𝑧₂−𝑎_1,2𝑧₁𝑧₂. (18) The coefficient of 𝐳^𝛾 in the Taylor expansion of [𝐴₂(𝐳)]^−𝜆, can be computed as follows 𝑐_𝛾𝜆,𝑃₂=(𝜆)_{max(𝛾1,𝛾2)}^{min(𝛾1,𝛾2)}_ℓ=0(𝜆+ℓ)_{min(𝛾1,𝛾2)−ℓ}𝛾₁−ℓ!𝛾₂−ℓ!ℓ!×𝑎^𝛾1₁^−ℓ𝑎^𝛾2₂^−ℓ𝑏^ℓ_1,2. (19)

Remark 2. The result (A.11) was mentioned in [7] without the factorization leading to (19). If 𝑎₁≠0 and 𝑎₂≠0, an equivalent formulation of (19) is 𝑐_𝛾𝜆,𝑃₂=𝑎^𝛾1₁𝑎^𝛾2₂𝛾₁!𝛾₂!(𝜆)_{max(𝛾1,𝛾2)}×^{min(𝛾1,𝛾2)}_ℓ=0(𝜆+ℓ)_{min(𝛾1,𝛾2)−ℓ}×𝛾₁ℓ𝛾₂ℓ𝑏ℓ!^1,2𝑎₁𝑎₂^ℓ. (20)

3.2. Trivariate NMDs

Theorem 6. Consider the affine polynomial with the three variables 𝐳=(𝑧₁,𝑧₂,𝑧₃) defined by 𝐴₃(𝐳)=1−_{𝑇∈𝔓∗3}𝑎_𝑇𝑧^𝑇. (21) The coefficient of 𝐳^𝛾 in the Taylor expansion of [𝐴₃(𝐳)]^−𝜆 are 𝑐_𝛾𝜆,𝑃₃=^𝛾1_𝛽1^𝛾₌₀²_𝛽2^𝛾₌₀³_𝛽3=0^{⌊|𝛾−𝛽|/2⌋}_{𝑣=𝛾−𝛽}(𝜆)_𝑣×𝑏^{𝑣−𝛾1+𝛽1}_2,3𝑏^{𝑣−𝛾2+𝛽2}_1,3𝑏^{𝑣−𝛾3+𝛽3}_1,2∏³_𝑖=1𝑣−𝛾_𝑖+𝛽_𝑖!𝑏^{|𝛾−𝛽|−2𝑣}_1,2,3!×𝛾−𝛽−2𝑣𝜆+𝛾₁−𝛽_1𝛽1𝛽₁!𝜆+𝛾₂−𝛽_2𝛽2𝛽₂!×𝜆+𝛾₃−𝛽_3𝛽3𝛽₃!𝑎^𝛽1₁𝑎^𝛽2₂𝑎^𝛽3₃. (22) When 𝑎_𝑖≠0, 𝑏_𝑖,𝑗≠0, 𝑖≠𝑗, 𝑖=1,2,3, 𝑗=1,2,3, and 𝑏_1,2,3≠0, an equivalent expression is 𝑐_𝛾𝜆,𝑃₃=𝑎^𝛾1₁𝑎^𝛾2₂𝑎^𝛾3₃^𝛾1_𝛽1^𝛾₌₀²_𝛽2^𝛾₌₀³_𝛽3=0^{⌊|𝛾−𝛽|/2⌋}_{𝑣=𝛾−𝛽}(𝜆)_𝑣×𝑏_2,3𝑏_1,3𝑏_1,2/𝑏²_1,2,3^𝑣𝛾−𝛽−2𝑣!∏³_𝑖=1𝑣−𝛾_𝑖+𝛽_𝑖!׳_𝑖=1𝜆+𝛾_𝑖−𝛽_𝑖𝛽𝑖𝛽_𝑖!𝑏_1,2,3𝑎₁𝑏_2,3^{𝛾1−𝛽1}×𝑏_1,2,3𝑎₂𝑏_1,3^{𝛾2−𝛽2}𝑏_1,2,3𝑎₃𝑏_1,2^{𝛾3−𝛽3}=(𝜆)_‖𝛾‖𝐚^𝛾𝛾!^𝛾1_𝛼1^𝛾₌₀²_𝛼2^𝛾₌₀³_𝛼3=0^{⌊|𝛼|/2⌋}_{‖𝛼‖𝑣=}(𝜆+‖𝛼‖)_{𝑣−‖𝛼‖}(|𝛼|−2𝑣)!∏³_𝑖=1𝑣−𝛼_𝑖!×𝑏_2,3𝑏_1,3𝑏_1,2𝑏²_1,2,3^𝑣𝑏𝛼!^1,2,3𝑎₁𝑏_2,3^𝛼1𝑏_1,2,3𝑎₂𝑏_1,3^𝛼2𝑏_1,2,3𝑎₃𝑏_1,2^𝛼3׳_𝑖=1𝜆+𝛼_{𝑖𝛾𝑖−𝛼𝑖}𝛾_𝑖𝛼_𝑖. (23)(24)

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 vector 𝐀=(𝐴_𝑋,𝐴_𝑌), called the Jones vector, whose covariance matrix, called the polarization matrix, is 𝔼𝐴Γ=^𝑋𝐴^∗_𝑋𝔼𝐴_𝑋𝐴^∗_𝑌𝔼𝐴_𝑌𝐴^∗_𝑋𝔼𝐴_𝑌𝐴^∗_𝑌≜𝑎₁𝑎₃+𝑖𝑎₄𝑎₃−𝑖𝑎₄𝑎₂, (25) where* denotes the complex conjugate. The covariance matrix Γ is a nonnegative Hermitian matrix whose diagonal terms are the intensity components in the 𝑋 and 𝑌 directions. The cross terms of Γ are the correlations between the Jones components. If we assume a fully developed speckle, the Jones vector 𝐀 is distributed according to a complex Gaussian distribution whose probability density function (pdf) is [15]: 𝑝(𝐀)=1𝜋²|Γ|exp−𝐀^†Γ⁻¹𝐀, (26) where |Γ| is the determinant of the matrix Γ and † denotes the conjugate transpose operator. As a consequence, the statistical properties of 𝐀 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): 𝐼₁=𝐴_𝑋²,𝐼₂=𝐴_𝑌²,𝐼₃=12𝐴_𝑋²+12𝐴_𝑌²+Re𝐴_𝑋𝐴^∗_𝑌,𝐼₄=12𝐴_𝑋²+12𝐴_𝑌²+Im𝐴_𝑋𝐴^∗_𝑌. (27) The state of the polarization of the light is classically characterized by the square DoP defined by [15, pages 134–136] 𝑃²|Γ|=1−4[trace(Γ)]²4𝑎=1−¹𝑎₂−𝑎²₃+𝑎²₄𝑎₁+𝑎₂², (28) where trace(Γ) is the trace of the matrix Γ. The light is totally depolarized for 𝑃=0, totally polarized for 𝑃=1, and partially polarized when 𝑃∈]0,1[. As a consequence, estimating the square DoP of a polarimetric image is important in many practical applications. Different estimation methods of 𝑃² using several combinations of intensity images were studied in [20]. Since only one realization of the random vector 𝐈=(𝐼₁,…,𝐼₄)^𝑇 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 vector 𝐍=(𝑁₁,…,𝑁₄) such that the conditional distributions of the random variables 𝑁_𝑙∣𝐼_𝑙, for 𝑙=1,…,4 are independent and distributed according to Poisson distributions with means 𝐼_𝑙, for 𝑙=1,…,4. The resulting joint distribution of 𝐍 is a multivariate mixed Poisson distribution [22]: 𝑃(𝐍=𝐤)=⋯_(ℝ+)4⁴_𝑙=1𝐼^𝑘𝑙_𝑙𝑘_𝑙!exp−𝐼_𝑙𝑓(𝐈)𝑑𝐈, (29) where 𝐤=(𝑘₁,…,𝑘₄), 𝑘_𝑖∈ℕ, and 𝑓(𝐈) is the joint pdf of the intensity vector. This section studies estimators of the square DoP 𝑃² defined in (28) based on several vectors 𝐍¹,…,𝐍^𝑛 belonging to the estimation window. These estimators are constructed from estimates of the covariance matrix elements 𝑎_𝑖, 𝑖=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 vector 𝐈 is known to be a multivariate gamma distribution whose Laplace transform is [20] 𝐸exp⁴_𝑘=1𝑧_𝑘𝐼_𝑘=1𝑃₄(𝐳), (30) where the affine polynomial 𝑃₄ is 𝑃₄(𝐳)=1+𝐳𝝁+𝑘_𝑎2𝑧₁𝑧₂+𝑧₃𝑧₄+𝑧₁+𝑧₂𝑧₃+𝑧₄ (31) with 𝐳=(𝑧₁,…,𝑧₄) and 𝑘_𝑎=12𝑎₁𝑎₂−𝑎²₃−𝑎²₄,𝝁=𝑎₁,𝑎₂,𝑎₃+𝑎₁+𝑎₂2,𝑎₄+𝑎₁+𝑎₂2^𝑇. (32) As a consequence, the distribution of 𝐍 is an NMD whose PGF can be written as (the interested reader is invited to consult [22] for more details)𝐺_𝐍(𝐳)=1𝑃₄𝑧₁−1,𝑧₂−1,𝑧₃−1,𝑧₄.−1 (33) The results of Section 2 allow us to compute the masses of 𝐍 that will be useful for studying the maximum likelihood estimator (MLE) of the square DoP.

4.2. MLE Using Three Polarimetric Images

The PGF of 𝐍=(𝑁₁,𝑁₂,𝑁₃) can be computed from (33) by setting 𝑧₄=1. The following result can be obtained: 𝐺_𝐍(𝐳)=1𝑃₃(𝐳) (34) with 𝑃₃(𝐳)=𝑃₃(𝟎)+𝑧₁𝜇₁−3𝑘_𝑎+𝑧₂𝜇₂−3𝑘_𝑎+𝑧₃𝜇₃−2𝑘_𝑎+𝑘_𝑎2𝑧₁𝑧₂+𝑧₁𝑧₃+𝑧₂𝑧₃, (35)𝐳=(𝑧₁,𝑧₂,𝑧₃), and 𝑃₃(𝟎)=1−∑³_𝑖=1𝜇_𝑖+4𝑘_𝑎. The results of Section 3.2 can then be used to express the masses of 𝐍 as a function of 𝜽=(𝑎₁,𝑎₂,𝑎₃,𝑎²₄)^𝑇.

The ML estimator of 𝜽 based on several vectors 𝐍^𝑘 belonging to the estimation window (where 𝑘=1,…,𝐾 and 𝐾 is the number of pixels of the observation window) is obtained by maximizing the log-likelihood 𝑙₃𝐍¹,…,𝐍^𝐾∣𝜽=^𝐾_𝑘=1log𝑃𝐍^𝑘 (36) with respect to 𝜽 (note again that 𝑃(𝐍^𝑘) 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 𝜽=(𝑎₁,𝑎₂,𝑎₃,𝑎²₄)^𝑇, are then plugged into (28) to provide the ML estimator of the square DoP based on three polarimetric images: 𝑃²4𝑎=1−¹𝑎₂−𝑎²₃+𝑎²₄𝑎₁+𝑎₂². (37)

4.3. MLE Using Two Polarimetric Images

The PGF of 𝐍=(𝑁₁,𝑁₂) can be computed from (34) by setting 𝑧₃=1. The following result can be obtained: 𝐺_𝐍(𝐳)=1𝑃₂(𝐳) (38) with 𝑃₂(𝐳)=𝑃₂(𝟎)+𝑧₁𝜇₁−2𝑘_𝑎+𝑧₂𝜇₂−2𝑘_𝑎+2𝑘_𝑎𝑧₁𝑧₂, (39)𝐳=(𝑧₁,𝑧₂) and 𝑃₂(𝟎)=1−∑²_𝑖=1𝜇_𝑖+2𝑘_𝑎. The results of Section 3.1 can then be used to express the masses of 𝐍 as a function of 𝜽=(𝑎₁,𝑎₂,𝑘_𝑎)^𝑇.

The MLE of 𝜽 based on several vectors 𝐍^𝑘 belonging to the estimation window is obtained by maximizing the log-likelihood 𝑙₂𝐍¹,…,𝐍^𝐾∣𝜽=^𝐾_𝑘=1log𝑃𝐍^𝑘 (40) with respect to 𝜽 (note that 𝑃(𝐍^𝑘) 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 𝜽=(𝑎₁,𝑎₂,𝑘_𝑎)^𝑇, are then plugged into (28) to provide the MLE of the square DoP based on two polarimetric images: 𝑃²=1−8𝑘_𝑎𝑎₁+𝑎₂². (41)

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 of 𝑎_𝑖, for 𝑖=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 Γ_𝑗, 𝑗=0,…,10, have been normalized so that the mean of the total intensity 𝐄[𝐼₁+𝐼₂]=𝐄[𝑁₁+𝑁₂]=𝑎₁+𝑎₂ 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.

Figure 1(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 𝑛=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 𝑃². The MLE is known to be asymptotically unbiased and efficient under mild regularity conditions (that are satisfied for 𝑃²). 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 Γ₀ and Γ₁₀. 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 𝑃² 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 of 𝑃² using 2 or 3 images, obtained from 1000 Monte-Carlo runs (𝑛=51×51).

(a) Low flux

(b) High flux

Figure 1 

log MSEs of the square DoP estimates using 2 and 3 images versus 𝑃2 for the set of polarization matrices defined in Table 1 under (a) low flux and (b) high flux assumptions (𝑛=51×51, ML: maximum likelihood estimators, and Asympt.: theoretical asymptotic value of the log MSE for a given estimator).

In 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 vector 𝐈 is assumed to be known. Thus, the high flux MLEs using two and three images can be derived from the intensity vectors 𝐈^𝑘=(𝐼^𝑘₁,𝐼^𝑘₂)^𝑇 and 𝐈^𝑘̃=(𝐼^𝑘₁,𝐼^𝑘₂,𝐼^𝑘₃)^𝑇 [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 𝑃²=1/3 and decreasing MSEs as 𝑃² 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 size 𝑛. 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 Γ₂ and Γ₈). 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, 𝑛>25×25. Moreover, the gain of a performance using 3 images instead of 2 is more significant for small values of 𝑃² (indeed, the difference between the different curves is more pronounced in the left figure corresponding to 𝑃²=0.2 than in the right figure corresponding to 𝑃²=0.8.)

Table 3 

Simulation results for the estimation of 𝑃² using 2 or 3 images, obtained from 1000 Monte-Carlo runs (matrix Γ₂).

Table 4 

Simulation results for the estimation of 𝑃² using 2 or 3 images, obtained from 1000 Monte-Carlo runs (matrix Γ₈).

(a) Matrix Γ2

(b) Matrix Γ8

Figure 2 

log MSE of the estimated square DoP 𝑃2 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).

Figures 3(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.

(a) Matrix Γ2

(b) Matrix Γ8

Figure 3 

log MSE of the estimated square DoP 𝑃2 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).

5.2. Application to Synthetic Polarimetric Images

In order to appreciate the estimation performance on polarimetric images, we consider a synthetic polarimetric image of size 512×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 𝐍=(𝑁₁,𝑁₂,𝑁₃,𝑁₄)^𝑇 for each pixel. The square DoP of each pixel 𝐱^(𝑖,𝑗) (for 𝑖,𝑗=1,…,512) has been estimated from vectors belonging to windows of size 𝑛=15×15 centered around the pixel of coordinates (𝑖,𝑗) 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 Figure 4.

(a) Scene

(b) Theoretical squared DoP

Figure 4 

Composition of the scene used to generate synthetic polarimetric low flux images and associated theoretical squared DoP.

(a) Low flux intensity 𝑁1

(b) Low flux intensity 𝑁2

(c) Low flux intensity 𝑁3

(d) Low flux intensity 𝑁4

Figure 5 

Synthetic intensity images (negative colors) for the scene depicted in Figure 4 and described in Table 5.

(a) MLE 2 images

(b) MLE 3 images

Figure 6 

Estimates of 𝑃2 using 2 or 3 low flux intensity images for the synthetic polarimetric images for an estimation window of size 𝑛=15×15. MLE: maximum likelihood estimator.

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 Figure 7). 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 𝑁₁, 𝑁₂, 𝑁₃, and 𝑁₄ 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 𝑁₁+𝑁₂ and can hardly be distinguished without a polarimetric processing.

(a) Low flux intensity 𝑁1

(b) Low flux intensity 𝑁2

(c) Low flux intensity 𝑁3

(d) Low flux intensity 𝑁4

Figure 7 

Real-world polarimetric intensity images of a scene composed of a plastic disk (left) and a steel disk (right).

Figure 8 shows the ML estimates of the square DoPs 𝑃² for 3 images and an estimation window of size 𝑛=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 of 𝑃2 using the 3 low flux intensity images 𝑁1, 𝑁2, and 𝑁3 for the real polarimetric images for an estimation window of size 𝑛=9×9. “MLE”: maximum likelihood estimator.

Appendix

Proofs of the Theorems

Proof of Theorem 1. The set of affine polynomials with real coefficients and variables (𝑧₁,…,𝑧_𝑛) is a vector space of dimension 2^𝑛 spanned by the basis (𝐳^𝑇)_{𝑇∈𝔓𝑛}. The set of polynomials (𝐳^𝑇∏_{𝑡∈[𝑛]⧵𝑇}(1−𝑎_𝑡𝑧_𝑡))_{𝑇∈𝔓𝑛} is another basis of this vector space. The proof of the theorem is obtained by expressing the coefficients of 𝐴_𝑛 in this latter basis. Considering the expansion 𝐴_𝑛𝑧₁,…,𝑧_𝑛=𝑑^𝑛_{∅[𝑛]𝑖∈}1−𝑎_𝑖𝑧_𝑖−_{𝑇∈𝔓∗𝑛}𝑑^𝑛_𝑇𝐳^𝑇_{[𝑛]𝑖∈⧵𝑇}1−𝑎_𝑖𝑧_𝑖, (A.1) the following results can be obtained.(1)Substituting 𝐳=𝟎 in (A.1) leads to 𝑑^𝑛_∅=1.(2)Substituting 𝐳_𝑖=(𝛿_𝑖(1),…,𝛿_𝑖(𝑛)) in (A.1), where 𝛿_𝑖(𝑗)=1 if 𝑖=𝑗 and 𝛿_𝑖(𝑗)=0 else leads to 1−𝑎_𝑖=1−𝑎_𝑖−𝑑^𝑛_{𝑖}, (A.2) or equivalently 𝑑^𝑛_{𝑖}=0,𝑖∈[𝑛]. (A.3)(3)Without loss of generality, if |𝑇|=𝑘>1, the coefficients 𝑑^𝑛_𝑇 can be computed for 𝑇={1,…,𝑘}. Indeed, consider a permutation 𝜎 such that 𝑇={𝜎(1),…,𝜎(𝑘)}. If 𝑑^𝑛_{1,…,𝑘}=𝑓({𝑎_𝑆}_𝑆∈𝔓^∗_𝑘), we have 𝑑^𝑛_𝑇=𝑓({𝑎_𝜎(𝑆)}_{𝑆∈𝔓∗𝑘}). Using the relation 1−𝐴_𝑛𝑧₁,…,𝑧_𝑛−1,0=_{𝑇∈𝔓∗𝑛−1}𝑎_𝑇𝑧^𝑇=1−𝐴_𝑛−1𝑧₁,…,𝑧_𝑛−1, (A.4) the expansion 𝐴_𝑛𝑧₁,…,𝑧_𝑛=_{[𝑛]𝑖∈}1−𝑎_𝑖𝑧_𝑖−_{𝑇∈𝔓∗𝑛,|𝑇|⩾2}𝑑^𝑛_𝑇𝐳^𝑇_{[𝑛]𝑖∈⧵𝑇}1−𝑎_𝑖𝑧_𝑖, (A.5) yields for any 𝑘<𝑛, 𝑑^𝑛_[𝑘]=𝑑_[𝑘]^𝑛−1. (A.6) In order to determine 𝑑^𝑛_[𝑛], we can assume 𝑎_𝑖≠0,  𝑖∈[𝑛] and substitute 𝑧_𝑖=1/𝑎_𝑖,  𝑖∈[𝑛] in (A.5). The following result is obtained 1−𝑛−_{𝑇∈𝔓∗𝑛,|𝑇|⩾2}𝑎_𝑇𝑎^𝑇𝑑=−^𝑛𝑇𝑎^[𝑛], (A.7) hence 𝑑^𝑛_[𝑛]=(𝑛−1)+_{𝑇∈𝔓∗𝑛,|𝑇|⩾2}𝑎_𝑇𝑎^𝑇𝑎^[𝑛]=_{𝑇∈𝔓∗𝑛,|𝑇|⩾2}𝑎_𝑇𝑎^{[𝑛]⧵𝑇}+(𝑛−1)𝑎^[𝑛]. (A.8)

Proof of Theorem 3. The relation (10) leads to 𝐴_𝑛(𝐳)^−𝜆=1−𝑄_𝑛𝑧₁1−𝑎₁𝑧₁𝑧,…,^𝑛1−𝑎_𝑛𝑧_𝑛^−𝜆×^𝑛_𝑖=11−𝑎_𝑖𝑧_𝑖^−𝜆=_{𝛼∈ℕ𝑛}𝑐_𝛼𝜆,1−𝑄_𝑛𝐳𝟏−𝐚𝐳^𝛼×^𝑛_𝑖=11−𝑎_𝑖𝑧_𝑖^−𝜆=_{𝛼∈ℕ𝑛}𝑐_𝛼𝜆,1−𝑄_𝑛𝐳^𝛼(𝟏−𝐚𝐳)^{−(𝜆𝟏+𝛼)}=_{𝛼∈ℕ𝑛}𝑐_𝛼𝜆,1−𝑄_𝑛𝐳^𝛼_{𝛽∈ℕ𝑛}(𝜆𝟏+𝛼)_𝛽𝐚^𝛽𝐳^𝛽𝛽!=_{𝛼∈ℕ𝑛𝛽∈ℕ𝑛}𝑐_𝛼𝜆,1−𝑄_𝑛(𝜆𝟏+𝛼)_𝛽𝐚^𝛽𝐳𝛽!^𝛼+𝛽=_{𝛾∈ℕ𝑛𝛼+𝛽=𝛾}𝑐_𝛼𝜆,1−𝑄_𝑛(𝜆𝟏+𝛼)_𝛽𝐚^𝛽𝐳𝛽!^𝛾 (A.9) which proves (15). Straightforward computations allow us to obtain the equalities (16) and (17) from (15).

Proof of Theorem 4. Denote ‖𝛼‖=max_{𝑖=1,…,𝑛}(𝛼_𝑖), |𝛼|=∑^𝑛_𝑖=1𝛼_𝑖 and introduce the notation of [6] 𝑐_𝛼𝜆,𝑃_𝑛=_{𝑘∈𝐾𝛼}(𝜆)_|𝑘|𝐚^𝑘,𝑘! (A.10) where (𝜆)_𝑘=𝜆(𝜆+1)⋯(𝜆+𝑘−1)=Γ(𝜆+𝑘)/Γ(𝜆). By using Theorem 3, the following results can be obtained 𝑐_𝛾𝜆,𝑃₂=_{𝛼+𝛽=𝛾}𝑐_𝛼𝜆,1−𝑄₂(𝜆𝟏+𝛼)_𝛽𝐚^𝛽=𝛽!_{min(𝛾1,𝛾2)ℓ=0}(𝜆)_ℓ(𝜆+ℓ)_𝛾1−ℓ(𝜆+ℓ)_𝛾2−ℓ×𝑎^𝛾1₁^−ℓ𝑎^𝛾2₂^−ℓ𝑏^ℓ_1,2𝛾₁−ℓ!𝛾₂=−ℓ!ℓ!_{min(𝛾1,𝛾2)ℓ=0}Γ(𝜆+ℓ)Γ(𝜆)Γ𝜆+𝛾₁Γ()Γ𝜆+ℓ𝜆+𝛾₂Γ()×𝑎𝜆+ℓ_𝛾11−ℓ𝑎^𝛾2₂^−ℓ𝑏^ℓ_1,2𝛾₁−ℓ!𝛾₂=Γ𝛾−ℓ!ℓ!𝜆+max₁,𝛾₂Γ(𝜆)^{min(𝛾1,𝛾2)}_ℓ=0Γ𝜆+min𝛾₁,𝛾₂Γ(𝜆)×Γ(𝜆)Γ()𝑏𝜆+ℓ𝑙!^ℓ1,22_𝑖=1𝑎^𝛾𝑖_𝑖^−ℓ𝛾_𝑖!=(𝜆)−ℓ_{max(𝛾1,𝛾2)min(𝛾1,𝛾2)ℓ=0}(𝜆+ℓ)_{min(𝛾1,𝛾2)−ℓ}𝛾₁−ℓ!𝛾₂−ℓ!ℓ!×𝑎^𝛾1₁^−ℓ𝑎^𝛾2₂^−ℓ𝑏^ℓ_1,2. (A.11)

Proof of Theorem 6. The relation (16) with 𝑄₃(𝐳)=𝑏_1,2𝑧₁𝑧₂+𝑏_1,3𝑧₁𝑧₃+𝑏_2,3𝑧₂𝑧₃+𝑏_1,2,3𝑧₁𝑧₂𝑧₃ leads to (22). By using the trivial equality 𝑏^{𝑣−𝛾1+𝛽1}_2,3𝑏^{𝑣−𝛾2+𝛽2}_1,3𝑏^{𝑣−𝛾3+𝛽3}_1,2𝑏^{|𝛾−𝛽|−2𝑣}_1,2,3𝑎^𝛽1₁𝑎^𝛽2₂𝑎^𝛽3₃=𝑏_2,3𝑏_1,3𝑏_1,2𝑏²_1,2,3^𝑣3_𝑖=1𝑎^𝛾𝑖_𝑖𝑏_1,2,3𝑎₁𝑏_2,3^{𝛾1−𝛽1}×𝑏_1,2,3𝑎₂𝑏_1,3^{𝛾2−𝛽2}𝑏_1,2,3𝑎₃𝑏_1,2^{𝛾3−𝛽3} (A.12) we easily obtain (23). Note that for ‖𝛼‖⩽𝑣⩽⌊|𝛼|/2⌋, we have (𝜆)_𝑣=(𝜆)_‖𝛼‖(𝜆+‖𝛼‖)_{𝑣−‖𝛼‖}. By substituting 𝛼_𝑖=𝛾_𝑖−𝛽_𝑖, 𝑖=1,2,3 in (23), the last result (24) can be obtained.

Acknowledgments

The authors would like to thank Gérard Letac for fruitful discussions regarding multivariate gamma distributions and negative multinomial distributions. The Authors are also grateful to Mehdi Alouini for acquiring and providing them the polarimetric images.