Monte Carlo simulation of polarized photon scattering in anisotropic media

Tianliang Yun; Nan Zeng; Wei Li; Dongzhi Li; Xiaoyu Jiang; Hui Ma

doi:10.1364/OE.17.016590

1. Introduction

Monte Carlo simulation has been an effective method for the investigation of light propagation in turbid media, such as biological tissues. The source code of a Monte Carlo program for simulating the behavior of photon transport in multilayered turbid media was made public [1]. In the program, the photons are unpolarized and the scattering phase function follows the Henyey-Greenstein function for isotropic media.

After the very early report on Monte Carlo simulation of polarized light scattering [2], many groups have studied the polarization characteristics of light propagation in turbid media, and reported different ways of tracking the polarization status of scattering photons in isotropic media [3–5]. Several authors calculated backscattering Mueller matrix of a solution of microspheres [6–8]. In addition, a three dimensional Monte Carlo has been developed to study optically active molecules in turbid media [9]. Most of these previous implementations describe polarization in terms of Stokes vector.

In the past few years there are increasingly interests in the light-scattering properties of tissues for therapeutic and diagnostic purposes. These applications motivate efforts in improving numerical simulations for real biological samples, which often contain complex tissue microstructure. Some groups have investigated light scattering from nonspherical particles, such as disc [10] and ellipsoid [11]. Many tissues contain microstructure of ordered elongated subunits, such as the myofibrils in muscles or the collagen fibers in skin, tendons, or ligaments. They are anisotropic and cannot be described by an isotropic scattering model, such as randomly distributed spherical or non spherical scatterers. Recently, Monte Carlo simulations of unpolarized photon propagation in anisotropic media, such as dentin [12] and softwood [13], have been reported. The anisotropic media were approximated to a mixture of aligned cylinders of infinite length and solid microspheres [14].

In this paper, we use a similar sphere-cylinder model for anisotropic scattering media. This model consists of both spherical and infinite long cylindrical scattering particles. A polarization-dependent Monte Carlo program is developed to simulate the behaviors of polarized photons as they propagate through the medium. Densities and sizes of the spherical and cylindrical scatterers, as well as the orientation of the cylinders are variables in the simulation of different anisotropic media. The good agreement between the simulation and the experimental results from polarization imaging confirms the validity of this anisotropic scattering model and of the polarization-dependent Monte Carlo program.

2. Methods: Monte Carlo simulation

2.1 Polarized photon scattering at an infinitely long cylinder

Using the continuity conditions on the boundary of a cylinder, the scalar wave equation $\nabla^{2} ψ + k^{2} ψ = 0$ has been solved analytically in cylindrical polar coordinate for light scattering by an infinitely long cylinder. The scattered fields in matrix form and relevant parameters are given as [15]:

(\begin{array}{l} E_{∥ s} \\ E_{⊥ s} \end{array}) = e^{i 3 π / 4} \sqrt{\frac{2}{π k r \sin ζ}} e^{i k (r \sin ζ - z \cos ζ)} (\begin{array}{c} \begin{array}{l} T_{1} \\ T_{3} \end{array} & \begin{array}{l} T_{4} \\ T_{2} \end{array} \end{array}) (\begin{array}{l} E_{∥ i} \\ E_{⊥ i} \end{array})

{\begin{cases} T_{1} = \sum_{- \infty}^{\infty} b_{n I} e^{- i n Θ} = b_{0 I} + 2 \sum_{n = 1}^{\infty} b_{n I} \cos (n Θ) \\ T_{2} = \sum_{- \infty}^{\infty} a_{n I I} e^{- i n Θ} = a_{0 I I} + 2 \sum_{n = 1}^{\infty} a_{n I I} \cos (n Θ) \\ T_{3} = \sum_{- \infty}^{\infty} a_{n I} e^{- i n Θ} = - 2 i \sum_{n = 1}^{\infty} a_{n I} \sin (n Θ) \\ T_{4} = \sum_{- \infty}^{\infty} b_{n I I} e^{- i n Θ} = - 2 i \sum_{n = 1}^{\infty} b_{n I I} \sin (n Θ) \end{cases}

{\begin{cases} \begin{array}{c} a_{n I} = \frac{C_{n} V_{n} - B_{n} D_{n}}{W_{n} V_{n} + i D_{n}^{2}}, & a_{n I I} = - \frac{A_{n} V_{n} - i C_{n} D_{n}}{W_{n} V_{n} + i D_{n}^{2}} \end{array} \\ \begin{array}{c} b_{n I} = \frac{W_{n} B_{n} + i D_{n} C_{n}}{W_{n} V_{n} + i D_{n}^{2}}, & b_{n I I} = - i \frac{C_{n} W_{n} + A_{n} D_{n}}{W_{n} V_{n} + i D_{n}^{2}} \end{array} \end{cases}

{\begin{cases} A_{n} = i ξ [ξ J_{n}^{'} (η) J_{n} (ξ) - η J_{n} (η) J_{n}^{'} (ξ)] \\ B_{n} = ξ [m^{2} ξ J_{n}^{'} (η) J_{n} (ξ) - η J_{n} (η) J_{n}^{'} (ξ)] \\ C_{n} = n \cos ζ η J_{n} (η) J_{n} (ξ) (\frac{ξ^{2}}{η^{2}} - 1) \\ D_{n} = n \cos ζ η J_{n} (η) H_{n}^{(1)} (ξ) (\frac{ξ^{2}}{η^{2}} - 1) \\ W_{n} = i ξ [η J_{n} (η) H_{n} {^{(1)}}^{'} (ξ) - ξ J_{n}^{'} (η) H_{n}^{(1)} (ξ)] \\ V_{n} = ξ [m^{2} ξ J_{n}^{'} (η) H_{n}^{(1)} (ξ) - η J_{n} (η) H_{n} {^{(1)}}^{'} (ξ)] \end{cases}

{\begin{cases} Θ = π - θ \\ ξ = x \sin ζ \\ η = x \sqrt{m^{2} - \cos^{2} ζ} \\ x = k a \end{cases}

where

E_{∥ i}

and

E_{⊥ i}

represent the incident electrical field, r and z are cylindrical polar coordinates, a is the radius of the cylinder, m is the refractive index of the cylinder relative to that of the surrounding medium, k is the wave number, ζ and θ are the incident angle and the angle between the incident and scattering planes as shown in Fig. 1 .

J_{n}

is the Bessel function and

H_{n}^{(1)} = J_{n} + i Y_{n}

is the Hankel function. The Mueller matrix

M (ζ, θ)

can be derived from Eq. (1) according to the definition of Stokes vector:

M (ζ, θ) = \frac{2}{π k r \sin ζ} (\begin{array}{c} m_{11} & m_{12} & m_{13} & m_{14} \\ m_{21} & m_{22} & m_{23} & m_{24} \\ m_{31} & m_{32} & m_{33} & m_{34} \\ m_{41} & m_{42} & m_{43} & m_{44} \end{array})

where

{\begin{cases} m_{11} = (| T_{1} |^{2} + | T_{2} |^{2} + | T_{3} |^{2} + | T_{4} |^{2}) / 2 \\ m_{12} = (| T_{1} |^{2} - | T_{2} |^{2} + | T_{3} |^{2} - | T_{4} |^{2}) / 2 \\ m_{13} = Re {T_{1} T_{4}^{*} + T_{2} T_{3}^{*}} \\ m_{14} = Im {T_{1} T_{4}^{*} - T_{2} T_{3}^{*}} \\ m_{21} = (| T_{1} |^{2} - | T_{2} |^{2} - | T_{3} |^{2} + | T_{4} |^{2}) / 2 \\ m_{22} = (| T_{1} |^{2} + | T_{2} |^{2} - | T_{3} |^{2} - | T_{4} |^{2}) / 2 \\ m_{23} = Re {T_{1} T_{4}^{*} - T_{2} T_{3}^{*}} \\ m_{24} = Im {T_{1} T_{4}^{*} + T_{2} T_{3}^{*}} \\ m_{31} = Re {T_{1} T_{3}^{*} + T_{2} T_{4}^{*}} \\ m_{32} = Re {T_{1} T_{3}^{*} - T_{2} T_{4}^{*}} \\ m_{33} = Re {T_{1}^{*} T_{2} + T_{3}^{*} T_{4}} \\ m_{34} = Im {T_{1} T_{2}^{*} + T_{3} T_{4}^{*}} \\ m_{41} = Im {T_{1}^{*} T_{3} + T_{2} T_{4}^{*}} \\ m_{42} = Im {T_{1}^{*} T_{3} - T_{2} T_{4}^{*}} \\ m_{43} = Im {T_{1}^{*} T_{2} - T_{3}^{*} T_{4}} \\ m_{44} = Re {T_{1}^{*} T_{2} - T_{3}^{*} T_{4}} \end{cases}

The matrix elements of

M (ζ, θ)

can be calculated using Eq. (2-7). In our simulation, the angle resolutions of ζ and θ of

M (ζ, θ)

are 1°.

Fig. 1 Scheme of light scattering characteristics by an infinitely long cylinder

Download Full Size | PDF

2.1.1 The phase function

The spatial distribution of scattered photons for oblique incidence at an infinitely long cylinder follows a conical shape around the cylinder with the half angle ζ which is the same as the incident angle [see Fig. 1]. The scattering probability around the cone depends on the scattering direction defined by θ. When a scattering occurs, a new reference frame is defined by the incident direction of the photon and the orientation of the cylinder. The Stokes vector must be transformed so that the polarization is properly represented with respect to the new reference frame. Then the new Stokes vector of scattered photon at any scattering direction θ can be calculated as follow:

S_{s} (I_{s}, Q_{s}, U_{s}, V_{s}) = M (ζ, θ) R S

where R is the rotational matrix transforming the initial polarization state S to the new incident reference frame,

M (ζ, θ)

is the Mueller matrix whose elements are calculated according to Eq. (2-7). We then obtain the scattering phase functions corresponding to different incident angle ζ and polarization state S by normalizing the Stokes vector

I_{s} (θ)

. For unpolarized incoming photons, the scattering phase functions at different incident angles (Fig. 2(a) ) are identical to the previously published data (Fig. 2(b) in Ref [16].) which was calculated using the unpolarized photon scattering model. Figure 2(b) shows phase functions corresponding to normal incidence but different polarizations. In both figures, we choose 2 µm as the diameter of cylindrical scatterers. The refractive indices of the cylinder and its surrounding are 1.33 and 1.52 respectively. The wavelength of the incident light is 633 nm.

Fig. 2 (a) Phase functions for unpolarized light with different incident angles (b) Phase functions for different polarizations of normally incident light

Download Full Size | PDF

2.1.2 The scattering coefficient

The scattering coefficient of the microspheres μ_s,sph is a constant once their size and refractive index are set. However, the scattering coefficient of cylinders, μ_s,cyl, varies with the angle ζ between the direction of incident photon and the cylinder [16]:

μ_{s, c y l} (ζ) = Q_{s c a} (ζ) \cdot d \cdot C_{A}

where d and

C_{A}

are the diameter and density of the cylinders.

Q_{s c a}

is the scattering efficiency of a single cylinder, which is a function of the angle ζ and the polarization state of the incident photon. In the program, we first determine ζ according to the spatial orientation of the cylinder and the propagation direction of the photon. If the electric field is parallel to the cylinder, efficiency

Q_{s c a, / /} (ζ)

is calculated using the analytic solution of the Maxwell equation [15]. For arbitrary polarization states

S = [I, Q, U, V]

,

Q_{s c a, S} (ζ)

is calculated as follow:

\begin{array}{l} Q_{s c a, S} (ζ) = \frac{W_{s, S}}{W_{s, / /}} \cdot Q_{s c a, / /} (ζ) = \frac{\sum_{θ = 0}^{2 π} I_{s, S} (θ)}{\sum_{θ = 0}^{2 π} I_{s, / /} (θ)} \cdot Q_{s c a, / /} (ζ) \\ = \frac{\sum_{θ = 0}^{2 π} M (ζ, θ) \cdot {(I, Q, U, V)}^{T} (1)}{\sum_{θ = 0}^{2 π} M (ζ, θ) \cdot {(1, 1, 0, 0)}^{T} (1)} \cdot Q_{s c a, / /} (ζ) \end{array}

2.2 The framework of program

The Monte Carlo program is based on a previous version for simulating polarized photon propagation in isotropic media, such as intralipid. The previous program follows the same structure as the public code [1], and takes reference of other published work for the calculation of polarized light scattering at microspheres [4,5].

To simulate anisotropic behavior of the media, the new program introduces infinitely long cylinders into the scattering media. The number of different cylinders as well as their sizes, concentrations and orientation functions are preset parameters in the input data file. New subroutines are added for calculations of cylindrical scatters, such as scattering probability, phase function, rotation of the Stokes vector reference plane etc.

A flow-chart of the program is shown in Fig. 3 . For polarized Monte Carlo, the polarization states of photons are tracked step by step in addition to their spatial positions and direction of motions. The detailed descriptions of the program are as follow:

Fig. 3 Flow chart of Monte Carlo program of polarized light scattering in medium consisting of spheres and cylinders

Download Full Size | PDF

1) The program first establishes the laboratory coordinate system: $X - Y$ plane represents the surface of semi-infinite medium and Z is the depth direction. Two unit vectors $\overset{⇀}{v}$ and $\overset{⇀}{w}$ define the E field reference frame. $\overset{⇀}{v} = [v_{x}, v_{y}, v_{z}] = [0, - 1, 0]$ . A photon of arbitrary polarization state is represented by a Stokes vector $S = [I, Q, U, V]$ in the reference frame. A photon is launched at an initial position and direction determined by the geometry and the incident angle of the beam. The propagation direction of the photon is represented by a unit vector $\overset{⇀}{u}$ , and $\overset{⇀}{w} = \overset{⇀}{u} \times \overset{⇀}{v}$ .
2) Before being scattered, the photon moves to a new position which is calculated based on the extinction coefficient $μ_{t}$ and a pseudo random number δ generated in the interval (0, 1).
$Δ s = - \frac{\ln (δ)}{μ_{t}}$
where $μ_{t} = μ_{a} + μ_{s, s p h} + μ_{s, c y l} (ζ)$ , $μ_{a}$ is the absorption coefficient, $μ_{s, s p h}$ and $μ_{s, c y l}$ are the scattering coefficient of spherical and cylindrical scatterers respectively. $μ_{s, c y l}$ is a function of ζ as mentioned in section 2.1.2. Thus, the new position of the photon is updated as follow:
${\begin{cases} x^{'} = x + u_{x} Δ s \\ y^{'} = y + u_{y} Δ s \\ z^{'} = z + u_{z} Δ s \end{cases}$
As in other Monte Carlo simulations [1,4], we use a weight W to describe the absorption of light. By updating this weight after every step according to $μ_{s} / (μ_{s} + μ_{a})$ , the weight of the photon Wwhich initially equals to 1 becomes
$W_{a f t e r} = W_{b e f o r e} \cdot (\frac{μ_{s}}{μ_{s} + μ_{a}})$
where $W_{b e f o r e}$ is the weight before this step.
3) The photon may hit either a spherical or cylindrical scatter,, so we need to make the choice before the next scattering event occurs. First, we determine ζ according to the spatial orientation distribution of the cylinders and the propagation direction of the photon, and calculate $μ_{s, c y l}$ using Eq. (9). Then considering $p_{c y l} = μ_{s, c y l} / (μ_{s, s p h} + μ_{s, c y l})$ for the probability of scattering at cylinders and $p_{s p h} = 1 - p_{c y l}$ for spheres, we make a statistical choice on what type of scatterer the photon hits.
4) The program carries out rotation of the Stokes vector reference frame and computes the phase function with the respective scattering matrix M. Details are significantly different for scattering at spheres and cylinders. The procedures for spherical scatterers are similar to a published work [4,5]. The procedure for scattering at a cylinder is as follow.
- (a) Once the spatial orientation of cylinder ( $\overset{⇀}{l}$ ) and the direction of photon propagation ( ${\overset{⇀}{u}}_{0}$ ) are known, the incident plane and the new reference frame ( ${\overset{⇀}{v}}_{i}$ ) are determined. ${\overset{⇀}{v}}_{i} = \overset{⇀}{l} \times {\overset{⇀}{u}}_{0}$ .
- (b) The Stokes vector is rotated by an angle β for its reference plane to coincide with the incident plane, as shown in Fig. 4. β is the angle between ${\overset{⇀}{v}}_{i}$ and ${\overset{⇀}{v}}_{0}$ . The rotational matrix is
  $R (β) = [\begin{array}{c} 1 & 0 & 0 & 0 \\ 0 & \cos (2 β) & \sin (2 β) & 0 \\ 0 & - \sin (2 β) & \cos (2 β) & 0 \\ 0 & 0 & 0 & 1 \end{array}]$

Fig. 4 Infinite cylinder obliquely illuminated by a plane wave

Download Full Size | PDF

(c) The new polarization state is calculated as follow,
$S^{'} (θ) = M (ζ, θ) R (β) S$
The phase function and a random number determine the scattering direction ${\overset{⇀}{u}}_{s}$ .
(d) Update the new reference frame of Stokes vector: ${\overset{⇀}{v}}_{s} = \overset{⇀}{l} \times {\overset{⇀}{u}}_{s}$ .
- 5) Repeat the above steps. If the weight W falls below a preset threshold, the photon is absorbed and the calculation terminates. If the photon moves beyond the boundary of the medium, the corresponding Stokes vector is calculated as
  $S_{o u t} = W \cdot S$
  Before the photon is collected by the detector, its polarization state needs to be transformed so that its reference plane coincides with the detector. The first rotation about the direction of photon propagation $\overset{⇀}{u}$ makes $v_{z} = 0$ . The angle of this rotation ε is given by
  $\begin{array}{l} ε = 0 when v_{z} = 0 and u_{z} = 0 \\ ε = \tan^{- 1} (\frac{v_{z}}{w_{z}}) in all other cases \end{array}$
  The second rotation by an angle φ about the Z-axis makes $\overset{⇀}{v}$ parallel to the Y-axis. φ is calculated as
  $\begin{array}{l} φ = - \tan^{- 1} (\frac{u_{y}}{u_{x}}) for backscattered photon \\ φ = \tan^{- 1} (\frac{u_{y}}{u_{x}}) for transmitted photon \end{array}$
  Thus the final Stokes vector is:
  $S_{f i n a l} = R (φ) R (ε) S_{o u t}$
- 6) Repeat above procedures for the next photon until the preset photon number is reached.

2.3 Model of the anisotropic medium

In this work, we approximate the anisotropic medium, such as the fibrous tissues, to a mixture of poly-dispersed microspheres and cylindrical scatterers of infinite length. The sizes and concentrations of the spheres and cylinders can be varied independently. The cylinders are aligned following a Gaussian distribution. The orientation and standard deviation $Δ η$ of the distribution function are also variable. Based on such a simple sphere-cylinder model, more complicated structure can be simulated.

3. Performance tests of the Monte Carlo program

3.1 Computation speed

In a typical test run using a personal computer with Intel Core2 Duo E6300 1.86 GHz processor and 2 G memory, it costs about 50 minutes to complete a simulation with the following parameters: thickness 2 cm, scattering coefficient $100 c m^{- 1}$ , absorption coefficient $1 c m^{- 1}$ , one million photons. With the same number of photons, the thickness and the total scattering coefficient of the medium significantly influence the computation time. The average scattering numbers of photons increase with the thickness or the total scattering coefficient. For polarization imaging simulation, the diffused photons quickly loose their polarization as they scatter around in the medium and may have little contribution in polarization sensitive detection. Thus, one may cut down the computation time by setting up a threshold for the maximum number of scattering a photon may experience.

3.2 Comparison with simulation of unpolarized light scattering in anisotropic medium

To test the validity of the program, we simulate polarized photon back scattering in an anisotropic scattering medium and compared the results with a previously published simulation for unpolarized photon back scattering using the same medium [14]. We use linearly polarized photons with random polarization angles to simulate the unpolarized situation. The scattering medium is semi-infinite and contains infinitely long cylindrical scatters only. The cylinders are aligned in the plane and around the x axis following a Gaussian distribution with standard deviation $Δ η =$ 10°. The diameter of the cylinders is $2 μ m$ and the absorption coefficient is $μ_{a} = 0.1$ cm ⁻¹. The scattering coefficient is chosen as 1420 cm ⁻¹ which was not given explicitly in Ref [14].

In the simulation, linearly polarized photons with random polarization angle incident perpendicularly onto the surface of the medium (the X-Y plane) at $x = y = 0$ . Stokes vectors of the backscattered photons are calculated to obtain the backward spatial distribution of $S = [I, Q, U, V]$ . The dimensionless spatially resolved reflectance $R (x, y)$ is then calculated by normalizing $I (x, y)$ to its maximum. The newly simulated result is shown in Fig. 5 using the same scale and format as Fig. 2(a) in Ref [14]. The plot scales with scattering coefficient. For 1420 cm ⁻¹, the two figures look similar.

Fig. 5 R(x, y) pattern obtained from Monte Carlo simulation of a semi-infinite turbid medium containing cylinders along the axis with Δη = 10°

Download Full Size | PDF

3.3 Comparison with polarization imaging experiments

More tests are carried out to compare the simulated results of reflective polarization imaging against those from experiments.

3.3.1 Experiment setup

The schematic diagram of the experimental setup is shown in Fig. 6 [17]. The 650 nm light from a 1 W LED with proper attenuation is collimated by a lens and passes through a linear polarizer P1 and then illuminates the sample at an angle of 20° to the normal direction of its surface. Backscattered photons from the sample first pass though a linear polarizer P2, then are recorded by a CCD.

Fig. 6 Scheme of experimental setup

Download Full Size | PDF

3.3.2 Imaging and data analysis

In the experiment, the sample surface is set as the X-Y plane and the fibrous structure is in the plane and parallel to the X axis. We record a series of images $I (θ_{i}, θ_{s})$ , where $θ_{i}$ and $θ_{s}$ correspond to the angles of the polarizer and the analyzer respectively, and calculate linear differential polarization (LDP) between orthogonal polarizations [17].

L D P (θ_{i}, θ_{s}) = I (θ_{i}, θ_{s}) - I (θ_{i}, θ_{s} + π / 2)

In an early study [17], we reported that LDP can be expressed analytically as:

L D P (θ_{i}, θ_{s}) = \sqrt{A \cos (4 θ_{s} - φ_{1}) + B} \cdot \cos (2 θ_{i} - φ_{2} (θ_{s})) + C \cos (2 θ_{s} - φ_{3})

Fitting the LDP images pixel by pixel to the above expression results in new images for each of the fitting parameters, which are related to the orders of alignment and the direction of the fibrous structure.

3.3.3 Result of the experiment

In the experiment, the sample is a three-layers medium which consists of a slab of well aligned silk fibers submerged in microsphere solution. The first and third layers are 1 cm thick solutions of $1.5 μ m$ diameter polystyrene microsphere in water (International Laboratory, USA). The refractive index of the microsphere and water are 1.59 and 1.33 respectively. The scattering coefficients are 19.5 cm ⁻¹, 13 cm ⁻¹, 6.5 cm ⁻¹ and 3.25 cm ⁻¹ determined by the concentration of the microspheres. The second layer contains only well aligned silk fibers (provided by Guangxi Institute of Supervision and Testing on Product Quality). The thickness of the silk fiber slab is 0.1 cm and the refractive index is 1.56. The diameter of the cylindrical scatters is taken as $1.5 μ m$ , since both literature and electron microscope images show that silk fibers contains micro fibers of such size..

We simulate the experimental results using both the sphere-only and sphere-cylinder models. The sphere-only model cannot give reasonable good match between the experimental and simulated results.

For the sphere-cylinder model, we first chose the experimental data with fixed microspheres scattering coefficient 3.3 cm ⁻¹, simulate the LDP curve using different parameters for the silk fiber slab until a good match between simulated and experimental results is reached. Then, we simulate the experimental results corresponding to different microsphere scattering coefficient. As shown in Fig. 7 , the simulations of different samples using the sphere-cylinder model agree well with the experiments. The scattering coefficient of silks fiber slab is 70 cm ⁻¹. The standard deviation of the Gaussian distribution for spatial alignment of the fibers is 5°.

Fig. 7 Dependence of Linear Differential Polarization to incident polarization angle

Download Full Size | PDF

4. Conclusion

In summary, we present a Monte Carlo program that simulates the propagation of polarized light in anisotropic media. The medium for the simulation is composed of spherical and infinitely long cylindrical scatterers. Its anisotropic characters can be adjusted by changing the orientation distribution of the cylindrical scatterers, and the composition and the diameter of the two types of scatterers. The program is capable of simulating the behavior of polarized light propagation in anisotropic media, obtaining the spatial distribution and polarization states of the scattered photons. To validate the program, we demonstrated that a simulation using photons of randomly linear polarization agrees well with a previously published simulation using unpolarized photons. We also carried on more tests using samples which contain spheres and cylinders. The good agreement between the simulated and experimental results for such sphere-cylinder samples presents further evidence to the validity of the program.

Acknowledgement

This work has been supported by National Natural Science Foundation of China (grants 60578003 and 60778044), and Ministry of Science and Technology (grant 2006CB70570).

References and links

1. L. H. Wang, S. L. Jacques, and L.-Q. Zheng, “MCML--Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Methods Programs Biomed. 47(2), 131–146 (1995). [PubMed]

2. G. W. Kattawar and G. N. Plass, “Radiance and polarization of multiple scattered light from haze and clouds,” Appl. Opt. 7(8), 1519–1527 (1968). [PubMed]

3. P. Bruscaglioni, G. Zaccanti, and Q. Wei, “Transmission of a pulsed polarized light beam through thick turbid media: numerical results,” Appl. Opt. 32(30), 6142–6150 (1993). [PubMed]

4. J. Ramella-Roman, S. Prahl, and S. L. Jacques, “Three Monte Carlo programs of polarized light transport into scattering media: part I,” Opt. Express 13(12), 4420–4438 (2005). [PubMed]

5. J. C. Ramella-Roman, S. A. Prahl, and S. L. Jacques, “Three Monte Carlo programs of polarized light transport into scattering media: part II,” Opt. Express 13(25), 10392–10405 (2005). [PubMed]

6. S. Bartel and A. H. Hielscher, “Monte carlo simulations of the diffuse backscattering mueller matrix for highly scattering media,” Appl. Opt. 39(10), 1580–1588 (2000).

7. B. D. Cameron, M. J. Rakovic, M. Mehrübeoglu, G. W. Kattawar, S. Rastegar, L. V. Wang, and G. L. Coté, “Measurement and calculation of the two-dimensional backscattering Mueller matrix of a turbid medium,” Opt. Lett. 23(7), 485–487 (1998).

8. M. J. Rakovi, G. W. Kattawar, M. B. Mehru Beo Lu, B. D. Cameron, L. V. Wang, S. Rastegar, and G. L. Coté, “Light backscattering polarization patterns from turbid media: theory and experiment,” Appl. Opt. 38(15), 3399–3408 (1999).

9. D. Côté and I. Vitkin, “Robust concentration determination of optically active molecules in turbid media with validated three-dimensional polarization sensitive Monte Carlo calculations,” Opt. Express 13(1), 148–163 (2005). [PubMed]

10. J. He, A. Karlsson, J. Swartling, and S. Andersson-Engels, “Light scattering by multiple red blood cells,” J. Opt. Soc. Am. A 21, 1953–1961 (2004).

11. J. D. Keener, K. J. Chalut, J. W. Pyhtila, and A. Wax, “Application of Mie theory to determine the structure of spheroidal scatterers in biological materials,” Opt. Lett. 32(10), 1326–1328 (2007). [PubMed]

12. A. Kienle and R. Hibst, “Light guiding in biological tissue due to scattering,” Phys. Rev. Lett. 97(1), 018104 (2006). [PubMed]

13. A. Kienle, C. D’Andrea, F. Foschum, P. Taroni, and A. Pifferi, “Light propagation in dry and wet softwood,” Opt. Express 16(13), 9895–9906 (2008). [PubMed]

14. A. Kienle, F. K. Forster, and R. Hibst, “Anisotropy of light propagation in biological tissue,” Opt. Lett. 29(22), 2617–2619 (2004). [PubMed]

15. C. F. Bohren, and D. R. Huffman, Absorption and Scattering of Light by Small Particles (John Wiley and Sons, New York, 1983).

16. A. Kienle, F. Forster, R. Diebolder, and R. Hibst, “Light propagation in dentin: influence of microstructure on anisotropy,” Phys. Med. Biol. 48(7–N), 14 (2003).

17. X. Y. Jiang, N. Zeng, Y. H. He, and H. Ma, “Investigation of Linear Polarization Difference Imaging Based on Rotation of Incident and Backscattered Polarization Angles,” Pro. Biochem. Biophys. 34, 659 (2007).

Monte Carlo simulation of polarized photon scattering in anisotropic media

Abstract

1. Introduction

2. Methods: Monte Carlo simulation

2.1 Polarized photon scattering at an infinitely long cylinder

2.1.1 The phase function

2.1.2 The scattering coefficient

2.2 The framework of program

2.3 Model of the anisotropic medium

3. Performance tests of the Monte Carlo program

3.1 Computation speed

3.2 Comparison with simulation of unpolarized light scattering in anisotropic medium

3.3 Comparison with polarization imaging experiments

3.3.1 Experiment setup

3.3.2 Imaging and data analysis

3.3.3 Result of the experiment

4. Conclusion

Acknowledgement

References and links

Cited By

Figures (7)

Equations (21)

Optics Express