The use of polarimetric techniques has received considerable recent attention for biomedical diagnosis.1 This is because the polarization properties of scattered light contain additional diagnostic information on tissue, which cannot be obtained from the use of unpolarized light. For example, collagen (a structural protein) and other fibrous tissues exhibit linear birefringence due to their oriented structures. In various kinds of tissue abnormalities, the structural (and functional) properties of collagen and other fibrous structures show distinct changes. Measurement of the linear birefringence of tissue may thus find useful applications in noninvasive diagnosis of various pathological diseases.2, 3, 4, 5 Glucose is another important tissue constituent that possesses circular birefringence due to its asymmetric chiral structure. Its presence in tissue leads to rotation of the plane of linearly polarized light about the axis of propagation (known as optical rotation). Measurements of optical rotation may offer an attractive approach for noninvasive monitoring of tissue glucose levels, and several studies have investigated the polarization properties of scattered light from chiral turbid medium. 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16 In addition to these applications, polarimetric measurements are also finding use for contrast enhancement in biomedical imaging.17, 18, 19, 20, 21
The use of the polarimetric techniques for biological tissue characterization is, however, severely hampered by multiple scattering effects. In a turbid medium such as tissue, the incident polarized light becomes strongly depolarized due to multiple scattering, the magnitude and rate of depolarization being determined by a large number of parameters such as the concentration, size, shape, refractive index of the scatterers, detection geometry, and by the incident light’s state of polarization.22, 23, 24, 25 Polarization modulation and synchronous detection methods, therefore, have been developed to extract the weak fraction of the polarization preserving component of scattered light from a turbid medium and analyze it for determination of the diagnostically relevant polarization properties such as linear birefringence and optical activity.7, 9, 10, 26 In particular, this approach is being investigated to determine of the concentration of optically active glucose in human tissue. However, studies have shown that even when depolarization due to multiple scattering is taken into account, the relatively small rotation of the polarization vector arising from optical activity in the medium is hidden by the much larger change in the orientation angle of the polarization vector due to scattering.7, 27 In addition to the scattering effects, interpretation of optical rotation using the conventional Stokes polarimetry-based approach is further confounded by the effects of linear birefringence and its complex orientation in real tissues.
Because the Mueller matrix (which transforms the Stokes vector describing the polarization properties of incident light to the Stokes vector of scattered light) description contains complete information about all the polarization properties of a medium,28, 29 measurements of full Mueller matrix elements may be used to address these issues. However, when many optical polarization effects are simultaneously occurring in the sample (as is the case for tissue), the resulting elements of the Mueller matrix reflect several “lumped” effects, thus hindering their unique interpretation. A methodology that allows the extraction of the individual “intrinsic” material polarimetry characteristics from the lumped Mueller matrix description is thus needed. The independent constituent polarization parameters—linear retardance ( , difference in phase between two orthogonal linear polarization, and its orientation angle ), circular retardance or optical rotation ( , difference in phase between right and left circularly polarized light), diattenuation ( , differential attenuation of orthogonal polarizations for both linear and circular polarization states), and depolarization coefficient ( , linear and circular)—can be efficiently extracted by decomposing the measured Mueller matrix into the “basis” Mueller matrices of a retarder, diattenuator, and depolarizer.30 The individual polarization parameters obtained through this process should provide more intrinsic information on the underlying morphological features and physiological state of tissue and may thus offer a more fundamental (less empirical) way to understand its structure and behavior.
Indeed, preliminary studies have shown that Mueller matrix decomposition can be used to extract the pure optical rotation component corresponding to the concentration of chiral substances by decoupling the contribution of scattering-induced rotation of the polarization vector that primarily arises due to the contribution of singly (or weakly) scattered photons in a scattering medium.31 In a highly scattering medium such as thick tissue, however, the contribution of singly scattered photons is expected to be small and the majority of the detected photons will be multiply scattered. Further, unlike the ideal situation, in actual tissue, all the polarization effects are exhibited simultaneously and not in an ordered, sequential manner. To realize the potential of the Mueller matrix decomposition approach for quantification and interpretation of the useful polarization parameters (linear retardance and optical rotation in particular) in real tissues, it is necessary to test the validity of the decomposition process on multiply scattering media exhibiting simultaneous linear birefringence, optical activity, and depolarization. A quantitative understanding of the confounding effects of multiple scattering on the values for the polarization parameters estimated through the decomposition process is also essential.
We therefore investigate these aspects both theoretically and experimentally. The effects of multiple scattering, the propagation path of scattered photons, and detection geometry on the values for and were studied systematically using polarization-sensitive Monte Carlo simulations of the propagation of polarized light in birefringent, chiral, multiply scattering media having varying optical properties. In the simulations, the simultaneous effects of linear birefringence and optical activity on the photons between successive scattering events were incorporated using the Jones N-matrix formalism.32 In the corresponding experimental studies, Mueller matrices were recorded from realistic tissue phantoms exhibiting controllable scattering, linear birefringence, and optical activity. The values for the polarization parameters estimated through the decomposition process were compared with the known optical properties of the phantoms. The results of these theoretical investigations and the experimental studies are reported in this communication.
The organization of this paper is as follows. Section 2 provides a description of the experimental systems and methods as well as the optical phantoms used in this study. The forward polarization-sensitive Monte Carlo simulations of propagation of polarized light in multiply scattering media and the inverse process for polar decomposition of the resulting Mueller matrix are described in Sec. 3. The results of decomposition of the experimental and the Monte Carlo–generated Mueller matrices are presented in Sec. 4 to validate the method and to elucidate the observed trends. Section 5 concludes the paper with a discussion on the implications of these results and potential use of the Mueller matrix decomposition methodology in diagnostic photomedicine.
Experimental Methods and Materials
The Mueller matrices were recorded from the tissue phantoms using a photoelastic modulation (PEM)–based polarimeter.7, 32 A schematic of the experimental setup is shown in Fig. 1 . Unpolarized light at from a He-Ne laser (LHRR-1200M, Research Electro-Optics, Boulder, Colorado) was first passed through a mechanical chopper (C) operating at a frequency and then through a linear polarizer (P1) and a removable quarter wave plate (WP1). The linear polarizer and quarter wave plate combination enabled generation of the four required input polarization states, (Stokes vector ), (Stokes vector ), (Stokes vector ) linear polarization, and circular polarization (Stokes vector ). The four set of Stokes vectors of light exiting the sample in selected detection directions were measured for each of the four input polarization states. For this study, we chose forward detection geometry; other detection geometries (including backscattering) are ongoing. The detection optics consists of another removable quarter wave plate (WP2) with its fast axis oriented at . When WP2 was in place, we could measure the Stokes parameters and (linear polarization),7 and when it was removed, we could measure the Stokes parameter (circular polarization).32 The forward scattered light from the sample then passed through a photoelastic modulator (PEM, IS-90, Hinds Instruments, Hillsboro, Oregon, operating at a frequency ) with its fast axis oriented at . The light finally passed through an analyzer (P2) oriented at . The resulting modulated intensity of the scattered light was collected using a pair of lenses (detection area of and acceptance angle ) and was relayed to an avalanche photodiode (APD) detector (C5460, Hamamatsu Photonics, Hamamatsu City, Japan). The detected signal was sent to a lock-in amplifier (SR 830, Stanford Research Systems, Sunnyvale, California) with reference input of the amplifier toggled between the chopper and the PEM controller ( and harmonics).
The Stokes vector of the output light for each of the four input polarization states were measured using the detected signal at the first harmonic of the chopper frequency and the first and the second harmonics of the PEM frequencies.7, 32 The Mueller matrix was generated using standard relationships between its 16 elements and the measured output Stokes parameters for each of the four input polarization states.28, 29., , , and correspond to the measured Stokes parameters for state of polarization of incident light of 0-, 45-, and linear polarization and circular polarization respectively. All the Mueller matrices were normalized with respect to the element .
The optical phantoms used in this study were polyacrylamide gels having strain-induced linear birefringence, sucrose-induced optical activity (magnitude of optical activity was , corresponding to concentration of sucrose) and polystyrene microspheres–induced scattering (mean diameter , refractive index ).32, 33 The refractive index of polyacrylamide was measured to be at the wavelength . The phantom preparation details are described in Ref. 32. The values for scattering coefficient of the phantoms at were calculated using the known concentration of the scatterers and Mie theory.34 The value for anisotropy parameter ( , average cosine of scattering angle) of the scattering phantom at was calculated34 to be .
To apply controllable strain to produce linear birefringence, one end of the polyacrylamide phantoms (dimension of ) was clamped to a mount and the other end to a linear translational stage. The phantoms were stretched along the vertical direction (along the long axis of the sample) to introduce varying linear birefringence with its axis along the direction of strain. Mueller matrices were recorded in the forward detection geometry from the phantoms with varying extension (varying degree of strain-induced birefringence).
Polarization-Sensitive Monte Carlo Simulations
Polarization-sensitive Monte Carlo simulation software, developed in our group,27, 32, 35 was used to generate the Mueller matrices of scattered light exiting the medium in the forward direction. In these simulations, the photons are propagated between scattering events, as determined by pseudorandom sampling of scattering mean free path, similar to conventional Monte Carlo model36 for light propagation in a multiply scattering medium. The Mie theory34 computed scattering phase function is used for sampling the scattering angles and computing the polarization changes of photons after successive scattering events. The polarization information is tracked in the form of a Stokes vector for each individual photon packet. The scattering histories of a large number of such photon packets are tracked as they propagate through the medium and are summed to yield the macroscopic polarization parameters of interest. To model the situation typically encountered in tissue, in addition to the multiple scattering effects, the simultaneous effects of linear birefringence and optical activity are also incorporated in the model. Note that this is not an obvious modeling step. Because matrix multiplication of the Mueller matrices are not commutative, different orders in which these effects are applied to the photon between scattering events will have different effects on the polarization. Therefore, this was accomplished by combining the effects into a single matrix describing them simultaneously, through the use of Jones N-matrix formalism37, 38 in the polarization-sensitive Monte Carlo simulation code. In this formalism, the matrix of the sample is represented as an exponential sum of differential matrices, known as N-matrices. Each matrix in this sum corresponds to a single optical property (e.g., linear birefringence or optical activity). These N-matrices are then used to express the M-matrix for the combined effect. The resulting M-matrix is applied to the photon between scattering events to model the simultaneous occurrence of linear birefringence and optical activity in a scattering medium. The details of this formalism are provided in Ref. 32.
Simulations were run for a given set of optical parameters of the scattering medium. The scattering phase function was computed using Mie theory for a set of input parameters, the diameter of scatterer , refractive index of scatterer , refractive index of the surrounding medium , and wavelength of light . The computed phase function, the scattering coefficient , and absorption coefficient of the medium were then used in the simulations. In the simulations, optical activity is specified as a rotation of in degrees per centimeter. Linear birefringence is modeled through the anisotropy in refractive indices , the difference in refractive index along the extraordinary axis (the axis along which the refractive index differs from other directions), and the ordinary axis. The anisotropy in refractive indices leads to a difference in the phase for light polarized parallel and perpendicular to the direction of the extraordinary axis. In this formalism, it is assumed that refractive index differs only along one axis (the extraordinary axis) and that the direction of the axis and the difference in refractive indices is constant throughout the scattering medium. In each simulation, a specific direction of this axis of linear birefringence was chosen. To generate the complete 16-element Mueller matrix, the state of polarization of incident light was changed sequentially to horizontal (Stokes vector ), vertical , linear polarization, and circular polarization ; the Stokes vectors of light exiting the sample were recorded for each given incident polarization state. The Mueller matrix was generated using Eq. 1, the standard relationships between its 16 elements and the generated output Stokes parameters for the four input polarization states.28, 29
The Mueller matrices were generated for a slab of scattering medium, having varying optical properties, for light exiting the medium in the forward direction. The average photon path length of light exiting the scattering medium was also recorded in the simulations. The photon collection geometry was chosen to have a detection area of and an acceptance angle of to match the experiment. The absorption was selected to be small, and accordingly, the value for absorption coefficient was kept at in all the simulations. To have satisfactory statistics, all the simulations were run with input photons.
Polar Decomposition of Mueller Matrices to Extract Individual Polarization Parameters
Polar decomposition of Mueller matrix is a robust mathematical tool for interpretation of the polarization characteristics of any medium. The method decomposes an arbitrary Mueller matrix M into the product of three elementary matrices representing a depolarizer , a retarder , and a diattenuator . The process and its validity were first demonstrated by Lu and Chipman30 in clear media. It is important to note that because the multiplication of the Mueller matrices is not commutative (as per the earlier N-matrix discussion), the results of the decomposition depend upon the order in which the three elementary matrices are multiplied. It has been shown previously that the six possible decompositions can be classified into two families, depending upon the order in which the depolarization and the diattenuator matrices are multiplied.39 Because the family in which the diattenuator matrix comes ahead of the depolarizer matrix always lead to a physically realizable Mueller matrix,39 this order of decomposition has been adopted in the present study.
The diattenuation matrix is defined asis a submatrix (the standard form of which is provided in reference 30) is the diattenuation vector and is defined as are elements of the Mueller matrix M. The coefficients and represent linear diattenuation for horizontal (vertical) and linear polarization, respectively, and the coefficient represents circular diattenuation.
The decomposed depolarization matrix has the formis the depolarization submatrix and the parameter depends on polarizance and diattenuation .30
The diagonal elements of the depolarization matrix can be used to calculate the depolarization coefficients [ , are depolarization coefficients for incident horizontal (or vertical) and linearly polarized light, and is the depolarization coefficient for incident circularly polarized light]. The net depolarization coefficient is defined as. The latter represents the value of degree of polarization resulting from several lumped polarization effects and also depends on the incident Stokes vector. In contrast, the depolarization coefficient defined by Eq. 7 represents the pure depolarizing transfer function of the medium.
The value for total retardance ( , a parameter that represents the combined effect of linear and circular birefringence) can be determined from the decomposed retardance matrix using the relationship30and orientation angle of the axis of linear retarder with respect to the horizontal ) and a circular retarder (optical rotation with magnitude of )31 , optical rotation , and linear retardance can be worked out as and linear retardance can be determined from the elements of the matrix as , , and determined through Eqs. 8, 11, 12 are not influenced by the order of multiplication of the two matrices in Eq. 9.
Once the value for is obtained, the optical rotation matrix formed using the value for can be inverse multiplied with the total retardance matrix to yield the linear retardance matrix . The linear retardance matrix thus obtained can be written in the formcan be obtained using the elements of the sub-matrix as30 is Levi-Civita permutation symbol.
This then allows us to determine the orientation angle of the axis of linear retardance (with respect to the horizontal axis) from the elements of the retardance vector as
The decomposition process just described was applied to the experimental Mueller matrices recorded from the polyacrylamide phantoms and to the Mueller matrices generated through Monte Carlo simulations. The values for diattenuation , depolarization coefficient , total retardance , optical rotation , and linear retardance ( and its orientation angle ) were determined from the decomposed matrices using Eqs. 5, 7, 8, 11, 12, 15, respectively. It should be mentioned that decomposition of a Mueller matrix from a multiply scattering medium exhibiting simultaneous linear birefringence and optical activity results in small but nonzero values for the off-diagonal elements of the depolarization submatrix [following Eq. 6]. In an ideal situation, should be a diagonal matrix with zero off-diagonal elements and its three diagonal elements represent depolarization coefficients for 0- and linear polarizations and circular polarization. The nonzero values for the off-diagonal elements of the decomposed depolarization submatrix arise because unlike the ordered multiplication of the three elementary matrices of a depolarizer, a retarder, and a diattenuator, in a birefringent chiral turbid medium, the individual constituent polarization effects are exhibited in an arbitrary order. To facilitate quantitative comparison of the depolarization strength of incident linearly [depolarization coefficients and ] and circularly [depolarization coefficient ] polarized light, the depolarization submatrix was presented in a diagonalized form for all the Mueller matrix decompositions. This way of representing the depolarizer matrix did not influence the polarization parameters estimated through the decomposition process.
Results and Discussion
Mueller matrices were experimentally recorded in the forward scattering geometry from the polyacrylamide phantoms having varying degrees of strain-induced birefringence ( applied along the vertical direction), varying scattering coefficients ( , 3, 4, and ), and sucrose-induced optical activity (concentration of corresponding to magnitude of optical activity ). Measurements were also performed on a clear (nonscattering) polyacrylamide phantom at different extensions to obtain the calibration curve for the strain-induced birefringence with extension. Measurements were performed three times from each of the experimental phantoms and the mean values were taken to construct the Mueller matrices. Standard deviations in the values for the normalized Mueller matrix elements were found to be in the range 0.0015 to 0.0050.
Table 1 shows the experimental Mueller matrix and the decomposed (following the procedure described in Sec. 3.2) depolarization , retardance , and diattenuation matrices for the nonscattering phantom ( , thickness ) with extension of and having a concentration of sucrose of . The values for all the retrieved polarization parameters [ ( and ), , and , determined from the equations in Secs. 3.2.1, 3.2.2, 3.2.3] are also listed in the table. As expected, the nonscattering phantom does not show any appreciable depolarization . The value for is for an extension of and the estimated orientation angle of the axis of linear retardance is reasonably close to the actual orientation of the linear birefringence axis ( for strain applied along the vertical direction). The value for increased from for varying extension of . Due to presence of optical rotation, the estimated value for the total retardance parameter [that represents the combined effect of linear retardance and optical rotation, see Eq. 10] is marginally larger than the corresponding value for . The estimate for , however, is, lower than that expected for concentration of sucrose of ( for a path length of ). The contraction of the phantom due to longitudinal stretching reduces the effective path length of transmitted light, contributing to the observed lower value of . The reduction in path length at extension is (path length of as compared with for the phantom with no strain) using the value for Poisson ratio of polyacrylamide.32, 40 However, even when accounting for this thickness reduction, the new expected is still higher than the one derived from polar decomposition. This difference is possibly due to an uncertainty in the concentration of sucrose during the fabrication of the phantom. It is also seen from Table 1 that decomposition of Mueller matrix yields a small but nonzero value for the diattenuation parameter . The presence of a small amount of dichroic absorption (at ) due to anisotropic alignment of the polymer molecules in the polyacrylamide phantom may contribute to this slight nonzero value of diattenuation.
(a) The experimentally recorded Mueller matrix and the decomposed basis matrices for a nonscattering phantom (μs=0mm−1) with strain-induced birefringence (extension of 4mm ) and a concentration of sucrose of 1M (χ=1.96degcm−1) . (b) The values for the different polarization parameters extracted from the decomposed matrices. M(1.0000−0.02290.00270.0058−0.01860.9956−0.03610.0318−0.01290.03920.2207−0.96560.00140.02800.97060.2231) MΔ(1.00000000.00410.996900−0.007000.99150−0.0019000.9966)MR(1.000000000.9988−0.03620.028100.03930.2227−0.974100.03200.97420.2239)MD(1.0000−0.02290.00270.0058−0.02291.0000−0.0000−0.00010.0027−0.00000.99970.00000.0058−0.00010.00000.9997)
Table 2 displays the experimental Mueller matrix, the decomposed basis matrices, and the retrieved polarization parameters for a nonbirefringent , achiral (concentration of , ), turbid phantom ( , , optical thickness , ). Ideally, the Mueller matrix of this phantom should model an isotropic depolarizer. The values for , , and obtained following decomposition are indeed found to be low for this nonbirefringent, achiral sample. The small nonzero value for (0.087) may arise due to the small amount of strain present during clamping of the sample to the mount and translational stage. Another possible cause, the scattering-induced linear retardance,31 which arises primarily from singly (or weakly) scattered photons, is not expected to contribute to the nonzero value for , because for this scattering sample having a value for optical thickness , multiply scattered photons are the dominant contributor to the detected photons in the forward detection geometry. Further, in agreement with previous reports,23, 24, 25 for this anisotropic scattering sample ( , prepared using large-sized scatterer with diameter ), depolarization of circularly polarized light is weaker than depolarization of linearly polarized light [ and ]. This was previously shown to arise due to the fact that in scattering media composed of large scatterers , the randomization of helicity (which is the major cause of depolarization of circularly polarized light) is considerably reduced due to predominantly forward scattering, leading to weaker depolarization of circularly polarized light.23, 24, 25 Note that the depolarization of linearly polarized light is almost independent of the orientation angle (horizontal, vertical, and ) of the incident linear polarization vector [depolarization coefficients ]. This makes sense because in an isotropic, multiply scattering medium, randomization of the field vector’s direction is the primary cause of depolarization of incident linearly polarized light. The linear depolarization coefficients of the resultant depolarizing Mueller matrix in such a situation will be the same for any orientation of the incident linear polarization vector.41
(a) The experimentally recorded Mueller matrix and the decomposed basis matrices for a nonbirefringent (extension=0mm) , achiral (concentration of sucrose=0M , χ=degcm−1 ), turbid phantom (optical thickness τ=30 , μs=3mm−1 , thickness t=10mm , g=0.95 ). (b) The values for the different polarization parameters extracted from the decomposed matrices. M(1.0000−0.01900.00300.0114−0.01830.7454−0.00250.00190.00130.00850.7378−0.06480.00340.00370.07470.8530) MΔ(1.0000000−0.00420.7468000.000000.73940−0.0065000.8566)MR(1.000000001.0000−0.0076−0.000500.00750.9961−0.087400.00120.08740.9962)MD(1.0000−0.01900.00300.0114−0.01900.9999−0.0000−0.00010.0030−0.00000.99980.00000.0114−0.00010.00000.9998)
The experimental Mueller matrix, the decomposed basis matrices, and the retrieved polarization parameters for the birefringent , chiral ( sucrose, ), and turbid phantom ( , ) are shown in Table 3 . Note the complicated nature of the Mueller matrix of this complex phantom and relatively unequivocal nature of the three basis matrices derived from the decomposition process. The elements of the depolarization matrix and the resulting value for the depolarization coefficient ( ) are similar to the corresponding values for the pure depolarizing phantom having the same scattering coefficient (Table 2, difference in the value for ). The decomposition process thus successfully decouples the depolarization effects (due to multiple scattering) and retrieves the values for the parameters (and ) and . However, the estimate for is found to be larger than that for the clear phantom having the same concentration of sucrose ( for the scattering phantom as compared with for the clear phantom). This is due to the increase in the propagation path of light in a multiply scattering medium. The average photon path length of light exiting the scattering sample in the forward direction was calculated to be from Monte Carlo simulation (optical parameters of the medium: , diameter of scatterer , , , ). For the calculation of average photon path length, the thickness of the scattering medium was taken to be as per the Poisson ratio discussion earlier. The value for calculated using the optical rotation of the nonscattering phantom ( from Table 1) and the Monte Carlo–generated average photon path length comes out to be reasonably close to that estimated through decomposition of the Mueller matrix recorded from the birefringent, chiral, turbid phantom (calculated value of as compared with estimated for the turbid phantom).
(a) The experimentally recorded Mueller matrix and the decomposed basis matrices for a birefringent (extension=4mm) , chiral (concentration of sucrose=1M , χ=1.96degcm−1 ), turbid ( μs=3mm−1 , g=0.95 ) phantom. (b) The values for the different polarization parameters extracted from the decomposed matrices. M(1.0000−0.03120.0029−0.0066−0.02140.7678−0.03700.0204−0.00550.02300.1043−0.77350.00140.03900.79720.1920) MΔ(1.00000000.00280.755200−0.010200.768900.0016000.8454)MR(1.000000000.9984−0.05130.024200.03330.1844−0.982300.04600.98150.1858)MD(1.0000−0.03120.0029−0.0066−0.03121.0000−0.00000.00010.0029−0.00000.9995−0.0000−0.00660.0001−0.00000.9995)
Although the estimated value for of the turbid phantom is larger than that for the clear phantom ( for the turbid phantom as compared with for the clear phantom having a path length of ), the value is lower than that one would expect for average photon path length of . This is illustrated in Fig. 2 where the estimates for and of the chiral phantoms having varying degrees of strain-induced birefringence (extension of ) are displayed. Results are shown for both the clear and the turbid phantoms. The increase in the value for as a result of increased average photon path length in the turbid phantom as compared with the clear phantom can be seen quite clearly. The gradual decrease in the value for with increasing longitudinal stretching is consistent with resulting lateral contraction of the phantom, reducing the effective path length. In contrast to , the expected increase in as a result of increased average photon path length in the turbid phantom is not seen. The results of studies on phantoms having other values of scattering coefficients ( and , data not shown) also yielded similar trends to that observed for the phantom with .
To understand these trends further and to investigate the effect of multiple scattering, propagation path, and detection geometry on and parameters, Monte Carlo simulations were carried out.
The Monte Carlo simulations were run with input optical parameters chosen to mimic the optical properties of the phantoms used in the experimental studies. The photon collection geometry was chosen to have a detection area of and an acceptance angle of around the forward directed ballistic beam. Standard deviations (originating from random statistical errors) in the different elements of the Monte Carlo–generated Mueller matrices were found to be in the range 0.0004 to 0.0025. Table 4 shows the Mueller matrix, the decomposed basis matrices, and the retrieved polarization parameters for a nonbirefringent, achiral, turbid medium (linear birefringence , , , , ) composed of monodisperse spherical scatterers ( , , , ). In agreement with the experimental results (Table 2), polarization preservation of circularly polarized light is found to be larger than linearly polarized light [ and ] for this forward scattering turbid medium. The elements of the decomposed depolarization matrix are also close to the experimental Mueller matrix results (difference in the value for ) for the phantom having similar scattering coefficient (small differences may arise due to slight mismatch in the collection geometry employed in the experiments and in the simulation, and possibly due to subtle errors in concentration of scatterers in the phantom). Further, as expected, the values for , , and of this nonbirefringent, achiral, turbid medium are also small. The small value of for this nonbirefringent turbid medium confirms that scattering-induced linear retardance should not contribute to the experimentally observed slightly larger value for for the phantom having the same scattering coefficient of .
(a) The Monte Carlo simulation–generated Mueller matrix and the decomposed basis matrices for a nonbirefringent, achiral, turbid medium ( Δn=0 , χ=0degcm−1 , τ=30 , μs=3mm−1 , t=10mm , g=0.95 ). (b) The values for the different polarization parameters extracted from the decomposed matrices. M(1.0000−0.0025−0.0026−0.0024−0.00080.77870.00010.0009−0.00200.00030.77670.00340.00020.0004−0.00040.8809) MΔ(1.00000000.00120.7788000.000000.776600.0024000.8810)MR(1.000000001.0000−0.00010.000300.00011.00000.00230−0.0003−0.00231.0000)MD(1.0000−0.0025−0.0026−0.0024−0.00251.00000.00000.0000−0.00260.00001.00000.0000−0.00240.00000.00001.0000)
Figure 3 shows the variation of the parameters , , , and of transmitted light as a function of optical thickness (varying and fixed path length of ) for the nonbirefringent, achiral, turbid medium. As expected, depolarization decreases with increasing value of . In the absence of chiral molecules, the values for are low at all the optical thickness values of the achiral, turbid medium. Further, both the parameters and are for the range of to 80. As discussed previously, the scattering-induced linear retardance and diattenuation arises primarily from singly (or weakly) scattered photons and the magnitude is greater at large scattering angles. Because for detection around the forward direction, the contribution of such photons is negligible, the scattering-induced retardance and diattenuation are weak in this geometry.
Table 5 presents the Monte Carlo–generated Mueller matrix, the decomposed basis matrices, and the retrieved polarization parameters for a birefringent, chiral, turbid medium ( , corresponding to for a path length of , optical activity , , , ). The axis of linear birefringence was along the vertical direction in the simulation. The values for and estimated following the decomposition process are reasonably close to those of the experimental phantom having similar properties (Table 3). Slightly larger disagreement in possibly arises due to an uncertainty in the concentration of sucrose during the fabrication of the phantom (as noted previously). Incorporation of the path length reduction due to longitudinal stretching ( as compared with ) in the simulations did not result in significant changes in the estimated values for and (reduction in and were ). As for the nonbirefringent, achiral, turbid medium (Table 4), this birefringent, chiral, turbid medium prepared using large-sized scatterers (diameter ) also exhibits weaker depolarization of circularly polarized light than linearly polarized light . The elements of and the resulting value for are also similar to the corresponding values for the nonbirefringent, achiral, turbid medium with the same scattering properties.
(a) The Monte Carlo simulation–generated Mueller matrix and the decomposed basis matrices for a birefringent, chiral, turbid medium ( Δn=1.36×10−5 , χ=1.96degcm−1 , τ=30 , μs=3mm−1 , t=10mm , g=0.95 ). (b) The values for the different polarization parameters extracted from the decomposed matrices. M(1.0000−0.00390.0022−0.00970.00120.77420.0310−0.0393−0.0061−0.04290.1234−0.7949−0.0009−0.02510.81460.1983) MΔ(1.00000000.00380.775800−0.014200.78360−0.0009000.8595)MR(1.000000000.99790.0412−0.04900−0.05610.1948−0.97920−0.03080.98000.1967)MD(1.0000−0.00390.0022−0.0097−0.00391.0000−0.00000.00000.0022−0.00000.9999−0.0000−0.00970.0000−0.00001.0000)
Note that at an even larger value of the input linear birefringence , the decomposed circular depolarization coefficient was found to be lower than the corresponding value for the nonbirefringent, turbid medium (the trend can be seen from Tables 4, 5 also). Because for this forward scattering medium , depolarization of circularly polarized light is weaker than depolarization of linearly polarized light,23, 24, 25 while propagating through such a birefringent (with axis of birefringence along the vertical direction), turbid medium, a greater fraction of incident circularly polarized light gets transferred to linear polarization as compared with transfer of linear polarization to circular polarization. This also led to a difference in the value for the depolarization coefficient for linearly polarized light and horizontally (or vertically) polarized light . This suggests that even in the presence of strong multiple scattering, for a turbid medium that exhibits linear birefringence, the depolarization of linearly polarized light may depend on the orientation angle of the incident linear polarization vector.
To gain quantitative understanding of the dependence of the and parameters on the propagation path of multiply scattered photons, Fig. 4 shows the Monte Carlo results for the variation of and as a function of the thickness of the same birefringent, chiral, scattering medium. The top axis shows the calculated average photon path length of light exiting the scattering medium. The value for increases with increasing average photon path length. However, the absolute values for were marginally lower compared to calculations using the linear relationship ( photon path length). It should be noted here that the average path length has contributions from both the polarization-preserving and the depolarized photons. The fact that the propagation path of the polarization-preserving photons (which would show experimentally detectable optical rotation) is shorter than the average photon path length42 should account for the observed lower value for .
It is interesting that although the value for increases with increasing average photon path length, the increase is much slower as compared with . A similar trend was also observed in the experimental studies. This is because in a turbid medium, the forward scattered light does not travel in a straight line but rather along a curved path, the curvature being controlled by the values for and . Though such paths do not have any influence on (as long as the average photon path length is the same), the net value for would be significantly affected because the linear birefringence has a specific orientation axis ( in the simulations and in the experiments).
To understand this, Mueller matrices were generated for transmitted light collected at the exit face of the birefringent ( and ), chiral , turbid medium ( , , ) at different spatial positions away from the position of the ballistic beam along the horizontal ( axis, perpendicular to the direction of linear birefringence) and vertical ( axis, parallel to the direction of linear birefringence) directions. The variations for the parameter as a function of the distance from the ballistic beam position (along both and axis) at the exit face are displayed in Fig. 5 . At spatial positions along the axis, increases with increasing distance from the center of the ballistic beam (i.e., with increasing average photon path length). In contrast, for detection at spatial positions along the axis, the value for shows gradual decrease with increasing distance from the center of the ballistic beam. This is because a larger component of the photon propagation path is along the axis of birefringence leading to a reduction in net linear retardance (because propagation along the direction of the birefringence axis does not yield any retardance). Because such differences in the photon propagation path for the two different detection geometries should have no influence on , the estimates for were found to be identical for similar detection positions either along the or the axis at the exit face of the medium (data not shown). It thus appears that the value for is influenced by the detection position and the orientation angle of the axis of linear birefringence in a birefringent, turbid medium.
To examine this further, in Fig. 6 , we show the variation of as a function of the orientation angle of the axis of linear birefringence for the birefringent, chiral, turbid medium ( , , , , ). Results are shown for three different detection positions at the exit face of the medium, around the position of the ballistic beam [coordinate (0,0)], at away from the position of the ballistic beam along the horizontal [coordinate (3,0)] and vertical [coordinate (0,3)] axes, respectively. The value for is observed to vary considerably with a change in the birefringence orientation angle for off-axis detection [at spatial positions (3,0) or (0,3)]. As noted previously, the minimum values for are found at spatial positions along the direction of the birefringence axis. Due to the symmetry of position coordinates with respect to the direction of birefringence axis, the two values of [for spatial positions (3,0) and (0,3)] are also observed to converge for the value of . In contrast, for detection around the position of the ballistic beam [coordinate (0, 0)], the estimates for do not show any significant variation with a change in the orientation angle of the axis of birefringence. Thus for simultaneous determination of the intrinsic values for the parameters and of a birefringent, chiral, turbid medium in the forward scattering geometry, detection around the direction of propagation of the ballistic beam may be preferable.
The presented results demonstrate that decomposition of the Mueller matrix can be used for simultaneous determination of the intrinsic parameters (and ) and of a birefringent, chiral, turbid medium. For conceptual and practical reasons, the extension of the developed methodology to backward detection geometry is warranted. However, this will be more complex, because in the backward detection geometry, the contribution of backscattered (large angle scattered) photons would be greater and thus the interplay of the scattering-induced linear retardance and diattenuation with the intrinsic values for and would be more coupled. Future studies will investigate this aspect in detail and will concentrate on developing an appropriate strategy to minimize the scattering-induced artifacts for simultaneous determination of and in the backscattering geometry using Mueller matrix decomposition.
We have investigated the validity of Mueller matrix decomposition methodology to extract the individual intrinsic polarimetry characteristics from a multiply scattering medium exhibiting simultaneous linear birefringence and optical activity. Decomposition of experimental Mueller matrices recorded from tissue phantoms exhibiting controllable scattering, linear birefringence, and optical activity yielded satisfactory estimates for linear retardance and optical rotation . Reliable estimates for linear and circular depolarization coefficients were also obtained. The effects of multiple scattering, the propagation path of scattered photons, and detection geometry on and were further investigated using polarization-sensitive Monte Carlo simulations. The simulation results were found to corroborate the experimental findings and demonstrate that the Mueller matrix decomposition methodology can successfully extract the individual polarization characteristics of a medium that exhibits simultaneous linear birefringence, optical activity, and multiple scattering effects.
The ability of this technique to determine the constituent contributions of individual processes should prove valuable in diagnostic photomedicine. The value for linear retardance estimated using this approach would represent its average value over the volume of the turbid medium (tissue) probed by the polarization-preserving photons43 and could monitor changes in tissue birefringence as a signature of tissue abnormality and/or response to therapy. Note that in actual tissues due to larger value of and due to presence of absorption, the sampling depth of the polarization-preserving photons will be considerably reduced (typically a few millimeters) as compared with that of the phantoms used in this study. As a representative application, we are currently investigating the use of the decomposition-derived linear retardance for monitoring myocardial tissue regeneration following stem cell treatments of myocardial infarction. Measurement (in the forward detection geometry) through thick ex vivo myocardial tissue samples are yielding encouraging results and show promise for the use of this approach for monitoring of stem-cell–based treatments of myocardial infarction. These results will be reported in a separate publication. Further, the decomposition method may also find useful applications in quantifying the optical rotations due to blood glucose in diabetic patients, potentially impacting the (currently unsolved) noninvasive glucose-monitoring problem. The results presented here demonstrate the ability to quantify optical activity (albeit at higher than physiological levels) despite depolarization and other effects in multiply scattering media, using the decomposition methodology. Ordinarily, the small chirality-caused optical rotation is swamped by much larger scattering-induced rotation, and it is the decomposition method that enables the glucose-related signal to be derived and quantified from the Mueller matrix results. The application of this approach to real tissues with lower blood glucose levels ( ), however, will, require further refinements of the highly sensitive Mueller matrix measurements to detect small changes in the matrix elements corresponding to the physiological glucose levels. This and the combined use of the Mueller matrix decomposition approach and spectroscopic polarimetry (coupled with chemometric data analysis to isolate the rotation due to glucose from that caused by other chiral constituent confounders in tissue) is currently under investigation in our laboratory.
The authors would like to thank Dr. Daniel Coté for his help in developing the polarization-sensitive Monte Carlo code. Natural Sciences and Engineering Research Council of Canada and the Ontario Graduate Scholarship program are gratefully acknowledged for their financial support.