Mueller matrix polarimetry with 3D integral imaging

In this paper, we introduce the Mueller matrix imaging concepts for 3D Integral Imaging Polarimetry. The Mueller matrix of a complex scene is measured and estimated with 3D integral imaging. This information can be used to analyze the complex polarimetric behavior of any 3D scene. In particular, we show that the degree of polarization can be estimated at any selected plane for any arbitrary synthetic illumination source which may be difficult to produce in practice. This tool might open new perspectives for polarimetric analysis in the 3D domain. Also, we illustrate that 2D polarimetric images are noisier than 3D reconstructed polarimetric integral imaging. To the best of our knowledge, this is the first report on Mueller matrix polarimetry in 3D Integral Imaging. © 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement


Introduction
Polarization imaging techniques are receiving increasing interest because they provide a way to overcome the intrinsic limitations of standard light intensity imaging [1,2]. Whereas conventional cameras record the intensity of the electric field, polarimetric measurements provide information related to the direction and phase difference of the electric field and the ratio between fully polarized and unpolarized light at each point of the scene. This is possible because of the recent availability of polarization image sensors (see, for instance [3][4][5]). Based on the well-known Mueller matrix theory, it is a relatively straightforward task to implement polarimetric Mueller imaging techniques using standard polarization equipment plus a CCD camera [6][7][8][9][10]. Polarization based imaging might be much more sensitive than conventional imaging. These techniques could have high impact in fields such as machine vision, target detection in turbid media or underwater imaging, provided they are adequately implemented according to the requirements of each particular problem [11][12][13]. For example, in machine vision applications, one of the required capabilities is the ability for plane depth discrimination within the scene. For this particular case, besides using polarization imaging, it may be appropriate to take advantage of other imaging techniques such as 3D Integral Imaging (InIm) [14][15][16][17][18][19][20][21][22][23][24][25]. In summary, to enhance the overall capabilities of the entire imaging process, it can be sometimes useful to have a 3D image acquisition step with a subsequent polarization-based imaging and processing stages. In this way, the specific advantages of the two techniques are jointly considered.
In the present work, we propose Mueller matrix imaging using 3D integral imaging. The Mueller matrix for the 3D scene is generated and used in several scenarios. For instance, we show that reconstructed 3D information displays less noise than the corresponding 2D images. Moreover, it is possible to select a particular plane where the in-focus polarimetric information is calculated. Also, provided that the 3D Mueller matrix InIm is known, we show that it is possible to generate polarimetric landscapes corresponding to any arbitrary synthetic polarimetric illuminations which may be difficult to produce in practice.
The paper is organized as follows. In section 2, we review basic concepts on Stokes imaging acquisition and Mueller matrix calculation. In section 3, we describe the optical setup and introduce the scene used in the experiments. The proposed imaging model is tested in section 4 and we compare the experimental Stokes images with those obtained with the proposed model. In next section 5, we use synthetic light in order to artificially produce polarimetric images that cannot be easily generated in practice. Finally, we present our conclusions in Section 6.

Review on Mueller matrix imaging
The Mueller matrix formalism [6] can be applied to 3D scenes illuminated by quasi monochromatic light. This assumption is based on two facts that hold in our experiments: i) the illuminating light has an almost parallel beam configuration, ii) the acquisition camera detects the light scattered by the scene in a certain well-defined direction. Thus, the directional character of both the illumination and the light detection, suggests that the light we collect at each pixel of the camera is the result of an overall interaction of the incident beam with the scene that can be modelled on a Mueller matrix basis. Therefore, we need to determine the Mueller matrix that fully represents each single pixel of the detected images. Eventually, the suitability of our central assumption can only be tested and confirmed by experiments. This is the main objective of the present work.
Let us briefly summarize the Mueller matrix formalism, adapted to our situation. The Mueller matrix of a sample illuminated by partially polarized quasi monochromatic light describes the changes in the state of polarization of the light when it is reflected from the sample, assuming the validity of a linear system scheme. The input and output data are the four Stokes parameters of light. Thus, the corresponding 4x4 Mueller matrix ( , ) i j M at pixel ( , ) Stokes vectors of the illumination source and the recorded light, respectively. In particular, note that we assume the scene is illuminated by a source with a constant polarization state, i.e. not dependent on the pixel. To determine the 16 parameters of matrix M, a suitable number of measurements to solve the system of linear equations in Eq. (1) is required. The Mueller matrix could be calculated using a maximum of four well-selected states of polarization. However, this approach requires the use of wave plates with tunable phase difference [26]. In order to find accurate results and due to the inherent experimental inaccuracies, overdetermination in the system of linear equations may be necessary. Accordingly, we illuminate the scene with 6 polarization states (L 1 to L 6 ) that correspond to the following set of Stokes parameters (Table 1): Input states of polarization L 1, L 2, L 4 and L 5 describe linearly fully polarized light in the 0°, 45°, 90° and 135° directions respectively. States L 3 and L 6 correspond to right and left circularly polarized sources. These six input polarization states are selected because they are simple to produce in the laboratory by using a polarizer and a quarter wave plate. Any polarization state on the sphere of Poincaré can be determined as a combination of states L 1 to L 6 .
For each illumination L 1 to L 6 , we take 6 snapshots (image captures) that enable the measurement of the 4 output Stokes parameters as follows. k I i j is the measured light intensity when an ideal linear polarizer (at an angle of 0° with the X axis) is placed before the detector, and L k stands for the k state of polarization (k = 1,…6, as in Table 1). 45  .
Similarly, we take 2 more snapshots (L ; , ), (L ; , ), RC k LC k I i j I i j whose intensities are transmitted through a quarter-wave plate plus a linear polarizer at angles suitable for allowing passing right (or left) circularly polarized light. Then, we define: .
Due to the experimental character of our work, we have considered more convenient to define: .
The degree of polarization (DoP) for illumination L k at each pixel ( , ) i j of the output is: For each pixel ( , ) i j , we define two 4x6 matrices as follows: and ( )   Each component of ( ) is a Stokes-parameter image obtained when the scene is illuminated with a light source with one of the input polarization states shown in Table 1. Then, the 36 intensity images give rise to an overdetermined linear equation system that can be written as follows: It is important to note that Eq. (9) is similar to Eq. (1), but there are several significant differences. First, the input and output quantities are arranged according to a 4x6 matrix. Second, it is worth to point out that Eq. (9) leads to a best fit estimation of the Mueller matrix ( , ) , that depends on the number of measurements carried out, i.e. the number of columns of ( , ) out i j V and .
inp V Incidentally, form the numerical point of view, modern computing platforms allow an optimal treatment of the mathematical equations involved in Eq. (1), so that finding matrix ( , ) est i j M is relatively straightforward.

Experimental setup
In order to test the validity of our approach, we arranged an experimental setup consisting of a scene with several elements, illuminated with a light source that provides a fairly parallel and monochromatic beam. The image capture can be done with a photographic camera. Integral imaging can be implemented with an array of conventional 2D cameras, or a single 2D camera with a lenslet array, or a single moving 2D camera. Or, we may use an available plenoptic camera. Figure 1(a) shows a sketch of the experimental arrangement, and the scene used in the experiments in Fig. 1(b). The scene combines several elements made of metal, plastic or glass and a very reflective stage; moreover, some objects are submerged in liquid.
In what follows, we refer to plane F (objects located around 90 cm far from the camera) and plane N (objects located around 70 cm far from the camera) for 3D experiments. Several details of the experiment are worth to discuss. First, the directional character of both the illumination and the detection processes is preserved. This is a requirement for the validity of the mathematical model we use. Note that the Mueller matrix we compute does not represent the scene in any general sense since it is only representative of a particular illumination and detection condition. Second, since we need to add or subtract experimental snapshots for computing the Stokes parameters, it is important that all the images we take have the same exposure conditions. In particular, we avoided overexposed areas. Finally, the images were carefully recorded in such a way that no motion between different frames was detected.
Regarding the detection of the light, besides the need for a quarter-wave plate and a linear polarizer, a photographic camera should be used. In our particular case, we are interested in obtaining focused sharp images at different object planes. In the experiments, we use a LYTRO ILLUM light-field camera for convenience. The specs of this particular device can be found in [27]. It is worth to point out that this particular device provides 15x15 elementary views with 434x625 color pixels. Since we use a green LED, the corresponding green channel of the light-field files is used in the calculations.
Subsequently, the required 3D image processing is carried out using an open code Light Field Toolbox that enables us to use the 3D imaging capabilities with moderate computational effort [28,29]. For example, with respect to the ability to focus only on a specific object plane at a given depth, it will be almost sufficient for us to use the LFFiltShiftSum filter. This basic filter is suitable for processing the images registered by the 3D camera so as to refocus it on a single plane, reducing the intrinsic depth of field of the photographic process. Thus, by selecting the object distance, we obtain quantities 0 (L ; , ),  6). We should point out that from the images given by the lenslets of the 3D camera, one can select the frame given by a single detector of each lenslet. Thus, this would lead to a conventional 2D photographic image.

Testing the proposed model
The proposed polarimetric 3D model is tested in this section. Assuming that the pixel to pixel Mueller matrix of the scene is known, we can compare the experimentally recorded images used to compute the Stokes values of the scene with those synthetically obtained using the estimated Mueller matrix. The agreement between experimental measurements and computed quantities should be the basic argument for the acceptance or rejection of our approach. Explicitly, we need to compare the estimated polarimetric images ( , ) with the experimental ones ( , ) out i j V . The evaluation of these differences provides an estimate for the validity of the proposed procedure. For comparison, we use the experimentally measured 3D DoP ( , ) out i j and the numerically estimated 3D DoP ( , ). est i j Figs.  Table 1 or Eq. (7).
Note that small differences on the polarimetric landscape can be detected depending on the type of illumination used such as the clock on the left side facing the scene or the reflections on the stage. In addition, objects away from the focused planes do not provide any polarimetric information. For example, the characters of the label wrapped in plastic (upper left part of the image) clearly appear in Fig. 2 but are not visible in Fig. 3. est i j In order to determine how close are DoP ( , ) out i j and DoP ( , ) est i j for the six illuminations considered, we used the Structural Similarity Index Measure (SSIM) [30] and the normalized correlation coefficient. Both metrics range from 0 (completely dissimilar images) to 1 (identical images). The results are presented in Table 2. Interestingly, the similarity between the experimental polarimetric images and the computational ones is very high. In both planes, SSIM and correlation values are around 0.8 and 0.9, respectively.
In our experiments, Mueller matrix ( , ) est i j M provides a reasonably accurate account of the polarimetric behavior of the 3D scene, according to the specific recording conditions sketched in Fig. 1.   As discussed earlier, it is possible to select the 2D image frame provided by each lenslet. Using a single 2D image provided by the central lens of the microlens array, we calculated the DoP. This would provide a comparison between the 2D and the 3D imaging technique.
The results are presented in Fig. 6 for the experimental, DoP ( , ) out i j ) and Fig. 7 for the estimated DoP ( , ) est i j ). As expected, DoP ( , ) out i j and DoP ( , ) est i j are very similar for the 2D imaging. However, the amount of noise present in the DoP of 2D image is very high, and the DoP values seem to be overestimated if compared with the results shown in Figs. 2-5. 3Ddimensional InIm reconstruction is obtained by weighted averaging of the different perspective view images, and thus noise is minimized. It can be shown that InIm reconstruction is optimum in maximum likelihood sense under certain conditions [31].
Note that 2D DoP estimation can be improved if a good quality conventional camera is used. Figure 8 shows the experimental DoP ( , ) out i j using a Canon EOS 400D. The results are less noisy; however, some areas in the background where no objects are present display DoP values close to 1 for certain illumination states (L 1 and L 4 ).  In order to provide a quantitative comparison among the DOP images, we calculated the Blind Referenceless Image Spatial Quality Evaluator (BRISQUE) for the images considered [32]. BRISQUE is a widely used estimator for assessing the quality of an image. This metric evaluates the statistics of the image in order to find possible distortions. Interestingly, this descriptor does not perform any comparison with another image that can be considered perfect. Moreover, the resulting values range from 0 (best quality) to 100 (worst quality). Other common quality metrics such as the peak signal-to-noise ratio requires a reference image for calculating the mean squared error. Table 3 shows the BRISQUE values obtained. It is apparent that (i) 3D InIm DOP distributions display better image quality than the corresponding single frame ones (columns #2 and #3 versus #4, and columns #6 and #7 versus #8). (ii) The image quality of the distributions generated using the conventional camera (columns #5 and #9) is better than the corresponding ones generated from a single frame (columns #4 and #8). Nevertheless, any comparison between these two set of values is not possible because of the different specifications of both devices. (iii) Finally, it is worth to point out that BRISQUE figures for numerically estimated DOPs (columns #6 to #9) are always lower that those experimentally estimated (columns #2 to #5).

Estimation of the DoP using synthetic light
The Mueller matrix ( , ) est i j M can be used to compute the expected output Stokes parameters ( , ) out i j S from any arbitrary input illumination inp S by simply using Eq. (1). The Mueller matrix formalism can be used for producing synthetic light with arbitrary states of polarization. This includes any state of polarization located on the surface of the Poincaré sphere (elliptical polarization) or inside the sphere (partially polarized light). Sometimes, these polarization states are difficult to be generated in real life experiment because special equipment is required. For instance, in order to produce partially polarized light it is necessary to combine non-polarized and totally polarized light. This means that the experimental illumination system might become more complicated. However, once the Mueller matrix of the scene is determined, it is possible to calculate the full polarimetric information corresponding to an illumination source with any arbitrary state of polarization. Note that this is a general result and may include even structured illumination, i.e. sources with a Stokes input vector ( , ) inp i j S that might change at every point of the wave-front.
In the examples presented in Fig. 9, we calculated the DoP ( , ) est i j at plane F for a set of The last of the image series displayed in Fig. 9 shows a resulting DoP ( , ) est i j image which is a combination of the six partially polarized sources inp S obtained by fusion. Each pixel of this distribution is calculated as the maximum value of the DoP ( , ) est i j set. This means that all the pixels of the scene with a high capability of polarizing light are shown in this image, irrespective of the values of the illumination state inp S . Image fusion algorithms might be an interesting approach for displaying polarimetric information within the context of object occlusion. Fig. 9. Numerical estimation of DoP ( , ) est i j at plane F using synthetic polarized light.

Conclusion
In this paper, we extended the Mueller matrix imaging concepts to the 3D integral imaging field. We have shown that the use of polarimetric light field information poses several advantages when compared with conventional 2D polarimetric imaging.
We have shown that it is possible to estimate polarimetric information for any single selected plane. We have illustrated the fact that 2D polarimetric images are noisier than 3D reconstructed integral imaging. This is due to the fact that integral imaging computations involve weighted averaging and may be statistically optimum in maximum likelihood sense, thus noise is substantially reduced. Also, the 3D Mueller Matrix is used to produce computational polarimetric information using simulated synthetic light sources. This is particularly useful if the required state of polarization is difficult to be generated. We have shown illustrative examples with natural, partially polarized, and synthetically-produced light. Interestingly, the results obtained in these conditions can be combined with image fusion techniques in order to produce polarimetric maps that cannot be obtained by using a single source. The proposed approach can have broad applications in micro and macro scale 3D integral imaging [33][34][35][36][37].