Theory of autocalibration feasibility and precision in full Stokes polarization imagers

Xiaobo Li; Xiaobo Li; Xiaobo Li; Benjamin Le Teurnier; Matthieu Boffety; Tiegen Liu; Tiegen Liu; Haofeng Hu; Haofeng Hu; Haofeng Hu; François Goudail; François Goudail

doi:10.1364/OE.390882

1. Introduction

Polarization imaging consists in measuring the state of polarization of the light backscattered by each point of a scene. It is a powerful tool with such applications as remote sensing, image restoration, material characterization [1–3]. The complete polarization state can be characterized by the Stokes vector [4]. Up to now, various setups for measuring an image of the Stokes vector have been proposed [5]. They all rely on acquiring a number $K$ of intensity measurements with different settings of the polarimeter and inverting these measurements to estimate the Stokes vector. Optimization of this sequence of measurements to minimize the estimation variance of the Stokes vector has long been investigated [6–8], and the minimal number of intensity measurements to estimate the four-dimensional Stokes vector is $K=4$. Thus, if more measurements are performed, they are redundant. This redundancy can be leveraged to improve the estimation precision [9,10] but also to jointly estimate the Stokes vector and one or several parameters of the polarimeter itself. Indeed, estimation of the Stokes vector relies on the precise knowledge of the optical parameters of the polarimeter, such as the directions of the polarizer and the retarder, the retardance of the retarder, etc. In practice, some of these parameters may vary after calibration due to environmental disturbance, such as temperature fluctuations [9,10]. Thus, they should be dynamically calibrated during the polarimetric measurements. Measurement redundancy can be leveraged for this purpose. The process of jointly estimating the Stokes vector and some parameters of the polarimeter will be called "autocalibration" in the following.

For example, Shibata et al. [9] showed that redundancy makes it possible to dynamically autocalibrate the retardance in imaging polarimeters based on a rotating retarder and a linear division-of-focal-plane (DoFP) camera. The validity domain and the precision of this method was characterized and optimized in Ref. [10], where it was shown that with this setup, autocalibration is possible only for some values of the input Stokes vector. The purpose of the present article is to build a general theory of the feasibility and of the precision of joint estimation of Stokes vector and polarimeter parameters. This theory is applicable to any polarization imaging setup defined by its measurement matrix. It will be illustrated on the example of retardance autocalibration in polarimetric imagers based on a rotating retarder, a rotating polarizer, and a standard or a DoFP camera [4,7,11–13]. These results are important to determine the best compromise between autocalibration capacity and polarimetric precision in practical applications.

The paper is organized as follows. In Section 2, we define the general theory of autocalibration based on the Cramer-Rao lower bound (CRLB). This theory makes it possible to determine the feasibility of autocalibration for any polarimeter architecture, any number of autocalibrated parameters and any value of the input Stokes vector. In Section 3, we specialize this theory to the autocalibration of a single parameter to get closed-form expressions of the CRLB that lend themselves to simpler physical interpretations. These results are applied in Sections 4 and 5 to retardance autocalibration in three of the most common setups of full Stokes imaging. We show that although these setups all yield optimal estimation precision of the Stokes vector, they do not have the same capacity and performance for autocalibration: They have to be specifically optimized to give the best possible autocalibration results. Finally, we summarize the obtained results, propose conclusion and draw perspectives in Section 6.

2. General theory for autocalibration

Let us denote $\textbf {S} = {({S_0},{S_1},{S_2},{S_3})^{T}}$ the four-dimension Stokes vector, where the superscript $T$ denotes transposition. $\textbf {S}$ can also be parameterized as [4]:

(1)$$\textbf{S} = {S_0}\left[ {\begin{array}{c} 1\\ {{\cal P}\cos 2\alpha \cos 2\varepsilon }\\ {{\cal P}\sin 2\alpha \cos 2\varepsilon }\\ {{\cal P}\sin 2\varepsilon } \end{array}} \right]$$

where ${S_0}$ represents the light intensity, ${\cal P} \in [0,1]$ the degree of polarization (DoP), $\alpha \in [ - 90^{\circ },90^{\circ }]$ the angle of polarization (AoP) and $\varepsilon \in [ - 45^{\circ },45^{\circ }]$ the ellipticity. Let us assume that $K$ intensity measurements are performed to estimate the Stokes vector in one pixel of the image. The $k^{th}$ intensity measurement can be written as [7]:

(2)$${i}_{k}=\frac{1}{2} \mathbf{T}_{k}^{T} \mathbf{S} + {n}_k, \quad k \in[1,K]$$

where ${n}_k$ denotes the measurement noise assumed additive, Gaussian, white, and of variance $\sigma ^{2}$, $\mathbf {T}_k$ denotes the $k^{th}$ measurement vector of the polarizer. Let us stack the measurements $i_k$ in a $K$-component vector $\mathbf {I}$ and define the matrix $\mathbb {W}$ which $k^{th}$ row is equal to the vector $\mathbf {T}_k^{T}$. With this notation, Eq. (2) can be written in matrix form as:

(3)$$\mathbf{I}= \mathbb{W}(\boldsymbol \eta)\mathbf{S} + \mathbf{N}.$$

where the symbol $\boldsymbol \eta = [\eta _1, \eta _2,\ldots , \eta _M]$ denotes a M-dimensional set of parameters on which the matrix $\mathbb {W}$ depends. These parameters $\eta _i$ may be the directions of the optical elements or the retardance of the retarder, etc. $\mathbf {N}$ is a white Gaussian noise vector. If $\boldsymbol \eta$ is known, the least mean square estimate of the input Stokes vector is $\hat {\mathbf {S}}=\mathbb {W}^{+} \mathbf {I}$, where $\mathbb {W}^{+}=\left (\mathbb {W}^{T} \mathbb {W}\right )^{-1} \mathbb {W}^{T}$ denotes the Moore Penrose pseudo-inverse matrix [14]. However, in practice, $\boldsymbol \eta$ may be unknown and need to be dynamically calibrated. The purpose for this paper is to address this issue by handling the simultaneous estimation of the Stokes vector $\mathbf {S}$ and of the set of parameters $\boldsymbol \eta$ as a joint estimation problem in the presence of noise.

We assume that the polarimetric image is formed of a number $P$ of pixels indexed by $p \in [1,P]$. This set can represent the pixels over the whole sensor, or over a restricted area. Each pixel $p$ may observe a different Stokes vector $\mathbf {S}_p$, but we assume that the unknown parameters $\boldsymbol \eta$ are common to all the pixels. Similarly, in most polarization imager architectures, the measurement matrix is the same for all pixels, and is thus independent of $p$. However, this may not always be the case. For example, in architectures based on DoFP cameras, the spatial variations of the micro-polarizer characteristics over the field due to manufacturing defects make the measurement matrix dependent of $p$ [15,16]. Thus our model takes into account this variation, by denoting $\mathbb {W}_p(\boldsymbol \eta )$ the measurement matrix at the $p^{th}$ pixel.

Since the noise model is Gaussian, white, and additive of variance $\sigma ^{2}$, the loglikelihood of the measurements has the following expression:

(4)$$\ell (\mathcal{S},\boldsymbol \eta) ={-} \frac{1}{2\sigma^{2}} \sum_{p=1}^{P} \| \textbf{I}_p - \mathbb{W}_p (\boldsymbol \eta) \textbf{S}_p \|^{2}$$

In this equation, the symbol $\mathcal {S} = \{ \textbf {S}_p \}, p\in [1,P]$ represents the set of the observed Stokes vectors, and $\textbf {I}_p$ represents the $K$-dimensional intensity measurement vector at the $p^{th}$ pixel.

Our goal is to jointly estimate the Stokes vectors $\textbf {S}_p$ in each pixel and the set of parameters $\boldsymbol \eta$ thanks to measurement redundancy. Let us compute the CRLB of this joint estimation problem. Using the expression of the loglikelihood in Eq. (4), it is easily shown that the Fisher matrix has the following expression:

(5)$$\mathbb{F} = \left[ \begin{array}{cc} \mathbb{A} & \mathbb{B}^{T} \\ \mathbb{B} & \mathbb{C} \end{array} \right]$$

where $\mathbb {A}$ is a matrix of dimension $M\times M$ with

(6)$$\mathbb{A}_{ij} ={-} \left\langle \frac{\partial ^{2} \ell}{\partial \eta_i \partial \eta _j} \right\rangle = \frac{1}{\sigma^{2}} \left( \sum_{p=1}^{P} \textbf{S}_{p}^{T} \frac{\partial \mathbb{W}_p^{T} }{\partial \eta _i} \frac{\partial \mathbb{W}_p}{\partial \eta _j} \textbf{S}_{p} \right)$$

where $\left \langle \cdot \right \rangle$ denotes ensemble averaging over the noise realizations. Please note that in the subsequent developments, we will denote the matrices $\mathbb {W}_p (\boldsymbol \eta )$ simply as $\mathbb {W}_p$, skipping the mention of the dependence on $\boldsymbol \eta$ for the sake of readability. The matrix $\mathbb {B}$ is a concatenation of $P$ matrices $\mathbb {B}_p$ of dimension $4\times M$ such that

(7)$$\mathbb{B}^{T} = \left[ \begin{array}{cccc} \mathbb{B}_1^{T} & \mathbb{B}_2^{T} & \ldots & \mathbb{B}_P^{T} \end{array} \right]^{T}$$

and each column of $\mathbb {B}_p$ is equal to:

(8)$$[\mathbb{B}_p]_{{\bullet} j} ={-} \left\langle \frac{\partial ^{2} \ell}{\partial \eta_j \partial \textbf{S}_p} \right\rangle = \frac{1}{\sigma^{2}} \mathbb{W}^{T}_p \frac{\partial \mathbb{W}_p}{\partial \eta _j} \textbf{S}_{p}$$

The matrix $\mathbb {C}$ is a block diagonal matrix composed of $P$ blocks $\mathbb {C}_p$ of size $4\times 4$ such that

(9)$$\mathbb{C}_p ={-} \left\langle \frac{\partial ^{2} \ell}{ \partial \textbf{S}_p^{2}} \right\rangle = \frac{1}{\sigma^{2}} \mathbb{W}_p^{T} \mathbb{W}_p$$

The expression of the inverse Fisher matrix is [17]:

(10)$$\mathbb{J}=\mathbb{F}^{{-}1}=\left[\begin{array}{cc} \mathbb{M}^{{-}1} & -\mathbb{M}^{{-}1} \mathbb{B}^{T} \mathbb{C}^{{-}1} \\ -\mathbb{C}^{{-}1} \mathbb{B} \mathbb{M}^{{-}1} & \mathbb{C}^{{-}1} + \mathbb{C}^{{-}1} \mathbb{B} {\mathbb{M}^{{-}1}} \mathbb{B}^{T} \mathbb{C}^{{-}1} \end{array}\right]$$

where $\mathbb {M}=\mathbb {A}-\mathbb {B}^{T} \mathbb {C}^{-1}\mathbb {B}$. Since the matrix $\mathbb {C}$ is block diagonal with blocks $\mathbb {C}_p$, the matrix $\mathbb {C}^{-1}$ is also block diagonal with blocks equal to $\mathbb {C}_p^{-1}$. Hence the matrix $\mathbb {M}$ has the following expression:

(11)$$\mathbb{M} = \mathbb{A} - \sum_{p=1}^{P} \mathbb{B}_p^{T} \mathbb{C}_p^{{-}1} \mathbb{B}_p$$

According to Eq. (10), the $i^{th}$ diagonal element of the matrix $\mathbb {M}^{-1}$, namely, $\mathbb {M}^{-1}_{ii}$, is the CRLB of the parameter $\eta _i$. The CRLB of a parameter represents a lower bound on its estimation variance with an unbiased estimator [18,19]. It is thus an indicator of the intrinsic difficulty of estimating this parameter. It is seen in Eq. (11) that the CRLB depends on the expression of the matrices $\mathbb {W}_p$ and on the observed Stokes vectors $\textbf {S}_p$. In particular, it is possible to estimate the parameter set $\boldsymbol \eta$ - i.e., to have a finite CRLB - only if the matrix $\mathbb {M}$ is non singular. The covariance matrix of estimation of the Stokes vector $\mathbf {S}_p$ is the $p^{th}$ diagonal block of the matrix in the bottom right side of Eq. (10). Its expression is

(12)$${\Gamma}_{_{\mathbf{S}_p}} = \mathbb{C}_p^{{-}1} + \mathbb{C}_p^{{-}1} \mathbb{B}_p {\mathbb{M}^{{-}1}} ( \mathbb{C}_p^{{-}1} \mathbb{B}_p)^{T}$$

where we have used the fact that the matrices $\mathbb {C}_p$ are symmetric. The CRLB of the Stokes parameter $[{S}_p]_i$, $i \in [0,3]$ is the $i^{th}$ diagonal element of this matrix.

The general expressions of the CRLB of $\boldsymbol \eta$ and $\mathbf {S}_p$ can be deduced from Eqs. 11 and (12), for any imaging polarimeter architecture defined by the measurement matrices $\mathbb {W}_p$ and any input Stokes vector. In the next section, we specialize these general expressions to a particular case of practical interest, in order to study more easily their physical meaning and their potential applications.

3. Autocalibration of a single parameter

Let us now assume that a single parameter is to be estimated, i.e., the parameter vector $\boldsymbol \eta$ is a scalar $\eta$. In this case, the matrix $\mathbb {M}^{-1}$ is also a scalar which represents the CRLB of the unique parameter $\eta$. According to Eqs. (6), (7), and (10), the matrix $\mathbb {A}$ is a scalar $a$ and the matrix $\mathbb {B}_p$ a 4 dimensional vector:

(13)$$a = \mathbb{A} = \frac{1}{\sigma^{2}}\sum_{p=1}^{P} \left\|\frac{\partial \mathbb{W}_p}{\partial \eta} \mathbf{S}_p\right\|^{2} \quad \hbox{and} \quad \mathbb{B}_p = \frac{1}{\sigma^{2}} \mathbb{W}_p^{T} \frac{\partial \mathbb{W}_p}{\partial \eta} \mathbf{S}_p$$

In consequence, using Eq. (9), the second additive term of $\mathbb {M}$ in Eq. (11) can be expressed as:

(14)$$x = \sum_{p=1}^{P} \mathbb{B}_p^{T} \mathbb{C}_p^{{-}1}\mathbb{B}_p = \frac{1}{\sigma^{2}} \sum_{p=1}^{P} \left\|\mathbb{P}_{W_p} \frac{\partial \mathbb{W}_p}{\partial \eta} \mathbf{S}_p\right\|^{2}$$

where $\mathbb {P}_{W_p}=\mathbb {W}_p \left (\mathbb {W}_p^{T} \mathbb {W}_p\right )^{-1} \mathbb {W}_p^{T}$ is the projector on the subspace spanned by the rows (or columns) of matrix $\mathbb {W}_p$. It can also be noted that:

(15)$$a-x = {\frac{1}{\sigma^{2}}} \sum_{p=1}^{P} \left( \left\| \frac{\partial \mathbb{W}_p}{\partial \eta} \textbf{S}_p \right \|^{2} - \left\| \mathbb{P}_{W_p} \frac{\partial \mathbb{W}_p}{\partial \eta} \textbf{S}_p \right \|^{2} \right)= {\frac{1}{\sigma^{2}}} \sum_{p=1}^{P} \left\| \mathbb{P}^{{\perp}}_{W_p} \frac{\partial \mathbb{W}_p}{\partial \eta} \textbf{S}_p \right\|^{2}$$

where $\mathbb {P}^{\perp }_{W_p}$ is the projector orthogonal the subspace spanned by the rows (or columns) of $\mathbb {W}_p$, that is:

(16)$$\mathbb{P}_{W_p}^{{\perp}} = \mathbb{I} - \mathbb{P}_{W_p} = \mathbb{I} - \mathbb{W}_p \left(\mathbb{W}_p^{T} \mathbb{W}_p\right)^{{-}1} \mathbb{W}_p^{T}$$

where $\mathbb {I}$ denotes the $K\times K$ identity matrix. Consequently, from Eqs. (11) and (15), one obtains the CRLB of the parameter $\eta$:

(17)$$\hbox{CRLB}[{\eta}] = \mathbb{M}^{{-}1} = \frac{{1}}{a-x} = \frac{\sigma^{2}}{\sum_{p=1}^{P} \left\| \mathbb{P}^{{\perp}}_{W_p} \frac{\partial \mathbb{W}_p}{\partial \eta} \textbf{S}_p \right \|^{2}}$$

The CRLB of the Stokes vector element $\{S_p\}_i~,~i\in [0,3]$, denoted CRLB$[\{S_{p}\}_i]$ in the following, is given by the diagonal elements of the matrix defined in Eq. (12). It is easily shown to be equal to:

(18)$$\hbox{CRLB}[\{S_{p}\}_i]_{\eta {\hbox{uk}}}~~=~~ \hbox{CRLB}[\{S_{p}\}_i]_{\eta \hbox{k}} ~+~ \frac{\left[ \mathbb{W}_p^{+} \frac{\partial \mathbb{W}_p}{\partial \eta} \textbf{S}{_p} \right]_{i}^{2}}{\sum_{p=1}^{P} \left\| \mathbb{P}^{{\perp}}_{W_p} \frac{\partial \mathbb{W}_p}{\partial \eta} \textbf{S}_p \right\|^{2}} \sigma^{2}$$

where the subscripts $\eta$k stands for "$\eta$ known" and $\eta$uk for "$\eta$ unknown". This expression is the sum of two terms. The first one:

(19)$$\hbox{CRLB}[\{S_{p}\}_i]_{\eta \hbox{k}} = \sigma^{2} \left[ (\mathbb{W}_p^{T} \mathbb{W}_p) ^{{-}1} \right]_{ii}$$

is the CRLB of $\{S_{p}\}_i$ when $\eta$ is known. It has to be noticed that this CRLB is exactly equal to the actual variance since we are in the presence of additive white Gaussian noise. The second term, which is non-negative, represents the increase of CRLB incurred by the fact that $\eta$ is unknown and estimated.

Let us now further simplify the problem by assuming that the measurement matrices at all pixels are equal, that is, $\forall p, \mathbb {W}_p = \mathbb {W}$, which is the case for most imaging polarimeters. Moreover, in order to facilitate the interpretation of the results, we also assume that the observed Stokes vectors are identical for all pixels, that is, $\forall p \in [1,P], \mathbf {S}_p = \mathbf {S}$. In this case, Eq. (17) simplifies to:

(20)$$\hbox{CRLB}[{\eta}] = \frac{\sigma^{2}}{P} \, \frac{1}{\left\| \mathbb{P}^{{\perp}}_W \frac{\partial \mathbb{W}}{\partial \eta} \textbf{S} \right \|^{2}}$$

As expected, the CRLB of $\eta$ is inversely proportional to the number $P$ of pixels used for calibration. Moreover, the condition for estimation of the parameter $\eta$ to be possible is clear: the input Stokes vector $\textbf {S}$ must not belong to the null space of the matrix $\mathbb {P}^{\perp }_W \frac {\partial \mathbb {W}}{\partial \eta }$. The capacity to autocalibrate $\eta$ thus depends both on the matrix $\mathbb {W}$ (and on the way it depends on $\eta$) and on the input Stokes vector. Similarly, from Eq. (18), the CRLBs of the Stokes parameters are equal to:

(21)$$\hbox{CRLB}[S_i ]_{\eta \hbox{uk}}~~=~~ \hbox{CRLB}[{S_i}]_{\eta \hbox{k}} ~ +~ \frac{\sigma^{2}}{P}\, \frac{\left[ \mathbb{W}^{+} \frac{\partial \mathbb{W}}{\partial \eta} \textbf{S} \right]_{i}^{2}}{\left\| \mathbb{P}^{{\perp}}_W \frac{\partial \mathbb{W}}{\partial \eta} \textbf{S} \right\|^{2}}$$

It is seen that the second additive term, which represents the increase of CRLB incurred by the fact of estimating $\eta$, is null if and only if the coefficient $\left [\mathbb {W}^{+} \frac {\partial \mathbb {W}}{\partial \eta } \mathbf {S}\right ]_{i}$ is equal to zero. In this case, estimating $\eta$ has no influence on the CRLB of $\mathbf {S}$. It is also seen that this term is multiplied by $1/P$: this means that the dependence of $\hbox {CRLB}[S_i ]_{\eta \hbox {uk}}$ on estimation of $\eta$ is reduced as the number of pixels used for this estimation increases.

In the next section, we illustrate the practical use of these results to characterize the autocalibration performance of a classical family of imaging polarimeters.

4. Autocalibration in the case of single-retarder-based polarimeters

In order to illustrate the previous results, let us now consider a class of polarimetric imagers composed of a polarizer and a retarder that may have variable orientations. The analysis vectors have the following expression:

(22)$$\mathbf{T}_{k}= [1, {a}_k \cos \delta + {b}_k, {c}_k \cos \delta \ + {d}_k, {e}_k \sin \delta]^{T},$$

where $\delta$ denotes the retardance of the retarder, and

(23)$$\begin{aligned} &{a}_k=\sin 2\theta_k \left[\sin 2 \left(\theta_k - \varphi_k\right)\right]; {b}_k=\cos 2\theta_k \left[\cos 2\left(\theta_k - \varphi_k\right)\right]; {c}_k = \cos 2\theta_k \left[\sin 2\left( \varphi_k-\theta_k \right)\right]\\ &{d}_k = \sin 2\theta_k \left[\cos 2\left( \varphi_k-\theta_k \right)\right]; {e}_k = \left[\sin 2\left(\varphi_{k}-\theta_{k} \right)\right] \end{aligned}$$

where $\boldsymbol \theta =(\theta _1,\ldots ,\theta _K)$ and $\boldsymbol \varphi =(\varphi _1,\ldots ,\varphi _K)$ denote the sets of $K$ directions of the retarder and of the polarizer, respectively [4,14].

In practice, the directions of the polarizer and of the retarder can be reasonably assumed stable during the measurement process, and not to depend on the environmental constraints such as the temperature. In contrast, the retardance $\delta$ of the retarder may be more sensitive to the environment and it makes sense to autocalibrate it [9]. Therefore, we consider in the following that $\delta$ is the only parameter to autocalibrate. Moreover, we assume that the measurement matrices and the observed Stokes vectors are common to all the pixels, so that Eqs. 20 and 21 can be used to compute the CRLB of $\delta$ and $\mathbf {S}$. The measurement matrix $\mathbb {W}$ can be computed by stacking the measurement vectors $\mathbf {T}_k$, and its derivative with respect to $\delta$ by using Eqs. (22) and 23. Moreover, it is shown in Appendix A. that for this type of setups, the first and fourth columns of the matrix $\mathbb {P}^{\perp }_W \frac {\partial \mathbb {W}}{\partial \eta }$ are zero. Using this result and Eq. (20), one obtains:

(24)$$\begin{aligned} \hbox{CRLB} [{\delta}] &=& \frac{1}{P \, \hbox{SNR}_\delta ^{2} }\frac{1}{\left\|\mathbb{Q}\mathbf{\bar{s}}\right\|^{2}} \quad \hbox{with} \quad \hbox{SNR}_\delta = \frac{S_0}{\sigma} \hbox{DoLP} \quad \hbox{and} \quad \hbox{DoLP} = \frac{\sqrt{S_1^{2}+S_2^{2}}}{S_0} \end{aligned}$$

where the "structure matrix" $\mathbb {Q}$ is a $4\times 2$ matrix composed of the second and third columns of the matrix $\mathbb {P}^{\perp }_W \frac {\partial \mathbb {W}}{\partial \eta }$, DoLP is the degree of linear polarization and

(25)$$\mathbf{\bar{s}} = \frac{(S_1,S_2)^{T}}{\sqrt{S_1^{2}+S_2^{2}}} = (\cos 2 \alpha,\sin 2 \alpha)^{T}$$

is the normalized linear Stokes vector. In Eq. (25), $\alpha$ denotes the AoP of the Stokes vector (see Eq. (1)).

From Eq. (24), it is immediately seen that autocalibration of $\delta$ is not possible for all input Stokes vectors. Indeed, the CRLB is infinite if the DoLP is zero. In consequence, autocalibration is impossible if the input Stokes vector is circular or totally depolarized. Moreover, autocalibration is all the more precise that the DoLP is large, that is, the input Stokes vector is polarized and close to linear. These properties had already been observed in Ref. [10] in the case of a DoFP-based full Stokes polarization imager. We show here that they are common to all the polarimetric imagers involving rotating polarizers and linear retarders. More generally, according to Eq. (24), it is obvious that the set of input polarization states for which autocalibration of the retardance is possible depends on the rank of the matrix $\mathbb {Q}$.

• If $\hbox {rank}(\mathbb {Q})=0$, $\mathbb {Q}$ is a zero matrix and $\left \|\mathbb {Q} \mathbf {\bar {s}} \right \| = 0$. Thus, autocalibration is impossible for any input Stokes vector.
• If $\hbox {rank}(\mathbb {Q})=1$, the null space of $\mathbb {Q}$ is of dimension 1. If $\mathbf {\bar {s}}$ belongs to this null space, CRLB$[\delta ]$ is infinite and autocalibration is impossible. In the Poincaré sphere, this null space is a disc passing by the center of the sphere and the poles (Fig. 1(a)). The orientation of this disc corresponds to a given value of the AoP. For all input Stokes vectors having this AoP, autocalibration is impossible.
• If $\hbox {rank}(\mathbb {Q})=2$, the null space is empty. Consequently, autocalibration is possible for all values of the AoP. It is impossible only for Stokes vectors which principal state of polarization is circular. These invalid input Stokes vectors form a line joining the two poles of the Poincaré Sphere (Fig. 1(b)). It is the largest possible validity domain for this family of polarization imagers.

In the next section, we consider different implementations of this family of polarization imagers and investigate, in each case, their validity domain for autocalibration by analyzing their structure matrix $\mathbb {Q}$.

Fig. 1. Representation on the Poincaré sphere of the set of input Stokes vectors for which retardance autocalibration is impossible when (a) $\operatorname {rank}\left (\mathbb {Q}\right )=1$; (b) $\operatorname {rank}\left (\mathbb {Q}\right )=2$. The red solid disk or line represent these sets.

Download Full Size | PDF

5. Autocalibration for different setups of Stokes imagers

In the following, we compare the autocalibration capacity and performance of the three most common setups of full Stokes polarimeters. They are represented in Fig. 2 and defined as follows:

• The setup RRFP includes a rotating retarder and a fixed polarizer, shown in Fig. 2(a).
• The setup RRRP includes a rotating retarder and a rotating polarizer, shown in Fig. 2(b).
• The setup RRDoFP includes a rotating retarder placed in front of a linear DoFP polarization camera having a grid of micro-polarizers in front of the sensor [15,16,20,21], shown in Fig. 2(c).

Fig. 2. Three most common setups of full Stokes polarimeter: (a) RRFP, (b) RRRP, (c) RRDoFP. R: retarder, P: polarizer, $S_{in}$ denotes the input Stokes vector.

Download Full Size | PDF

In the following, we investigate the validity domain for retardance autocalibration of these three setups by analyzing their structure matrix $\mathbb {Q}$.

5.1 Autocalibration of the RRFP setup

The RRFP setup has been used as a component of numerous polarimeters [11,12]. We assume that the direction of the fixed polarizer is set to a given angle $\varphi$ and the direction of the retarder in the $k^{th}$ intensity measurement is $\theta _k$. At least $K=4$ measurements are necessary to estimate the Stokes vector, but since redundancy is necessary for autocalibration, $K$ has to be greater than or equal to 5.

We show in Appendix B. that for all the possible measurement strategies, that is, whatever the number $K$ of measurements, the polarizer angle $\varphi$ and the set of retarder angles $\boldsymbol \theta$, one has $\operatorname {rank}(\mathbb {Q})=0$. Hence with this setup, autocalibration is impossible for any input Stokes vector. In other words, the type of redundancy introduced by using more than 4 measurements in this setup is not adapted to estimation of the retardance.

5.2 Autocalibration of the RRRP setup

This is the most common setup for Stokes polarimeters [4,7]. Since the polarizer can be rotated, it has a larger number of degrees of freedom than the RRFP setup. Let us first assume that $K=5$, which corresponds to the least level of redundancy. We first consider a simple example of setting in which the nominal value of the retardance is $\delta = 90^{\circ }$ and the rotation angles of the retarder and the polarizer are:

(26)$$\boldsymbol \theta = (0^{{\circ}},0^{{\circ}},0^{{\circ}},90^{{\circ}},45^{{\circ}}) ~~~~;~~~~ \boldsymbol \varphi = (0^{{\circ}},90^{{\circ}},45^{{\circ}},45^{{\circ}},45^{{\circ}}).$$

It is easily shown that the first column of the corresponding matrix $\mathbb {Q}$ is null, hence $\operatorname {rank}(\mathbb {Q})=1$. This means that when the input Stokes vector satisfies $s_2 = 0$ and thus AoP$= 90^{\circ }$, the denominator $\left \|\mathbb {Q} \mathbf {\bar {s}}\right \|^{2}$ in Eq. (24) is equal to zero: Autocalibration for $\delta$ is thus impossible. We have also checked that the result is the same for all the possible nominal values of the retardance.

In practice, the setting above is not optimal for Stokes vector estimation since it does not minimize the Equally Weighted Variance (EWV) defined as EWV $= \operatorname {trace} \left [ \left ( \mathbb {W}^{T} \mathbb {W} \right )^{ - 1} \right ]$ where trace $[\cdot ]$ denotes the trace of a matrix [8,22]. One of the optimal configurations that minimizes the EWV is $\delta =90^{\circ }$ and

(27)$$\boldsymbol \theta = (100.2^{{\circ}},68.6^{{\circ}},140.2^{{\circ}},19.1^{{\circ}},163.5^{{\circ}}) ~~~~;~~~~ \boldsymbol \varphi = (92.5^{{\circ}},130.3^{{\circ}},70.4^{{\circ}},121.5^{{\circ}},54.8^{{\circ}}).$$

The corresponding matrix $\mathbb {Q}$ is also of rank 1. The null space consists of the Stokes vectors with $\mathrm {AoP}\approx 31.8^{\circ }$. We have checked numerous settings by numerical simulation, all the results showed that the rank of the matrix $\mathbb {Q}$ is always equal to 1. Therefore, there exists no strategy based on $K=5$ intensity measurements can lead to $\operatorname {rank}(\mathbb {Q})=2$.

The only way to expand the validity domain of autocalibration is thus to increase the number of intensity measurements. Let us first consider the optimal setting based on $K=6$ measurements with $\delta =90^{\circ }$ and

(28)$$\boldsymbol \theta = (0^{{\circ}},0^{{\circ}},0^{{\circ}},45^{{\circ}},45^{{\circ}},45^{{\circ}}) ~~~~;~~~~ \boldsymbol \varphi = (0^{{\circ}},90^{{\circ}},45^{{\circ}},0^{{\circ}},45^{{\circ}},135^{{\circ}})$$

This setting reaches the minimal possible value EWV$=20/3$ [14,23]. Its measurement frame defines a regular octahedron (Platonic polyhedron) inscribed in the Poincaré Sphere (shown in Fig. 3(a)), which is an example of spherical 3-design [24,25]. However, the rank of the corresponding structure matrix is $\operatorname {rank}(\mathbb {Q})=1$. Its null space is spanned by the Stokes vectors of AoP $= 45^{\circ }$.

Fig. 3. Measurement frames defining polyhedra inscribed on the Poincaré Sphere and corresponding to the optimal settings in (a) Eq. (28): Platonic (octahedron), and (b) Eq. (32): Non-Platonic.

Download Full Size | PDF

An interesting question is whether there exists any configuration with 6 intensity measurements and $\operatorname {rank}(\mathbb {Q})=2$. The answer is yes, there exists an infinity of such configurations. However, each of them has a different matrix $\mathbb {Q}$ which corresponds to a different value of CRLB$[\delta ]$. Thus, the question is how to find the best one? A natural optimization criterion is to determine the configuration that minimizes the largest value of CRLB$[\delta ]$ as $\mathbf {\bar {s}}$ varies, for a given value $\delta$ of the nominal retardance. However, a configuration optimizing this criterion could lead to poor values of the EWV. We thus chose to determine the configurations that optimize this criterion while reaching the optimal value of EWV$=20/3$. In other words, we solve following constrained minmax problem:

(29)$$(\boldsymbol \theta;\boldsymbol \varphi ) = \arg \mathop {\min }_{\boldsymbol \theta;\boldsymbol \varphi} \left\{\max _{\|\bar{s}\|=1} \frac{\sigma^{2}}{\left\|\mathbb{Q}\mathbf{\bar{s}}\right\|^{2}}\right\} ~~~~ {s.t ~~ \hbox{EWV} = 20/3}$$

One has [26]:

(30)$$\max _{\|\bar{s}\|=1} \frac{1}{\left\|\mathbb{Q}\mathbf{\bar{s}}\right\|^{2}} = \frac{1}{\min\limits _{\|\bar{s}\|=1} \mathbf{\bar{s}}^{T} \mathbb{Q}^{T} \mathbb{Q} \mathbf{\bar{s}}} = \frac{1}{\lambda_{min}^{2}}$$

where $\lambda _{\min }$ is the minimal singular value of $\mathbb {Q}$. Consequently, the minmax optimization problem expressed in Eq. (29) is equivalent to the following constrained optimization problem:

(31)$$(\boldsymbol \theta;\boldsymbol \varphi ) = \arg \mathop {\max }_{\boldsymbol \theta;\boldsymbol \varphi} \left\{ {{\lambda_{\min}}} \right\} ~~~~ {s.t} ~~ {\hbox{EWV} = 20/3}$$

This constrained optimization problem is numerically solved by employing sequential quadratic programming (SQP) algorithm. It turns large nonlinear problems into a sequence of small quadratic problems to reduce the computational load, and thus can handle any degree of non-linearity including non-linearity in the constraints [27]. It has been implemented with the nonlinear programming solver "fmincon" in MATLAB. When $\delta =90^{\circ }$, one of the optimal settings found with this method is:

(32)$$\begin{aligned} \boldsymbol \theta & = (57.0^{ }\circ,42.6^{{\circ}},177.0^{ }\circ,102.6^{{\circ}},117.0^{{\circ}},162.6^{{\circ}});\\ \boldsymbol \varphi & = (129.4^{{\circ}},150.3^{{\circ}},69.4^{{\circ}},30.3^{{\circ}},9.4^{{\circ}},90.3^{{\circ}}). \end{aligned}$$

It correspond to a maximal value of $\lambda _{min} = 1/2$. For this optimal configuration, the matrix $\mathbb {Q}^{T}\mathbb {Q}$ is diagonal:

(33)$${{\mathbb{Q}}^{T}}{\mathbb{Q}} = \frac{1}{4}\left[ {\begin{array}{cc} 1 & {0}\\ {0} & 1 \end{array}} \right]$$

and thus CRLB$[\delta ]$ computed with Eq. (24) is independent of the AoP of the input Stokes vector. Moreover, the CRLB of the Stokes parameters can be computed with Eq. (21). One obtains:

(34)$$\hbox{CRLB}[\delta] = \frac{4}{P \, \hbox{SNR}_\delta^{2}} \quad \hbox{and} \quad \hbox{CRLB}[\mathbf{S}] = {\sigma ^{2}} \left[\frac{2}{3},2,2,2 \right]$$

It is observed that for this configuration, the CRLB of the four Stokes parameters have the same value as that when $\delta$ is known, which can be computed with Eq. (19). Indeed, the second term of Eq. (21) is equal to zero.

We have represented this optimal configuration on the Poincaré sphere in Fig. 3(b). It is interesting to compare it with the configuration defined in Eq. (28) for which autocalibration was not possible for all values of the AoP since $\operatorname {rank}(\mathbb {Q})=1$ (see Fig. 3(a)). We can observe that, contrary to the configuration in Eq. (28), it does not form a regular octahedron structure, but corresponds to a rotation of the measurement states on the lower yellow circle relative to those on the upper yellow circle as depicted in Fig. 3(b). It is very interesting to remember the results in Ref. [23], in which this type of structure corresponds to spherical 2-designs.

As a summary, for $K=6$, among all the configurations that optimize the EWV, some are more attractive than other in terms autocalibration capacity. We have determined the configurations that jointly maximize this capacity and optimize the EWV in the presence of additive noise.

5.3 Autocalibration of a rotating retarder and a DoFP polarization camera (RRDoFP)

Let us now consider the RRDoFP setup. Such polarization imagers have recently aroused strong interest due to the integrated structure of DoFP cameras and to their ability to perform real-time polarimetric measurements [10,14–16,20,21]. In contrast with the two setups considered above, each pixel of a DoFP cameras has a micro-polarizer of fixed orientation in front of it. The linear Stokes vector is estimated from the intensity measurements of a set of 4 neighboring pixels having 4 different micro-polarizer orientations. These sets of pixels are called "superpixels". The advantage of these cameras is to obtain the linear Stokes vector from one single image acquisition, at the price of a reduction of the spatial sampling rate and apparition of spatial imaging artifacts [28–30]. Note that in this article, we only consider estimation within isolated superpixels. We do not investigate the influence of autocalibration on image resolution and artifacts, which is an interesting perspective for future work.

Moreover, it has been shown that the full Stokes vector can be estimated by placing a retarder in front of the DoFP camera and performing at least two image acquisitions with different orientations of the retarder [13,31]. Importantly, it has to be noted that in this setup, one image acquisition corresponds to $4$ intensity measurements since the Stokes vector is estimated from a whole superpixel. In the following, we will denote $N$ the number of image acquisitions, which corresponds to $K=4N$ intensity measurements.

If $N=2$ and the rotation angles of the retarder during the two acquisitions are $\boldsymbol \theta = (\theta _1,\theta _2)$, we determined after cumbersome but elementary calculations that the structure matrix $\mathbb {Q}$ is always of rank one and its null space is spanned by the Stokes vector of AoP$= -(\theta _1 + \theta _2)/2$: For this value of the AoP, autocalibration is not possible. This result is independent of the nominal value of the retardance $\delta$. It generalizes the results of Ref. [9], where it was shown, on a specific example, that autocalibration independent of the AoP requires at least three acquisitions.

Let us now consider that $N\geq 3$ and that the retarder angles used at each image acquisition are evenly spaced. The precision reachable with these configurations has been studied in Ref. [14]. The variances (and thus the CRLB) of the four Stokes parameters, and the EWV are:

(35)$$\hbox{CRLB}[\mathbf{S}]_{\delta \hbox{k}} = \frac{\sigma^{2}}{N} \left\{ 1, \frac{4}{1+c}, \frac{4}{1+c}, \frac{2}{1-c}\right\} \quad \hbox{and} \quad \hbox{EWV}_{\delta \hbox{k}} = \frac{\sigma^{2}}{N} \frac{11-6 c- c^{2}}{1-c^{2}}$$

where $c=\cos ^{2} \delta$ and the subscript $\delta$k stands for "$\delta$ known". In this expression, the number of superpixels used to estimate $\delta$ is equal to $P_s=1$. The value of the EWV is minimal when $c=1/3$, that is, $\delta =54.74^{\circ }$ or $125.26^{\circ }$. In this case, it takes the value $10 \sigma ^{2}/N$. If $\delta =90^{\circ }$, it is slightly larger and equal to $11 \sigma ^{2}/N$.

When $\delta$ is unknown and has to be autocalibrated, it is easily shown that the structure matrix satisfies $\operatorname {rank} [\mathbb {Q}]=2$ for any value of $\delta$. By applying Eq. (20), one obtains:

(36)$$\hbox{CRLB}[\delta] = \frac{4}{N} \frac{1}{P_s \, \hbox{SNR}_\delta^{2}} \frac{1+c}{1-c}$$

where $P_s$ is the number of superpixels used to estimate $\delta$. This result is consistent with Ref. [10], where a similar result has been derived with a different method. It is seen in Eq. (36) that $\hbox {CRLB}[\delta ]$ is independent of the AoP and depends only on the SNR$_\delta$ defined in Eq. (24). It is also seen that the value of $\hbox {CRLB}[\delta ]$ is minimal for $\delta =90^{\circ }$. The CRLBs of the Stokes parameters can be found from Eq. (21):

(37)$$\hbox{CRLB} [\mathbf{S}]_{\delta \hbox{uk}} = \hbox{CRLB} [\mathbf{S}]_{\delta \hbox{k}} + \frac{1}{P_s} \frac{4 \sigma^{2}}{N} \left\{ 0 , \frac{c}{1+c} \cos^{2} 2 \alpha, \frac{c}{1+c} \sin^{2} 2 \alpha, \frac{c(1 + c)}{(1-c)^{2}} \tan^{2} 2\varepsilon \right\}$$

where the subscript $\delta \hbox {uk}$ stands for "$\delta$ unknown". It is seen that the fact of estimating $\delta$ adds a non-negative value to the CRLB, which is consistent with Eq. (21). This value is zero for the parameter $S_0$, thus the CRLB of intensity estimation is not modified. Moreover, the CRLB of $S_1$ and $S_2$ are increased by an amount that depends on the AoP, but their values remain bounded. In contrast, the CRLB of $S_3$ diverges as the ellipticity $\varepsilon$ (defined in Eq. (1)) tends to $\pm 45^{\circ }$, that is, when the input Stokes vector becomes circular. However, it has to be noted that if $c=0$, that is, $\delta =90^{\circ }$, the CRLB of $\mathbf {S}$ is independent of the AoP and the ellipticity: one has $\hbox {CRLB} [\mathbf {S}]_{\delta \hbox {uk}} = \hbox {CRLB} [\mathbf {S}]_{\delta \hbox {k}}$.

The EWV, defined as the sum of the CRLBs of the four Stokes parameters defined in Eq. (37), has the following value:

(38)$$\hbox{EWV}_{\delta \hbox{uk}} = \hbox{EWV}_{\delta \hbox{k}} + \frac{1}{P_s} \frac{4\sigma^{2}}{N} \left( \frac{c}{1+c}+ \frac{c(1 + c)}{(1 - c)^{2}} \tan^{2} 2\varepsilon \right)$$

It does not depend on the AoP of the input Stokes vector, but depends in general on its ellipticity $\varepsilon$. There is however an exception when $\delta = 90^{\circ }$, since in this case, the EWV is totally independent of the ellipticity and always equal to $11 \sigma ^{2}/N$. We have plotted in Fig. 4(a) the values of the EWV when $\delta$ is known (see Eq. (35)) and when it is unknown (see Eq. (38)) for different values of the ellipticity $\varepsilon$ as a function of the nominal value of $\delta$. We assumed that $P_s=1$. It is seen that the value of $\delta$ for which the EWV is minimal depends on $\varepsilon$. It is equal to $65.55^{\circ }$ or $114.45^{\circ }$ when $\varepsilon =0$, and converges to $90^{\circ }$ as $\varepsilon$ increases. It is easily shown to be exactly equal to $90^{\circ }$ as soon as $\varepsilon \geq \arcsin (1/\sqrt {3})$.

Fig. 4. Value of EWV$_{\delta \hbox {uk}}$ as a function of the nominal value of the retardance $\delta$ for different values (a) of the ellipticity $\varepsilon$, and (b) of number of superpixels $P_s$. In (a), $P_s=1$, and in (b), $\varepsilon = 30^{\circ }$.

Download Full Size | PDF

Moreover, it is seen in Eq. (38) that the second additive term due to estimating $\delta$ is inversely proportional to the number $P_s$ of superpixels used for estimation. We have plotted in Fig. 4(b) the value of EWV$_{\delta \hbox {uk}}$ as a function of the nominal value of $\delta$ for $\varepsilon =30^{\circ }$ and for different values of $P_s$. It can be seen that, with the increase of $P_s$, the values of EWV tends to that obtained when $\delta$ is known. That means that, in practice, the increase of CRLB incurred by the fact of estimating $\delta$ can be reduced by using enough different super pixels for estimation

As a summary, the optimal value of the nominal retardance $\delta$ for autocalibration depends on $\varepsilon$ and $P_s$. However, the value of $\hbox {EWV}_{\delta \hbox {uk}}$ obtained with $\delta =90^{\circ }$ is always $11 \sigma ^{2}/N$, whereas it varies from $10 \sigma ^{2}/N$ to infinity when $\delta =125.26^{\circ }$. Hence it can be considered that the preferable nominal value of retardance is $\delta = 90^{\circ }$.

5.4 Experimental validation

The results presented in this article are based on the CRLB, which is a lower bound on the estimation variance that can be reached with an unbiased estimator. They thus represent potential performance. It is important to check that this performance can be actually reached on data from real-world polarization cameras. In order to perform this experimental validation, we used a DoFP camera made by LUCID Vision Labs equipped with the polarization sensor Sony IMX250MZR CMOS. The input Stokes vector is generated by a halogen lamp followed by a spectral bandpass filer at 650nm with 10nm width and a polarizer to produce a narrow band, fully linearly polarized light oriented at $\alpha = 0^{\circ }$. We place in front of the DoFP camera a rotatable quarter wave plate (QWP), and we perform $N=3$ acquisitions with 3 different orientations of the QWP at $0^{\circ }, 60^{\circ }$ and $120^{\circ }$ in order to estimate the input Stokes vector.

We have plotted in Fig. 5 the experimental standard deviation (STD) of estimation of $\delta$, $S_3$ and $S_1$ as the function of the SNR, defined as SNR$=S_0/\sigma$. It is to be noted that in this setup, the main source of noise is not the additive noise introduced by the camera but the photon noise. It can been shown that if the Stokes vector intensity $S_0$ on the camera (expressed in number of photo-electrons) is fixed, photon noise is equivalent to an additive (albeit non Gaussian) noise of variance $\sigma ^{2}=S_0/2$ [32]. One thus has SNR$=\sqrt {2S_0}$. In consequence, one just has to change the intensity of illumination in order to change the SNR in experiments. For each value of SNR, we perform $10^{4}$ acquisitions of the Stokes vector and estimate the STD of the different parameters on each single pixel from this set of data.

Fig. 5. Estimation standard deviation (dotted lines with markers) and CRLB$^{1/2}$ (solid lines) as function of the SNR for the following parameters: (a) $\delta$, (b) $S_3$, and (c) $S_1$. The input Stokes vector is such that $\alpha = 0^{\circ }$ , $\varepsilon = 0^{\circ }$, and the measurement configuration is "0-60-120". In (b) and (c), the standard deviations of $S_3$ and $S_1$ are divided by $S_0 \times \hbox {DoLP}$.

Download Full Size | PDF

Figure 5(a) represents the experimental STD of ${\delta }$ (red dotted line with markers "I") and CRLB$[\delta ]^{1/2}$ (solid blue line). It is observed that for high values of the SNR, the STD is very close to the CRLB, which means that the CRLB is an accurate representation of the actual estimation variance. For the lowest SNR values, one observes, however, a slight deviation. This deviation happens when the SNR value is around 6, and the corresponding STD of $\delta$ is around $7^{\circ }$.

Figure 5(b) represents the estimation STD of $S_3$ using the estimated value of $\delta$ (red dotted lines with markers "I") and using the known value $\delta = 90^{\circ }$ (green dotted lines with markers "o"). The values of the STD are divided by a factor $S_0 \times \hbox {DoLP}$. Both $S_0$ and DoLP are estimated and averaged over all the $10^{4}$ acquisitions and over 10 pixels in order to get reliable estimation. The value of $\hbox {CRLB} [S_3]_{\delta \hbox {uk}}^{1/2}$ (see Eq. (37)) is also plotted for comparison. It is noted that in our case, $\hbox {CRLB} [\mathbf {S}]_{\delta \hbox {uk}}=\hbox {CRLB} [\mathbf {S}]_{\delta \hbox {k}}$ since $\delta = 90^{\circ }$ (see Section 5.3). We observe that for high values of the SNR, the experimental value of the STD is equal to the CRLB. However, from a value of the SNR around 6, one observes a sharp departure of the STD with respect to the CRLB. By looking back at Fig. 5(a), one can note that this SNR value corresponds to the already observed small divergence of the STD of ${\delta }$ with respect to the CRLB, and to a value of this STD around $7^{\circ }$. The same behavior was observed in Ref. [10] on Monte Carlo simulations. The reason for this phenomenon is that when the estimation of the parameter $\delta$ is too noisy, the Stokes vector $S_3$, that relies on it, cannot be properly estimated. It is known that the CRLB accurately represents the actual STD only if the estimated value of the parameter is not too far from the true value [18]. This condition is less and less respected as the SNR decreases. We have also represented in Fig. 5(c) the estimation STD of the Stokes parameter $S_1$, divided by a factor $S_0\times \hbox {DoLP}$. We observe that the STD fits the CRLB for any value of the SNR. Thus, in contrast with $S_3$, the estimation accuracy of $\delta$ has no influence on the STD of $S_1$. We have checked that this is also the case with $S_2$.

As a conclusion, these experimental results validate the use of the CRLB a surrogate of actual STD for studying autocalibration feasibility and performance of polarization imagers. One just has to bear in mind that the CRLB expressions given in this article are valid only for sufficient values of the SNR, or, which is equivalent, when the estimation STD of $\delta$ is not too large (less than $7^{\circ }$ in our case).

6. Summary and conclusion

We have proposed in this article a general theory to characterize the feasibility of joint estimation of Stokes vector and parameters of polarimetic imagers. It gives a simple way to determine the feasibility and the precision of autocalibration of any imager. This theory has been illustrated on the problem of retardance autocalibration in three of the most common polarimetric imaging architectures. The obtained results are summarized in Table 1. The main conclusion is that although these three setups can achieve the similar values of the EWV when they are properly optimized, their properties regarding autocalibration are very different. With the RRFP setup, autocalibration is impossible whatever the number of intensity measurements. For the other setups, autocalibration is possible, and increasing the number $K$ of acquisitions widens its domain of feasibility: to get the maximal domain of feasibility, one needs at least 6 intensity measurements for the RRRP setup and 3 image acquisitions for the RRDoFP setup.

Table 1. Feasibility of retardance autocalibration for different setups. $\alpha _0$ denotes a given value of the AoP. $N$ denotes the number of image acquisitions for RRDoFP setup and $K$ denotes the number of intensity measurements for RRFP or RRRP setup.

View Table

These results are important for choosing a polarimetric setup adapted to the precision and autocalibration requirements of the application at hand. It should be noted that the approach developed in this article can be applied to characterize the autocalibration feasibility and performance of any imaging setup, and of any other parameter than retardance. Examples of such parameters are the direction and diattenuation of the polarizer or of the micro-polarizers, and the direction of the retarder, be it a waveplate or a liquid crystal variable retarder (LCVR). Of course, for some parameters or combination of parameters, our approach could lead to the conclusion that autocalibration is not feasible. Other interesting perspectives are to evaluate the influence of the arrangement of the micro-polarizer array (see Ref. [33] and references therein) on the feasibility and precision of autocalibration, and to generalize this approach to multi-retarder-based polarimeters, for example the Mueller matrix polarimeters [34,35].

A. Rank of the structure matrix $\mathbb {Q}$

According to Eq. (22), one has:

(39)$$\forall k\in [1,K], ~~ \left[\frac{\partial \mathbb{W}}{\partial \delta}\right]_k = \frac{\partial \mathbf{T}_k}{\partial \delta} = \left[0, - {a}_k \, \sin \delta , - {c}_k \, \sin \delta , {e}_k \, \cos \delta \right]^{T}.$$

It is thus obvious that the first column of the matrix $\mathbb {P}^{\perp }_{W} \frac {\partial \mathbb {W}}{\partial \eta }$ is zero. Moreover, the last column of $\frac {\partial \mathbb {W}}{\partial \eta }$ is equal to $(\cos \delta ) \, \mathbf {e}$, where $\mathbf {e}=\left [e_1, e_2,\ldots e _K\right ]^{T}$. It is thus collinear to the fourth column of $\mathbb {W}$, which is equal to $(\sin \delta ) \, \mathbf {e}$. Since $\mathbb {P}_{W}^{\perp }$ is the projector orthogonal to the subspace spanned by the rows (or columns) of the matrix $\mathbb {W}$, one has $\mathbb {P}_{W}^{\perp } \mathbf {e} = \mathbf {0}$. Hence the fourth column of $\mathbb {P}^{\perp }_{W} \frac {\partial \mathbb {W}}{\partial \eta }$ is zero.

B. Demonstration of $\bf {\operatorname {rank}(\mathbb {Q}) =0}$ for the RRFP setup

Let us consider the input Stokes vector $\mathbf {S}_a = \left [1,1,0,0\right ]^{T}$ and $\mathbf {S}_b = \left [1,0,1,0\right ]^{T}$. For any orientation $\varphi$ of the fixed polarimeter, elementary calculations show that:

(40)$$\frac{\partial \mathbb{W}}{\partial \delta} \mathbf{S}_a = \alpha\mathbb{W}_{{\cdot} 1} + \beta \mathbb{W}_{{\cdot} 2} ~~~\hbox{and}~~~ \frac{\partial \mathbb{W}}{\partial \delta} \mathbf{S}_b = \lambda \mathbb{W}_{{\cdot} 1} + \mu \mathbb{W}_{{\cdot} 3},$$

where $\mathbb {W}_{\cdot k}$ denotes the $k^{th}$ column of the matrix $\mathbb {W}$, $\alpha = - \beta \cos 2\varphi$, $\beta = \frac {{\sin \delta }}{{1 - \cos \delta }}$, $\lambda = - \mu \sin 2\varphi$, and $\mu = \frac {{\sin \delta }}{{1 - \cos \delta }}$. Since the two values in Eq. (40) are linear combinations of columns of $\mathbb {W}$, one has, for the same reason as in Appendix A.,

(41)$$\mathbb{P}^{{\perp}}_{W} \frac{\partial \mathbb{W}}{\partial \eta} \mathbf{S}_a = \mathbb{P}^{{\perp}}_{W} \frac{\partial \mathbb{W}}{\partial \eta} \mathbf{S}_b = \mathbf{0} \Longrightarrow \mathbb{Q} \mathbf{\bar{s}}_a = \mathbb{Q} \mathbf{\bar{s}}_b = \mathbf{0}$$

where the vector $\mathbf {\bar {s}}$ is defined in Eq. (25). Since the rank of the matrix $\mathbb {Q}$ is less than or equal to 2, and the vectors $\mathbf {\bar {s}}_a$ and $\mathbf {\bar {s}}_b$ are linearly independent, it means that $\mathbb {Q} = \mathbb {O}$ and $\operatorname {rank}(\mathbb {Q}) =0$.

Funding

National Natural Science Foundation of China (61775163); Direction Générale de l’Armement; Mission pour la Recherche et l’Innovation Scientifique (MRIS).

Disclosures

The authors declare no conflicts of interest.

References

1. V. Thilak, D. G. Voelz, and C. D. Creusere, “Polarization-based index of refraction and reflection angle estimation for remote sensing applications,” Appl. Opt. 46(30), 7527–7536 (2007). [CrossRef]

2. X. Li, H. Hu, L. Zhao, H. Wang, Y. Yu, L. Wu, and T. Liu, “Polarimetric image recovery method combining histogram stretching for underwater imaging,” Sci. Rep. 8(1), 12430 (2018). [CrossRef]

3. M. Foldyna, A. De Martino, R. Ossikovski, E. Garcia-Caurel, and C. Licitra, “Characterization of grating structures by mueller polarimetry in presence of strong depolarization due to finite spot size,” Opt. Commun. 282(5), 735–741 (2009). [CrossRef]

4. D. H. Goldstein, Polarized light (CRC, 2016).

5. J. S. Tyo, D. L. Goldstein, D. B. Chenault, and J. A. Shaw, “Review of passive imaging polarimetry for remote sensing applications,” Appl. Opt. 45(22), 5453–5469 (2006). [CrossRef]

6. R. Azzam, I. Elminyawi, and A. El-Saba, “General analysis and optimization of the four-detector photopolarimeter,” J. Opt. Soc. Am. A 5(5), 681–689 (1988). [CrossRef]

7. X. Li, T. Liu, B. Huang, Z. Song, and H. Hu, “Optimal distribution of eintegration time for intensity measurements in stokes polarimetry,” Opt. Express 23(21), 27690–27699 (2015). [CrossRef]

8. D. Sabatke, M. Descour, E. Dereniak, W. Sweatt, S. Kemme, and G. Phipps, “Optimization of retardance for a complete stokes polarimeter,” Opt. Lett. 25(11), 802–804 (2000). [CrossRef]

9. S. Shibata, N. Hagen, and Y. Otani, “Robust full stokes imaging polarimeter with dynamic calibration,” Opt. Lett. 44(4), 891–894 (2019). [CrossRef]

10. F. Goudail, X. Li, M. Boffety, S. Roussel, T. Liu, and H. Hu, “Precision of retardance autocalibration in full-stokes division-of-focal-plane imaging polarimeters,” Opt. Lett. 44(22), 5410–5413 (2019). [CrossRef]

11. Z. Sekera, “Light scattering in the atmosphere and the polarization of sky light,” J. Opt. Soc. Am. 47(6), 484–490 (1957). [CrossRef]

12. D. Sabatke, M. Descour, E. Dereniak, W. Sweatt, S. Kemme, and G. Phipps, “Optimization of retardance for a complete stokes polarimeter,” Opt. Lett. 25(11), 802–804 (2000). [CrossRef]

13. S. Roussel, M. Boffety, and F. Goudail, “On the optimal ways to perform full stokes measurements with a linear division-of-focal-plane polarimetric imager and a retarder,” Opt. Lett. 44(11), 2927–2930 (2019). [CrossRef]

14. X. Li, H. Hu, F. Goudail, and T. Liu, “Fundamental precision limits of full stokes polarimeters based on dofp polarization cameras for an arbitrary number of acquisitions,” Opt. Express 27(22), 31261–31272 (2019). [CrossRef]

15. S. B. Powell and V. Gruev, “Calibration methods for division-of-focal-plane polarimeters,” Opt. Express 21(18), 21039–21055 (2013). [CrossRef]

16. S. Roussel, M. Boffety, and F. Goudail, “Polarimetric precision of micropolarizer grid-based camera in the presence of additive and poisson shot noise,” Opt. Express 26(23), 29968–29982 (2018). [CrossRef]

17. C. D. Meyer, Matrix analysis and applied linear algebra, vol. 71 (Siam, 2000).

18. S. M. Kay, Fundamentals of statistical signal processing (Prentice Hall PTR, 1993).

19. P. Réfrégier, Noise theory and application to physics: from fluctuations to information (Springer Science & Business Media, 2004).

20. Z. Chen, X. Wang, and R. Liang, “Calibration method of microgrid polarimeters with image interpolation,” Appl. Opt. 54(5), 995–1001 (2015). [CrossRef]

21. B. Feng, Z. Shi, H. Liu, L. Liu, Y. Zhao, and J. Zhang, “Polarized-pixel performance model for dofp polarimeter,” J. Opt. 20(6), 065703 (2018). [CrossRef]

22. T. Mu, Z. Chen, C. Zhang, and R. Liang, “Optimal design and performance metric of broadband full-stokes polarimeters with immunity to poisson and gaussian noise,” Opt. Express 24(26), 29691–29704 (2016). [CrossRef]

23. M. R. Foreman, A. Favaro, and A. Aiello, “Optimal frames for polarization state reconstruction,” Phys. Rev. Lett. 115(26), 263901 (2015). [CrossRef]

24. N. J. A. Sloane and R. H. Hardin, Spherical designs in 3d, 2002, http://neilsloane.com/sphdesigns/dim3.

25. F. Goudail, “Equalized estimation of stokes parameters in the presence of poisson noise for any number of polarization analysis states,” Opt. Lett. 41(24), 5772–5775 (2016). [CrossRef]

26. R. Bhatia, Matrix analysis, vol. 169 (Springer Science & Business Media, 2013).

27. J. Nocedal and S. Wright, Numerical optimization (Springer Science & Business Media, 2006).

28. B. M. Ratliff, C. F. LaCasse, and J. S. Tyo, “Interpolation strategies for reducing ifov artifacts in microgrid polarimeter imagery,” Opt. Express 17(11), 9112–9125 (2009). [CrossRef]

29. J. S. Tyo, C. F. LaCasse, and B. M. Ratliff, “Total elimination of sampling errors in polarization imagery obtained with integrated microgrid polarimeters,” Opt. Lett. 34(20), 3187–3189 (2009). [CrossRef]

30. S. Gao and V. Gruev, “Gradient-based interpolation method for division-of-focal-plane polarimeters,” Opt. Express 21(1), 1137–1151 (2013). [CrossRef]

31. J. Qi, C. He, and D. S. Elson, “Real time complete stokes polarimetric imager based on a linear polarizer array camera for tissue polarimetric imaging,” Biomed. Opt. Express 8(11), 4933–4946 (2017). [CrossRef]

32. J. Dai and F. Goudail, “Precision analysis of arbitrary full-stokes polarimeters in the presence of additive and poisson noise,” J. Opt. Soc. Am. A 36(7), 1229–1240 (2019). [CrossRef]

33. A. S. Alenin, I. J. Vaughn, and J. S. Tyo, “Optimal bandwidth and systematic error of full-stokes micropolarizer arrays,” Appl. Opt. 57(9), 2327–2336 (2018). [CrossRef]

34. J. J. G. Perez and R. Ossikovski, Polarized light and the Mueller matrix approach (CRC, 2017).

35. R. Ossikovski and O. Arteaga, “Instrument-dependent method for obtaining a nondepolarizing estimate from an experimental mueller matrix,” Opt. Eng. 58(08), 1 (2019). [CrossRef]

Setup	$K$ (or $N$ )	rank( $Q$ )	Invalid $S$
RRFP	$\geq 5$	0	All
RRRP	5	1	AoP = $α_{0}$
	$\geq 6$	1	AoP = $α_{0}$
		2	DoLP = 0
RRDoFP	2	1	AoP = $α_{0}$
	$\geq 3$	2	DoLP = 0

Theory of autocalibration feasibility and precision in full Stokes polarization imagers

Abstract

1. Introduction

2. General theory for autocalibration

3. Autocalibration of a single parameter

4. Autocalibration in the case of single-retarder-based polarimeters

5. Autocalibration for different setups of Stokes imagers

5.1 Autocalibration of the RRFP setup

5.2 Autocalibration of the RRRP setup

5.3 Autocalibration of a rotating retarder and a DoFP polarization camera (RRDoFP)

5.4 Experimental validation

6. Summary and conclusion

A. Rank of the structure matrix $\mathbb {Q}$

B. Demonstration of $\bf {\operatorname {rank}(\mathbb {Q}) =0}$ for the RRFP setup

Funding

Disclosures

References

Cited By

Figures (5)

Tables (1)

Equations (41)

Optics Express