1. INTRODUCTION

In various application domains such as biomedical diagnosis, defense, or remote sensing, standard-intensity imaging techniques sometimes fail to reveal relevant contrasts or gather sufficient information. In these domains, polarimetric approaches have proved efficient at enhancing the estimation or detection capabilities of the imaging systems [1–7]. Most often, the polarimetric information is provided by a scalar polarization contrast image, which offers complementary contrast information with respect to the conventional intensity image. This polarization contrast image usually corresponds to the map of the degree of polarization (DOP) of the light backscattered at each location of the scene. Four polarimetric measurements are theoretically needed to determine such a DOP image [8]. However, to reduce both cost and acquisition time, simplified two-channel imaging modalities are usually preferred in active polarimetric imaging, such as in orthogonal states contrast (OSC) imaging. This approach consists of computing a contrast map $\mathrm{OSC}=({\mathbf{x}}_S-{\mathbf{x}}_P)/({\mathbf{x}}_S+{\mathbf{x}}_P)$ between two polarimetric images (${\mathbf{x}}_S$ and ${\mathbf{x}}_P$) of the scene, acquired through a linear polarizer oriented along two orthogonal directions; these are denoted $S$ and $P$ throughout this paper, where $S$ denotes the polarization direction of the illumination source. Owing to its instrumental simplicity, which can be further improved using voltage-controlled electro-optic devices [9,10] or appropriate optical design [11], and due to the fact that OSC provides a good estimate of the DOP (under the assumption that the imaged objects are purely depolarizing [4,12]), OSC imaging is today widely used in many applications [12].

On the other hand, compressive sensing (CS) and CS-derived original imaging concepts such as the single-pixel camera (SPC) have attracted much attention in recent years [13,14]. With this approach, the measurement process relies on the spatial sampling of the image of interest with a digital micromirror device (DMD), and on numerical reconstruction of the image from intensity measurements on a single photodetector for different sampling patterns on the DMD, allowing a compressed version of the image to be recovered from the photocurrent signal acquired. More recently, the concept of SPC has been applied to a number of domains including, among others, multi/hyperspectral imaging [15–18], THz imaging [14], or random media-assisted CS [19]. However, despite the burgeoning interest in CS, few attempts to perform polarimetric CS imaging have been reported so far [16,20–23]. The imaging setups proposed in these works are all directly based on the SPC concept, where polarimetric sensitivity was simply gained by detecting the optical signals through appropriate polarization-analyzing devices during sequential acquisitions, or with a unique acquisition on several detectors after appropriate beam splitting of the light reflected by the DMD. In these works, the reconstruction process consisted of solving as many CS reconstruction problems as polarimetric channels were considered (two or four). More precisely, in Refs. [20,21], the SPC scheme was readily improved by adding a rotating polarizer in front of the detector. These techniques can operate as polarimetric imaging systems only at the expense of a two-fold (respectively four-fold) increase in the measurement time, while at the same time suffering from a loss in intensity due to the use of a polarization analyzer. In Refs. [22,23], the measurement time was limited, but the complexity of the imaging system was largely increased, to the expense of important losses in the imaging setup.

In this paper, we revisit the problem of two-channel polarimetric CS by proposing an original polarimetric imaging architecture using two single-pixel detectors. The proposed setup is inspired from the initial concept of SPC, but does not require any polarization-analyzing element as it relies on the tiny polarization dependence of reflection coefficients of the DMD itself. Contrary to previous attempts in polarimetric CS, the polarimetric information is obtained through a single temporal data acquisition step on the two photodetectors, and the polarimetric channels are recovered simultaneously from a single reconstruction step. It also offers in principle the best detectivity tradeoff, as all the light directed towards the DMD is involved in the imaging process without passing through any polarization analysis component. As with the SPC concept, one interesting aspect of this approach is its potential of being applied to any spectral domain, even at wavelengths where neither focal plane arrays nor polarization analyzing components (polarizers, waveplates) are available. It can be noted that recent works in compressive imaging have proposed to simultaneously detect the light reflected in both directions of the micromirrors, either to optimize measurement dynamics [24], light throughput [25], or to combine two types of detectors [26], but never, to our best knowledge, to gain polarimetric sensitivity as proposed in this paper.

This paper is organized as follows: in Section 2, the principle and the optical setup proposed to achieve CS polarimetric imaging are detailed. Then in Section 3, various algorithms are presented to tackle this original CS inverse problem, along with several possible algorithmic optimizations such as reweighted approaches or constrained algorithms that can be implemented to improve the reconstruction results. The performance of the algorithms and of the two-pixel polarimetric CS imaging approach are finally discussed in Section 4 through numerical simulations on 1D and 2D test signals. The influence of experimental imperfections on the reconstruction quality is also analyzed in this section, before conclusions and perspectives are provided in Section 5.

2. PRINCIPLE OF TWO-PIXEL COMPRESSIVE SENSING POLARIMETRIC IMAGING

A. Two-Pixel CS Polarimetric Imaging Setup

The polarimetric CS imaging approach proposed relies on the SPC imaging architecture. As will be demonstrated, this makes it possible to perform standard intensity and polarimetric contrast imaging using CS, without requiring any polarization analysis component. The corresponding experimental setup is sketched in Fig. 1. In the context of active imaging, we assume that the scene or object of interest is enlightened by a narrow-band light source (such as a laser or LED). An image of the scene is formed through a lens $L$ onto the surface of the DMD, which spatially samples it by applying a controlled binary pattern to the micromirrors. The detection setup is strictly similar to the SPC architecture, with a first photodetector ($P_1$) used to detect light reflected in the first reflection direction through a lens ($L^{\prime }$). However, the light reflected in the second direction of the tilted mirrors is directed towards a second photodetector ($P_2$) via a lens ($L^{\prime }$), instead of being discarded as in the original SPC scheme. By applying a series of different patterns on the DMD, the detection of the total light intensity reflected in directions 1 and 2 by photodetectors $P_1$ and $P_2$ provides two temporal signals that are sampled and digitized on an analog to digital converter before they can be used for numerical inversion of the intensity and polarimetric contrast images.

 

Fig. 1. Sketch of the proposed two-pixel CS polarimetric imaging setup. It is inspired from the concept of SPC, where the image is first spatially sampled by a DMD which reflects light in two directions, and where the total photon flux is detected on a single photodetector in each direction. The illuminating light source can be either polarized (as represented) or not.

Download Full Size | PDF

Throughout this paper, we will denote the total intensity image of the scene by a single dimensional row vector ${\mathbf{x}}_T={\{x_{T_i}\}}_{i=1,\dots ,N}$ of length $N$. We assume that the total intensity image can be written as the sum of two polarimetric components, i.e., ${\mathbf{x}}_T={\mathbf{x}}_S+{\mathbf{x}}_P$. Subscripts $S$ and $P$ denote two orthogonal linear polarization directions w.r.t. the orientation of the DMD surface and to the direction of linear polarization of the illuminating beam ($S$), as sketched in Fig. 1. We assume that these two components are compressible in the same sparse representation ${\{{\mathrm{\Psi }}_k\}}_{k=1,\dots ,N^{\prime }}$ (with ${\mathrm{\Psi }}_k$ row vectors), i.e., they can be written as ${\mathbf{x}}_{S,P}={\mathbf{s}}_{S,P}\mathbf{\Psi }$, where $\mathbf{\Psi }=[{\mathrm{\Psi }}_1,\dots ,{\mathrm{\Psi }}_{N^{\prime }}]$ is a $N^{\prime }\times N$ matrix with $N^{\prime }\ge N$ and ${\mathbf{s}}_{S,P}$ are $N^{\prime }\times 1$ row vectors containing the expansion coefficients. In the compressed sensing framework, compressibility means that most of these expansion coefficients have small amplitude. Only the few large-amplitude coefficients code for the salient information of the polarimetric signals. This assumption generally refers to the approximately sparse signal model in the CS literature [27].

Upon reflection on the DMD surface, the polarimetric components of the original image ${\mathbf{x}}_T$ are altered by the Fresnel’s reflection coefficients (in intensity) of the mirror, depending on the reflection angle. Let us denote by $r_i^S$ and $r_i^P$ the Fresnel’s reflection coefficients (in intensity) of the mirror for each tilt direction $i=1,2$ and respectively for the $S$ and $P$ polarimetric components of the image formed on the DMD surface. Neglecting absorption effects, these coefficients are real and verify $r_{1,2}^{S,P}\in [0,1]$. With such notations, the images respectively reflected towards detectors $P_1$ and $P_2$ read ${\tilde{\mathbf{y}}}_1^{°}=r_1^S{\mathbf{x}}_S+r_1^P{\mathbf{x}}_P$ and ${\tilde{\mathbf{y}}}_2^{°}=r_2^S{\mathbf{x}}_S+r_2^P{\mathbf{x}}_P$, which can be rewritten in a compact form as

When a given pattern indexed by $k$ is applied on the DMD to spatially sample the image, the total intensity reflected towards direction 1 and integrated on detector $P_1$ can be denoted by ${\tilde{y}}_1^{(k)}={\tilde{\mathbf{y}}}_1^{°}{\varphi }^{(k)}$, where ${\varphi }^{(k)}$ is a binary valued $N$-dimensional column vector encoding the set of orientations of the individual mirrors (DMD pattern). Similarly, ${\tilde{y}}_2^{(k)}={\tilde{\mathbf{y}}}_2^{°}{\overline{\varphi }}^{(k)}$, where ${\overline{\varphi }}^{(k)}={\mathbf{1}}_N-{\varphi }^{(k)}$ is the complement of vector ${\varphi }^{(k)}$. Then, when $M$ measurements are accumulated with various sets of pseudo-random configurations of the DMD, the detected intensities ${\{{\tilde{y}}_i^{(1)}\dots {\tilde{y}}_i^{(M)}\}}_{i=1,2}$ form a $2\times M$ measurement matrix

Under these assumptions, we will show that such a simple setup suffices to retrieve a compressive measurement of the polarimetric components ${\mathbf{x}}_S$ and ${\mathbf{x}}_P$, and thus of the intensity image ${\mathbf{x}}_T={\mathbf{x}}_S+{\mathbf{x}}_P$ and OSC image $\mathrm{OSC}=({\mathbf{x}}_{\mathbf{S}}-{\mathbf{x}}_{\mathbf{P}})/{\mathbf{x}}_T$. It is worth noting here that all the available light incoming from the scene is detected, thereby offering optimal energy balance, and that no polarimetric optical component has been inserted in the setup. Indeed, to achieve polarimetric sensitivity, this CS imaging setup relies on the slight variation of the Fresnel coefficients for the $S$ and $P$-polarized components of light as a function of the angle of incidence of the incoming light beam. As sketched in Fig. 1, the angle of incidence on the mirror is denoted by ${\theta }_1=\theta -t_1$ (respectively ${\theta }_2=\theta -t_2$) for tilt direction 1 (respectively direction 2), where $\theta$ is the angle of incidence with respect to the DMD surface, and where $t_1$ and $t_2$ denote the tilt angle of the mirrors (typically $t_1=12°$ and $t_2=-12°$ on most DMDs [28]). This dependency with $\theta$ is illustrated in Fig. 2(a), where we plotted the evolution of the four reflection coefficients $r_{1,2}^{S,P}$ at wavelength 780 nm for aluminum mirrors when $\theta$ varies from 17° to 65°. The influence of the value of these coefficients on the CS reconstruction problem is discussed in Section 2.C. At this level, it is interesting to note that polarization vision mechanisms in some animal species rely on such polarization sensitivity of the Fresnel reflection or refraction coefficients [29].

 

Fig. 2. (a) Evolution of the reflection coefficients in intensity for two tilt directions and two polarimetric components $S$ and $P$, as a function of incidence angle $\theta$ on the DMD surface at wavelength 780 nm. (b) Evolution of the condition number $\kappa (\mathbf{A})$ of matrix $\mathbf{A}$ as a function of $\theta$ at 780 nm. (c) Contour plot of $\kappa (\mathbf{A})$ as a function of $\theta$ and wavelength. The black cross indicates the physical situation addressed in the numerical experiments.

Download Full Size | PDF

B. Description of the CS Inverse Problem

With the preceding notations, we will now describe this imaging process as a CS inverse problem. We first rewrite the binary sensing matrix $\mathit{\varphi }$ as $\mathit{\varphi }=({\mathbf{1}}_{N,M}+\mathbf{\varphi })/2$ and consequently $\overline{\mathit{\varphi }}=({\mathbf{1}}_{N,M}-\mathbf{\varphi })/2$, where the elements of $\mathbf{\varphi }$ take on $-1$ and $+1$ values. We impose that, for all DMD configurations, 50% of the micromirrors are oriented towards direction 1 such that $\sum _{i=1}^N{\mathrm{\varphi }}_i^{(k)}=0$. In such a case, the measurement matrix can then be rewritten as $\tilde{\mathbf{Y}}=\overline{\mathbf{Y}}+\mathbf{Y}$, with $\overline{\mathbf{Y}}=\tilde{\mathbf{Y}}°{\mathbf{1}}_{N,M}/2$ and

For the sake of simplicity, we will assume in the following that the constant term $\overline{\mathbf{Y}}$ can be easily estimated and subtracted out from the measured data, e.g., by averaging the acquisitions over all $M$ DMD pattern realizations considered. In this case, the CS problem that must be solved is given in Eq. (3) and reads $\mathbf{Y}=\mathbf{Q}\tilde{\mathbf{Y}}°\mathbf{\varphi }/2=\mathbf{Q}\tilde{\mathbf{A}}\mathbf{X}\mathbf{\varphi }/2$. For the sake of clarity, and taking into account an additive noise contribution $\mathbf{b}$ on the detected intensities, we propose to rewrite the problem with simplified notations as

We also introduce $\mathbf{Y°}=\mathbf{AX}$, such that $\mathbf{Y}=\mathbf{Y°}\mathbf{\varphi }$.

As a result, the $N$-dimensional polarimetric components of the image contained in $\mathbf{X}$ can be in principle recovered from a number of measurements $M\ll N$ provided the problem described in the preceding equation can be solved. Contrarily to most CS inverse problems that have been considered so far, we are facing an additional difficulty in this particular situation, as the signals to recover are strongly mixed in the measurement process via the mixing matrix $\mathbf{A}$. Indeed, as illustrated in Fig. 2(a), the reflection coefficients on metals are usually quite similar for polarization directions $S$ and $P$, causing a strong crosstalk between the two components of interest. As a consequence, the signals detected at photodetectors $P_1$ and $P_2$ are almost perfectly anticorrelated, the polarimetric information lying in the tiny discrepancies between these two signals. This is illustrated in Fig. 3, where simulated intensity signals are plotted. As will be shown in Section 3, several approaches can be used to tackle this unmixing/CS reconstruction problem, either by considering the two problems independently or by solving them simultaneously in the recovery process.

 

Fig. 3. Inset (a): synthetic 1D polarimetric test signal used to assess reconstruction algorithm performance. An example of reconstructed signal is also given (see text for details). Inset (b): example of measured intensities on photodetectors $P_1$ and $P_2$ for $M=306$ (compression 40%) different binary patterns (Hadamard) applied on the DMD. Main figure: plot of intensity ${\mathbf{y}}_2$ as a function of ${\mathbf{y}}_1$, revealing strong anticorrelation between the two detected signals.

Download Full Size | PDF

C. Experimental Parameters and Imperfections

Before we detail the reconstruction algorithms used to achieve polarimetric CS imaging with the two-pixel camera setup proposed, let us analyze the possible influence of some experimental parameters on the reconstruction quality, and how these parameters can be optimized.

Obviously, an important parameter that will control the difficulty of the unmixing problem is the mixing matrix $\mathbf{A}$, which depends on the wavelength and bandwidth of the illuminating source, on the incidence angles ${\theta }_1$ and ${\theta }_2$ on the two tilt positions, and on the optical coating of the reflective surface of the micromirrors. Sticking to the specifications of commercially available DMDs [28], we have simulated the values of the $r_{1,2}^{S,P}$ coefficients for aluminum-coated metallic micromirrors with tilt angles $\pm 12°$, for varying wavelengths over the spectral bandwidth of commercially available DMDs (450–850 nm) and for varying incidence angle $\theta$. Of course, efficient polarimetric imagers based on focal plane arrays can be found over this (visible) spectral range, which the proposed technique could hardly compete with. However, the two-pixel polarimetric approach could be operated at any other wavelengths, with a similar study and using the same algorithms. The expression of the Fresnel’s intensity reflection coefficients on metals is readily obtained from optics textbooks [30], and with refractive index of aluminum tabulated in Ref. [31]. The evolution of parameters $r_{1,2}^{S,P}$ with $\theta$ is plotted in Fig. 2(a) for a wavelength of 780 nm, showing that the four reflection coefficients considered do not differ much (typically, less that 15% difference for reasonable incidence angles). As the unmixing step in the reconstruction problem consists basically of an “inversion” of matrix $\mathbf{A}$, the condition number $\kappa (\mathbf{A})={\Vert {\mathbf{A}}^{-1}\Vert }_2·{\Vert \mathbf{A}\Vert }_2$ naturally gives an indication of the difficulty of such an inversion procedure. It can be noted here that the aforementioned implementations of SPC setups with polarization analyzers [16,20–23] would mathematically correspond to the case of a mixing matrix with perfect conditioning. This would naturally yield easier reconstruction but at the expense of doubling the measurement time, decreasing the SNR due to the intensity losses incurred by the polarizers, and reducing spectral tolerance due to the presence of polarization analysis elements.

To assess the behavior of the condition number of our two-pixel setup with experimental parameters, we have plotted in Fig. 2(b) the evolution of $\kappa (\mathbf{A})$ as a function of $\theta$, which confirms that high incidence angles may be favored. We further analyzed the evolution of $\kappa (\mathbf{A})$ with $\theta$ and with the illumination wavelength. The contour plot in Fig. 2(c) shows that the condition number can vary by a factor of 5 across the range of wavelength and incidence angle considered, higher wavelengths being best adapted to maximize the polarimetric sensitivity of the reflection on aluminum mirrors. As a result, owing to the availability of laser and LED sources at such wavelength, and to keep reasonable incidence angles, we will consider throughout the remainder of this paper the situation of a monochromatic illumination at 780 nm with incidence angle $\theta =50°$ (black cross in Fig. 2), yielding a reasonably low condition number of $\kappa (\mathbf{A})=26$. It must be noted here that operating the same approach with a wide-band light source might lead to higher values of $\kappa (\mathbf{A})$, as the Fresnel’s coefficients of the mirrors might be averaged out over the considered light spectrum. It can be noted, however, that dielectric coated micromirrors could offer stronger polarization dependence of the reflection coefficients, but at the expense of totally revisiting the fabrication process of DMDs. For the sake of simplicity, we will neglect in this paper the influence of the anti-reflection coated cover slit that protects the DMD surface [28]. We also neglect the dispersion of the values of coefficients $r_{1,2}^{S,P}$ with the source spectral bandwidth and with the slight variation of incidence angle when the image of the scene is formed on the DMD surface. All these possible sources of imperfections can be neglected here to study the principle of the two-pixel polarimetric camera, but may be tackled in further work addressing its experimental validation.

Nevertheless, we shall analyze in our simulations how a bias in the estimated incidence angle $\theta$ can have an impact on reconstruction quality. Indeed, it is quite obvious that a bias in $\theta$ will lead to the use of an incorrect mixing matrix in the inversion process, hence hindering the recovery of a satisfactory polarization contrast image. Moreover, we also consider the influence of possible individual random errors in the tilt angles of each micromirror of the DMD. Indeed, the typical individual tilt error is about $\pm 1°$ to $\pm 2°$ according to standard DMD specifications [28], and thus its consequences on the reconstruction process may not be negligible. If global angle bias on $\theta$ and individual tilt errors can affect the reconstruction, their influence on image quality is likely to be very different. Global bias on $\theta$ can be basically treated as a calibration error, whereas individual tilt errors would rather behave as additional random noise in the inversion process.

Lastly, we will consider that the only source of noise is the photodetector electronic noise, and we assume statistical independence between the noise at photodetectors $P_1$ and $P_2$, and between noise realizations as the DMD patterns are varied. In the preceding inverse problem, noise vector $\mathbf{b}$ can thus be modeled by a centered Gaussian distribution of variance ${\sigma }^2$, i.e., $\forall i\in \{1,2\},\forall j\in [1,M],{\mathbf{b}}_{ij}\sim \mathcal{N}(0,{\sigma }^2)$.

3. RECONSTRUCTION ALGORITHMS

In this section, we describe different strategies to reconstruct the polarimetric data $\mathbf{X}$ from the compressed measurements $\mathbf{Y}=\mathbf{AX}\mathbf{\varphi }+\mathbf{b}$. In the framework of compressed sensing, the signal-recovery problem is generally described as a “large $p$, small $n$” problem, where the number of unknowns, i.e., the number of samples in the polarimetric components, can be much larger than the number of measurements in $\mathbf{Y}$. Hence, it is essential to enforce additional constraints on the signal to be recovered, which eventually boils down to solving a minimization problem of the form

In the context of CS [32,33], it is customary to enforce the sparsity of the unknown variable $\mathbf{X}$ in some signal representation $\mathbf{\Psi }$ chosen a priori. In the following sections, we describe and compare several strategies that are precisely dedicated to solving the two-pixel polarimetric compressed sensing recovery problem.

A. Two-Stage Reconstruction Approach

The recovery problem in Eq. (5) can be described as the combination of two classical inverse problems: a compressed sensing problem and a source separation problem. A first straightforward approach consists of tackling alternately both problems. Then, recovering the polarimetric components can be performed with the following two-stage approach.

Since $\mathbf{A}$ is invertible, the solution of this problem is given by $\widehat{\mathbf{X}}={\mathbf{A}}^{-1}\widehat{\mathbf{Y}°}$.

B. Combined Sparse Reconstruction

Despite its simplicity, this two-stage approach suffers from a major drawback: the mixed components $\mathbf{Y}°$ won’t be perfectly estimated, especially when only a few measurements in $\mathbf{Y}$ are available and when noise contaminates the data. These estimation errors will be amplified in the unmixing stage. Since the mixing matrix $\mathbf{A}$ is likely to be ill-conditioned, these errors will largely impact the accuracy of the reconstruction process.

A more effective strategy consists of jointly tackling both the compressed sensing recovery and unmixing problems. Extending standard reconstruction procedures yields the following optimization problem:

C. Constrained Sparse Reconstruction

Further improving the accuracy of the components’ recovery requires imposing additional, more physical, constraints on $\mathbf{X}$. In the context of polarimetric data, each component ${\mathbf{x}}_S$ and ${\mathbf{x}}_P$ has naturally non-negative values. Additionally, under active polarized illumination (as sketched in Fig. 1) and assuming purely depolarizing samples, the components must verify the inequality ${\mathbf{x}}_S⪰{\mathbf{x}}_P$. Indeed, most natural objects or scenes mostly depolarize the received light and, under such polarized illumination, lead to a decrease of ${\mathbf{x}}_S$ and an increase of ${\mathbf{x}}_P$, reaching ${\mathbf{x}}_S={\mathbf{x}}_P$ for a fully depolarizing object. In this section, we propose extending the reweighted FISTA to enforce these additional constraints. The optimization problem to tackle is as follows:

In contrast to the standard problem in Eq. (8), the problem in Eq. (9) is composed of a sum of convex penalizations that cannot be tackled with the FISTA. For that end, the generalized forward backward (GFB) algorithm [34] is the perfect algorithm to solve such an optimization problem. The application of the GFB to Eq. (9) is described in Algorithm 1, where $\gamma$ stands for the gradient path length used in the algorithm, and $L$ denotes the maximum iteration number.

In this algorithm, the application of each penalization or constraint is performed independently on distinct intermediate variables ${\{{\mathbf{U}}_i\}}_{i=1,\dots ,3}$. Updating each of these variables only requires the current estimate of the polarimetric components ${\mathbf{X}}^{(t)}$, the gradient of the quadratic data fidelity term $\mathbf{G}=-{\mathbf{A}}^T(\mathbf{Y}-{\mathbf{AX}}^{(t)}\mathbf{\varphi }){\mathbf{\varphi }}^T$, as well as the so-called proximal operator of the penalization or constraint.

The proximal operators required in the preceding algorithms are defined in Appendix B. The convergence of the GFB algorithm is guaranteed as long as the gradient path length verifies $\gamma <\frac{1}{{\Vert \mathbf{A}\Vert }_2^2}$, where the spectral norm of $\mathbf{A}$ is defined as its largest singular value. The scalar weights ${\{{\mu }_i\}}_{i=1,\dots ,3}$ must be strictly positive and have to add up to 1. Hence, the final update of the polarimetric components is a weighted average of the different intermediate variables. In the remainder of this paper, these weights will all be set equal to 1/3. The proposed GFB-based algorithm is initialized with the polarimetric signals provided by the reweighted FISTA described in the preceding section. Since problem 9 is convex, this initialization does not change the solution but dramatically reduces its computational cost.

Algorithm 1. Combined GFB

View Table | View all tables in this article

D. Optimization of Algorithm Parameters

The three preceding signal recovery approaches require the tuning of a certain number of parameters, whose setting is essential for an accurate estimation:

4. NUMERICAL RESULTS

In this section, we analyze the performance of the reconstruction algorithms and regularization procedures just described. These algorithms will also be compared, in terms of robustness, to some of the experimental imperfections mentioned in Section 2. This analysis will be conducted on a 1D test signal for the sake of computation speed. The quality of the reconstructed signals will be evaluated throughout this section by computing the peak signal to noise ratio (PSNR) of a concatenated vector of the reconstructed polarimetric components, i.e., $[{\widehat{\mathbf{x}}}_S{\widehat{\mathbf{x}}}_P]$. Then we will present some polarimetric imaging numerical results on 2D signals that demonstrate the ability of a two-pixel polarimetric camera to provide satisfactory compressed polarization contrast images at low cost and low complexity.

A. Comparison of Algorithm Performance and Robustness

1. Description of the 1D Test Signal

In order to optimize computation resources, we generated a synthetic polarimetric 1D data sample of length $N=512$ that will be used throughout this subsection to compare the performance of the just-given algorithms. The corresponding signals ${\mathbf{x}}_S$ and ${\mathbf{x}}_P$ are plotted in inset (a) of Fig. 3 with blue and red solid lines, respectively. The simulated polarimetric components verify the positivity constraint ${\mathbf{x}}_S⪰{\mathbf{x}}_P$, and it can be noted that their supports are not joint for the sake of generality. In inset (b) of Fig. 3, we plot a set of simulated detected intensity signals ${\mathbf{y}}_1$ and ${\mathbf{y}}_2$, with 40% compression rate ($M=306$) and $\mathrm{SNR}=40\mathrm{dB}$. The mixing matrix A used to generate the data has been calculated assuming an incidence angle $\theta =50°$ and wavelength of 780 nm (i.e., the optimal conditions identified in Section 2.C). The patterns used on the DMD to sample the image were generated from a randomly scrambled Hadamard transform of the signal $\mathbf{X}$ of size $N=512$. Unless otherwise specified, the same experimental conditions will be assumed for all numerical results presented in this paper. In Fig. 3, ${\mathbf{y}}_2$ is plotted as a function of ${\mathbf{y}}_1$, showing the strong anticorrelation existing between the two detected signals. Lastly, in inset (a) of Fig. 3, we show an example of reconstructed signals ${\widehat{\mathbf{x}}}_S$ and ${\widehat{\mathbf{x}}}_P$ obtained with the reweighted FISTA and yielding $\mathrm{PSNR}=37.7\mathrm{dB}$. As for all reconstructions of the 1D signal presented subsequently, the reconstruction process makes use of the unidimensional undecimated wavelet transform with the Haar filter, which is well suited to sparsely represent piecewise constant signals. It can be checked in inset (a) of Fig. 3 that the algorithm has not been constrained to ensure positivity of either ${\widehat{\mathbf{x}}}_S$, ${\widehat{\mathbf{x}}}_P$, or ${\widehat{\mathbf{x}}}_S-{\widehat{\mathbf{x}}}_P$.

2. Algorithm Performance

Using this 1D test signal, we are now able to compare the performance of the various algorithms as a function of the different parameters involved in the imaging process. Let us first analyze the influence of the data SNR on the reconstruction quality for an intermediate compression rate of 40%. The evolution of the PSNR of $\widehat{\mathbf{X}}$ as a function of the SNR is given in Fig. 4. All the data points in Fig. 4 have been obtained from 30 realizations of the numerical experiments, with the error bars indicating the interquartile range, i.e., distance between the first to third quartile of the data series.

 

Fig. 4. Evolution of the PSNR of the signal $\widehat{\mathbf{X}}$ reconstructed with the four compared algorithms as a function of detected signal SNR for a compression rate of 40%. The symbols represent the mean PSNR over 30 realizations, with error bars indicating the interquartile range. The lines are guides for the eyes.

Download Full Size | PDF

It can first be seen that all the algorithms asymptotically exhibit a linear evolution of their PSNR as a function of the SNR. Then, it is interesting to note that for intermediate values of SNR ($10\mathrm{dB}<\mathrm{SNR}<50\mathrm{dB}$), the two-step approach underperforms with respect to the simplest implementation of the combined reconstruction approach (denoted by combined FISTA). However, as soon as a reweighted procedure is implemented, solving the CS and the unmixing problems simultaneously (algorithm denoted as combined rFISTA) provides an asymptotical gain of about 12 dB in PSNR with respect to the two-step algorithm. Lastly, imposing physical positivity constraints on ${\widehat{\mathbf{x}}}_S$, ${\widehat{\mathbf{x}}}_P$, and ${\widehat{\mathbf{x}}}_S-{\widehat{\mathbf{x}}}_P$ through the implementation of the GFB algorithm (denoted as combined GFB) does not bring any additional gain in performance for highest values of SNR. However, in noisy situations, for $\mathrm{SNR}<50\mathrm{dB}$, the positivity constraints prove efficient at improving the reconstruction quality. A maximum gain of almost 10 dB is obtained for $\mathrm{SNR}=0\mathrm{dB}$, a situation where the unconstrained rFISTA approach fails to surpass the reconstruction quality obtained with the two-step procedure.

We now analyze the influence of the compression rate and incidence angle $\theta$ on the PSNR of the reconstructed signals. As discussed in Section 2.3, the incidence angle controls the condition number $\kappa (\mathbf{A})$ of the mixing matrix, and hence the difficulty of the unmixing problem. For this numerical experiment, we consider only the combined-rFISTA (plotted in blue solid line in the previous figure), with fixed number of iterations ($L={10}^4$) and with SNR of 60 dB. A 2D plot of the average PSNR obtained over 10 realizations of each numerical experiment is provided in Fig. 5 for 40 values of $\theta$ ranging between 25° and 65°, and 49 values of compression rate between 0% (no compression) to 96%. It can be observed that the reconstruction quality obeys a classical “phase transition” behavior frequently observed in CS problems, with three main domains which can be identified. First, for compression rates below 50%, the reconstruction is almost perfect (PSNR $\ge 55\mathrm{dB}$), whatever be the value of $\theta$. Only for the highest values of $\kappa (\mathbf{A})$ (i.e., for $\theta \le 30°$) does the algorithm fails at reaching a PSNR of 50 dB, but it remains above 45 dB. Then, for intermediate values of compression rate between 50% and 80%, the reconstruction quality gradually decreases while remaining above 35 dB. In that regime, the influence of the condition number $\kappa (\mathbf{A})$ appears clearly on the reconstruction quality. Lastly, a sharp transition occurs around 85% compression rate, above which no satisfactory reconstruction can be obtained regardless of $\theta$. This 2D map is rather encouraging towards possible experimental implementations of the two-pixel polarimetric CS camera. Indeed, it shows that, provided a good SNR can be ensured (here 60 dB), polarimetric signals can be retrieved at relatively high compression rates ($\le 80\%$) and for reasonable incidence angles on the DMD.

 

Fig. 5. 2D map of the PSNR (averaged over 10 realizations) of $\widehat{\mathbf{X}}$ reconstructed with the combined-reweighting FISTA, as a function of incidence angle $\theta$ and compression rate. The corresponding condition number $\kappa (\mathbf{A})$ is indicated as a function of $\theta$.

Download Full Size | PDF

3. Robustness to Incidence Angle Bias and Tilt Errors

Before presenting results on a polarimetric imaging scenario with 2D test signals, we report a last numerical experiment conducted on the 1D test signal to assess the robustness of the various algorithms to experimental imperfections, such as a bias on the incidence angle $\theta$, and random errors on the micromirror tilt angles.

The influence of a bias in $\theta$ has been analyzed as follows. Measurement vector $\mathbf{Y}$ was generated assuming a SNR of 80 dB and an incidence angle $\theta =50°$, as in Fig. 4. However, the reconstruction procedure was run with an incorrect mixing matrix $\mathbf{A}$, assuming a wrong incidence angle of $\theta +\delta \theta$. This way, we mimicked an experimental bias $\delta \theta$ comprised between 0.05° and 10° on the incidence angle. The reconstruction PSNR obtained with different algorithms is plotted in Fig. 6 in solid lines for a compression rate of 0%. It can be seen that the presence of a bias in reconstruction angle leads to a linear degradation of the PSNR in a log–log scale for all three algorithms tested. Even for a very small bias ($\delta \theta =0.05°$), a significant drop of more than 20 dB in reconstruction quality can be observed for all three algorithms with respect to unbiased reconstruction, but they still offer satisfactory recovery quality (PNSR $\ge 55\mathrm{dB}$). Applying physical constraints in the reconstruction procedure with the GFB algorithm seems to alleviate the degradation for any magnitude of $\delta \theta$, as far as simple angular bias is considered. However, the PSNR value of the reconstructed vector $\widehat{\mathbf{X}}$, which is used to gauge reconstruction quality, is not satisfactory here: it has indeed been observed on reconstructed signals that applying constraints with an erroneous mixing matrix often leads to singular results where the second polarimetric component ${\mathbf{x}}_P$ is forced to zero, thus yielding null polarimetric contrast. This is interpreted by the fact that the physical constraints applied can be no longer valid for the measured data with an incorrect matrix $\mathbf{A}$.

 

Fig. 6. Influence of the bias on incidence angle $\theta$ in the reconstruction quality (PSNR) for 0% compression (solid lines). PSNR of the reconstructed signal in the presence of bias and random tilt error on individual micromirrors (dashed line).

Download Full Size | PDF

Concerning the influence of random tilt errors on the micromirrors, we simulated this effect by computing individual mixing matrices A for all $N$ micromirrors simulated, assuming that the tilt angles $t_1$ and $t_2$ were affected by a random error. For the sake of computational speed, we assumed a uniform distribution over 11 values of the tilt error between $\pm 1°$. The PSNR of the reconstructed signal with angular bias and uniform tilt error is also plotted in Fig. 6 in dashed lines for 0% compression. On the one hand, the additional random tilts do not modify significantly the results for highest values of angular bias ($\delta \theta \ge 1°$), where the incorrect “average” matrix A is responsible for most of the quality degradation. On the other hand, for very low values of $\delta \theta$, the PSNR reaches a limit upper value of about 35–40 dB due to the presence of random tilt errors, which seems to affect more strongly the constrained version of the algorithm. Despite such degradation of reconstruction quality, these numerical experiments show that, with imperfect experimental configurations, using a combined recovery approach instead of a two-step reconstruction procedure can be advantageous. Correct PSNR values ($\ge 35\mathrm{dB}$) can be reached with the rFISTA or GFB algorithms with small angular biases ($\le 1°$) and in the presence of random tilt errors.

B. Numerical Polarimetric Imaging Results on 2D Signals

In this last section, we analyze the ability of the two-pixel polarimetric CS camera to retrieve polarimetric contrast images from simulated 2D data. Owing to its performance and simple implementation with respect to constrained GFB, we only consider in this section the reweighted FISTA implementing a combined reconstruction procedure. We first consider in the next subsection a simple imaging scenario to study the influence of polarimetric contrast on the reconstruction quality, before a more realistic example of polarimetric image reconstruction is given in Section 4.B.2.

1. Influence of Polarimetric Contrast

For this first imaging scenario, we consider a square object with high total intensity value over a dark background, forming an image of $N=128\times 128$ pixels. This intensity image ${\mathbf{x}}_T$, plotted in Fig. 7, is supposed to be the sum of two polarimetric image components ${\mathbf{x}}_S$ and ${\mathbf{x}}_P$, yielding a true OSC map also plotted in Fig. 7. In this scenario, we assume that a first object (smaller square) cannot be distinguished from the second object (bigger square) on the intensity image ${\mathbf{x}}_T$, while the OSC map allows the two objects to be clearly identified. The background is assumed to be totally depolarizing ($\mathrm{OSC}=0$). The smaller square object is always supposed to be slightly depolarizing ($\mathrm{OSC}=0.8$), whereas the OSC of the second object (bigger square) is varied between 0 and 1 in the following numerical experiments. Figure 7 shows an example of reconstruction of the total intensity and the OSC map with the combined-rFISTA (with Haar wavelets) for $\mathrm{SNR}=40\mathrm{dB}$, compression rate 40%, incidence angle $\theta =50°$, and no bias or uncertainty on the tilt directions. The reconstruction quality is visually very good ($\mathrm{PSNR}=52.8\mathrm{dB}$), as evidenced by the reconstruction error map of the total intensity and OSC shown in Fig. 7. This simple example demonstrates the possibility of providing relevant polarimetric information from the proposed two-pixel polarimetric CS camera concept.

 

Fig. 7. Simulated polarimetric imaging scenario containing two objects over constant background. First row: total intensity image ${\mathbf{x}}_T$ (true), ${\widehat{\mathbf{x}}}_T$ (reconstruction), and reconstruction error map. The two objects are indistinguishable in the total intensity image. Second row: true and reconstructed OSC image, and reconstruction error map.

Download Full Size | PDF

In the first row (a) of Fig. 8, we plotted the reconstructed OSC maps for three other values of the OSC of the second object (bigger square). The reconstruction error map is given in the second row (b), while the evolution of the PSNR with the OSC of the second object is plotted in Fig. 8(c). These results show the fact that the reconstruction quality naturally decreases when the polarimetric contrast between the two objects is reduced, i.e., when the reconstruction problem becomes more difficult. This can be easily understood as the smallest variations of contrast are likely to be buried in the noise and totally filtered out by the regularization process. This is clearly seen for $\mathrm{OSC}=0.6$, where the sharpest features of the first object disappear in the reconstructed OSC map. Obviously, reconstruction quality is maximized when the OSC of the second object reaches 0.8, i.e., a single object is to be identified in the image, thus yielding a simpler reconstruction problem.

 

Fig. 8. (a) Reconstructed OSC image and (b) reconstruction error map for various values of the OSC of the second object. (c) Evolution of the reconstruction PSNR of $\mathbf{X}$ when the OSC of the second object is varied between 0 and 1.

Download Full Size | PDF

2. Example of Reconstruction on Realistic Image Data

Lastly, we present an example of reconstructed polarimetric image on a more realistic imaging scenario, considering an approximately sparse natural image. For this purpose, we consider a true intensity image ${\mathbf{x}}_T$ of the cameraman with size $N=512\times 512$, as plotted in Fig. 9. Appropriate polarimetric components ${\mathbf{x}}_S$ and ${\mathbf{x}}_P$ were generated so that a true OSC map would reveal four hidden objects (three in the grass, one in the buildings) over a depolarizing background, as can be seen in Fig. 9. The reconstruction results with Symmlet wavelet transform and a compression rate of 60% are also displayed in Fig. 9, along with reconstruction error maps. The total intensity image is almost perfectly reconstructed, as would be the case with a SPC imaging system. However, the four hidden objects remain of course invisible in the reconstructed image ${\widehat{\mathbf{x}}}_T$. Contrarily, the reconstructed OSC map makes it possible to identify the presence of the four hidden objects by revealing their polarimetric contrast over the background. The analysis of the reconstruction error map of OSC shows that the polarimetric information about the four hidden objects is fairly retrieved. However one can notice significant reconstruction errors in the darkest regions of the image (cameraman and tripod), where the impact of reconstruction errors on the polarimetric components ${\mathbf{x}}_S$ and ${\mathbf{x}}_S$ are strongly amplified by the computation of the OSC (division by very low intensity values). These imperfections could be lowered in the future by refining regularization parameters and constraints in the algorithm implementation, or in further algorithmic developments on this subject, by applying reconstruction constraints directly on the OSC map itself.

 

Fig. 9. Example of numerical polarimetric imaging experiment on realistic $512\times 512$ pixel image data. First row: total intensity image ${\mathbf{x}}_T$, reconstruction ${\widehat{\mathbf{x}}}_T$, and error map with 60% compression rate. The total intensity image is the standard cameraman test image, where the four hidden objects are totally hidden. Second row: true OSC image, reconstructed OSC image, and error map. The four hidden objects are revealed by their polarimetric contrast over depolarizing background on the OSC images.

Download Full Size | PDF

5. CONCLUSION AND PERSPECTIVES

In this paper, we have proposed a new concept of CS polarimetric imaging inspired by the SPC principle. Relying on the tiny differences in reflection coefficients of mirrors with incidence angle and polarization direction, the setup proposed allows intensity and polarimetric contrast information to be recovered from the temporal acquisition on two single-pixel detectors, and without requiring any polarization analysis optical component. We have shown that this recovery problem could be analyzed as joint CS and source separation tasks that can be tackled independently and successively using standard approaches (respectively ${\ell }_1$ minimization and direct matrix inversion), or optimally treated in an original combined reconstruction approach. For this purpose, we have presented different versions of a combined algorithm, including FISTA implementation of the combined optimization problem, with potential reweighted steps that have proved efficient at increasing the reconstruction quality. Moreover, to enforce additional physical constraints on the measured data, we have proposed a constrained sparse reconstruction method based on the GFB algorithm, allowing the reconstruction quality to be improved in low-SNR conditions.

Numerical simulations have permitted us to validate these approaches and to analyze the influence of experimental conditions such as incidence angle, SNR, micromirror tilt errors, etc. Generally speaking, these results are encouraging towards the experimental implementation of such an imaging setup, as good reconstruction quality was obtained for reasonable incidence angles on the mirror, even in the presence of a small bias or randomness in the mirror tilt direction. Significant compression rates could be achieved while still offering sufficient reconstruction quality, as illustrated in the 2D simulation results presented.

Experimental validation of this computational imaging approach appears as a natural perspective to this work. Such a proof of concept will require us to clear a number of experimental issues with regard to the previous SPC implementations, such as limited measurement dynamics, diffraction on the micromirror edges that could affect the intensity balance between the two photodetectors, and relatively important incidence angles that can cause vignetting and foreshortening effects. Another interesting research track is the design of a blind calibration method so as to compensate for possible errors in the estimation of the mixing matrix $\mathbf{A}$ involved in the reconstruction process. More generally, applying CS concepts to more sophisticated multi-channel polarimetric imaging techniques which are characterized by very specific algebraic constraints is likely to raise interesting reconstruction issues.

APPENDIX A: FORWARD–BACKWARD SPLITTING ALGORITHM

Whether in the two-step or the combined sparse reconstruction approach, recovering the polarimetric signal $\mathbf{X}$ requires solving optimization problems of the form

$$\underset{\mathbf{X}}{\min }\mathcal{P}(\mathbf{X})+\mathcal{L}(\mathbf{X}),$$

where the data fidelity term $\mathcal{L}(\mathbf{X})$ is a quadratic norm that is differentiable and whose gradient is Lipschitz with some constant $L_i$, and $\mathcal{P}(\mathbf{X})$ is a convex but non-smooth penalization. The minimization problem can be carried out efficiently with recent proximal algorithms [37], and more precisely with the forward–backward splitting (FBS) algorithm [38]. The FBS algorithm is an iterative procedure that can be described with the following update rule at iteration $t$:

$${\mathbf{X}}^{(t+1)}={\mathrm{prox}}_{\gamma \mathcal{P}}({\mathbf{X}}^{(t)}-\gamma \nabla \mathcal{L}({\mathbf{X}}^{(t)})),$$

where $\nabla \mathcal{L}({\mathbf{X}}^{(t)})$ stands for the derivative of $\mathcal{L}$ in ${\mathbf{X}}^{(t)}$ and $\gamma$ is the gradient path length. The proximal operator ${\mathrm{prox}}_{\gamma \mathcal{P}}$ is defined by the solution to the problem

$${\mathrm{prox}}_{\gamma \mathcal{P}}(\mathbf{z})={\mathrm{Argmin}}_{\mathbf{X}}\gamma \mathcal{P}(\mathbf{X})+\frac{1}{2}{\Vert \mathbf{z}-\mathbf{X}\Vert }_F^2.$$

While the proximal operator of the penalization function is defined as the solution of an optimization problem, the standard proximal operators admit a closed-form expression (see Appendix B).

In this paper, the term $\mathcal{P}$ will penalize non-sparse solutions and will be based on the ${\ell }_1$ norm. Introduced in Ref. [36], the re-weighted ${\ell }_1$ norm further introduces weights $\mathbf{\Lambda }$ that aim at reducing the bias induced by the standard ${\ell }_1$ norm. Consequently, the penalization term used in this paper will take the generic form

$$\mathcal{P}(\mathbf{X})={\Vert \mathbf{\Lambda }\odot (\mathbf{X}{\mathbf{\Psi }}^T)\Vert }_{{\ell }_1},$$

where $\mathbf{\Psi }$ is the signal representation where sparsity is modeled. The final optimization procedure then alternates between updates of $\mathbf{X}$ for fixed weights $\mathbf{\Lambda }$ and updates of these weights, which has been showed to dramatically improve the accuracy of the reconstruction.

In the proposed reweighted algorithm, updates of the polarimetric signals are carried out using an accelerated version of the FBS algorithm coined FISTA [39]. The generic description of the reweighted FISTA is given in Algorithm 2.

This procedure generally increases the accuracy of the reconstruction process with few updates of the weights $\mathbf{\Lambda }$, typically $L_{\text{outer}}$ between 3 and 5. Details about practical parameter tuning are given in Section 3.D.

Algorithm 2. Reweighted FISTA

View Table | View all tables in this article

APPENDIX B: USEFUL PROXIMAL OPERATORS

Hereafter, we describe different proximal operators that are used in the proposed reconstruction algorithms.

Reweighted ${\ell }_1$

Assuming that the sparse signal representation $\mathbf{\Psi }$ is an orthogonal matrix, the proximal operator of ${\Vert \mathbf{\Lambda }\odot (\mathbf{X}{\mathbf{\Psi }}^T)\Vert }_{{\ell }_1}$ is defined as

(B1)$${\mathrm{prox}}_{{\Vert \mathbf{\Lambda }\odot (.{\mathbf{\Psi }}^T)\Vert }_{{\ell }_1}}(\mathbf{Z})=\mathbf{\Psi }{\mathcal{S}}_{\mathbf{\Lambda }}(\mathbf{Z}{\mathbf{\Psi }}^T),$$

where the weighted soft-thresholding operator ${\mathcal{S}}_{\mathbf{\Lambda }}$ is defined as

(B2)$$\begin{gathered} \forall i,j;{\mathcal{S}}_{\mathbf{\Lambda }}({\mathbf{Z}}_{ij}) \\ =\{\begin{array}{ll}{\mathbf{Z}}_{ij}-{\mathbf{\Lambda }}_{ij}\mathrm{sign}({\mathbf{Z}}_{ij})& \mathrm{if}|{\mathbf{Z}}_{ij}|>{\mathbf{\Lambda }}_{ij}\\ 0& \mathrm{otherwise}\end{array}. \end{gathered}$$

Positivity Constraint

The proximal operator of the positivity constraint is defined as the orthogonal projector onto the non-negative orthant:

(B3)$$\begin{gathered} \forall i,j;{\mathrm{prox}}_{i_{\mathbf{X}⪰0}(·)}({\mathbf{Z}}_{ij}) \\ =\{\begin{array}{cc}{\mathbf{Z}}_{ij}& \mathrm{if}{\mathbf{Z}}_{ij}>0\\ 0& \mathrm{otherwise}\end{array}. \end{gathered}$$

Inequality Constraint

The inequality constraint $x_s⪰x_p$ can be recast as $\mathbf{DX}⪰0$, where $\mathbf{D}=[1,-1]$. Its proximal operator is defined as the orthogonal projector onto the convex set $\{\mathbf{Z};\forall i,j,{[\mathbf{DZ}]}_{ij}\ge 0\}$, reading

$${\mathrm{prox}}_{i_{\mathbf{D}·⪰0}}(\mathbf{Z})={\mathrm{Argmin}}_{\mathbf{DX}⪰0}\frac{1}{2}{\Vert \mathbf{Z}-\mathbf{X}\Vert }_F^2.$$

This expression can be more conveniently recast in a Lagrangian formulation by introducing the Lagrange multipliers $\mathit{\pi }$

$${\mathrm{prox}}_{i_{\mathbf{D}·⪰0}}(\mathbf{Z})={\mathrm{Argmin}}_{\mathbf{X}}⟨\mathbf{DX},\mathit{\pi }⟩+\frac{1}{2}{\Vert \mathbf{Z}-\mathbf{X}\Vert }_F^2.$$

Its optimum is obtained for $\mathbf{X}=\mathbf{Z}-{\mathbf{D}}^T\mathit{\pi }$, which should verify the constraint $\mathbf{DX}⪰0$. This entails

$$\mathbf{DZ}-2\mathit{\pi }⪰0.$$

Consequently, the Lagrange multipliers must take the values

(B4)$$\forall j;{\mathit{\pi }}_j=\{\begin{array}{ll}\frac{1}{2}{[\mathbf{DZ}]}_j& \mathrm{if}{[\mathbf{DZ}]}_j<0\\ 0& \mathrm{otherwise}\end{array}.$$

From this expression of the Lagrange multipliers, the proximal operator is then defined as

(B5)$${\mathrm{prox}}_{i_{\mathbf{D}·⪰0}}(\mathbf{Z})=\mathbf{Z}-{\mathbf{D}}^T\mathit{\pi }.$$

Funding

H2020 European Research Council (ERC) (LENA Grant 678282).

REFERENCES

1. W. Groner, J. W. Winkelman, A. G. Harris, C. Ince, G. J. Bouma, K. Messmer, and R. G. Nadeau, “Orthogonal polarization spectral imaging: a new method for study of the microcirculation,” Nat. Med. 5, 1209–1212 (1999). [CrossRef]  

2. F. Boulvert, B. Boulbry, G. Le Brun, B. Le Jeune, S. Rivet, and J. Cariou, “Analysis of the depolarizing properties of irradiated pig skin,” J. Opt. A 7, 21–28 (2005). [CrossRef]  

3. M. Anastasiadou, A. D. Martino, D. Clement, F. Liège, B. Laude-Boulesteix, N. Quang, J. Dreyfuss, B. Huynh, A. Nazac, L. Schwartz, and H. Cohen, “Polarimetric imaging for the diagnosis of cervical cancer,” Phys. Stat. Solidi C 5, 1423–1426 (2008). [CrossRef]  

4. S. Breugnot and P. Clémenceau, “Modeling and performances of a polarization active imager at lambda = 806 nm,” Proc. SPIE 3707, 449–460 (1999). [CrossRef]  

5. M. Alouini, F. Goudail, A. Grisard, J. Bourderionnet, D. Dolfi, A. Bénière, I. Baarstad, T. Løke, P. Kaspersen, X. Normandin, and G. Berginc, “Near-infrared active polarimetric and multispectral laboratory demonstrator for target detection,” Appl. Opt. 48, 1610–1618 (2009). [CrossRef]  

6. J. Fade, S. Panigrahi, A. Carré, L. Frein, C. Hamel, F. Bretenaker, H. Ramachandran, and M. Alouini, “Long-range polarimetric imaging through fog,” Appl. Opt. 53, 3854–3865 (2014). [CrossRef]  

7. F. Snik, J. Craven-Jones, M. Escuti, S. Fineschi, D. Harrington, A. De Martino, D. Mawet, J. Riedi, and J. S. Tyo, “An overview of polarimetric sensing techniques and technology with applications to different research fields,” Proc. SPIE 9099, 90990B (2014). [CrossRef]  

8. C. Brosseau, Fundamentals of Polarized Light—A Statistical Optics Approach (Wiley, 1998).

9. N. Gupta, R. Dahmani, and S. Choy, “Acousto-optic tunable filter based visible- to near-infrared spectropolarimetric imager,” Opt. Eng. 41, 1033–1038 (2002). [CrossRef]  

10. G. Anna, H. Sauer, F. Goudail, and D. Dolfi, “Fully tunable active polarization imager for contrast enhancement and partial polarimetry,” Appl. Opt. 51, 5302–5309 (2012). [CrossRef]  

11. A. Bénière, M. Alouini, F. Goudail, and D. Dolfi, “Design and experimental validation of a snapshot polarization contrast imager,” Appl. Opt. 48, 5764–5773 (2009). [CrossRef]  

12. F. Goudail and P. Réfrégier, Statistical Image Processing Techniques for Noisy Images: an Application Oriented Approach (Kluwer, 2004).

13. R. G. Baraniuk, “Single-pixel imaging via compressive sampling,” IEEE Signal Process. Mag. 25(2), 83–91 (2008). [CrossRef]  

14. W. L. Chan, K. Charan, D. Takhar, K. F. Kelly, R. G. Baraniuk, and D. M. Mittleman, “A single-pixel terahertz imaging system based on compressed sensing,” Appl. Phys. Lett. 93, 121105 (2008). [CrossRef]  

15. A. Wagadarikar, R. John, R. Willett, and D. Brady, “Single disperser design for coded aperture snapshot spectral imaging,” Appl. Opt. 47, B44–B51 (2008). [CrossRef]  

16. A. Asensio Ramos and A. López Ariste, “Compressive sensing for spectroscopy and polarimetry,” arXiv:0909.4439 (2009).

17. V. Studer, J. Bobin, M. Chahid, H. S. Mousavi, E. Candes, and M. Dahan, “Compressive fluorescence microscopy for biological and hyperspectral imaging,” Proc. Natl. Acad. Sci. U.S.A. 109, E1679–E1687 (2012). [CrossRef]  

18. Y. August and A. Stern, “Compressive sensing spectrometry based on liquid crystal devices,” Opt. Lett. 38, 4996–4999 (2013). [CrossRef]  

19. A. Liutkus, D. Martina, S. Popoff, G. Chardon, O. Katz, G. Lerosey, S. Gigan, L. Daudet, and I. Carron, “Imaging with nature: a universal analog compressive imager using a multiply scattering medium,” arXiv:1309.0425 (2013).

20. V. Durán, P. Clemente, M. Fernández-Alonso, E. Tajahuerce, and J. Lancis, “Single-pixel polarimetric imaging,” Opt. Lett. 37, 824–826 (2012). [CrossRef]  

21. F. Soldevila, E. Irles, V. Durán, P. Clemente, M. Fernández-Alonso, E. Tajahuerce, and J. Lancis, “Single-pixel polarimetric imaging spectrometer by compressive sensing,” Appl. Phys. B 113, 551–558 (2013). [CrossRef]  

22. S. S. Welsh, M. P. Edgar, R. Bowman, B. Sun, and M. J. Padgett, “Near video-rate linear stokes imaging with single-pixel detectors,” J. Opt. 17, 025705 (2015). [CrossRef]  

23. C. Fu, H. Arguello, B. M. Sadler, and G. R. Arce, “Compressive spectral polarization imaging by a pixelized polarizer and colored patterned detector,” J. Opt. Soc. Am. A 32, 2178–2188 (2015). [CrossRef]  

24. F. Soldevila, P. Clemente, E. Tajahuerce, N. Uribe-Patarroyo, P. Andrés, and J. Lancis, “Computational imaging with a balanced detector,” Sci. Rep. 6, 29181 (2016). [CrossRef]  

25. J. Liang, C. Ma, L. Zhu, Y. Chen, L. Gao, and L. V. Wang, “Single-shot real-time video recording of a photonic mach cone induced by a scattered light pulse,” Sci. Adv. 3, e1601814 (2017). [CrossRef]  

26. N. A. Riza, J. P. La Torre, and M. J. Amin, “CAOS–CMOS camera,” Opt. Express 24, 13444–13458 (2016). [CrossRef]  

27. Y. C. Eldar and G. Kutyniok, Compressed Sensing: Theory and Applications (Cambridge University, 2012).

28. Texas Instruments, “DLP6500-0.65-1080p datasheet,” Rev. B (2016).

29. I. N. Flamarique, C. W. Hawryshyn, and F. I. Hárosi, “Double-cone internal reflection as a basis for polarization detection in fish,” J. Opt. Soc. Am. A 15, 349–358 (1998). [CrossRef]  

30. M. Born and E. Wolf, Principles of Optics, 7th ed. (Cambridge University, 1999), pp. 752–754.

31. A. D. Rakić, A. B. Djurišić, J. M. Elazar, and M. L. Majewski, “Optical properties of metallic films for vertical-cavity optoelectronic devices,” Appl. Opt. 37, 5271–5283 (1998). [CrossRef]  

32. E. Candès and M. Wakin, “People hearing without listening: an introduction to compressive sampling,” IEEE Signal Process. Mag. 25(2), 21–30 (2008). [CrossRef]  

33. E. Candès, “Compressive sampling,” in International Congress of Mathematics, Madrid, Spain (2006).

34. H. Raguet, J. Fadili, and G. Peyre, “A generalized forward-backward splitting,” SIAM J. Imaging Sci. 6, 1199–1226 (2013). [CrossRef]  

35. J.-L. Starck, J. Fadili, and F. Murtagh, “The undecimated wavelet decomposition and its reconstruction,” IEEE Trans. Image Process. 16, 297–309 (2007). [CrossRef]  

36. E. J. Candes, M. B. Wakin, and S. P. Boyd, “Enhancing sparsity by reweighted L1 minimization,” J. Fourier Anal. Appl. 14, 877–905 (2008). [CrossRef]  

37. N. Parikh and S. Boyd, “Proximal algorithms,” Found. Trends Optim. 1, 123–231 (2014).

38. P. L. Combettes and V. R. Wajs, “Signal recovery by proximal forward-backward splitting,” Multiscale Model. Simul. 4, 1168–1200 (2005). [CrossRef]  

39. A. Beck and M. Teboulle, “Fast iterative shrinkage-thresholding algorithm for linear inverse problems,” SIAM J. Imaging Sci. 2, 183–202 (2009). [CrossRef]