Calibration of a collinear backscattering Mueller matrix imaging system

A collinear backscattering Mueller matrix (CBMM) imaging system has clear advantages in the detection of bulk biological tissues, which are highly scattering and depolarizing. Due to the double-pass configuration and noise in the system, the calibration of a collinear backscattering Mueller matrix imaging system is usually complex and of poor accuracy. In this work, we propose an alternative modified eigenvalue calibration method (ECM) based on the equivalent standard sample. For better noise suppression and higher calibration accuracy, we design the distribution of polarization states over the Poincaré sphere and solve for the parameters of equivalent standard samples by means of an optimization. Compared to other variants of the eigenvalue calibration method used in the double-pass system, the accuracy of the proposed method is improved by more than 40 times. The comparison results with the error model-based calibration methods indicate that the modified eigenvalue calibration method generally gives the best accuracy and precision, as well as the best reliability.

easy due to its simple optical system. One of the most popular calibration methods used to calibrate the TMM system is the analysis calibration method (ACM), which builds an error model with five main systemic errors and obtains the error values by the Fourier analysis method [12][13][14]. The calibration of a CBMM system is more complex due to the variety of error sources, for example, the parasitic polarization of the beam splitter also affects the SOP of light. In our previous work, to calibrate the CBMM system we proposed a numerical calibration method (NCM) that takes into account the parasitic polarization of the beam splitter and solves the system errors numerically by Levenberg-Marquardt (LM) algorithm [15]. However, the systematic error model-based calibration methods have the following drawbacks, they require a priori knowledge of the system and are prone to over-calibration when the systematic error model is overly complex.
Different from the error model-based calibration methods, the eigenvalue calibration method (ECM) does not require any priori information such as the polarization states produced in the system [16,17]. As a general method for the calibration of Mueller matrix imaging systems, the ECM has shown an increase in popularity [18]. However, due to the complex double-pass configuration, the original ECM is not suitable for the calibration of the CBMM system and some modifications of the original ECM are necessary [19]. Besides, the ECM is inefficient in the case of low signal-tonoise ratio (SNR) since it assumes that there is no noise in the measurement system [20].
In this paper, we propose an alternative modified ECM for the CBMM system. Unlike other variants of the ECM, we simplify the double-pass calibration to a single-pass calibration by the mirror symmetry and the concept of the equivalent standard sample; it is much simpler in theory as it does not need to address issues such as matrix reciprocity and eigenvalue squaring. For better calibration performance, we optimize the distribution of analysis states over the Poincaré sphere taking into account the inherent constraints of polarization modulators, and calibrate the parameters of equivalent standard samples by means of an optimization. To compare the performance of the modified ECM and the error model-based calibration methods, a series of experiments with various samples are performed.

Experimental setup
The optical path diagram of the CBMM system is shown in Figure 1. A LED (center wavelength = 610 nm and FWHM = 20 nm) is used as the light source in the system. The output beam from the LED is collimated by an objective lens (L1) and a convex lens (L2). P1 and P2 are linear polarizers (LPVISE100-A, 400-700 nm, Thorlabs, Inc., U.S.). R1 and R2 are zero-order quarter-wave plates (GCL-060402, 633 nm, Daheng, Optic, China). The polarization state generator (PSG) consists of a fixed polarizer (P1) and a rotating quarter-wave plate (R1), while the polarization state analyzer (PSA) consists of a rotating polarizer (P2) and a rotating quarter-wave plate (R2). Wave plates R1, R2 and polarizer P2 are driven by three servo motors (PRM1Z8E, Thorlabs, Inc., U.S.), respectively. The polarization-modulated beam is reflected by the sample at normal incidence after passing through a beam splitter (GCC-403102, Daheng, Optic, China), and then the backscattered light is imaged by a charged couple device (CCD) camera (PCO. panda 4.2, 16bit, Inc., Germany) after passing through the PSA and a low numerical aperture (NA) objective lens (L3).

The modified eigenvalue calibration method
Previous studies on the eigenvalue calibration method in the double-pass system needed to consider the complex reciprocal relationship between the Mueller matrix of the reflection mirror and the Mueller matrix of the sample [19,21]. In this work, we propose an alternative modified ECM that effectively simplifies the calibration theory by utilizing the mirror symmetry and equivalent standard samples.
In the absence of depolarization, the Mueller matrices of polarizers and wave-plates can be written uniformly as follows [22]: where τ is the transmission coefficient, φ and Δ are ellipsometric parameters, θ stands for the azimuth angle, and R(θ) is the rotation matrix. Assuming that the observer is at the detector and looks in the direction from which the light is incoming, the counterclockwise direction of rotation is chosen as the positive

FIGURE 1
Optical path diagram of the collinear backscattering Mueller matrix imaging system.
Frontiers in Physics frontiersin.org direction of the azimuth angle θ. In addition, it is necessary to assume that the light is reflected by a perfect flat mirror at normal incidence and the polarization properties of optical elements are independent of the direction of propagation of the light. The mirror Mueller matrix M m is actually the mirror transformation matrix satisfying M m M m E, where E stands for the identity matrix. It is obvious that any left rotation transformation can be equated to a right rotation transformation in mirror space, and vice versa. According to this simple property, the mirror transformation of an arbitrary Mueller matrix described by Eq. 1 can be written as follows, Given the properties of mirror symmetry, the double-pass process for an arbitrary standard sample can be simplified to the left multiplication of Mueller matrix M m with an equivalent Mueller matrix M′, where M f and M b represent the Mueller matrix of the forward and backward propagating path of a standard sample, respectively. The parameters of the equivalent Mueller matrix satisfy the following relationships: τ′ τ 2 (1 + cos 2 2φ), θ′ θ, φ′ 1/2asin[(1 − cos 2 2φ)/(1 + cos 2 2φ)], and Δ′ 2Δ, which means that the equivalent standard sample has the same azimuth angle as the standard sample and the double phase retardance. It is clear that Eq. 3 transforms the double-pass process into a single-pass process by using the equivalent Mueller matrix, hence in the next theoretical analysis only the equivalent standard sample needs to be considered. Similar to the calibration process for single-pass system, air, ± 45°l inear polarizers, and a 30°λ/8 wave plate whose equivalent Mueller matrix is a quarter-wave plate are chosen as the four standard samples. The detected intensity matrix I i of the i th equivalent standard sample can be described by the following equation, where W and A are the polarization generation matrix and the polarization analysis matrix, respectively, and M ′ i stands for the equivalent Mueller matrix of the i th standard sample. Obviously, the equivalent Mueller matrix of air is still an identity matrix, therefore I 1 can be simplified as, Define matrix C i as follows: where superscript ' †' denotes the Moore-Penrose pseudo-inverse. Notice that C i and M ′ i are similar matrices in this modified ECM, which means that they share the same eigenvalues. Hence, it is possible to partially reconstruct the four equivalent Mueller matrices using the following eigenvalue relationships [16], where λ i stands for the eigenvalue of matrix C i . Nevertheless, it is found that the reconstruction of the equivalent Mueller matrices by Eq. 7 may be inefficient in the actual experiment, because the process of computing ratios of eigenvalues is highly noise sensitive. An optimization algorithm is introduced in the next section to exactly reconstruct each equivalent Mueller matrix, and the results given by Eq. 7 can be used as a starting point for the optimization algorithm. Left multiplying Eq. 6 on both sides by matrix W, one obtains: where the superscript 'T' represents the transpose. To uniquely determine the vec(W), a positive semidefinite Hermitian matrix L is constructed by using the H i matrix of each equivalent standard sample, Theoretically, the matrix L should only have one null eigenvalue and the corresponding eigenvector is vec(W). Once polarization generation matrix W is determined, polarization analysis matrix A can be easily obtained by the pseudo-inverse matrix of W as follows: After calibration, the backscattering Mueller matrix of an arbitrary sample can be calculated directly from the detected intensity matrix of sample (I sample ) as shown below:

Optimization of the modified ECM
The ECM assumes that all the Mueller matrices of standard samples are perfectly known. However, due to noise in the system, the reconstruction of the equivalent Mueller matrices by calculating the ratios of eigenvalues becomes inefficient; in fact, even if the whole system is perfectly known, the noise still causes the SOP of light to deviate from the theoretical value, thus decreasing the calibration performance. For higher calibration accuracy, two optimizations for the modified ECM are introduced: designing a reasonable polarization state analyzer (PSA) to restrain noise and calculating the parameters of the equivalent Mueller matrices by a noise insensitive method.  [27]. Nevertheless, the polarization modulation module of a CBMM system based on dual-rotating retarders is usually composed of a fixed polarizer and a rotating quarter-wave plate, which results that the modulated polarization states cannot cover the entire sphere but only form an 8-shaped curve on the surface of Poincaré sphere. In addition, a non-ideal quarter wave plate causes the 8-shaped curve incapably crossing the two poles of Poincaré sphere. In the experiment, the phase retardance of the quarter-wave plate used in the PSA is approximately 96°. Taking into account the inherent limitations of polarization modulators, we design the distribution of polarization states on the surface of Poincaré sphere for a fixed polarizer (FP) based PSA and a rotating polarizer (RP) based PSA, respectively.
Two quantitative indicators the condition number (CN) and the equally weighted variance (EWV) are used to estimate the performance of PSA [28,29], where the μ j is the singular value of a matrix A. Theoretically, the optimal distribution of polarization states with constraints on the surface of Poincaré sphere can still be obtained by minimizing the CN or EWV with the genetic algorithm (GA) [30]. Considering the two constraints of a fixed polarizer and a non-ideal waveplate, part of the optimal modulation frames for a FP-based PSA and a RPbased PSA are shown in Figure 2. It is found that the optimal results given by the CN and EWV are consistent with each other. For a FPbased PSA, the optimal frames are irregular polyhedra with small volumes since the modulable polarization states are restricted to the 8-shaped curve on the surface of Poincaré sphere, as shown in Figures 2A-C. And for a RP-based PSA, the polarization states in the polar regions on Poincaré sphere become unmodifiable as the phase retardance of the quarter-wave plate deviates from 90°, as shown in Figures 2D-F. Generally, the influence of a non-ideal quarter wave plate on the optimal results is insignificant, the optimal modulation frames are still inscribed regular polyhedra with only the optimal frame for RP-6 being slightly tilted.
To compare the performance of different optimal modulation frames in the Mueller matrix measurements, a homogeneous wave plate is used as the sample and its Mueller matrix image is measured under different polarization analysis modes. Then the phase retardance image of the sample for each of optimal modulation frames is calculated by Mueller matrix polar decomposition (MMPD) method, respectively [31]. The SNR calculated for each phase retardance image is shown in Table 1. The CNs of the RP mode can reach 3 √ , which is considered as the theoretical minimum CN for an optimal PSA [32], while the CNs of the FP mode cannot reach this limit due to the inherent constraint of the FP mode. For the same N, the optimal modulation frame of the RP mode has smaller optimization indicators, i.e., the CN and the EWV, and a higher SNR, thanks to the higher degree of freedom of , it can be seen that optimizing the EWV is effective in improving the SNR of the measurement results despite the fact that the CN remains the same, which indicates that the EWV is a better indicator. More importantly, Table 1 proves that optimizing the distribution of polarization states on the surface of Poincaré sphere can effectively suppress noise and improve the SNR of measurement results, even under an optimally constrained condition.
In the previous calibration study, the ellipsometric parameters of standard samples are obtained by calculating the ratios of the eigenvalues of the matrix C i according to Eq. 7. Nevertheless, due to noise in the system, it is found that such a Mueller matrix reconstruction method for equivalent standard samples is inaccurate; especially the transmission coefficient τ ′ i and the retardance Δ ′ i are calculated with poor accuracy. Assuming that the transmission coefficients of the ± 45°linear polarizers are the same, i.e., τ 2 ′ τ 3 ′ , and their phase retardances can be ignored, i.e., Δ 2 ′ Δ 3 ′ 0, there are only 6 parameters of equivalent standard samples that need to be accurately calibrated: one phase retardance Δ 4 ′ , two transmission coefficients τ 2 ′ and τ 4 ′ , and three azimuth angles θ 2 , θ 3 and θ 4 . In the absence of noise, if all parameters of equivalent standard samples are correct, all the eigenvalues of the Hermitian matrix L are non-zero except one; on the other hand, if the reconstruction of the equivalent standard samples is inaccurate, the system has only a trivial solution and all the eigenvalues of L are nonzero [16]. Owing to noise in the system, the minimum eigenvalue of the system is no longer zero, however the eigenvalues still satisfy the following relationship that the more accurately the equivalent standard samples are reconstructed, the smaller the minimum eigenvalue of L is with respect to the 15 other eigenvalues [16]. Therefore, it is possible to calculate the 6 parameters of equivalent standard samples by means of an optimization with the result given by Eq. 7 as a starting point. Similar study has shown that the reconstruction of matrices by means of an optimization has better noise resistance than the direct calculation of the ratios of eigenvalues [20]. In the modified ECM, we adopt the following fitness function, where λ 16 is the smallest eigenvalue and λ 15 is the second smallest eigenvalue of matrix L, respectively. The fitness function can be easily solved by various optimization algorithms such as the genetic algorithm. Hence, the Mueller matrix reconstruction for equivalent standard samples is transformed from calculating the ratios of the eigenvalues of the matrix C i to an overall optimization for the matrix L.
For the same set of experimental data, we compared the accuracy and relative error (RE) of the parameters of equivalent standard samples calculated by different methods, as shown in Table 2. It can be seen that the RE calculated by the previous method is as high as 5%, while the RE calculated by means of an optimization is less than 1.6%, only 1/3 of the former. In other words, the accuracy of the proposed method to reconstruct the Mueller matrix of the equivalent standard sample is three times that of the previous method.
The accuracy of the equivalent standard sample has a direct impact on the final accuracy of the calibration. The comparison results of the eigenvalues of Matrix L calculated with different methods are shown in Figure 3. Thanks to the more accurate Mueller matrices of equivalent standard samples, the minimum eigenvalue of the proposed method is 1.45E-05, which is two orders of magnitude lower than that of the previous method. In ECM, the error estimator ε λ 16 /λ 15 is commonly used to assess the actual relative error of the entire system after calibration [31]. The error estimator ε of the previous method is 0.72, while that of the new method is 0.016, which also means that the calibration accuracy of the proposed method is more than 40 times higher than that of the previous method.

Results and discussion
To compare the calibration performance of the modified ECM with the error model-based methods such as the ACM and NCM, we measured several samples including a polarizer, a wave plate, and a well-aligned silk sample. For the modified ECM, we adopt the RP-8 polarization state analysis mode and reconstruct the Mueller matrices of equivalent standard samples by means of an optimization with the genetic algorithm. For the ACM and NCM, the quarter-wave plates in PSG and PSA rotate synchronously with a fixed ratio of angular velocity 1:5 and a total of 30 measurements are required to calculate the Mueller matrix image of the sample. During calibration with the ACM and the NCM, the two polarizers in the PSG and PSA are fixed, in other words the PSA is in the FP mode. Since the transmission Mueller matrices of transparent samples can be accurately measured in a calibrated TMM system, their true equivalent Mueller matrices which can be calculated by Eq. 3 are in fact well-known, namely M true . Therefore, it is simple to compare the differences of the To estimate the performance of different calibration methods, four accuracy indicators the average root mean square error (RMSE), the average maximum error (MAX ), the average relative diattenuation error (δD), and the average relative retardance error (δR) are defined as follows, where the "average" refers to the average over all pixels in the selected area of the Mueller matrix image, D and R represent the diattenuation function and the phase retardance function defined in the MMPD, respectively. RMSE and MAX are used to estimate the average and maximum deviation of the Mueller matrix measurement results given by different calibration methods from the true Mueller matrix M true , respectively. δD and δR are used to estimate the relative diattenuation error and the relative retardance error for different calibration methods, respectively. The Mueller Matrix images of a quarter-wave plate (GCL-060403, 488 nm, Daheng, Optic, China) and a polarizer (LPVISE100-A, 400-700 nm, Thorlabs, Inc., U.S.) are measured in the CBMM system based on three different calibration methods. The comparison results are shown in Figure 4. It is clear that the modified ECM has the best calibration performance, with the NCM the second best and the ACM the worst. For the wave plate sample, the RMSE and the MAX of the modified ECM are 0.0054 and 0.034, which are less than half and about 3/5 of the results of the NCM, respectively. Due to the absence of the beam splitter in the error model, the MAX of the ACM is close to 0.15, which indicates that the ACM is not suitable for the calibration of the CBMM system. Although the NCM takes into account the effect of the beam splitter, it is still less effective than the modified ECM, due to the deviation of the error model from the actual system and overcalibration. Calibration results for the polarizer sample with different methods are similar to the former results. In addition, the error bars of the modified ECM are also minimum, benefiting from the efficient noise suppression of the optimal RP-based PSA.
In the next experiment, a well-aligned silk sample is selected as the experimental sample, as shown in Figure 5. Because of its highly scattering property and typical birefringent nature, the silk sample is an ideal object of study for the Mueller matrix imaging system [34]. As the illumination light and detection light are collinear in the CBMM system, several polarization parameters independent of the azimuth angle of the sample, i.e., rotation invariant indicators can be derived from the Mueller matrix [35]. Obviously, if a CBMM system is well calibrated, the rotation invariant indicators separated from different azimuthal regions of a circularly symmetric sample should be constant; on the other hand, when a CBMM system is not well calibrated, these rotation invariant indicators may still be orientation dependent. Considering the circular symmetry of the silk sample, it is possible to estimate the calibration performance of different calibration methods in the view of symmetry.
In the experiment, we select three different azimuthal regions on the silk sample, i.e., ± 45°and 90°, as shown in the red

FIGURE 3
Comparison of the eigenvalues calculated by different methods.
Frontiers in Physics frontiersin.org

FIGURE 4
The calibration performance comparison of three calibration methods with different samples, (A) a wave plate sample, (B) a polarizer sample.

FIGURE 5
A well-aligned silk sample (A) and its Mueller matrix image (B). The rotation invariant indicator D L represents the linear diattenuation of the sample and r L relates to the ability of the sample to convert circularly polarized light into linearly polarized light [35]. Frequency distributions images of D L and r L obtained by different calibration methods are shown in Figure 6 [36]. It is clear that the rotation invariant curves obtained by the modified ECM for different azimuthal regions are perfectly coincident with each other, while the results of the error model-based calibration methods are partially non-overlap. In other words, the rotation invariant indicators measured by the modified ECM do not depend on the azimuth angle of samples, whereas the results given by the NCM and the ACM are still azimuth dependent. As a result, it can be deduced that only the calibration result of the modified ECM is reliable from the view of symmetry. In fact, the calibration performance of the NCM and the ACM depends on their error models, however, it is difficult for a theoretical error model to accurately describe the actual measurement system. In addition, solving an error model that is overly complex may lead to over-calibration, which affects the accuracy of the calibration. The major advantage of the modified ECM is that it does not require modelling the errors in the system, hence it is more suitable for the complex CBMM system. Thanks to the optimal modulation frames for the PSA and the optimization-based Mueller matrix reconstruction method, the modified ECM gives the best calibration performance in the CBMM system.

Conclusion
In this paper, based on the eigenvalue calibration method we propose a simple modified ECM to calibrate the CBMM system configured in double-pass. In the modified ECM, we transform the complex double-pass calibration problem into a simple single-pass form by using mirror symmetry and the concept of the equivalent standard sample. To suppress the effect of noise on calibration accuracy, we design the distribution of polarization states over Poincaré sphere considering the inherent constraints of polarization modulators, and reconstruct the Mueller matrix of equivalent standard samples by means of an optimization. It is found that an optimal PSA can effectively improve the SNR of the measurement data and the optimization-based Mueller matrix reconstruction method is more robust against noise than the direct calculation method. After optimization, the accuracy of the modified method is 40 times better than the previous variants of the ECM. The comparison results with the error model-based calibration methods indicate that the modified ECM has clear advantages in the calibration of the CBMM system with higher accuracy and superior reliability.

Data availability statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.