Survey of Demosaicking Methods for Polarization Filter Array Images

Snapshot polarization imaging has gained interest in the last few decades. Recent research and technology achievements defined the polarization Filter Array (PFA). It is dedicated to division-of-focal plane polarimeters, which permits to analyze the direction of light electric field oscillation. Its filters form a mosaicked pattern, in which each pixel only senses a fraction of the total polarization states, so the other missing polarization states have to be interpolated. As for Color or Spectral Filter Arrays (CFA or SFA), several dedicated demosaicking methods exist in the PFA literature. Such methods are mainly based on spatial correlation disregarding inter-channel correlation. We show that polarization channels are strongly correlated in images. We therefore propose to extend some demosaicking methods from CFA/SFA to PFA, and compare them with those that are PFA-oriented. Objective and subjective analysis show that the pseudo panchromatic image difference method provides the best results and can be used as benchmark for PFA demosaicking.


Introduction
Polarization imaging is a way to analyze the particular direction of oscillation of the electric field described by the light. In opposition with conventional color or multispectral imaging that sample the spectral information, polarization imaging considers the electric field as a vector. Such a vector field is contained in a plane perpendicular to the direction of propagation. As the wave travels, it can oscillate in one particular direction (linear polarization), or in an ellipse (elliptic or circular polarization). Values of polarization images depend on the polarization properties of both the light source and the objects that compose the observed scene. The light can be partially polarized or unpolarized, resulting from either a rapidly changing state of polarization, or an interference effect of polarization.
Several polarization imaging systems, called polarimeters, have been developed in the last past few decades for recovering the polarization state of a lighted scene from few acquisitions. Such systems combine a standard panchromatic imaging device with polarizing optics, e.g., polarization filter, liquid crystal modulator, or prism. Reviews of recent polarimeters have been achieved in the literature [1,2]. The most simple optical setup consists in the rotation of a linear polarization filter at several polarization angles in front of a camera. After a preliminary calibration step (radiometric and polarimetric), the polarization states of the incoming irradiance that reaches the sensor can be estimated. However, this setup is sequential and slow, since several image acquisitions at different filter orientations are required to recover the polarization states of a single scene. To overcome this limitation, Polarization Filter Array (PFA) imaging provides a way for snapshot acquisition that could be useful for many Let us consider the intensities of light measured I 0 , I 45 , I 90 , I 135 after the light is filtered by linear polarization filters (oriented at 0°, 45°, 90°, and 135°). In the literature, the choice of these 4 angles forms an optimal system for polarization imaging in the presence of noise, as described in [29]. The mathematical formulations for Stokes parameters and descriptors are as follows: The total incoming irradiance is represented by S 0 , S 1 is the horizontal/vertical polarization difference, whereas S 2 is the +45/−45°polarization difference. If we consider that channels I k , k ∈ {0, 45 It is useful to note that a radiometric calibration is very important in case of polarimetric imaging, even more than for color imaging, as the different channel errors are coupled, and thus it can invalidate greatly the parameter estimation [1]. An example of a complete 2D Stokes processing starting from a PFA image is given in Figure 2.
The purpose of this paper is to first study the correlation properties of polarization channels and their similarities with those of spectral channels. Then to review some existing interpolation strategies dedicated to filter array imaging, i.e., CFA, SFA, and PFA. Finally, we propose to evaluate objectively the methods and those we have adapted to the PFA case, in the special context of PFA. A diagram of the proposed analysis is shown in Figure 3. We organize the paper as follows. First, a data correlation study across the polarization channels is presented in Section 2. Next, different CFA, SFA, and PFA interpolation techniques are presented in Section 3. Results and discussion of the surveyed methods is proposed in Section 4. The paper ends with several conclusions in Section 5.   (Table 3) Demosaiced image Table 1  Table 2 St okes descript ors comput at ion Subject ive analysis Object ive analysis Figure 3. Overview of the analysis which is done in this paper.

Polarimetric Channel Correlation Study
All demosaicking methods estimate missing values using spatial (intra-channel) (i) and/or inter-channel (ii) correlation assumptions. (i) The spatial correlation assumes that; if a pixel p and its neighborhood belong to the same homogeneous area, the value of p is strongly correlated with the values in its neighborhood. Thus, assuming that a channel is composed of homogeneous areas separated by edges, the value of a pixel can be estimated by using its neighbors within the same homogeneous area. Spatial gradients are often used as indicators to determine whether two pixels belong to the same homogeneous area. Indeed, gradient considers the difference between values of two spatially close pixels. We can therefore assume that these pixels belong to the same homogeneous area if the gradient is low, and that they belong to different homogeneous area otherwise. (ii) The inter-channel correlation (also called spectral correlation in CFA and SFA imaging) assumes that the high frequencies (textures or edges) of the different channels are strongly correlated. If the filter array contains a spatially dominant band, demosaicking generally estimates the associated channel whose high frequencies can be faithfully reconstructed, then uses it as a guide to estimate other channels. The faithfully reconstructed image can be used to guide the high frequency estimation within the different channels [30].
Usual PFA demosaicking methods assume only spatial correlation, thus disregarding correlation among polarization channels. In order to extend CFA and SFA demosaicking methods that also use the inter-channel correlation to PFA images, we propose to compare the spatial and inter-channel correlations in multispectral images with those of polarization images. For this purpose, we use the database proposed in [31]. Images were acquired by the single-band near-infrared sensor from the JAI AD080 GE camera, coupled with a linear polarizer from Thorlabs (LPNIRE100-B). A precision motorized rotation stages (Agilis™ Piezo Motor Driven Rotation Stages) allowed to take the four images at four orientation angles ([I 0 , I 45 , I 90 , I 135 ] T ). A registration procedure aligned the images [32] pixel-to-pixel. The images were also calibrated with respect to the spatial deviation of the illuminant and the non-linearities. There are ten multispectral images, each one being provided with four different polarization angles k ∈ {0, 45, 90, 135}. Scenes imply different objects with materials like fabrics, plastics, food, color checkers, glass, and metal. Conditions of acquisition are constant for all scenes, i.e., constant illuminant source (tungsten halogen source) and location, constant field of view and constant lens aperture. Multispectral recoverred images are composed of six spectral channels: Five channels are associated with the visible domain, whereas one channel is associated with the Near-InfraRed domain (NIR). The six spectral channels C u are arranged so that their associated spectral band wavelengths increase with respect to u ∈ {1, . . . , 6}.
Let us first study the properties of multispectral images with respect to the polarization angle of analysis. For this purpose we assess the spatial correlation within a given channel C u using the Pearson correlation coefficient (PCC) between the value C u p of each pixel p and that of its right neighbor C u q at spatial distance 2. This coefficient is defined as [33] PCC where µ u is the mean value of channel C u . We also assess the inter-channel correlation using the PCC between any pair of spectral channels C u and C v , (u, v) ∈ {1, . . . , 6} 2 as Note that in Equations (4) and (5), we select a centered area excluding the 16 pixels on the image borders to avoid border effects, that are induced by the registration step used on raw images (described in [31]). Moreover the choice of 16 border pixels is done to match with the experimental assessment (see Section 4) of demosaicking methods presented in Section 3. Table 1 is the spatial correlation within each spectral channel and the inter-channel correlation between the six spectral channels according to each of the four polarization angles. Table 1 shows that the spatial correlation is relatively high (0.9504 on average over all channels and polarization angles), which validates the use of the spatial correlation assumption for both SFA and PFA demosaicking.
According to Table 1a,d, the spatial correlation has the same behavior for the four polarization angles. It also highlights that the channel C 4 has low spatial correlation. We believe that it is due to the database acquisition setup, which uses the dual-RGB method leading to unbalanced channel sensitivities. In this configuration, the spectral sensitivity function associated with the channel C 4 is lower than other channels over the spectrum. Thus, all channels don't share the same noise level, and poor information recovery (especially for C 4 ) could lead to low correlation values.
Regarding the spectral inter-channel correlation, the usual behavior is that close spectral channels in term of wavelength band are more correlated than distant ones, and channels in the visible are weakly correlated with the near-infrared channel [11]. Except the channel C 4 that exhibit low correlation values, this behavior is observed in Table 1.
Moreover the correlation between the NIR channel C 6 and other channels is low (ranges on average between 0.7953 and 0.8787), while the correlation between channels in the visible domain reaches up to 0.9596 (correlation between C 2 and C 3 ).  Table 1c. We can therefore expect that the polarization channels at 0°are more correlated with those at 135°than those at 45°or 90°. Table 1. Inter-channel correlation between the six spectral channels relatively to the polarization angle of analysis. Last line of each subtable is the spatial correlation within each channel. Values are averaged over the ten multispectral images from [31]. Last subtable (e) is the average over the four polarization channels.   Now, let us consider the polarization images composed of four polarization angles for a given spectral band. The spatial and inter-channel correlations are assessed using the PCC applied respectively to channels I k , k ∈ {0, 45, 90, 135} (see Equation (4)), and to any pair of polarization channels I k and I l , (k, l) ∈ {0, 45, 90, 135} 2 (see Equation (5)). Table 2 is the average polarization correlation between the four channels of polarization images, according to each of the six spectral bands. Results highlight that the spatial correlation is high and does not depend on the considered spectral band (except for channel C 4 ). Results also confirm that channel I 0 is highly correlated with channel I 135 and channel I 45 is highly correlated with channel I 90 .
In general terms, inter-channel correlation between polarization channels is higher than inter-channel correlation between spectral channels (see Table 1). Indeed, if the incoming irradiance is not polarized, the associated pixel has only the information of the total intensity divided by two, that is the same from one channel to another. Table 2. Inter-channel correlation between the four polarization channels according to the spectral band. Last line of each subtable is the spatial correlation within each channel. Values are averaged over the ten multispectral images from [31]. Last subtable is the average over the six spectral channels.  Since the inter-channel correlation is high in polarization images, we propose to apply SFA demosaicking schemes based on inter-channel correlation assumption on PFA images. For this purpose, we can choose the four polarization channels associated to any spectral band but not the one associated to C 4 . Since dual-RGB method is not applied for the channel C 6 , we selected it for the experimental assessment in Section 4.

Demosaicking Problem and Properties
A PFA camera provides a raw image I raw with X × Y pixels, in which a single polarization angle k ∈ {0, 45, 90, 135} is available at each pixel p according to the PFA arrangement. Let S be the set of all pixels (with a cardinal of |S| = X × Y) and S k be the pixel subset where the PFA samples the angle k, such that S = k∈{0,45,90,135} S k . A PFA can be defined as a function PFA: S → {0, 45, 90, 135} that associates to each pixel p its polarization angle. Therefore the pixel subset where the PFA samples the polarization angle k can be defined as S k = {p ∈ S, PFA(p) = k}. The raw image I raw can then be seen as a sampled version of the fully-defined reference image I = {I k } k∈{0,45,90,135} (that is unavailable in practice) according to the PFA: ∀p ∈ S, I raw The raw image can also be seen as the direct sum of four sparse channelsĨ k , k ∈ {0, 45, 90, 135} that contains the available values at pixel positions in S k and zero elsewhere: where denotes the element-wise product and m k is a binary mask defined at each pixel p as: Demosaicking is performed on each sparse channelĨ k to obtain an estimated imageÎ = {Î k } K k=1 with four fully-defined channels, among which three are estimated at each pixel p. For illustration purpose, Figure 4 shows the demosaicking problem formulation for a PFA raw image.
PSfrag replacements  In the following, we review the demosaicking methods dedicated to PFA. We also review those dedicated to CFA/SFA that can be used or adapted to our considered PFA. All these methods were either re-coded, adapted to PFA, or kindly provided by authors (online or in private). See Table 3 for an overview of all methods. Table 3. Summary of the color, spectral, and polarization filter array (CFA/SFA/PFA) interpolation methods. R, A and P abbreviations mean that the algorithms were Re-coded, Adapted, or Provided by the authors of the original work.
Year Code

PFA Demosaicking
Among PFA demosaicking methods, we exclude the learning-based methods [46], since they require well-adapted dictionaries, and methods that exploit multiple sampling of the raw data [47]. We also exclude the gradient-based method [48], since a SFA method has a very close behavior (BTES).

Bilinear with 5 Different Kernels (B)
Bilinear interpolation dedicated to PFA was firstly investigated by Ratliff et al. [34]. They employ three different bilinear and two weighted bilinear kernels (see Figure 5). Bilinear interpolation is one of the most commonly used technique due to its low computational complexity. It is based on space-invariant linear filtering. Two kinds of bilinear interpolations exist. One uses a linear combination of neighboring pixel values using equal weights. Another employs non-equal weights in accordance to the Euclidean distance between the interpolated pixel location and centers of neighboring pixels. The subtractive nature of the Stokes vector processing results in strong edge artifacts in the reconstructed images. Based on this assumption, authors define the term of Instantaneous Field-Of-View (IFOV), which is the local deviation between an ideal full-resolution polarimeter and the interpolated PFA pixel responses at each position. They evaluate the interpolation performance of the methods using synthetic data, in the frequency domain of the reconstructed Stokes images. As evaluation metrics, they computed the modulation and intermodulation transfer functions in the descriptor images, along with the DOLP Mean Squared Error (MSE). It is found that the larger the size of the kernels becomes, the more DOLP artifacts are reduced, at the cost of loosing the high spatial frequency features. They found that the 12-pixel neighborhood kernel (B 4 in the Figure 5) gives the best performance in term of IFOV removal. For algorithm implementations, we used the same weights as the original paper for the two weighted bilinear kernels. Figure 5. The five demosaicking kernels of the five bilinear methods B 1−5 from the work by Ratliff et al. [34]. It refers to the neighborhood used for interpolation. (a-c) Are simple bilinear kernels, whereas (d,e) are weighted bilinear kernels.

Linear System (LS)
Tyo et al. [35], in 2009, elaborates a new method to reconstruct the first three Stokes parameters directly from the mosaicked image, without estimatingÎ. The four polarization imagesÎ 0 ,Î 45 ,Î 90 , andÎ 135 are thus not available with this method. The philosophy starts from the analysis of a raw PFA image in the frequency domain. By doing the discrete 2D Fourier transform, they define the spatial low-pass and high-pass filters. They assume that S 0 , S 1 , and S 2 are spatially band limited in the frequency domain. The centering baseband of the Fourier transform represents S 0 , whereas the horizontal and vertical sidebands represent S 1 + S 2 and S 1 − S 2 respectively. They could then reconstruct S 1 and S 2 after applying the filters in the Fourier domain, and by doing the inverse Fourier transform of images.

Adaptive (A)
An extension of the bilinear interpolation was proposed by Ratliff et al. [36]. In this work, the principle is inspired by Ramanath et al. [49]. The loss of high frequency features in bilinear interpolation techniques is compensated by local computation using a 3 × 3 filtering. First, S 0 is approximated using not only I 0 and I 90 , but the four available neighboring intensities, as it is suggested in the literature [50]. A 2 × 2 low-pass filtering of the raw PFA image is performed with the kernel as follows: Then, intensity similarity masks and Euclidian distance masks are computed, in such a way that the weights are higher for pixels that have similar intensity values within a close neighborhood. These local interpolation weights are computed at each position in the image, and avoid interpolation across edges, and thus preserve high frequency contents. Results show that IFOV artifacts and false edges are minimized in the DOLP image, while high spatial frequencies are preserved. Only a subjective evaluation of their algorithm is performed in the article. The parameter ρ i in the paper was selected to be equal to 0.3 in our implementation.
In our work, the cubic and bicubic interpolation methods have been implemented using built-in functions from Matlab software (The MathWorks, Inc., Natick, MA, USA). Cubic interpolation uses the third order polynomial fitting to interpolate an area delimited by four corners, and uses three directional derivatives (horizontally, vertically and diagonally) as input. The cubic-spline method is a sequence of an horizontal interpolation and a vertical interpolation. Polynomial fitting (third order) is also used to reconstruct missing values from adjacent pixels, but with the additional constraint that the first and second derivative at the interpolation points are continuous. A modulation transfer function study is done to investigate on how the high spatial frequencies are preserved. A visual and objective evaluation (using MSE) are done on real data. Main results show that the cubic-spline methods performed the best in terms of visual artifacts removal and MSE. It appears that bilinear and weighted bilinear give the worst results.

Intensity Correlation among Polarization Channels (ICPC)
Another method by Zhang et al. [38] takes advantage of the correlations in PFA to enhance the spatial resolution of images. Spatial and polarization correlations between channels are investigated in a close pixel neighborhood, directly in the raw PFA image. Edges can not be accurately distinguished if the incoming light is polarized at some degree. Thus, in their work, edge orientations are estimated using the intensity correlation. They start by computing correlation errors from the assumption that edges have poor correlation within the pixel neighborhood. The correlation error magnitude reflects the presence of a homogeneous zone, or of a horizontal, vertical, or diagonal edge. For the interpolation in itself, a diagonal interpolation is firstly done by applying a bicubic-spline interpolation. Then, horizontal and vertical interpolation are performed by bicubic-spline interpolations, according to the correlation errors previously computed. Evaluation of their method is done by constructing a set of four ground truth polarization images using a linear polarizer rotated at four different angles. They found that their method performs better visual results, and has better RMSE compared to bilinear, bicubic, bicubic-spline, and gradient-based methods.

CFA Demosaicking
Bayer CFA has a dominant green band that represents half of pixels, and is used as a guide containing the high spatial frequencies. Therefore, CFA demosaicking methods generally estimate the green channel first in order to use the spectral correlation by considering that green channel is strongly correlated with blue and red channels. Here, we extend residual interpolation methods [39,40] from the CFA to the PFA pattern by considering the intensity image S 0 as a guide instead of the estimated green channel.

Residual Interpolation (RI)
Kiku et al. [39] propose a demosaicking scheme based on the residual interpolation. Their method requires a well estimated guide image, i.e., the estimated green channel that is dominant in the Bayer CFA raw image. Since there is no dominant band in our considered PFA, we adapt their method by using the intensity image S 0 as a guide. It is well estimated from a simple 2 × 2 bilinear kernel (see Equation (9)). Each channelÎ k is then recovered by following these successive steps:

1.
It computes a tentative estimated channelǏ k by applying the guided filter [51] and the guide image to the sparse channelĨ k . Note that such process modifies the raw values in the tentative estimated channelǏ k .

2.
It computes the residuals defined by a difference betweenĨ k and tentatively estimated channelǏ k at pixels in S k .

3.
It performs a bilinear interpolation of the residuals by using B 3 filter. 4.
The finally estimated channelÎ k is given by the summation of the tentative estimated channelǏ k and the interpolated residuals.

Adaptive Residual Interpolation (ARI)
Monno et al. [40] improve the RI by applying a Laplacian filter onĨ k and the guide before using the guided filter. The parameters for RI and ARI implementations are h = 5, v = 5, and = 0.

Binary-Three Edge Sensing (BTES)
In our knowledge, BTES interpolation [41] is the first SFA demosaicking method applicable on PFA raw images. This method improves the bilinear interpolation by considering weights inversely proportional to the directional gradient. It follows two steps, in which only a subset of pixels are estimated as shown in Figure 6. In a first step, a quarter of pixels are estimated using their four closest neighbors weighted with respect to the diagonal gradients. In a second step, the rest of the pixels ( card(S) 2 ) are estimated using their four closest neighbors weighted with respect to horizontal (for an horizontal neighbor) or vertical (for a vertical neighbor) gradients.  As bilinear interpolation, this method is only based on spatial correlation since there is no dominant channel. Other SFA demosaicking methods also consider the inter-channel correlation to estimate the missing channels.

Spectral Difference (SD)
Brauers and Aach [42] estimate missing values of a channel using the inter-channel correlation. They consider the available value in the raw image at the missing position, i.e., a pixel p of a channel I k is estimated using the information of channel I PFA(p) as follows:

1.
It computes the sparse difference channel∆ k,PFA(p) between channel I k and the channelÎ PFA(p) B 3 estimated by bilinear interpolation (using filter B 3 ) at pixels in S k .

2.
It estimates the fully-defined difference channel∆ k,PFA(p) B 3 by bilinear interpolation.

3.
The value ofÎ k p is given by the sum between the difference channel∆ k,PFA(p) B 3 and the available value at p in the raw image.
Mizutani et al. [57] further improve this method using an additional assumption: Spectrally close channels are more correlated than distant ones. Since this assumption is not validated for polarization images, we cannot use it in this context.

Vector Median (VM)
Wang et al. [43] consider that each pixel of an image as a vector with four dimensions. For each pixel p, the method defines many pseudo-pixels by column vectors ([I 0 p , I 45 p , I 90 p , I 135 p ] T in our case) according to the mosaic, and it affects the median pseudo-pixel to p. The pseudo-pixels at p represents all the possible combinations of the four channels in a 5 × 5 neighborhood around p. The four values of a pseudo-pixel are taken from spatially connected pixels. To preserve value discontinuities and color artifacts, authors also propose a post-processing in a 3D-spheric space.

Discrete Wavelet Transform (DWT)
Wang et al. [44] extend the DWT-based CFA demosaicking to SFA demosaicking. By considering an image as low-frequency (homogeneous areas) and high-frequency contents (edges), This approach assumes two things: The low-frequency content is well estimated by bilinear interpolation, and the high-frequency contents have to be determined more accurately and have to be the same in different channels. The algorithm first estimates a fully-defined multispectral imageÎ B 3 by bilinear interpolation, then applies five successive steps to each channelÎ k B 3 as follows: 1. It decomposesÎ k B 3 into K Down-Sampled (DS) images, so that each DS image is composed of pixels in S l , l ∈ {0, 45, 90, 135}. Note that one DS image is only composed of raw values.

2.
It decomposes each DS image into spatial frequency sub-bands by DWT using Haar wavelet D2.

3.
It replaces the spatial high-frequency sub-bands of all estimated DS images by those of the corresponding DS images of the mid-spectrum channel, assuming this is the sharpest one. In PFA images, there is no mid-spectrum channel, we therefore propose to use arbitrarily theÎ 90 DS images are transformed by inverse DWT.

5.
It recomposes the full-resolution channelÎ k from the four transformed DS images.

Multi-Local Directional Interpolation (MLDI)
Shinoda et al. [45] combine BTES and SD approaches into the MLDI method that follows the two steps of BTES. Each pixel is estimated using the difference planes, as in SD scheme. Moreover, instead of simply estimating the fully-defined difference planes by bilinear interpolation, the authors use directional gradients (following the step in BTES scheme), which improves the estimation. Shinoda et al. [45] also propose a post-processing that updates each estimated channel by taking into account the previous estimated values.

Pseudo-Panchromatic Image Difference (PPID)
The Pseudo-Panchromatic Image (PPI) is defined in each pixel as the average of all channels. By assuming that PPI values of neighboring pixels are strongly correlated, Mihoubi et al. [11] estimate the PPI from the PFA image by applying an averaging filter M proposed in [58]. Such filter estimates the PPI as the average value of all channels in a given neighborhood of each pixel. For this purpose, it takes all channels into account, while being as small as possible to avoid estimation errors. For our considered PFA arrangement, the filter M is adapted as: In the case of strong spectral correlation (≥ 0.9), authors propose to restore the edges of the estimated PPI using directional gradients. However, the condition is not validated for PFA images. The estimated PPI is thereafter used in a PPI difference scheme that estimates each channel k as follows: 1.
It computes the sparse difference channel∆ k,PPI between channel I k and the PPI at pixels in S k .

2.
It estimates the fully-defined difference channel∆ k,PPI by weighted bilinear interpolation in which the weights are inversely proportional to the gradients computed from the raw image. 3.
The finally estimated channelÎ k is the sum between the estimated PPI and the difference plane.

Pseudo-Panchromatic Image based Discrete Wavelet Transform (PPIDWT)
The PPI has similar information to the mid-spectrum channel, and it is better estimated. Mihoubi et al. [11] therefore propose to replace the spatial high-frequency sub-bands by those of the PPI instead of I 90 b 3 channel in the DWT scheme.

Experimental Setup
PFA image simulation is employed to assess the interpolation strategies. As for the correlation study in Section 2, the polarimetric images from the database of Lapray et al. [31] was used as references.
All methods of Table 3 were either re-coded (R), adapted to PFA (A), or provided by authors in Matlab/ImageJ language software (P). They are further integrated into the framework presented in Figure 4 in order to assess and compare the performances of demosaicking. Stokes descriptors are then computed for both reference and estimated images, according to Equations (1)-(3). To avoid precision errors during image conversions, all considered images and processing are using 32-bit float data representation.
We consider the Peak Signal-to-Noise Ratio (PSNR) as quality metric. Between each couple of reference (R) and estimated (E) channel/descriptor, the PSNR is computed as follows: where MSE(R, E) denotes the mean squared error between R and E. Because max p R can differ from a channel (or a descriptor) to another, Equation (11) takes into account this actual maximal level rather than the theoretical one to avoid misleading PSNR values. In PSNR computation, as for the previous correlation study, we exclude the 16 pixels in each of the four borders of the image to avoid inherent border effect related to either registration or demosaicking processing. Table 4 displays the PSNR values provided by the demosaicking methods on average over the ten database scenes. Results show that among bilinear filters, B 3 provides the best results for I 0 , I 45 , I 90 , I 135 , S 0 , S 1 , and DOLP, while B 4 slightly exceeds it for S2 and AOLP. Among PFA-oriented methods, CB and CBSP generally provide the best results. Our proposition to adapt RI and ARI CFA demosaicking methods to the PFA case (with S 0 as guide) provides better results than classical PFA-oriented methods. We remark that RI and ARI are very close together in the PSNR results. RI also provides the best results among all tested methods regarding S 2 parameter and DOLP descriptor.

Results and Discussion
For PFA methods, it is important to note that the output interpolated pixel is shifted by half pixel when using bilinear kernels B 1 , B 4 , and B 5 , compared to other bilinear kernels B 2 , B 3 . The output pixel is either aligned to the original interpolated pixel position center, or in the pixel boundaries. We did not correct for this misalignment because applying an image translation by half pixel needs an additional cubic or linear interpolation. So such a registration process cannot be used as a pre-processing for an acceptable comparison methodology over the bilinear demosaicking methods. Thus, the results for B 1 , B 4 , and B 5 should be interpreted with care.
For tested SFA-oriented methods, the use of spectral correlation generally provides better performance than simple bilinear interpolations. Moreover, methods based on gradient computation (BTES, MLDI, and PPID) exhibit the best demosaicking performances. By considering the PPI as a guide for demosaicking, PPID shows the best demosaicking performances among all methods for all polarization channels, also for S 0 , S 1 parameters and AOLP descriptors.
To visually compare the results provided by demosaicking methods on S 0 , AOLP, and DOLP descriptors, we select a zoomed area from the "macbeth_enhancement" scene of the database. Among demosaicking methods, we show the results of the most intuitive method (bilinear interpolation using B 3 filter), and the pseudo panchromatic image difference (PPID) that globally provides the best PSNR results. Figure 7 shows that there is no significant difference regarding the S 0 parameter, except that the two highlight dots are more apparent in PPID demosaicked image. Computing AOLP and DOLP parameter from a bilinearly interpolated image generates many artifacts that are fairly reduced using PPID demosaicking method. Generally speaking, we found that demosaicking that are dedicated to PFA don't necessary give better PSNR result. Thus, it was not obvious that considering color and spectral demosaicking techniques applied to PFA arrangement could be beneficial. The results highlights that this can benefit the pre-processing of PFA.
However, we can express some reservations about the results obtained. First, we limited our study on a relatively small database. Other polarization database in the literature [59] furbish only the Stokes parameters and polarization descriptors, but no fully-defined reference image I = {I k } k∈{0,45,90,135} are available. Natural scene samples could also be beneficial for a complementary algorithm classification. Secondly, the database used in this work was made with the same experimental conditions, i.e., constant angle/plan of incidence and a unique non-polarized illuminant. We think that supplementary tests of the best identified algorithms in an extended database containing a better statistical distribution of data could be valuable. Thirdly, the noise associated with reference images is not quantified, and is slightly different from a PFA acquisition system. We thus disregarded recent denoising-demosaicking algorithms that estimate sensor noise to improve the accuracy of demosaicking [60][61][62][63].
The arrangement of the filter array investigated consists in a 2 × 2 square periodic basic pattern with no dominant band. Our goal was to stay general and to apply the evaluation on a well-used pattern. But some other tentatives of designing extended pattern have been proposed in the literature [64], for a better spatial distribution of linear polarization filters. An extensive evaluation of best demosaicking algorithms on different arrangements would be considered in a future work.
We found that the acquisition setup may induces correlation between some polarized channels that could be exploited for demosaicking. Since these properties are data-dependent, we have chosen to not use them in our study, despite that they are used in few SFA demosaicking methods.
We remarked that some algorithms need more computation time than others, without necessary giving better results. No computational complexity consideration has been reported in this work. We think that there is a lake of information about these aspects in the original articles. Moreover, Matlab or ImageJ can not provide a consistent evaluation of the complexity of the selected algorithms, e.g., for their potential ability to be parallelized in a hardware acceleration for real-time computing.
(a) Reference S 0 (b) Reference AOLP  Figure 7. Zoom in of the "macbeth_enhancement" scene from the database of Lapray et al. [31].
(a-c) Images resulting S 0 , angle of linear polarization (AOLP), and degree of linear polarization (DOLP) processed using the full-resolution reference.

Conclusions
By considering the inter-channel correlation, CFA and SFA schemes aim to improve the spatial reconstruction of channels from the information of other channels. Experiments on the only available polarization image database have shown that such methods provides better results in term of PSNR than PFA-oriented methods. More particularly, we proposed to adapt two CFA demosaicking algorithms based on residual interpolation to the PFA case, and showed that they provide better results than classical PFA-oriented methods. Moreover, the SFA PPID method provides the overall best results, and largely reduces visual artifacts in the reconstructed polarization descriptors in comparison with bilinear method.
Correlation study has shown that the spectral band considered in the acquisition of polarization channels has no influence on the correlation between polarization channels. The correlation results from this study could be an input and provide assumptions for the design of new demosaicking algorithms applied on cameras that mix spectral and polarization filters.
All algorithms were tested on a small database of images. As future work, we hope that an extensive database of polarization and spectral images will be available soon in the research community. Thus, further evaluations on more various materials and image statistics would validate more deeply our conclusions. Furthermore, we believe that a real-time pipelined implementation of the PPID method using GPU or FPGA needs to be investigated, that would be a valuable tool for machine vision applications.