We used the dual rotating retarder polarimeter to measure the Mueller matrix in this study [33]. Here, a commercial transmission microscope (L2050, Liss Optical Instrument Factory, Guangzhou, China) was upgraded to a Mueller matrix microscope, whose schematic and photo are shown in Figure 1. We added a polarization state generator (PSG) and a polarization state analyzer (PSA) to the optical path [14]. The key components of PSG and PSA are the fixed polarizers (P1 and P2, extinction ratio 1000:1, Daheng Optics, China) whose axes are placed along the horizontal direction and the rotatable quarter-wave plates (R1 and R2, Daheng Optics, China). An LED lamp (3W, 633 nm, Δλ = 20 nm, Cree, China) was adopted as the illuminant for the microscope. The PSG generates light with different polarization states to transmit through the tissue section and objective lens (Lens, Olympus, Japan). After modulated by the tissue sample, the polarization states of transmitted light are analyzed by the PSA, and the results are received and imaged by a CMOS camera (MVCA023-10UM, 1920 × 1200 pixels, Hikvision, China).

Before the Mueller matrix measurement, we calibrated the microscope using the model based on Ref. [34] and the eigenvalue calibration method to ensure the maximum errors of the measured Mueller matrix elements are less than 1%. We used polarizers and quarter-wave plates in different directions of 0°, 45°, 90°, and 135° as the reference samples for measurement and calibrations. More details of the calibration procedure and results can be found in Ref. [35].

During the measurement, as shown in Figure 1, the two retarders (R1 and R2) rotate at constant rates, and the rotation of R2 is five times faster than that of R1. In this configuration, the two retarders are rotated harmonically and obtain 30 intensity images at discrete intervals in rotation angle. The images can be expanded in the Fourier series shown in Eq. 1, where I is the light intensity captured by the CCD each time, α_n and β_n are the Fourier coefficients, and ω is the rotating angle rate and set to 6° in our experiment [36]. As the rotation ratio is set to be 5:1, the Fourier expansion series are limited to 12 harmonics. There is a coefficient of 2 before nω in Eq. 1 because of the repetition of polarization elements during one complete revolution. I = α 0 + ∑ n = 1 12 ( α n cos ⁡ 2 n ω + β n sin ⁡ 2 n ω ) , (1)

After getting the Fourier coefficients, the Mueller matrix elements of the sample can be obtained through the detailed procedure described in Ref. [34].

The Mueller matrix is a comprehensive description of polarization properties of turbid media and an important tool for revealing the microstructural and even nanostructural information of biomedical samples [37–40]. However, the Mueller matrix elements can be affected by multiple structures simultaneously. Thus, the physical meaning and associations with certain structural information of individual Mueller matrix elements are unclear [11–13,28]. To deal with this problem, several Mueller matrix analyzing methods have been proposed to derive groups of polarization parameters, which can reflect the specific structural characteristics of tissues. Among those methods, the Mueller matrix polar decomposition (MMPD) proposed by Lu and Chipman is widely acknowledged [23,29,41,42], while the Mueller matrix transformation (MMT) method has also shown broad biomedical application prospects [30,41,42].

Previous studies demonstrated that the MMPD and corresponding MMT parameters have a good correlation [43–45]. Thus, they can be used to reveal characteristic structural features, especially linear birefringent structures of tissue samples [25,27,28]. In this study, for quantitative analysis of skin tissue samples, we used the MMPD technique, whose basic idea is to decompose a Mueller matrix into the product of three basic 4 × 4 matrices to represent basic light-medium interaction processes, namely, depolarization (Δ), retardance (R,) and diattenuation (D) shown in Eq. 2 [29]. In Eq. 2, M _D represents the diattenuation submatrix, M _R represents the retardance submatrix, and M _Δ represents the depolarization submatrix. In addition, D in submatrix M _D is a 3 × 1 vector and can be calculated directly from the Mueller matrix: D ₁ equals M ₁₂ , D ₂ equals M _13, and D ₃ equals M ₁₄ (D _i is the element of D vector, M _ij is the element in i row j column of Mueller matrix). m _D is 2 × 2 submatrix of M _D , m _R is 2 × 2 submatrix of M _R and m _Δ is 2 × 2 submatrix of M _Δ . 0 is the 3 × 1 zero vector and P _Δ is a 3 × 1 vector in M _Δ . More detailed calculation processes can be referred to in [29]. M = M Δ M R M D , M D = [ 1 D → T D → m D ] , M R = [ 1 0 → T 0 → m R ] , M ∆ = [ 1 0 → T P ∆ → m ∆ ] , (2)

Considering that the main polarization property induced by skin follicle structures is linear birefringence, we adopt the linear retardance parameter δ and its orientation parameter θ for further analysis. The two parameters can be calculated from the total retardance M _R matrix as shown in Eq. 3, where a ₁ and a ₂ are the vectorial elements of retardance [31]. m _R is the lower right 3 × 3 sub-matrix of M _R , ε _ijk is the Levi-Civita permutation symbol, and M _R (i,j) are the elements of M _R . R = cos − 1 [ t r ( M R ) 2 − 1 ] ,   a i = 1 2 s i n R ∑ j , k = 1 3 ε i j k ( m R ) j k δ = cos − 1 { [ ( M R ( 2,2 ) + M R ( 3,3 ) ) 2 + ( M R ( 2,3 ) − M R ( 3,2 ) ) 2 ] 1 2 − 1 } , (3) θ = 1 2 tan − 1 [ a 2 a 1 ]

In a grayscale image, the regular gray value distributions of all pixels compose the image texture. Each pixel has some correlation with other pixels with a certain distance. Due to grayscale discrepancy, the spatial correlation characteristics between pixels are generated. We can obtain a series of image texture characteristic values by analyzing the discrepancy of grayscale values among the pixels. By combining the sliding window method and some prevalent image texture analysis techniques, such as the Gray Level Co-occurrence Matrix (GLCM) [46], we can obtain texture information from the Mueller matrix parameter images. Specifically, we use a sliding window with certain size on the Mueller matrix-derived parameter images to calculate the image texture parameter images, which retain both the position and image texture information [47]. The sliding window method is similar to convolution at a conceptual level. By sliding the image with a certain length and width core, we can turn the image into a series of small units for further processing. For instance, the GLCM is a square matrix, but the images are usually rectangular arrays. The sliding window method can transfer the rectangular image into a series of square images with changeable size and step size for further texture analysis.

The GLCM method proposed by Haralick et al has been widely used for texture feature extraction and image texture analysis [46,48–50]. The texture of a grayscale map is formed by repeated occurrence of different grayscale values in a certain pattern. Thus, there is a certain gray relationship between the two pixels that are separated by a distance in the image space, which is the spatial correlation between the gray values in the image. The GLCM describes texture features by calculating the spatial correlation of the image gray level. We select four representative parameters, namely, contrast, energy, correlation, and inverse different moment (IDM) for follicle structure detection.

(1) Contrast is defined as Eq. 4:

The energy parameter can be obtained by summing the square values of each element in the gray level co-occurrence moment. The parameter reflects the image's gray distribution uniformity and texture fineness. If the texture in the original image changes equably, that is, the gray scale changes equably, the calculated energy has a large value.

(3) Correlation is defined as Eq. 6:

The correlation reflects the consistency of the image texture and the gray correlation between pixels. When the grayscale difference between pixels is small, the calculated correlation value is large. On the contrary, if the gray values of adjacent elements have a huge difference, the calculated correlation value is small.

(4) IDM is defined as Eq. 7:

The IDM reflects the degree of local texture change. When the gray values between pixels are equal, the minimum value of the denominator in the formula mentioned above is 1, and the calculated inverse moment is large. Therefore, a small IDM value means a large change in the image texture, while a large value shows a small change in the image's local texture.

Another image texture analysis method used in this study is the Tamura Image Processing Method (TIPM) proposed by Tamura [51]. Compared with the GLCM, the Tamura texture features are more in line with human texture visual perception. The TIPM mainly includes six texture features: coarseness, contrast, directionality, line-likeness, regularity, and roughness. Here, we mainly used the contrast feature, whose calculation procedure is shown in Eq. 8, where μ ₄ is the fourth-order moment about the mean and σ ² is the variance. F c o n = σ α 4 1 / 4 , α 4 = μ 4 σ 4 , (8)

The skin tissue samples used in this study are the 6 μm thick, dewaxed, unstained common, and healthy mouse dorsal skin tissue sections provided by the Experimental Research Center, China Academy of Chinese Medical Sciences. We selected sagittal section samples of mouse skin tissue to obtain complete skin hair follicle information. The sagittal section surface is parallel to the skin surface, enabling better visualization of the hair follicle structure. We prepared three groups of mouse skin tissue samples as shown in Figure 2. Before Mueller matrix imaging, the skin tissue sections were evaluated by an experienced pathologist, and the locations of hair follicles were identified. The study was approved by the Ethics Committee of the China Academy of Chinese Medical Sciences.

Figure 2 shows the imaging results of three mouse skin tissue samples. We can observe from Figures 2A–C of grayscale light-intensity images that the skin hair follicles have annulus structures as marked by the red circles. However, there are some other structures and dust that may also show similar imaging results to those of the hair follicles, for instance, in the lower left corner of (B), which can hardly be distinguished using the unpolarized intensity images, whereas it can be noticed from Figures 2D–F that the Mueller matrix elements, especially the lower right 3 × 3 elements, show characteristic intensity changes in the areas of hair follicles, suggesting the prominent linear retardance generated by the hair follicle micro- and nanostructures in the skin tissues [23,37,41].

To observe the birefringent hair follicle structures, we then calculated the MMPD parameters δ and θ of the three mouse skin tissue sections as shown in Figure 3 to represent the magnitude and orientation distribution of the linear birefringence of the sample, respectively [31]. Specifically, the MMPD δ parameter images shown in Figures 3A–C demonstrate that the hair follicles have circularly aligned structures around the center with larger linear retardance values than other skin structures. It should be mentioned that, as shown in Figure 3C, the hairs also represent strong linear retardance properties. Meanwhile, in the MMPD θ parameter images shown in Figures 3D–F, the hair follicles represent periodically changed colors, confirming the continuous orientation variations of the birefringent structures around the center of skin hair follicles. These periodic birefringent structures encircle the center of the follicle in circular distribution, and each periodicity is filled with a whole circle. Figures 3G–J also provide the zoom-in images of parameters δ and θ of two selected hair follicle regions, marked by the white squares in (A) and (B). Obviously, the values and distribution behaviors of the linear retardance induced by the characteristic microstructures can be used for hair follicle location identification. In addition, we can see that the structural integrity also contributes to retardance values, indicating that the MMPD parameters can be used as indicators for quantitative hair follicle health evaluations.

In Polarization staining images of skin hair follicles, we showed that the location and structure integrity information can be obtained by the MMPD parameter δ and θ images. Previously, we proposed a polarization staining method to combine different polarization parameter images into a single one to provide multi-dimensional information conveniently [52]. Here, in this study, the parameters δ and θ are projected into the HSV color space to generate polarization staining images shown in Figure 4, in which the brightness intensity encoded by parameter δ describes the linear retardance value, and the hue encoded by parameter θ shows the fast axis orientation of the linear birefringent structures. We can see from the results of all three samples shown in Figures 4A–C that the skin hair follicle's locations are clearly revealed: The larger the value of linear birefringence, the brighter the polarization staining images. Compared to the parameter θ images as Figure 3 shows, the hue weighted with the intensity images can more clearly represent the circularly distributed orientation of the hair follicles as shown in Figure 4. The results confirm that we can observe the structural features more directly using polarization staining images.

In previous sections, we first demonstrated that the MMPD parameter δ and θ images can be used to characterize the mouse skin hair follicle microstructures, and then we visualized hair follicle structural features by the polarization staining images based on MMPD parameters. In this section, for hair follicle location identification and quantitative structural evaluations, we apply the image segmentation method, sliding window algorithm, GLCM, and TIPM parameters [46,51] to process and analyze the polarization images.

To identify the hair follicle locations, we extract their structural features based on polarization parameter images first. Using the traditional image segmentation method of double-down sampling–binarization–etching–expanding–closing–hole filling, we can obtain the corresponding hair follicle area from the MMPD δ parameter image as Figure 5 shows. Figure 5 (D1)–(F1) are the segmentation results using original MMPD δ parameter images of three mouse skin samples shown as Figure 5 (A1)–(C1) at the reduction rate of four-fold down-sampling. The white areas are the segmented hair follicle regions whose pixel values are 1, while other regions are with pixel values of 0. We can notice that by eliminating the influence of the background noise in the polarization parameter images using a relatively high reduction rate, the location and number of the skin hair follicles can be identified using the image segmentation method. More specifically, it is expected that the size and number of the continuous areas in the images can be used to locate and count the hair follicles in the skin tissue.

For better identification results, here, we also use the sliding window method proposed in our previous study to show the local texture feature of the polarization parameter images [47]. We select TIPM and GLCM parameters as introduced in Polarization staining images of skin hair follicles, which present prominent differences between hair follicles and surrounding structures in MMPD δ and θ images, as the typical texture feature parameters to show the sliding window method results. The sliding window size was set to 30 × 30 pixels, and the step size was 10 pixels. Figure 5 shows the sliding window images of TIPM contrast (A2)–(C2), GLCM 1-energy (A3)–(C3), GLCM contrast (A4)–(C4), GLCM correlation (A5)–(C5), GLCM 1-IDM (A6)–(C6), and their corresponding segmentation results after etching–expanding–closing–hole filling processes. The sliding window method used here is based on a threshold of 0.5. As shown in Figure 5, the textural parameter values in each local region were calculated to build a new 2D image according to their original spatial position. It can be noticed that compared to the segmentation results using original MMPD δ parameter images shown in Figure 5D1, F1), the sliding window method based on different textual images can help suppress the influence of background noise. Thus, the locations of hair follicle regions can be more accurately identified as shown in Figure 5. From Figure 5D1, F1, we can see that the segmentation results based on the original MMPD δ parameter images are not suitable for hair follicle location, counting, and size calculating. It is supposed that since the contrast of the original MMPD δ parameter images is low, it does not fit the threshold. Thus, using different texture analysis methods to extract image features and improve image contrast is necessary. More specifically, among all the textural parameters, the segmentation results based on TIPM contrast (A2)–(C2), GLCM 1-energy (A3)–(C3), GLCM contrast (A4)–(C4), and GLCM 1-IDM (A6)–(C6) show good identifications, and all the hair follicles can be correctly recognized through these parameters, while those based on GLCM correlation (A5)–(C5) have low accuracy. It is supposed that the GLCM correlation is too sensitive to the texture changes; thus, the results mixed with noise are unrelated to hair follicles. This leads to lower accuracy here.

After obtaining the regions of hair follicles in the mouse skin tissue samples using polarization texture parameter image segmentation, we can then quantitatively evaluate the size of the follicle structures. In this section, all the hair follicles in the original images were confirmed and manually measured by experienced researchers from China Academy of Chinese Medical Sciences. It should be pointed out that, in this study, we are trying to provide an automatic polarimetric method for the researchers, who currently use manual measurement to detect and evaluate skin hair follicles, as a screening tool. Hence, here we compare our method to that of the experienced researchers’ measurement as the reference. Figure 6 shows the comparison between manual size measurement (long axis diameter) and calculation based on different polarization texture images of 23 randomly selected hair follicles via correlation scatter plots. Here, we chose the TIPM Contrast, GLCM Contrast, GLCM 1-DIM, and GLCM 1-Energy, which show good structure identification abilities in Figure 5, to acquire the size information (long axis diameter) for hair follicle structures. Specifically, in Figure 6, the manually measured hair follicle size values were obtained from the original unpolarized intensity images, and the calculated size values were extracted based on the segmentation results of corresponding sliding window polarization texture parameter images. Figure 6 also gives the coefficient of determination (R ² ) and the p value between the manually measured and calculated values for each polarization parameter. From Figure 6, we can observe that the sliding window images of TIPM Contrast and GLCM 1-IDM based on the MMPD δ parameter can provide relatively accurate size estimation results, whose R ² are more than 0.25, which means their Pearson’s correlation coefficients are larger than 0.5, showing moderate correlation. Also, the p values of TIPM Contrast and GLCM 1-IDM are less than 0.01, which strengthens the correlation significance. However, the results derived from sliding window images of GLCM Contrast and GLCM 1-Energy can reach a determination coefficient of less than 0.2, indicating they may not be precise size indicators. The relatively low R ² values may have resulted from the manual measurement. In the manual measurement process, the artificial experience has a prominent impact on the final measured values. The error resulting from manual measurement would cause the measured values of hair follicles of similar size to be distributed in a large range, while the calculated values from the algorithm are relatively stable. Thus, the results obtained from TIPM Contrast and GLCM 1-IDM based on the MMPD δ parameter have a lower standard deviation than those of the measured value. This can also be seen in Figure 6 as the data points seem to be aligned along horizontal lines. Moreover, in the case of a large distribution difference, there is still a moderate correlation between measured values and calculated TIPM Contrast values and GLCM 1-IDM values, which shows the effectiveness of the two methods.

Table 1 provides the mean relative error between manually measured size and calculated size values of 23 hair follicles for polarization texture parameters TIPM Contrast, GLCM Contrast, GLCM 1-DIM, and GLCM 1-Energy. The average relative error can be used to determine accurate parameters for segmentation and quantitative evaluation of hair follicle structure. The values shown in Table 1 confirm that among the four parameters, the TIPM Contrast and GLCM 1-DIM parameters have better performance than GLCM Contrast and GLCM 1-Energy, whose mean relative errors are 25 and 24.2%, respectively. It should be noted that to achieve good hair follicle structure identification and evaluation using the TIPM or GLCM parameters, the Mueller matrix microscope should be calibrated and samples should be prepared carefully. The results provided in Table 1; Figure 6 indicate that the Mueller matrix imaging method combined with polarization image segmentation, texture analysis, and calculation is expected to realize automatically hair follicle identification, detection, and quantitative evaluation. The information provided by Mueller matrix polarimetry shows great potential in dermatological studies and applications.