# Two-pixel polarimetric camera by compressive sensing

## Abstract

We propose an original concept of compressive sensing (CS) polarimetric imaging based on a digital micro-mirror (DMD) array and two single-pixel detectors. The polarimetric sensitivity of the proposed setup is due to an experimental imperfection of reflecting mirrors which is exploited here to form an original reconstruction problem, including a CS problem and a source separation task. We show that a two-step approach tackling each problem successively is outperformed by a dedicated combined reconstruction method, which is explicited in this article and preferably implemented through a reweighted FISTA algorithm. The combined reconstruction approach is then further improved by including physical constraints specific to the polarimetric imaging context considered, which are implemented in an original constrained GFB algorithm. Numerical simulations demonstrate the efficiency of the 2-pixel CS polarimetric imaging setup to retrieve polarimetric contrast data with significant compression rate and good reconstruction quality. The influence of experimental imperfections of the DMD are also analyzed through numerical simulations, and 2D polarimetric imaging reconstruction results are finally presented.

## 1Introduction

In various application domains such as biomedical diagnosis, defence or remote sensing, standard intensity imaging techniques sometimes fail to reveal relevant contrasts or to gather sufficient information. In these domains, polarimetric approaches have proved efficient to enhance the estimation or detection capabilities of the imaging systems [1]. Mostly often, the polarimetric information is provided by a scalar polarization contrast image, which offers complementary constrast 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 2-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 between two polarimetric images ( and ) of the scene, acquired through a linear polarizer oriented along two orthogonal directions, denoted and throughout the article where denotes the polarization direction of the illumination source. Owing to its instrumental simplicity, which can be further improved using voltage-controlled electro-optics devices [9] 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]), 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 these past years [13]. 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], THz imaging [14], or random media-assisted CS [19]. However, despite the swarming interest in CS, only few attempts were reported so far to perform polarimetric CS imaging [16]. The imaging setups proposed in these references 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 references, the reconstruction process consisted of solving as many CS reconstruction problems as polarimetric channels were considered (2 or 4). More precisely, in [20], 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 the loss in intensity due to the use of a polarization analyzer. In references [22], 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 article, we revisit the problem of 2-channel polarimetric CS by proposing an original polarimetric imaging architecture using two single-pixel detectors. The proposed setup is still inspired from the initial concept of SPC, but does not require any polarization analyzing element as it relies on imperfections of the DMD itself. Contrarily to previous attempts in polarimetric CS, the polarimetric information is obtained through a single temporal data acquisition 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.

This article 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 2-pixel polarimetric CS imaging approach are finally discussed in Section 4 through numerical simulations on 1D and 2D test signals. In this section, the influence of experimental imperfections on the reconstruction quality is also analyzed, before conclusions and perspectives of the article are provided in Section Section 5.

## 2Principle of 2-pixel CS polarimetric imaging

### 2.12-pixel CS polarimetric imaging setup

The polarimetric CS imaging approach proposed relies on the SPC imaging architecture. As will be demonstrated, it 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. ?. In the context of active imaging, we assume that the scene or object of interest is enlightened by an horizontally polarized light source. An image of the scene is formed through a lens onto the surface of the DMD, which spatially samples it by applying a controlled binary pattern on the micromirrors. The detection setup is strictly similar to the SPC architecture, with a first photodetector () used to detect light reflected in the first reflection direction through a lens (). However, the light reflected in the second direction of the tilted mirrors is directed towards a second photodetector () via a lens (), 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 and provides two temporal signals that are sampled and digitized on a Analog to Digital Converter (ADC), before they can be used for numerical inversion of the intensity and polarimetric contrast images.

Throughout this article, we will denote the total intensity image of the scene by a single dimensional row vector of length . We assume that the total intensity image can be written as the sum of two polarimetric components, i.e., . Subscripts and 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 (), as sketched in Fig. ?. We assume that these two components are compressible in the same sparse representation , i.e., they can be written as , where is a matrix with , and are column vectors containing the expansion coefficients. In the compressed sensing framework, compressibility means that most of these expansion coefficients have a 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 [24].

Upon reflection on the DMD surface, the polarimetric components of the original image are altered by the Fresnel’s reflection coefficients (in intensity) of the mirror, depending on the reflection angle. Let us denote by and the Fresnel’s reflection coefficients (in intensity) of the mirror for each tilt direction and respectively for the and polarimetric components of the image formed on the DMD surface. Neglecting absorption effects, these coefficients are real and verify . With such notations, the images respectively reflected towards detectors and read and , which can be rewritten in a compact form as

When a given pattern indexed by is applied on the DMD to spatially sample the image, the total intensity reflected towards direction 1 and integrated on detector can be denoted by , where is a binary valued -dimensional column vector encoding the set of orientations of the individual mirrors (DMD pattern). Similarly, , where is the complement of vector . Then, when measurements are accumulated with various sets of pseudo-random configurations of the DMD, the detected intensities organized in a measurement matrix read

with sensing matrix and .

Under the above assumptions, we will show that such a simple setup suffices to retrieve a compressive measurement of the polarimetric components and , and thus of the intensity image , and OSC image . 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 and -polarized components of light as a function of the angle of incidence of the incoming light beam. As sketched in Fig. ?, the angle of incidence on the mirror is denoted by (respectively ) for tilt direction (respectively direction ), where is the angle of incidence with respect to the DMD surface, and where and denote the tilt angle of the mirrors (typically and on most DMDs [25]). This dependency with is illustrated in Fig. ?.a, where we plotted the evolution of the four reflection coefficients at wavelength 780 nm for aluminum mirrors when varies from to . The influence of the value of these coefficients on the CS reconstruction problem is discussed below in Section 2.3 and their calculation is recalled in Appendix Section 7. 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[26].

### 2.2Description of the CS inverse problem

With the above notations, we will now describe this imaging process as a CS inverse problem. We first rewrite the binary sensing matrix as , and consequently, , where the elemets of take on and values. We impose that, for all DMD configurations, 50 of the micromirrors are oriented towards direction 1 such that . In such a case, the measurement matrix can then be rewritten as , with , and

For the sake of simplicity, we will assume in the following that the constant term can be easily estimated and subtracted out from the measured data, e.g. by averaging the acquisitions over all DMD pattern realizations considered. In this case, the CS problem that must be solved is given in Eq. (Equation 1) and reads . For the sake of clarity, and taking into account an additive noise contribution on the detected intensities, we propose to rewrite it with simplified notations as

We also introduce , such that .

As a result, the -dimensional polarimetric components of the image contained in can be in principle recovered from a number of measurements provided the problem described in the above 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 . Indeed, as illustrated in Fig. ?.a, the reflection coefficients on metals are usually quite similar for polarization directions and , causing a strong crosstalk between the two components of interest. As a consequence, the signals detected at photodetectors and are almost perfectly anticorrelated, the polarimetric information lying in the tiny discrepancies between these two signals. This is illustrated in Fig. ?, 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.

### 2.3Experimental parameters and imperfections

Before we detail the reconstruction algorithms used to achieve polarimetric CS imaging with the 2-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 , which depends on the wavelength and bandwidth of the illuminating source, on the incidence angles and 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 [25], we have simulated the values of the coefficients for aluminum-coated metallic micromirrors with tilt angles , for varying wavelengths over the spectral bandwidth of commercially available DMDs (450-850 nm), and for varying incidence angle . The expression of the Fresnel’s intensity reflection coefficients on metals is recalled in Appendix Section 7. The evolution of parameters with is plotted in Fig. ?.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). The unmixing step in the reconstruction problem consisting basically of an “inversion” of matrix , the condition number naturally gives an indication about the difficulty of such inversion procedure. We have thus plotted in Fig. ?.b the evolution of as a function of , which confirms that high incidence angles may be favored. We further analyzed the evolution of with and with the illumination wavelength. The contour plot in Fig. ?.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 article the situation of a monochromatic illumination at 780 nm, with incidence angle (black cross in Fig. ?), yielding a reasonably low condition number of .

It can be noted here that dielectric coated micromirrors could offer stronger polarization dependence of the reflection coefficients, but to the expense of totally revisiting the fabrication process of DMDs. For the sake of simplicity, we will neglect in this article the influence of the anti-reflection coated cover slit that protects the DMD surface [25]. We also neglect the dispersion of the values of coefficients 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 2-pixel polarimetric camera, but may be addressed in further work to achieve the experimental validation of this CS imaging scheme.

Nevertheless, we shall analyze in our simulations how a bias in the estimated incidence angle can have an impact on reconstruction quality. Indeed, it is quite obvious that a bias on will lead to using 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 to according to standard DMD’s specifications [25], thus its consequences on the reconstruction process may not be negligible. If global angle bias on and individual tilt errors can affect the reconstruction, their influence on image quality is likely to be very different. Global bias on can be basically treated as a calibration error, whereas individual tilt errors would rather behave as additional random noise on the inversion process.

Lastly, we will consider that the only source of noise is the photodetectors electronic noise, and we assume statistical independence between the noise at photodetectors and , and between noise realizations as the DMD patterns are varied. In the above inverse problem, noise vector can thus be modeled by a centered Gaussian distribution of variance , i.e., .

## 3Reconstruction algorithms

In this section, we describe different strategies to reconstruct the polarimetric data from the compressed measurements . In the framework of compressed sensing, the signal recovery problem is generally described as a “large , small ” problem, where the number of unknown, i.e. the number of samples in the polarimetric components, can be much larger than the number of measurements in . 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

where the first term is a penalization term that favors solutions with certain desired properties, and the second term is a data fidelity term that measures the discrepancy between the data and the model . The Frobenius norm is defined as .

In the context of CS [27], it is customary to enforce the sparsity of the unknown variable in some signal representation that is chosen a priori. In the following sections, we describe and compare several strategies that are precisely dedicated to solve the two-pixel polarimetric compressed sensing recovery problem.

### 3.12-stage reconstruction approach

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

**Compressed sensing:**denoting the non-compressed mixed polarimetric components by , the actual measurements can be defined as . Recovering then boils down to a standard compressed sensing recovery problem. This step is customarily solved by finding the minimum of the problem

where the -norm enforces the sparsity of in and stands for the regularization parameters, which is composed of strictly positive entries (see Section 3.4 for details about the parameters’ selection). This optimization problem is solved using the reweighting FISTA algorithm. This algorithm is referred to as Algorithm ? and is detailed in Appendix Section 8.

**Source separation:**once the mixed polarimetric components are estimated, retrieving from is equivalent to a source separation or unmixing problem, which can be tackled by minimizing the Euclidean distance between and the model as followsSince is invertible, the solution of this problem is given by .

### 3.2Combined sparse reconstruction

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

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

where the matrix stands for the same weight matrix we introduced in the 2-stage approach (see Section 3.4). In the next, we make use of the reweighted FISTA algorithm (see Appendix Section 8) to solve .

### 3.3Constrained sparse reconstruction

Further improving the accuracy of the components’ recovery requires imposing additional, more physical, constraints on . In the context of polarimetric data, each component and has naturally non-negative samples. As well, under active polarized illumination, and assuming purely depolarizing samples, the components must verify the following inequality: . In this section, we propose extending the reweighted FISTA algorithm to enforce these additional constraints. The optimization problem to tackle is described as follows

where stands for the characteristic function of the positive orthant and for the characteristic function of the convex set where .

In contrast to the standard problem in Eq. (Equation 6), the problem in Eq. is composed of a sum of convex penalizations that cannot be tackled with the FISTA algorithm. For that end, the Generalized Forward Backward (GFB) algorithm [29] is the perfect algorithm to solve such an optimization problem. The application of the GFB to Eq. is described in Alg. ?, where stands for the gradient path length used in the algorithm, and denotes the maximum iteration number.

In this algorithm, the application of each penalization or constraint is performed independently on distinct intermediate variables . Updating each of these variables only requires the current estimate of the polarimetric components , the gradient of the quadratic data fidelity term as well as the so-called proximal operator of the penalization or constraint. The proximal operators required in the above algorithms are defined in Appendix Section 8.

The convergence of the GFB algorithm is guaranteed as long as the gradient path length verifies , where the spectral norm of is defined as its largest singular value. The scalar weights must be strictly positive and have to sum to . 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 . The proposed GFB-based algorithm is initialized with the polarimetric signals provided by the reweighted FISTA algorithm described in the above section. Since the problem Equation 7 is convex, this initialization does not change the solution but it dramatically reduces its computational cost.

### 3.4Optimization of algorithm parameters

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

**Sparse signal representation:**the choice of the sparse signal representation highly depends on the geometrical content of the components . For instance, if the signal is assumed to be composed of oscillatory structures, it is customary to choose as a discrete cosine transform or its localized variant. In case mainly contains local anisotropic contour-like features, the curvelet transform is a good fit. In this framework, redundant wavelets are a generic choice that generally provides good recovery results for most reconstruction problems that involve natural images. In this article, will be chosen as an undecimated wavelet frame [30]. Strictly speaking, undecimated wavelet frames are not orthogonal representations, which entails that the proximal operator defined in Eq. of Appendix Section 8 is an approximation. Nevertheless, it is customary to use Eq. along with undecimated wavelets. Indeed, in that specific case, the Gram matrix of the representation is diagonally dominant, which makes it close to an isometry.**Regularization parameters:**whether it is in the 2-step reconstruction approach, or in the combined reconstruction approaches including the reweighted FISTA or the proposed constrained GFB, the regularization parameters contained in matrix aim at balancing between the sparsity constraint and the data fidelity term. These parameters define thresholds that are applied to the expansion coefficients of in the sparse representation , which eventually act as a denoising procedure. Therefore, in practice, these parameters are chosen so as to reject noise-dominated coefficients in . For that purpose, the weight matrix is built so that each of its elements are the product of two terms: where and .

The first term defines the global threshold per polarimetric component and its value is derived straight from the derivative of the data fidelity term:where the median absolute deviation () is a robust empirical estimator of a Gaussian noise standard deviation, is the gradient of the data fidelity term, and is a scalar that is generally chosen between and . This choice holds true for the three reconstruction approaches.

The extra parameters are the standard parameters that are defined in reweighted techniques. Following [31], these parameters are chosen based on some estimate of the polarimetric componentswhere is a small scalar. This procedure is applied to the reweighted FISTA algorithm as well as the constrained GFB.

**Number of iterations, reweighted steps and stopping criterion :**In the next, the maximum number of iterations is fixed to , except specified otherwise. For the reweighted algorithms, two reweighted steps () taking place respectively after and iterations were sufficient to maximize reconstruction quality. Lastly, for all reconstruction algorithms, the stopping criterion is based on the relative variation of the solution between two consecutive iterations: , where .

## 4Numerical results

In this section, we analyze the performance of the reconstruction algorithms and regularization procedures described above. 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 standardly evaluated throughout this section by computing the Peak Signal to Noise Ratio (PSNR) of a concatenated vector of the reconstructed polarimetric components, i.e., . Then, we will present some polarimetric imaging numerical results on 2D signals that demonstrate the ability of this concept of 2-pixel polarimetric camera to provide satisfactory compressed polarization contrast images at low cost and low complexity.

### 4.1Comparison of algorithms performance and robustness

#### Description of 1D test signal

In order to optimize computation resources, we generated a synthetic polarimetric 1D data sample of length , that will be used throughout this subsection to compare the performance of the above algorithms. The corresponding signals and are plotted in the inset (a) of Fig. ?, respectively with blue and red solid lines. The simulated polarimetric components verify the positivity constraint , and it can be noted that their supports are not joint, for the sake of generality. In the inset (b) of Fig. ?, we plot a set of simulated detected intensity signals and , with 40 compression rate (), and SNR = 40 dB. The mixing matrix used to generate the data has been calculated as detailed in Appendix Section 7, assuming an incidence angle and wavelength of 780 nm (i.e., the optimal conditions identified in Section 2.3). The patterns used on the DMD to sample the image were generated from a randomly scrambled Hadamard transform of the signal of size . Unless otherwise specified, the same experimental conditions will be assumed for all numerical results presented in this article. In Fig. ?, is plotted as a function of , evidencing the strong anticorrelation existing between the two detected signals. Lastly, in the inset (a) of Fig. ?, we show an example of reconstructed signals and obtained with the reweighted FISTA algorithm, and yielding PSNR = 37.7 dB. As for all reconstructions of the 1D signal presented below, the reconstruction process makes use of the unidimensional undecimated wavelet transform with the Haar filter, which is well suited to sparesely represent piecewise constant signals. It can be checked in the inset (a) of Fig. ? that the algorithm has not been constrained to ensure positivity of either , or .

#### Algorithms 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 as a function of the SNR is given in Fig. ?. All the data points in Fig. ? 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 the third quartile of the data series.

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 dBSNR50 dB), the 2-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 2-step algorithm. Lastly, imposing physical positivity constraints on , and 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 SNR 50 dB, the positivity constraints prove efficient to improve the reconstruction quality. A maximum gain of almost dB is obtained for SNR = dB, a situation where the unconstrained rFISTA approach fails to surpass the reconstruction quality obtained with the 2-step procedure.

We now analyze the influence of the compression rate and of the incidence angle on the PSNR of the reconstructed signals. As discussed in Section 2.3, the incidence angle controls the condition number of the mixing matrix, and hence the difficulty of the unmixing problem. For this numerical experiment, we consider only the combined-rFISTA algorithm (plotted in blue solid line in the previous figure), with fixed number of iterations (), and with a SNR of 60 dB. A 2D-plot of the average PSNR obtained over 10 realizations of each numerical experiment is provided in Fig. ?, for 40 values of ranging between to , 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” behaviour frequenty observed in CS problems, with three main domains which can be identified. Firstly, for compression rates below 50 , the reconstruction is almost perfect (PSNR dB), whatever be the value of . Only for highest values of (i.e., for ), the algorithm fails is reaching a PSNR of 50 dB, but remains above 45 dB. Then, for intermediate values of compression rate between 50 and 80 , the reconstruction quality gradually decreases while remaining above dB. In that regime, the influence of the condition number 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 . This 2D-map is rather encouraging towards possible experimental implementations of the 2-pixel polarimetric CS camera. Indeed, it shows that provided a good SNR can be ensured (here dB), polarimetric signals can be retrieved at relatively high compression rates (), and for reasonable incidence angles on the DMD.

#### 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 , and random errors on the micromirrors tilt angles.

The influence of a bias in has been analyzed as follows. Measurement vector was generated assuming a SNR of 80 dB, and an incidence angle , as in Fig. ?. However, the reconstruction procedure was run with an incorrect mixing matrix , assuming a wrong incidence angle of . This way, we mimicked an experimental bias comprised between and on the incidence angle. The reconstruction PSNR obtained with different algorithms is plotted in Fig. ? in solid lines for a compression rate of . 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 (), a significant drop of more than 20 dB in reconstruction quality can be observed for all three algorithms with respect to unbiased reconstruction, but still offering satisfactory recovery quality (PNSR dB). Applying physical constraints in the reconstruction procedure with the GFB algorithm seems to alleviate the degradation, for any magnitude of , as far as simple angular bias is considered. However, the PSNR value of the reconstructed vector 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 second polarimetric component 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 .

Concerning the influence of random tilt errors on the micromirrors, we simulated this effect by computing individual mixing matrices for all micromirrors simulated, assuming that the tilt angles and were affected by a random error. For the sake of computational speed, we assumed a uniform distribution over values of the tilt error between . The PSNR of the reconstructed signal with angular bias and uniform tilt error is also plotted in Fig. ? 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 () where the incorrect “average” matrix is responsible for most of the quality degradation. On the other hand, for very low values of , 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 evidence that with imperfect experimental configurations, using a combined recovery approach instead of a 2-step reconstruction procedure can be advantageous. Correct PSNR values ( 35 dB) can be reached with the rFISTA or GFB algorithms with small angular biases () and in the presence of random tilt errors.

### 4.2Numerical polarimetric imaging results on 2D signals

In this last section, we analyze the ability of the 2-pixel polarimetric CS camera to retrieve polarimetric contrast images from simulated 2D data. Owing to its performance and its simple implementation with respect to constrained GFB, we only consider in this section the reweighted FISTA algorithm 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 reconstuction quality, before a more realistic example of polarimetric image reconstruction is given in Section ?.

#### 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 pixels. This intensity image , plotted in Fig. ? is supposed to be the sum of two polarimetric image components and , yielding a true OSC map also plotted in Fig. ?. In this scenario, we assume that a first object (smaller square) cannot be distinguished from the second object (bigger square) on the intensity image , while OSC map allows the two objects to be clearly identified. The background is assumed totally depolarizing (OSC = ). The smaller square object is always supposed slightly depolarizing (OSC = ) whereas the OSC of the second object (bigger square) is varied between 0 and 1 in the following numerical experiments. Fig. ? shows an example of reconstruction of the total intensity and the OSC map with the combined-rFISTA algorithm (with Haar wavelets) for SNR = dB, compression rate , incidence angle and without bias or uncertainty on the tilt directions. The reconstruction quality is visually very good (PSNR = dB), as evidenced by the reconstruction error map of the total intensity and OSC shown in Fig. ?. This simple example demonstrates the possibility of providing relevant polarimetric information from the proposed concept of 2-pixel polarimetric CS camera.

In the first row (a) of Fig. ?, 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. ?.c. These results evidence 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 burried in the noise and totally filtered out by the regularization process. This is clearly seen for OSC = , where the sharpest features of the first object disappear in the reconstructed OSC map. Obviously, reconstruction quality is maximimed 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.

#### Example of reconstruction on realistic image data

Laslty, we present an example of reconstructed polarimetric image on a more realistic imaging scenario. For that purpose, we considered a true intensity image of the *cameraman* with size , as plotted in Fig. ?. Appropriate polarimetric components and were generated so that a true OSC map would reveal 4 hidden objects (3 in the grass, 1 in the buildings) over a depolarizing background, as can be seen in Fig. ?. The reconstruction results with Symmlet wavelet transform are also displayed in Fig. ?, 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 4 hidden objects remain of course invisible in the reconstructed image . Contrarily, the reconstructed OSC map makes it possible to identify the presence of the 4 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 4 hidden objects is fairly retrieved. However one can notice significant reconstruction errors in the darkest regions of the image (cameraman and tripod). These imperfections could be lowered in the future by refining regularization parameters and constraints in the algorithm implementation.

## 5Conclusion and perspectives

In this article, we have proposed a new concept of CS polarimetric imaging inspired from 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 informations 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 a joint CS and source separation tasks, that can be tackled independantly and successively using standard approaches (respectively minimization and direct matrix inversion), or optimally treated in an original combined reconstruction approach. For that 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 to increase 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 to validate these approaches and to analyze the influence of experimental conditions, such as incidence angle, SNR, micromirrors 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 mirrors tilt direction. Significant compression rates could be achieved while still offering sufficient reconstruction quality, as illustrated on the 2D simulation results presented.

Experimental implementation of this computational imaging approach appears as a natural perspective to this work. Another interesting research track is the design of a blind calibration method so as to compensate possible mis-estimation of the mixing matrix 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.

## Appendices

## 7Computation of the reflection coefficients

We consider a mirror with complex refraction index , which can be conveniently rewritten [32]. We consider the reflection of a beam propagating in an input medium consisting of air (). For an incidence angle on the mirror, one has from the generalized Snell’s law , with the refraction angle. Setting and , one has .

It can be shown [32] that with such notations the reflection coefficient in intensity for TE waves (S-polarized) can be written [32]

whereas for the TM (P-polarized) waves, one has

with and , and and .

From these equations, we were able to compute the micromirrors reflection coefficients for any incidence and any wavelength. The complex refraction index was obtained from a recent evaluation of aluminum reflection coefficients in Reference [33].

## 8Reconstruction algorithms

### 8.1Forward-backward splitting algorithm

Whether it is in the 2-step or in the combined sparse reconstruction approach, recovering the polarimetric signal requires solving optimization problems of the form

where the data fidelity term is a quadratic norm that is differentiable and whose gradient is Lipschitz with some constant , and is a convex but non-smooth penalization. Minimization problem can be carried out efficiently with recent proximal algorithms [34], and more precisely with the Forward-Backward Splitting (FBS) algorithm [35]. The FBS algorithm is an iterative procedure that can be described with the following update rule at iteration

where stands for the derivative of in and is the gradient path length. The proximal operator is defined by the solution to the problem

While the proximal operator of the penalization function is defined as the solution of an optimization problem, standardly used proximal operators admit a closed form expression (see Appendix Section 8.2).

In this article, the term will penalize non sparse solutions and will be based on the norm. Introduced in [31], the re-weighted norm further introduces weights that aim at reducing the bias induced by the standard norm. Consequently, the penalization term used in this article will take the generic form

where is the signal representation where sparsity is modeled. The final optimization procedure then alternates between updates of for fixed weights and updates of these weights, which has been showed to dramatically improves 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 [36]. The generic description of the reweighted FISTA algorithm is given in Algorithm ?

This procedure generally increases the accuracy of the reconstruction process with few updates of the weights , typically between to . Details about practical parameter tuning are given in Section 3.4.

### 8.2Useful proximal operators

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

#### Reweighted

Assuming that the sparse signal representation is an orthogonal matrix, the proximal operator of is defined as

where the weighted soft-thresholding operator is defined as

#### Positivity constraint

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

#### Inequality constraint

The inequality constraint can be recast as , where . Its proximal operator is defined as the orthogonal projector onto the convex set , which is defined as

This expression can be more conveniently recast in a Lagragian formulation by introducing the Lagrange multipliers

Its optimum is obtained for , which should verify the constraint . This entails

Consequently, the Lagrange multipliers must take the values

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

## Acknowledgments

The authors would like to thank Anthony Carré for his help in the 3D rendering of Fig. ?. This work is supported by the European Research Council through the grant LENA (contract no. 678282).

### Footnotes

- Institut FOTON, University of Rennes 1, CNRS, Campus de Beaulieu, Rennes, France
- CEA, IRFU, Service d’Astrophysique-SEDI, Gif-sur-Yvette, France

### References

- 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,”
*Nature Medicine*, vol. 5, no. 10, pp. 1209–1213, 1999. - F. Boulvert, B. Boulbry, G. Le Brun, B. Le Jeune, S. Rivet, and J. Cariou, “Analysis of the depolarizing properties of irradiated pig skin,”
*Journal of Optics A: Pure and Applied Optics*, vol. 7, no. 1, p. 21, 2005. - 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,”
*physica status solidi (c)*, vol. 5, no. 5, pp. 1423–1426, 2008. - S. Breugnot and P. Clémenceau, “Modeling and performances of a polarization active imager at lambda=806 nm,” in
*Laser Radar technology and applications IV*, G. W. K. C. Werner, Ed., vol. 3707.1em plus 0.5em minus 0.4emProc. SPIE, 1999, pp. 449–460. - 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.*, vol. 48, no. 8, pp. 1610–1618, 2009. - J. Fade, S. Panigrahi, A. Carré, L. Frein, C. Hamel, F. Bretenaker, H. Ramachandran, and M. Alouini, “Long-range polarimetric imaging through fog,”
*Applied optics*, vol. 53, no. 18, pp. 3854–3865, 2014. - 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,” in
*SPIE Sensing Technology+ Applications*.1em plus 0.5em minus 0.4emInternational Society for Optics and Photonics, 2014, pp. 90 990B–90 990B. - C. Brosseau,
*Fundamentals of polarized light - A statistical Optics approach*.1em plus 0.5em minus 0.4emJohn Wiley, 1998. - N. Gupta, R. Dahmani, and S. Choy, “Acousto-optic tunable filter based visible- to near-infrared spectropolarimetric imager,”
*Optical Engineering*, vol. 41, no. 5, pp. 1033–1038, 2002. - G. Anna, H. Sauer, F. Goudail, and D. Dolfi, “Fully tunable active polarization imager for contrast enhancement and partial polarimetry,”
*Applied optics*, vol. 51, no. 21, pp. 5302–5309, 2012. - A. Bénière, M. Alouini, F. Goudail, and D. Dolfi, “Design and experimental validation of a snapshot polarization contrast imager,”
*Appl. Opt.*, vol. 48, no. 30, pp. 5764–5773, 2009. - F. Goudail and P. Réfrégier,
*Statistical image processing techniques for noisy images: an application oriented approach*.1em plus 0.5em minus 0.4emNew York: Kluwer, 2004. - R. G. Baraniuk, “Single-pixel imaging via compressive sampling,”
*IEEE Signal Processing Magazine*, 2008. - 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,”
*Applied Physics Letters*, vol. 93, no. 12, p. 121105, 2008. - A. Wagadarikar, R. John, R. Willett, and D. Brady, “Single disperser design for coded aperture snapshot spectral imaging,”
*Applied optics*, vol. 47, no. 10, pp. B44–B51, 2008. - A. Asensio Ramos and A. López Ariste, “Compressive sensing for spectroscopy and polarimetry,”
*arXiv*, vol. 509, p. A49, Jan. 2010. - V. Studer, J. Bobin, M. Chahid, H. S. Mousavi, E. Candes, and M. Dahan, “Compressive fluorescence microscopy for biological and hyperspectral imaging,”
*Proceedings of the National Academy of Sciences*, vol. 109, no. 26, pp. E1679–E1687, 2012. - Y. August and A. Stern, “Compressive sensing spectrometry based on liquid crystal devices,”
*Optics letters*, vol. 38, no. 23, pp. 4996–4999, 2013. - 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 e-prints*, Sep. 2013. - V. Durán, P. Clemente, M. Fernández-Alonso, E. Tajahuerce, and J. Lancis, “Single-pixel polarimetric imaging,”
*Optics letters*, vol. 37, no. 5, pp. 824–826, 2012. - 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,”
*Applied Physics B*, vol. 113, no. 4, pp. 551–558, 2013. - S. S. Welsh, M. P. Edgar, R. Bowman, B. Sun, and M. J. Padgett, “Near video-rate linear stokes imaging with single-pixel detectors,”
*Journal of Optics*, vol. 17, no. 2, p. 025705, 2015. - 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*, vol. 32, no. 11, pp. 2178–2188, Nov 2015. - Y. C. Eldar and G. Kutyniok,
*Compressed sensing: theory and applications*.1em plus 0.5em minus 0.4emCambridge University Press, 2012. - =2plus 43minus 4
*DLP6500-0.65-1080p Datasheet*, Texas Instruments, october 2016, rev. B. [Online]. Available: http://www.ti.com/product/DLP6500FYE/datasheet =0pt - I. N. Flamarique, C. W. Hawryshyn, and F. I. Hárosi, “Double-cone internal reflection as a basis for polarization detection in fish,”
*JOSA A*, vol. 15, no. 2, pp. 349–358, 1998. - E. Candès and M. Wakin, “People hearing without listening: An introduction to compressive sampling,”
*IEEE Signal Processing Magazine*, March 2008, 21–30. - E. Candès, “Compressive sampling,”
*International Congress of Mathematics, Madrid*, 2006. - =2plus 43minus 4 H. Raguet, J. Fadili, and G. Peyre, “A Generalized Forward-Backward Splitting,”
*SIAM Journal on Imaging Sciences*, vol. 6, no. 3, pp. 1199–1226, 2013. [Online]. Available: http://hal.archives-ouvertes.fr/hal-00613637/ =0pt - J.-L. Starck, J. Fadili, and F. Murtagh, “The undecimated wavelet decomposition and its reconstruction,”
*IEEE Transactions on Image Processing*, vol. 16, 2007, 297–309. - E. J. Candes, M. B. Wakin, and S. P. Boyd, “Enhancing sparsity by reweighted L1 minimization,”
*Journal of Fourier Analysis and Applications*, vol. 14, no. 5, 2008, 877–905. - M. Born and E. Wolf, “Principles of optics, 7-th ed,”
*Cambridge University, Cambridge*, pp. 752–754, 1999. - A. D. Rakić, A. B. Djurišić, J. M. Elazar, and M. L. Majewski, “Optical properties of metallic films for vertical-cavity optoelectronic devices,”
*Applied optics*, vol. 37, no. 22, pp. 5271–5283, 1998. - N. Parikh and S. Boyd, “Proximal algorithms,”
*Foundations and Trends in Optimization*, vol. 1, no. 3, 2014. - P. L. Combettes and V. R. Wajs, “Signal recovery by proximal forward-backward splitting,”
*Multiscale Modeling & Simulation*, vol. 4, no. 4, 2005, 1168–1200. - A. Beck and M. Teboulle, “Fast iterative shrinkage-thresholding algorithm for linear inverse problems,”
*SIAM journal on imaging science*, vol. 2, 2009, 183–202.