Open Access
Add to BibSonomyAbstract
A time-resolved Monte Carlo technique was used to simulate the propagation of polarized light in turbid media. Calculated quantities include the reflection Mueller matrices, the transmission Mueller matrices, and the degree of polarization (DOP). The effects of the polarization state of the incident light and of the size of scatterers on the propagation of DOP were studied. Results are shown in animation sequences.
©2000 Optical Society of America
1. Introduction
Tissue optics has become an active area of research primarily because light is non-ionizing and it can furnish physiological information [1]. The primary challenge, however, is due to the strong scattering of light in biological tissues: multiply scattered photons degrade the imaging contrast and resolution. Several techniques have been studied to differentiate between weakly scattered and multiply scattered photons. Because multiply scattered photons usually have greater path lengths, they can be rejected with time gating [2, 3]. It is also widely recognized that the original polarization state is lost in multiply scattered light, but is partially preserved in weakly scattered light. Polarization techniques can thus be employed to discriminate weakly scattered light from multiply scattered light [4–6]. Of course, it is also possible to combine these two techniques [7].
The propagation of polarized light in turbid media is a complex process. Parameters, such as the size, shape, and density of the scatterers as well as the polarization state of the incident light, all play important roles [8, 9]. A good understanding of this process is essential for improving the polarization-based techniques. Because the number of scattering events is related to the optical path length and the time of propagation, a time-resolved study is needed to understand the evolution of polarization in turbid media.
We used a time-resolved Monte Carlo technique to simulate the propagation of polarized light in turbid media. In particular, we calculated the reflection Mueller matrices, the transmission Mueller matrices, and the evolution of the degree of polarization (DOP) in turbid media. We also studied the effects of the size of the scatterers and the polarization state of the source. Results are presented as animation sequences.
2. Methods
Several groups have used Monte Carlo techniques to simulate the steady-state backscattering Mueller matrix of a turbid medium [10, 11]. Whereas an indirect method utilizing the symmetry of the backscattering Mueller matrix was used in Ref. 10, the direct tracing method [11] was used in our simulation. The turbid medium was assumed to have a slab structure, on which a laboratory coordinate system was defined (Fig. 1). A pencil beam was incident upon the origin of the coordinate system at time zero along the Z axis.
The basic Stokes-Mueller formalism and the simulation of propagation of polarized light in turbid media have been described earlier [10, 11]. Briefly, the Stokes vector and the local coordinates of each incident photon packet were traced statistically. At each scattering event, the incoming Stokes vector of the photon packet was first transformed into the scattering plane through a rotation operator and then converted by
where S is the Stokes vector before scattering, but it is re-defined in the scattering plane; S’ is the Stokes vector of the scattered photon; θ is the polar scattering angle; and M is the single-scattering Mueller matrix, given by the Mie theory as [13]
The element m 11 satisfies the following normalization requirement:
The joint probability density function (pdf) of the polar angle θ and the azimuth angle ϕ is a function of the incident Stokes vector S={S0, S1, S2, S3}:
In our method, the polar angle θ is sampled according to m 11(θ) and the azimuth angle ϕ is sampled with the following function:
It is worth noting that a biased sampling technique was used in Ref. 10. The Stokes vectors of all the output-photon packets were transformed to the laboratory coordinate system and then accumulated to obtain the final Stokes vector. The Mueller matrix of the scattering media can be calculated algebraically from the Stokes vectors of four different incident polarization states [12]. The degree of polarization (DOP) was calculated by
To accelerate the computation we calculated the single-scattering Mueller matrix and the probability density functions of the scattering angles and stored them in arrays before tracing the photon packets. The path length of the photon packets was recorded to provide time-resolved information. For purposes of illustration, the scatterers were assumed to be spherical; the thickness of the scattering slab was taken to be 2 cm; the temporal resolution was 1.33 ps, corresponding to 0.4 mm in real space; the wavelength of light was 543 nm; the absorption coefficient was 0.01 cm-1; the index of refraction of the turbid medium was unity, matching that of the ambient. The dimensions of the pseudo-color images in the following section are 4, 4, and 2 cm along the X, Y, and Z axes, respectively.
3. Results
Figure 2 shows the reflection and the transmission Mueller matrices of a turbid medium with a scattering coefficient of 4 cm-1 and a radius of scatterers of 0.102 µm. The calculated Mueller-matrix elements were normalized to the m 11 element to compensate for the radial decay of intensity. Each of the images is displayed with its own color map to enhance the image contrast. The size of each image is 4×4 cm2.
The patterns of the reflection Mueller matrix are identical to those reported previously [10, 11]. The symmetries in the patterns can be explained by the symmetries in the single-scattering Mueller matrix and the medium [10]. The transmission Mueller matrix has different patterns from the reflection Mueller matrix. One of the noticeable differences lies in elements m 31 and m 13, which are anti-symmetric in the reflection Mueller matrix but symmetric in the transmission Mueller matrix. This difference is caused by the mirror effect in the reflection process of the scattered light.
Figure 3 shows the time-resolved DOP propagation in the turbid medium with right-circularly (R) and horizontal-linearly (H) polarized incident light. The scattering coefficient was 1.5 cm-1, and the radius of the scatterers was 0.051 µm. The anisotropic factor <cos(θ)> was 0.11. The transport mean free path was calculated to be 0.74 cm. In the simulation the Stokes vectors of the forward propagating photons were accumulated to calculate the DOP. As shown in the movies, the DOPs at the expanding edges of the distributed light remain near unity because these photons experience few scattering events. As the light propagates in the medium, the DOP in some regions decreases significantly. The DOP patterns are dependent on the single-scattering Mueller matrix and the density of scattering particles.
Fig.3. (787 KB) Movie of the DOP propagation in the slab. The X axis is along the horizontal direction, and the Z axis is along the vertical direction. R: right-circularly polarized incident light. H: horizontal-linearly polarized incident light.
From Fig. 3 the DOP images have different patterns for the R- and the H-polarized incident light. As expected, such a difference appears in the DOP images of transmitted light as well. As shown in Fig. 4, the DOP images of the transmitted light are rotationally symmetric for circularly polarized incident light, whereas such symmetry does not exist for linearly polarized incident light. This difference is related to the dependence of the scattering probability on the incident Stokes vector (Eqs. 4 and 5). The Stokes vectors of the R- and the H-polarized light are [1, 0, 0, 1]’ and [1, 1, 0, 0]’, respectively. According to Eq. 5, the R-polarized light has a uniform pdf for the azimuth angle, whereas the H-polarized light has a non-uniform one. Hence, the single-scattering pattern depends on the polarization state of the incident light. Because the change of DOP is related strongly to the single scattering events, it is understandable that the DOP images have different features for different incident Stokes vectors.
Fig. 4. (950 KB) Movie of the DOP of the transmitted light with (a) R- and (b) H-polarized incident light. The X axis is along the horizontal direction, and the Y axis is along the vertical direction.
To demonstrate the dependence of the evolution of DOP on the number of scattering events, we recorded the time-resolved images of the average number of scattering events of the transmitted light (Fig. 5). The simulation parameters are the same as those for Figs. 3 and 4. As it can be seen, the transmitted photons experience more and more scattering events as time elapses. If we compare Figs. 5 and 4, it is clear that the patterns of DOP are related directly to the patterns of the scattering counts: the DOP decreases as the number of scattering events increases. Nevertheless, the number of scattering events does not solely determine the change of DOP. Two dark regions are clearly visible in Fig. 4(b), but they are inseparable in Fig. 5(b). The change of DOP must also depend on the nature of each scattering event, determined by the single-scattering Mueller matrix.
Fig.5. (950 KB) Movie of the weighted-averaged numbers of scattering events for (a) R- and (b) H-polarized incident light. The numbers of scattering events are normalized to a maximum value of 7 for the plots. The X axis is along the horizontal direction, and the Y axis is along the vertical direction.
To study the effect of the size of the scatterers, we simulated the evolution of the DOP in a scattering medium with a different radius of 1.02 µm. The scattering coefficient of the medium was 14 cm-1, and the anisotropic factor was 0.91. The transport mean free path was 0.76 cm, which was similar to the value for Figs. 3–5. The time-resolved propagation of the DOP in the medium is shown in Fig. 6. The DOP movie of the transmitted light is shown in Fig. 7.
Fig.6. (787 KB) Movie of the DOP propagation in the slab. The X axis is along the horizontal direction, and the Z axis is along the vertical direction. R: right-circularly polarized incident light. H: horizontal-linearly polarized incident light.
Fig. 7. (950 KB) Movie of the DOP of the transmitted light with (a) R- and (b) H-polarized incident light. The X axis is along the horizontal direction, and the Y axis is along the vertical direction.
Note the different patterns in Figs. 6–7 and Figs. 3–4. The key difference is that the DOP of R-polarized incident light was preserved much better than the DOP of H-polarized incident light for the large scatterers, as shown in Figs. 6–7. In contrast, the DOP of H-polarized incident light was preserved better than the DOP of R-polarized incident light for the small scatterers, as shown in Figs. 3–4. This observation is consistent with previous experimental and theoretical findings [8, 9].
Another significant difference is that the DOP patterns for the large scatterers become rotationally symmetric even when the incident light is H-polarized. This phenomenon can be easily understood if we examine the probability distribution functions of the scattering angle for different particle sizes. The scattering angle θ is determined by m 11, and its probability density function ρ(θ) is 2πm 11sin(θ). The probability density function of the azimuth angle ϕ is a function of both ϕ and the incident Stokes vector, as defined in Eq. 5. The contribution of the ϕ-dependent term is proportional to |m 12/m 11|. The curves of ρ(θ) and |m 12/m 11| are shown in Fig. 8. When the scatterer size is small, ρ(θ) is approximately homogeneous and the photon is likely to be scattered into 60°–120° [Fig. 8(a)]. At these angles, the |m 12/m 11| ratio has large values and Eq. 5 depends strongly on ϕ. When the scatterer size is large, most of the photons are scattered into smaller angles [Fig. 8(b)]. The |m 12/m 11| ratio is small at small scattering angles, which means that the homogeneous-distribution term is dominant in the probability distribution function of the ϕ angle. As a consequence, the scattering process becomes rotationally symmetric for the larger particle sizes.
Fig. 8. Probability density function ρ(θ) and |m 12/m 11| at a particle radius of (a) 0.051 µm and (b) 1.02 µm.
As observed in Figs. 4 and 7, the DOP at the center of the 2D time-resolved images is smaller than that in the outer area. However, the DOP decreases radially in the time-integrated images (Fig. 9), which is in agreement with previous experimental results [9]. This is because the photons exiting at early times have a limited span in space and hence have a dominant weight in the central area. These early exiting photons better preserve the DOP because they experience fewer scattering events than the photons exiting at later times. Consequently, a combination of time-gating and polarization discrimination [7] has a better performance in rejecting multiply scattered photons than either technique alone. Figure 9 also reveals that the DOP decreases as the scattering coefficient increases. The radial distribution of the DOP becomes flat at a large scattering coefficient because of the increasing number of multiply scattered photons [9].
Fig. 9. Radial distribution of the DOP of the transmitted light for different scattering coefficients. The particle radius was 1.02 µm. The incident light was H polarized.
4. Conclusion
A Monte Carlo technique was employed to simulate the time-resolved propagation of polarized light in scattering media. Results are consistent with prior experimental findings. Hence, time-resolved simulation is a useful tool for understanding better the essential physical processes of polarization propagation in turbid media. Because of the nature of the Monte Carlo simulation, coherent phenomena, such as laser speckles, are not modeled. Nevertheless, the simulation method can be applied in the non-coherent regime or in the cases where the coherent effect is removed, such as ensemble-averaged measurements [10].
Acknowledgments
This project was sponsored in part by National Institutes of Health grants R29 CA68562, R01 CA71980, and R21 CA83760, National Science Foundation grant BES-9734491, and Texas Higher Education Coordinating Board grant 000512-0123-1999.
References and links
1. R. R. Alfano and J. G. Fujimoto, eds., Advances in Optical Imaging and Photon Migration, Vol. 2 of Topics in Optics and Photonics Series (Optical Society of America, Washington, D. C., 1996).
2. B. Das, K. Yoo, and R. R. Alfano, “Ultrafast time gated imaging,” Opt. Lett. 18, 1092–1094(1993). [CrossRef] [PubMed]
3. S. Marengo, C. Pepin, T. Goulet, and D. Houde, “Time-gated transillumination of objects in highly scattering media using a subpicosecond optical amplifier,” IEEE J. Sel. Top. Quant. 5, 895–901(1999). [CrossRef]
4. J. M. Schmitt, A. H. Gandjbakhche, and R. F. Bonner, “Use of polarized light to discriminate short-path photons in a multiply scattering medium,” Appl. Opt. 31, 6535–6546(1992). [CrossRef] [PubMed]
5. S. P. Morgan, M. P. Khong, and M. G. Somekh, “Effects of polarization state and scatterer concentration on optical imaging through scattering media,” Appl. Opt. 36, 1560–1565(1997). [CrossRef] [PubMed]
6. S. L. Jacques, J. R. Roman, and K. Lee, “Imaging superficial tissues with polarized light,” Lasers in Surg. & Med. 26, 119–129(2000). [CrossRef]
7. X. Liang, L. Wang, P. P. Ho, and R. R. Alfano, “Time-resolved polarization shadowgrams in turbid media,” Appl. Opt. 36, 2984–2989(1997). [CrossRef] [PubMed]
8. D. Bicout, C. Brosseau, A. S. Martinez, and J. M. Schmitt, “Depolarization of multiply scattered waves by spherical diffusers: influence of the size parameter,” Phy. Rev. E 49, 1767–1770(1994). [CrossRef]
9. V. Sankaran, K. Schonenberger, J. T. Walsh Jr., and D. J. Maitland, “Polarization discrimination of coherently propagation light in turbid media,” Appl. Opt. 38, 4252–4261(1999). [CrossRef]
10. M. J. Rakovic, G. W. Kattawar, M. Mehrubeoglu, B. D. Cameron, L. V. Wang, S. Rastegar, and G. L. Coté, “Light backscattering polarization patterns from turbid media: theory and experiments,” Appl. Opt. 38, 3399–3408(1999). [CrossRef]
11. S. Bartel and A. H. Hielscher, “Monte Carlo simulation of the diffuse backscattering Mueller matrix for highly scattering media,” Appl. Opt. 39, 1580–1588(2000). [CrossRef]
12. G. Yao and L. V. Wang, “Two dimensional depth resolved Mueller matrix measurement in biological tissue with optical coherence tomography,” Opt. Lett. 24, 537–539(1999). [CrossRef]
13. H. C. van de Hulst, Light Scattering by Small Particles (Dover, New York, 1981).
