Development of perturbation Monte Carlo methods for polarized light transport in a discrete particle scattering model.

We present a polarization-sensitive, transport-rigorous perturbation Monte Carlo (pMC) method to model the impact of optical property changes on reflectance measurements within a discrete particle scattering model. The model consists of three log-normally distributed populations of Mie scatterers that approximate biologically relevant cervical tissue properties. Our method provides reflectance estimates for perturbations across wavelength and/or scattering model parameters. We test our pMC model performance by perturbing across number densities and mean particle radii, and compare pMC reflectance estimates with those obtained from conventional Monte Carlo simulations. These tests allow us to explore different factors that control pMC performance and to evaluate the gains in computational efficiency that our pMC method provides.

1. Introduction

The study of in vivo and in vitro cell and tissue microstructure, which yields information about underlying cellular processes and responses, has been historically addressed by optical biopsy techniques. Indeed, there appears to be growing interest in the use of polarization-sensitive optical probes as a tool to non-invasively obtain information regarding subcellular microstructure. Mourant and co-workers have integrated polarized wavelength-dependent measurement capabilities into reflectance based optical probes to increase sensitivity to scatterer size [1]. Backman and co-workers have employed polarization-gating techniques to characterize and control the depth of the probed tissue volume [2]. Sokolov and co-workers have used a probe with angled illumination collection geometry that detects light polarization to obtain depth-specific information of scatterer size and size distribution [3]. These methods, however, all represent instrumentation-based approaches to extract information about the underlying subcellular structure. Here we propose a versatile and rigorous computational method to analyze polarized signals in order to provide insight into subcellular morphology.

The perturbation Monte Carlo (pMC) method is well-established and has been applied to a wide range of biomedical optics problems including: determination of layered tissue optical properties [4, 5], tissue image reconstruction [6], optical tomography [7], and time-resolved functional imaging [8, 9]. Here we provide a novel extension to pMC that builds upon our previous work [10] to characterize the utility of pMC in providing estimates of polarization-sensitive reflectance measurements. Our exploration in this study also serves as a demonstration of one of many applications that can benefit from the use of pMC and as a mechanism for discovering the factors that affect pMC performance for a class of similar problems.

2. Model problem

2.1. Light scattering model

We develop a pMC computational framework within a tissue representation that models optical scattering using three log-normal distributions to represent the polydispersity of the scattering medium [1,11-13]. The pMC method introduced here can accommodate any number (m) of scatterer groups with arbitrary distributions of scatterer size. The scattering coefficient of the ith distribution, μ, is calculated using where Q is the scattering efficiency, N is the product of number density and the scattering cross-section of the particle, t is a particle radius, and L is the ith log-normal distribution characterized by its mean radius and standard deviation, . The scattering coefficient of the medium can be calculated by summing the contribution of each distribution to the scattering coefficient [14]. The phase function of each scatterer population is calculated by integrating across radii values and using the log-normal distribution to properly weight the phase functions that Mie Theory provides: f(θ) = 1/μ(t, θ)Q(t)N(t)dt where f is the phase function of the ith distribution of scatterers, f is the phase function of a mono-dispersed distribution of scatterers, and θ is the scattering angle. We approximate these integrals for μ and f using trapezoid rule integration with uniform grid spacing with the width of the bin set at 0.001 µm for all three distributions.

Table 1 provides scattering parameters that approximate the contributions of protein complexes, organelles, and nuclei present in epithelial cervical tissues based on previous experimental results [1,11]. These parameter values serve as baseline values in the pMC results presented this paper and are meant to model a class of problems (similar to the problem in [11]) in which pMC application could potentially benefit forward and inverse solutions.

Table 1

Distributions of scatter sizes in the tissue model

i	Average Size, r¯_i±σ(μm)	Number Density, Ns,i (cm−3)	nmedium	nscatterer
1	0.03 ± 0.4	4 × 1013	1.33	1.46
2	0.45 ± 0.3	5 × 1010	1.35	1.40
3	4.8 ± 0.1	5 × 107	1.37	1.49

2.2. Measurement geometry

In this study we consider a probe with two detectors placed equidistant from a central source fiber which detect either unpolarized or polarized reflectance. A probe schematic is shown in Fig. 1(c). All detector fibers are angled at ω =20° relative to the outward pointing surface normal and towards the source fiber. The distance from the center of any detector fiber to the center of the source fiber is 550 µm and all the fibers are 480 µm in diameter. All fibers have a light cone detection half angle of 21.7°. The close source-detector separation and the angling of the detector fibers towards the source fibers promotes photon trajectories that interact with scatterers located at shallow depths.

Fig. 1

A conceptual diagram of the pMC application to the light scattering model used in this study. (a) Depiction of the light scattering model composed of three distinct distributions of scatterers. The dotted blue and red plots represent perturbations in mean radius and number density while solid green, blue and red plots represent “baseline” parameter values. (b) A schematic diagram of the source-detector configuration. (c) The top-down view of the probe. The yellow circle represents the source and the circles with numbers inside represent the two detectors. See Section 2.2 for probe details.

3. Theory

The radiative transport equation (RTE) is an integro-differential equation describing the transport of photons in a turbid medium that can be transformed to an integral equation: where , Ψ is the radiance or transport solution, the function Q models some photon source, R is the set of all possible positions in space, and S2 is the set of all solid angles. In this paper, is an integral operator, K is the RTE kernel that characterizes light transport from a location r′ and unit direction Ω′ to another location r and unit direction Ω. To predict the detection of photonic signals, we are typically interested in estimating one or more integrals I with the general form where g is a characteristic detector function and I is some value of interest, often reflectance or transmittance. A rigorous theory that exhibits the equivalence between an RTE-based analytic model, as described by Eqs. (1) and (2), and a Monte Carlo-based stochastic model is treated in detail in [10,15,16]. That equivalence produces the string of equations that relates the stochastic model with the physical model where is a random variable (estimator) on the sample space of biographies and M is a probability measure on . Roughly speaking, the measure M describes a probability density function defined on that can be interpreted as a “likelihood function” for each . This likelihood function is uniquely induced by the specific functions used to generate the biographies in (details can be found in [15]). In Eq. (3)E is the expected value of ξ with respect to the measure induced on by the random walk process used to generate biographies, the second equality in Eq. (3) states that ξ is an unbiased estimator of E[ξ], and the third equality in Eq. (3) demonstrates the equivalence between the analytic model and the stochastic model. This equivalence is necessary to support the derivation of the pMC model.

3.1. Perturbation as change of measure

Suppose that and are two different probability measures that are induced on by two different methods of generating photon biographies . For example, could represent the baseline measure (which for us is generated by the continuous ∈ absorption weighting (CAW) process [16]), and is a different measure representing any perturbation of . If we denote as the CAW-induced measure, the measure in our case would arise by replacing the baseline optical properties by any perturbation of them. The Radon-Nikodym Derivative, , may be used to estimate the reflectance obtained when replaces . The existence of the Radon-Nikodym Derivative can be guaranteed only if certain technical conditions are satisfied, namely, the sets of biographies that are “small” in the measure are also small in the measure. A detected photon measure can also be expressed as an expected value with respect to the modified measure of a modified random variable , where

We may interpret the application of the Radon-Nikodym derivative in the context of pMC as re-weighting the baseline collected photon weights to account for the change in the density functions used to generate the random walks.

3.2. Perturbation Monte Carlo

The pMC method is effective for estimating detected photon weights from small changes in parameters that characterize the forward RTE. When using pMC, a set of Monte Carlo biographies is generated to produce a baseline dataset consisting of the location, weight and incoming unit direction at every collision point for each collected photon. These biographies are then post-processed using pMC equations to determine the effect of perturbations in tissue composition on the detected photon weights. The advantage of this method is that it requires only a single set of Monte Carlo biographies, namely those in the baseline set, to perform simulations for tissues with perturbed properties using only a post-processing module. The computational expense of the post-processing module is 1-3 orders of magnitude smaller compared to that needed to generate a complete set of independent MC biographies.

The CAW random walk process is employed in all MC simulations in this study. Using CAW the photon is deweighted by the factor exp(−μ), where s is the pathlength between collisions and intercollision distances are sampled from . In addition to propagation and absorption, a photon may scatter through a given scattering angle as specified by the phase function. Our study uses a Distribution Selection Scattering Method (DSSM) where at each scattering interaction a photon has a probability of to interact with particles of the ith scattering distribution. Our previous paper [10] uses the Composite Phase Function Scattering Method (CPFSM). A proof of the equivalence between CPFSM and DSSM is provided in the Appendix. Once the scattering distribution has been selected, the phase function corresponding to that distribution, f(θ), will be used to sample the scattering angle, θ. This results in the following RTE kernel [16]:

In Eq. (5), is the probability of interacting with the ith distribution through the DSSM, f(Ω′→ Ω) is the probability density that takes the photon from a direction of Ω′ to Ω for the ith distribution of scatterers, and accounts for scattering probability density and the absorption attenuation associated with a photon pathlength s. This equation assumes a homogeneous absorption coefficient throughout the medium. Application of pMC through the Radon-Nikodym derivative using the kernel in Eq. (5) produces the perturbed estimator : where j1 is the total number of collisions experienced by the photon with respect to the first scattering distribution, j2 is the total number of collisions experienced by the photon with respect to the second scattering distribution, …, j is the total number of collisions experienced by the photon with respect to the last (the mth) scattering distribution of the light scattering model, L is the total pathlength traveled by the photon in the medium, and and are polar and azimuthal scattering angles, respectively, of the lth collision that result from interacting with the 1st distribution. The estimator ξ refers to the collected photon weight from the baseline cMC simulation, which can be polarized parallel or perpendicular to the source. All unhatted quantities refer to baseline parameter values for the quantity and hatted quantities refer to perturbed parameter values. is the total scattering coefficient of the perturbed medium, is the absorption coefficient of the perturbed medium, and is the perturbed phase function of the ith distribution. The sample mean of based on N photon biographies can be used as an estimator for the reflectance and is denoted as . The standard deviation of the sample mean is denoted by .

4. Methods

To test pMC performance, pMC reflectance estimates were compared to those obtained using cMC. cMC reflectance estimates were obtained by launching 20 million photons in a tissue representation characterized by a combination of the parameter values in Table 2 and parameter baseline values in Table 1. The ranges in Table 2 were chosen with the intent of modeling the length scales over which the index of refraction changes in epithelial cells. pMC estimates were obtained by launching 20 million photons in a single Monte Carlo simulation at the baseline parameter values in Table 1 and calculating the sample averages of the estimator to obtain perturbed reflectance values. Error bars shown on the plots represent the standard deviation of the sample means. All photon biographies were obtained from a modified version of the conventional Monte Carlo simulation code outlined in [17].

Table 2

Range of parameters examined in pMC calculations. Baseline and .

Baseline Value, α	Perturbed Parameter, α^	g ranges	μ¯_sˆ ranges (cm−1)
λ = 620 nm	500 – 720 nm	g¯=0.952−0.953	94.5 – 183
Ns,1= 4.00 × 1013 cm−3	2.00 – 6.00 × 1013cm−3	g1 = 0.183	123 – 127
Ns,2= 5.00 × 1010 cm−3	2.5 – 7.5 × 1010cm−3	g2 = 0.958	93.9 – 156
Ns,3 = 5.00 × 107 cm−3	2.5 – 7.5 × 107cm−3	g3 = 0.999	95.5 – 154
r¯_1=0.030μm	0.015–0.045 µm	g^_1=0.0502−0.345	121 – 152
r¯_2=0.45μm	0.38 – 0.52 µm	g^_2=0.945−0.966	94.2 – 174
r¯_3=4.8μm	3.9 – 5.7 µm	g^_3=0.998−0.999	93.8 – 173

cMC and pMC calculations were performed on the High Performance Cluster at University of California, Irvine and a single private computer, respectively. The cluster runs CentOS 6.6 and code on the cluster was compiled using gcc−4.8.2. The private computer has 2 × 1.33 GHz Quad-Core AMD processors and 8 GB 1333 MHz DDR2 RAM running Ubuntu 10.04.4 LTS and our code was compiled using gcc−4.1.3. Benchmark tests were performed to ensure similar performance across the two settings.

Both polarized and unpolarized simulations were performed using the same optical properties and probe geometry. Each cMC unpolarized light simulation launching 20 million photons took 22 – 45 minutes to run while the cMC simulations for polarized light transport took 6.7–18.9 hours. The post-processing module generating pMC reflectance estimates requires only ∼7 seconds and ∼13 seconds for unpolarized and polarized cases, respectively. Table 2 provides the perturbations used for the pMC tests we performed by applying perturbations in mean radii and number densities . Our model has three distributions of scatterers and each distribution has an associated anisotropy factor. We refer to the anisotropy factor of the ith distribution as g and the ensemble anisotropy factor of the medium as . Figure 2 shows how μ and g change as these parameters are changed with the parameter perturbation.

Fig. 2

(a) Changes in the mean radii and how it relates to g. (b) Changes in the ensemble scattering coefficient for all perturbations.

In both cMC and pMC simulations, we treat the source fiber center as the origin of our coordinate system. The center of Detector 1 in Fig. 1(c) resides along the y-axis and the center of Detector 2 resides along the x-axis of the coordinate system. Our pMC implementation requires the biographies of detected photons generated by an initial baseline cMC simulation.

In the case of polarized light transport, polarization is tracked by rotating the Stokes vector as outlined in [18]. All photons initially enter the medium linearly polarized and with polarization parallel to the x-axis. Each scattering event is characterized by polar and azimuthal angles, θ and ϕ, which are sampled from the phase function through a rejection sampling method used by Bartel and Hielscher [19]. Once the photon exits the medium, the Stokes vector is used to calculate the components of the photon weight that are, respectively, parallel and perpendicular to the x-axis.

5. Results

5.1. Scattering model parameter perturbations

First, we investigate the ability of our pMC method to predict changes in reflectance produced by perturbations in the mean radii and/or number density of one of the three particle populations. Figure 3 shows a comparison between cMC and pMC reflectance estimates for Detector 1, which tallies photon weights with parallel or perpendicular polarization relative to the source. Figs. 3(a), 3(c), and 3(e) provide a comparison of pMC and cMC reflectance estimates resulting from perturbations in number density for each of the three scattering particle populations, whereas Figs. 3(b), 3(d), and 3(f) compare pMC and cMC reflectance estimates for perturbations in mean radii.

Fig. 3

pMC and cMC estimates of reflectance of parallel and perpendicular polarization in Detector 1 for perturbations in number density ((a), (c), and (e)) and mean radii ((b), (d), and (f)). Solid lines display the trend for parallel reflectance estimates and dashed lines display the trend for perpendicular reflectance estimates.

The agreement between pMC and cMC estimates in Fig. 3 is linked to the range in which is perturbed and the degree in which is perturbed. For perturbations of the scattering model parameters, the two largest ranges in occur when the second and third mean radii are perturbed, as seen in Table 2. Although perturbations in the second and third mean radii produce similar ranges, the perturbations in the second mean radii coincides with a larger perturbation in anisotropy factor which causes a more pronounced disagreement between pMC and cMC. As shown in Table 2, we explored a perturbation range for mean radius that spans a larger range of ensemble values as compared to perturbations in the number density. For example, perturbations we examined for number density of the first particle distribution result in a range of ensemble scattering coefficients , whereas perturbations we examined for first mean radius in the range result in a broader range of . Similar observations can be made for the second and third populations. Because number density perturbations do not result in perturbation of the phase function and because of the smaller changes in , these pMC reflectance estimates have smaller variance than when perturbing mean radii. From this, we conclude that pMC performance is a function of the magnitude of changes as well as the changes in the phase function which can be characterized by the change in g.

Interestingly, Fig. 3 reveals that reflectance is highly sensitive to changes in the first distribution since the slopes of the plots in Figs. 3(a) and 3(b) are much greater than those in Fig. 3(c)–(f). The sensitivity to the parameters of the first distribution can be explained by the low anisotropy factor of the associated phase function, which heavily promotes backscattering (See Table 2). This probe’s bias toward detecting photons that have undergone backscattering is connected to the geometry of the probe which features angled detector fibers and a small source-detector separation. Based on ratios, we would expect that 3.0%, 49.8%, and 47.2% of all scattering interactions would be with the first, second, and third distributions, respectively. However, analysis of the collected photons reveal that 4.7%, 48.8%, 46.4% of the scattering interactions are with the particles from the first, second, and third distribution, respectively. This represents a 54% increase in the number of interactions that the detected photons have with particles from the first population over the expected number of interactions. This is evidence that collected photons have enriched sensitivity toward changes to parameters associated with the first particle distribution.

Another important observation from the results in Fig. 3 is that the polarization-sensitive measurements provide sensitivity to parameters of the second particle distribution. This result differs from findings in [10] where we saw sensitivity to only parameters of the first distribution with unpolarized reflectance estimates. Changes in the first and second scattering populations, which represent protein complexes and organelles, are associated with dysplasia [11] and sensitivity to parameters of both first and second scattering populations may be helpful. Figure 4 provides a comparison between pMC reflectance estimates for perturbations in the second mean radii for the cases of unpolarized and polarized light propagation and is intended to supplement the results in Fig. 3(d). In Fig. 4(a) and Fig. 4(b), the change in reflectance across values are observable where as in Fig. 4(c), the changes in reflectance are not discernible from the noise in pMC and cMC reflectance estimates.

Fig. 4

pMC estimates for perturbation in second mean radius in cases of (a) parallel polarization (b) perpendicular polarization and (c) unpolarized detection.

Although pMC’s reflectance estimates are particularly poor for Detector 2 with perpendicular polarization, pMC accurately captures the 6% change in reflectance at for Detector 2 with perpendicular polarization. The differences between pMC reflectance estimates in the two detectors and polarizations will be further addressed in §5.2.

In Fig. 5, we explore the ability of pMC to reproduce spectral reflectance measurements. The baseline wavelength is λ = 620 nm and the wavelength perturbations are performed in the range of λ = 500–720 nm. Alteration of the wavelength requires modification of both the scattering coefficients and phase functions for all three scatterer distributions. The spectral asymmetry in Fig. 5 of pMC’s performance relative to the baseline wavelength can be explained by the more rapid increase of the scattering coefficient as the wavelength decreases, as shown in Table 2. We also see that larger wavelength perturbations away from the baseline value of λ = 620nm result in larger standard deviations of the resulting pMC estimates. Perturbations from the baseline wavelength result in perturbations in both the ensemble scattering coefficient and the anisotropy factors, both of which drive increases in the variance of the estimator .

Fig. 5

pMC and cMC estimates of reflectance across wavelength perturbations for polarized light propagation.

To summarize, we find that pMC’s performance is negatively affected by: (1) increasing perturbations in the scattering coefficient and (2) extreme changes in the phase function which can be characterized by anisotropy factor perturbation. We have also demonstrated the utility of our method by showing that turning on polarization can add sensitivity to some parameters and by showing the accuracy of our method in the problem of perturbing across wavelength.

5.2. PMC performance and conventional Monte Carlo relative error

Second, we examine the relative performance of pMC reflectance estimates for parallel and perpendicular polarization as shown in Figs. 3 – 5. The top part of Table 3 shows the parameter ranges in which pMC estimates are within 5% of cMC reflectance estimates. The elements in the top part of Table 3 have a lower bound followed by an upper bound for the parameter values. These bounds are recorded as percentage changes from the baseline value for that parameter. For example, in the column Detector 1, ║ and in the row r1, the parameter range is listed as −50%, +40%. This means that pMC reflectance estimates for Detector 1 with parallel polarization were within 5% of cMC reflectance estimates for the range r1 = 0.015 – 0.042 µm since the baseline value for r1 is 0.03 µm. In the bottom part of Table 3, we compare the relative errors, which is defined as the standard deviation of the estimate divided by the mean, of pMC estimates at baseline parameter values.

Table 3

Relative error and parameter ranges for polarized pMC reflectance estimates that are within 5% of cMC reflectance estimates for each detector.

Parameter	Detector 1, ║	Detector 1, ⊥	Detector 2, ║	Detector 2, ⊥
λ	−12.9%, +16.1%	−6.4%, +12.9%	−6.4%, +12.9%	−6.4%, +6.4%
r1	−50%, +40%	−30%, +30%	30%, +10%	−20%, +10%
r2	−15.6%, +8.9%	−15.6%, +8.9%	−11.1%, +6.7%	−13.3%, 6.7%
r3	−14.6%, +16.7%	−6.25%, +4.2%	−14.6%, +10.4%	−8.3%, +8.3%
Ns,1	−50%, +50%	−50%, +50%	−50%, +50%	−50%, +50%
Ns,2	−40%, +50%	−30%, +40%	−50%, +20%	−50%, +30%
Ns,3	−50%, +50%	−50%, +30%	−40%, +40%	−10%, +20%
	Detector 1, ║	Detector, 1 ⊥	Detector 2, ║	Detector 2, ⊥
rel. err.	6.74 × 10−3	9.39 × 10−3	8.53 × 10−3	1.11 × 10−2

The parameter ranges for accurate pMC estimates are governed by (a) the inherent noise of the baseline cMC simulation and (b) the size of the perturbation characterized by the changes in , , and . The effect of the inherent noise of the baseline cMC simulation can be determined through the comparison of the pMC reflectance estimates between Detector 1, ║ and Detector 1, ⊥. This is because reflectance estimates use the exact same photon biographies and the differences between pMC performance can be attributed to the final weights. The fact that Detector 1, ⊥ provides accurate pMC reflectance estimates over smaller parameter ranges than Detector 1, ║ can be explained by the larger relative error in the baseline cMC simulation for Detector 1, ⊥. Similar arguments can be made for the range of accuracy of the pMC estimates for Detector 2, ║ and Detector 2, ⊥. The accuracy range reported for N1 in Table 3 is the full range of parameter values for N1 in this study; the true limit of pMC’s performance lies beyond the range listed for N1.

The observed disparity between the relative errors of reflectance estimates for Detectors 1 and 2 is consistent with prior experimental work [11, 20, 21] and is consistent with the underlying physics of our model problem, which causes Detector 1 to detect a larger reflectance polarized parallel to the x-axis than Detector 2. This asymmetry in photon collection can be explained by a simple model outlined in [21] for a similar optical probe geometry. This model uses three assumptions: (1) photons undergo only 2 scattering interactions before entering a detector, (2) photons that are collected by Detector 1 have scattering events in the Y-Z plane and photons that are collected by Detector 2 have scattering events in the X-Z plane and (3) photons enter the medium with a trajectory parallel to the z-axis. Figure 6 shows plots of |s1(θ)|2 and (s2θ)|2, where s1 and s2 are elements of the scattering amplitude matrix yielded by Mie Theory. The plot of |s1(θ)|2 (blue) is the scattered irradiance per unit incident irradiance given that the incident light is polarized perpendicular to the scattering plane and the plot of |s2(θ)|2 (red) is the scattered irradiance per unit incident irradiance given that the incident light is polarized parallel to the scattering plane [22]. Since each photon initially starts with polarization parallel to the x-axis, |s1(θ)|2 can be interpreted as the phase function of scattering events in the Y-Z plane and |s2(θ)|2 can be interpreted as the phase function of scattering events in the X-Z plane in the context of this three assumption model. The greatest difference between |s1(θ)|2 and |s2(θ)|2 occurs at a scattering angle of around 90°. In the event that a photon scatters at two 90° angles before entering the detector, then the implication of Fig. 6 is that a photon will end up in the Y-Z plane more frequently.

Fig. 6

A plot of the s1 and s2 components produced by the Mie Scattering Method.

In this section, we have shown that there is a link between the noise of the baseline cMC simulation and the subsequent pMC performance. This link is helpful in understanding the difference in performance between the pMC reflectance estimates for parallel and perpendicular polarizations. Furthermore, to ensure better pMC performance, it is critical to keep the relative error of the baseline cMC simulation as low as possible.

5.3. PMC and cMC computational efficiency comparisons

To examine the computational advantage provided by pMC over cMC, we examine their relative computational efficiencies, η, defined as where is the mean of the estimate, σ2 is the variance, and T is the time required for the calculation [23]. This figure of merit allows for both characterization of the performance of the estimator and the amount of computational resources used by combining the relative error of the estimator and the computational time. Furthermore, this figure of merit is dimensionless since relative error is proportional to where N is the number of photons launched and T is proportional to N; this lack of dimensionality then allows for comparison across different algorithms or estimators. Figs. 7(a) and 7(b) plots the computational efficiency for perturbations in wavelength, mean radii, and number densities. Both plots feature a relatively constant computational efficiency value for cMC estimates, while pMC estimates have the highest computational efficiency when the change in the ensemble scattering coefficient is zero. The gain in computational efficiency for pMC estimates degrades as the perturbation in the ensemble scattering coefficient increases, which results from an increase in standard deviation of the estimates. In Fig. 7(a), plots of the computational efficiency are shown for perturbations in the mean radii whereas in Fig. 7(b), plots of the computational efficiency are shown for perturbations in the number density. A comparison of the pMC computational efficiency plots in Fig. 7(a) and 7(b), reveals that changes in degrade the computational efficiency even though the range in is much smaller than perturbations for and . This degradation is the result of increases in the size of the standard deviation of the estimate, which is caused by the perturbation in anisotropy factor of the phase function of the first population.

Fig. 7

A comparison of the computational efficiency of pMC and cMC estimates for (a) perturbations in and (b) perturbations in . The pink and red symbols show the computational efficiency of pMC estimates. The blue and light blue symbols in the plot above are for the computational efficiencies of cMC estimates.

Recall that pMC is a true perturbation method; that is it works best if the perturbations are “small”. The pMC estimator is unbiased for ALL perturbations, but the larger the perturbation, the larger the statistical error it incurs. (When we say that the pMC estimator is unbiased we mean that on average the estimator will not tend to over-estimate or under-estimate the true value of the parameter.) Thus, in the limit of very large perturbations, it is more efficient computationally to examine independent Monte Carlo simulations of the baseline pMC and at the perturbed values of parameters of interest. As the size of the perturbation shrinks to zero, the advantage of pMC over independent conventional simulations tends to increase.

6. Summary and conclusion

We have developed a new pMC algorithm that provides estimates of linearly polarized reflectance resulting from linearly polarized incident light. Unlike most other pMC algorithms, this algorithm takes into account phase function perturbations due to varying scatterer parameters values or incident wavelength.

The performance of this implementation of the pMC method depends primarily on the magnitude of the perturbation of the anisotropy factor and the ensemble scattering coefficient of the medium. We have applied the algorithm to the problem of accurately solving a forward problem for a scattering model that simulates realistic tissue optical properties and have shown that the pMC method is able to produce estimates that agree in a probabilistic sense with cMC estimates for a limited range of s and g perturbations.

The ability to model polarized light transport is important because polarized light transport measurements can sometimes provide information not available from unpolarized measurements. We have previously shown that ratios of measurements made by this type of probe can differentiate changes in size from changes in concentration in monodisperse populations of scatterers [17] using single wavelength measurements. The scattering media modeled in this paper are more complex. We demonstrate that polarized measurements provide sensitivity to specific morphological parameters (r2, N2), which is not possible with unpolarized measurements. The sensitivity provides the opportunities to resolve inverse problems that may be out of reach when acquiring measurements of unpolarized light. Consequently, light transport parameters may be obtainable from more complicated systems.