2.1.

Diffusion Model

The diffusion equation for the fluence rate (strictly speaking, the diffuse photon density wave, but we use the term fluence interchangeably) in the frequency domain reads:

2.2.

Single Monte Carlo Simulation

For the basics of Monte Carlo simulation, refer to the paper by Jacques and Wang.¹¹ This section shows how the fluences resulting from different OPs are related each other and what scored quantities are needed.

To exploit the spherical symmetry of the geometry, the scoring volumes are concentric shells centered at the origin. The fluence as a function of time and radius then is proportional to the photon weights deposited in the shell, as expressed by the equation:

In Eq. 9, the fluence rate is proportional to the total number of interactions modified by the weight related to the albedo. The accuracy of this algorithm depends on the replacement of a series of numbers $n_k$ for different photons by one number $\overline{n}$ . MC simulation is a statistical technique, and the probability of one photon arriving at $(t,r)$ is obtained by sampling $N$ identical photons, each of them contributing to the sought probability with equal importance. Mathematically,

Equation 9 indicates that the fluence rate is related to the photon interaction density (P/N) and the scattering history $\overline{n}$ . The photon interaction density is a measure of how often a photon visits a particular differential volume in a particular time interval. The proportionality of this quantity to the fluence rate is reasonable, since we expect that the fluence rate at a position will be larger if photons have more chances to “hit” that position. Scattering history is a measure of how many scatterings a photon undergoes just before it reaches a particular space-time point. It is clear that more scatterings will certainly reduce the fluence rate at a position, because fewer photons will arrive at that position. The quantitative relationship between the fluence rate and the scattering history is, however, not linear but is given by a power law, with the albedo coming into play as the base. The previous analysis shows that determining the fluence rate is equivalent to determining the photon interaction density and scattering history. This equivalence can also be applied to other geometries such as semi-infinite and slab geometries.

Furthermore, we point out that $\overline{n}$ is independent of the position for a given value of time and is proportional to time, with the proportionality constant being the product of the scattering coefficient $\left({\mu }_s\right)$ and the speed of light $\left(v\right)$ in the medium. The reason is straightforward. The distance traversed by a photon between two successive scatterings is $v\bullet \mathrm{\Delta }t$ , where $\mathrm{\Delta }t$ is the time elapsed between the two scatterings. The total distance traveled is $vt$ after $\overline{n}$ scatterings. On average, that distance is $\overline{n}$ divided by ${\mu }_s$ . Combining the two equations results in the statement at the beginning of this paragraph. This observation offers convenience in calculating $\overline{n}$ in inverse computation, since one simple equation allows us to get the scattering history for arbitrary position and time, and it is not necessary to parameterize the SMC data.

An alternative method follows the paper by Kienle and Patterson.¹⁵ From the results of a MC simulation performed for a conservative medium characterized by ${\mu }_{sr}$ , fluence can be obtained for arbitrary ${\mu }_s$ and ${\mu }_a$ using Beer’s law and scaling relations:

The proposed method, which has two parameters to fit, looks more complicated than the previous method,¹⁵ which has one parameter to fit the fluence rate data. Actually, the computation burden for both methods is essentially the same, because the proposed method does not need a look-up table for calculating scattering history, which is just the product of light speed and travel time. Look-up tables must be built for interaction density in the proposed method, and fluence rate in conservative medium in the previous method.¹⁵

Time-resolved data are transformed to frequency domain data using Fourier transform as shown in Appendix A in Sec. 5.

2.3.

Delta- $P_{\mathit{1}}$ Model

This section first briefly reviews the derivation following the paper by Venugopalan, You, and Tromberg.⁷ This is extended to obtain the exact solution for a finite spherical source, and it is shown that the previously obtained analytical solution⁷ is approximate and valid only when the dimension of the source is small. The collimation approximation in the derivation is also discussed. The solution has been generalized to the frequency domain (FD), where the functional form resembles its steady-state solution.

2.4.

Isotropic Similarity Model

We assume that the fluence distribution is invariant when a triplet $({\mu }_a,{\mu }_s,g)$ is replaced by another triplet $({\mu }_a^n,{\mu }_s^n,g^n)$ , provided the first-order similarity relation holds:²⁶

Consider the real medium with OPs $({\mu }_a,{\mu }_s,g)$ replaced by the equivalent medium with OPs $({\mu }_a,{\mu }_s^{\prime },0)$ , where ${\mu }_s^{\prime }=(1-g){\mu }_s$ . The transport equation for the equivalent medium in the frequency domain is given by:

The solution for the frequency domain is more complicated and has not been published to our knowledge. We derive the solution by making use of reasonable approximations for the mathematics to be tractable. We first obtain a particular solution ${\phi }_p\left(\vec{r}\right)$ of Eq. 29.

The approximations employed in deriving Eqs. 34, 35 need further comment. In Eq. 34 second- and higher-order terms are lumped into $O\left({\omega }^2\right)$ , which really is $O\left(T{\omega }^2\right)$ , where $T$ is a function of spatial frequency, optical properties, and speed of light. It can be shown that this coefficient $T$ is substantially less than the coefficient in the first-order term for all $k$ , provided ${\mu }_t^{\prime }\gg \omega ∕v$ . Thus the approximation used here is reasonable. The validity of the approximation used in Eq. 35 becomes less predictable, since other coefficients, which are functions of spatial frequency and optical properties, contribute to the coefficient of ${\omega }^2$ , and due to the unknown magnitude of that coefficient, ${\omega }^2$ may not be small compared to the first-order term. Hence Eq. 36, resulting from the prior approximations, may not be correct for certain values of OPs. An alternative method to obtain the fluence and phase is to evaluate the Fourier-type integral [Eq. 33] directly in the complex domain using a proper numerical method. We performed evaluations using the double exponential formula.³² These two methods are compared in Sec. 3.

2.5.

Inversion Scheme

The fluence and phase data were generated using our conventional MC program in steady state and in the frequency domain with modulation frequency of $100\mathrm{MHz}$ . The geometry in the simulation program was the same as that described in the Introduction in Sec. 1. This program allowed for anisotropic light scattering by utilizing the Henyey-Greenstein (HG) phase function. Isotropic light scattering can also be simulated using the same HG phase function by setting the anisotropy factor to zero. For the source, light was emitted in all directions with equal probability (isotropic). It was a point object except for the delta- ${\mathrm{P}}_1$ model, where the light emitting tip was spherical, and the index of refraction (1.33) of the tip was assumed to be matched to the medium. The radial distance ranged from 0.1 to $100\mathrm{mm}$ with an increment of $0.1\mathrm{mm}$ for scoring.

Figure 1 shows the generic inversion scheme. Fluence and/or phase data and trial OPs were used to calculate ${\chi }^2$ using the selected forward model. The trial OPs were adjusted until $\mid ({\chi }_{prev}^2-{\chi }_{cur}^2)∕{\chi }_{cur}^2\mid <{10}^{-15}$ , where ${\chi }_{cur}^2$ and ${\chi }_{prev}^2$ refer to the values for the current and previous iterations. The value of ${10}^{-15}$ was selected because further reduction yielded insignificant change in estimated OPs. ${\chi }^2$ is defined as:

Fig. 1

Schematic representation of generic nonlinear optimization algorithm.

In SMC, the photon interaction density $(P∕N)$ and scattering history $\left(\overline{n}\right)$ were obtained by MC code with OPs $({\mu }_a,{\mu }_s,g)=(0.01{\mathrm{mm}}^{-1},10{\mathrm{mm}}^{-1},0.9)$ . The raw MC data are discrete in position and time. The discreteness has to be extrapolated to continuity before the inverse algorithm can be applied, otherwise the inversion will terminate prematurely with totally unacceptable estimated optical properties. In other words, data in both $\overline{n}$ and $P$ have to be parameterized. The scattering history is independent of position and is linearly related to time. The proportionality constant, although in principle equal to the scattering coefficient times the light speed in the medium as discussed in Sec. 2 was found by using linear regression to be $2.1457∕\sec$ . The calculation of interaction density is more subtle in that it does not appear to be a simple mathematical function of space and time. A practical way to solve this problem is to fit the MC generated interaction density as a function of time to a $20\text{-}\mathrm{deg}$ polynomial at each fixed position. The data for other points can be interpolated using the MatLab function polyval. Although in principle SMC can use any optical properties to calculate phase and fluence at any position, the reduction of dynamic range due to the scaling effects limits its use for large optical coefficients at large distances. Hence a limit of ${\mu }_t=25{\mathrm{mm}}^{-1}$ was set to prevent the algorithm from searching the parameter space where large coefficients lie. For low total interaction coefficients (around $1.00{\mathrm{mm}}^{-1}$ ), the data developed previously were not used, since the scaling factor of 10 severely reduces the valid range of the fitting algorithm. Hence a dataset using $1.00{\mathrm{mm}}^{-1}$ as a reference total interaction coefficient was developed, following the same procedure described before.

3.1.

Verification of Single Monte Carlo

The verification of the theory of SMC was performed in two steps: validation of interaction history and that of scaling relations. To validate the replacement of a series $n_k$ by one universal number $\overline{n}$ , we chose four space-time points at different distances from the source, as shown in Table 1 . We tracked 30,000 photons for each point to get sufficient statistics. The arithmetic mean number of interactions before a photon reaches the chosen point is virtually identical to the fitted Gaussian mean. The spread measured by ${\mathrm{\sigma }}_{\mathrm{gau}}$ decreases steadily to less than 5%. The reduced ${\chi }^2$ is between 0.64 and 2, indicating it is highly likely that $n_k$ is normally distributed, although a small deviation occurs for a position very close to the source. Satisfactory agreement was also seen when the corresponding frequency histograms were drawn compared to the fitted curves (not shown here in the interest of space).

Table 1

Data for testing universality of scattering history n¯ . θari and θgau denotes the arithmetic and Gaussian mean, respectively. X2∕v denotes the reduced chi squared.

To test the scaling relationships, we chose three representative points at 1, 10, and $40\mathrm{mm}$ from the source. Figure 2 shows the comparison between the fluence rates in the time domain with ${\mu }_a=0.01{\mathrm{mm}}^{-1}$ and ${\mu }_s^{\prime }=1.00{\mathrm{mm}}^{-1}$ using traditional Monte Carlo (TMC) simulation and SMC at all three distances. Also shown is the fluence rate calculated by diffusion theory. In Figs. 2a and 2b, where the points are not far from the source, the diffusion theory differs from SMC and TMC at early times, while SMC and TMC are essentially identical. Far from the source, as shown in Fig. 2c, all three methods converge as expected. Various combinations of optical properties as summarized in Table 2 were also used to test the validity of the scaling relations in SMC. The four combinations were chosen to span the range typical of soft tissues in the red or near infrared. All the results from SMC agreed very well with TMC.

Fig. 2

Comparison of the fluence rates in time domain obtained from diffusion model (dotted line), TMC (solid line), and SMC (dashed line) at distances of (a) 1, (b) 10, and (c) $40\mathrm{mm}$ from the source for ${\mu }_a=0.01{\mathrm{mm}}^{-1}$ , ${\mu }_s^{\prime }=1.00{\mathrm{mm}}^{-1}$ . Note that time scale for short distance is reduced to emphasize the difference in fluence rates at an early time.

Table 2

Combinations of optical properties.

In Fig. 3 , the frequency domain fluence and phase were obtained for OPs: ${\mu }_a=0.20{\mathrm{mm}}^{-1}$ and ${\mu }_s^{\prime }=0.50{\mathrm{mm}}^{-1}$ . Since the reduced scattering coefficient is comparable to the absorption coefficient, it is expected that the diffusion theory may not work very well even at large distances. This is evident in Fig. 3, as is the good agreement between SMC and TMC.

Fig. 3

Comparison of the fluence and phase in frequency domain obtained from diffusion mode (dotted line), TMC (solid line), and SMC (dashed line) versus distances for ${\mu }_a=0.20{\mathrm{mm}}^{-1}$ , ${\mu }_s^{\prime }=0.50{\mathrm{mm}}^{-1}$ .

3.2.

Effect of Finite Dimension of the Source on Fluence in Delta- ${\mathrm{P}}_1$ Model

First we investigate whether the delta- ${\mathrm{P}}_1$ model describes fluence better than the SDA. It has been shown ^{5, 6, 7, 8, 9} that the delta- ${\mathrm{P}}_1$ model works better than the SDA in the region proximal to the source and in a transition region for a wide range of albedos. This may not be true far from the source, where the uncollided beam is negligible and diffusion dominates. The delta- ${\mathrm{P}}_1$ model recovers Beer’s law in the far region only when absorption is dominant, but when scattering is comparable to or larger than absorption, it resembles the SDA more than the true MC result. In this respect, the delta- ${\mathrm{P}}_1$ model is most valuable in dealing with the fluence and current near the source. Figures 4a and 4c show the fluence and phase, respectively, as functions of distance using the delta- ${\mathrm{P}}_1$ model, MC simulation, and the SDA for OPs: ${\mu }_a=0.2{\mathrm{mm}}^{-1}$ and ${\mu }_s^{\prime }=0.5{\mathrm{mm}}^{-1}$ . The fluence plot shows that the delta- ${\mathrm{P}}_1$ model is closer to MC than the SDA for distance less than $0.5\mathrm{mm}$ [see Fig. 4b]. Beyond this range, the delta- ${\mathrm{P}}_1$ model deviates from MC and starts to resemble SDA. For the phase, the delta- ${\mathrm{P}}_1$ model provides a worse description than SDA at all distances (including the region close to the source).

Fig. 4

(a) and (b) Fluence and (c) phase versus distance from the source. Plot (b) shows the behavior close to the source. Fluence (phase) obtained from delta- ${\mathrm{P}}_1$ model (solid line), MC (dashed line), and diffusion model (dotted red line) as the result of OPs: ${\mu }_a=0.2{\mathrm{mm}}^{-1}$ , ${\mu }_s^{\prime }=0.5{\mathrm{mm}}^{-1}$ (color online only).

The second term in Eq. 23 represents the effect of source dimension on the fluence distribution. Two sets of OPs[ $({\mu }_a=0.01{\mathrm{mm}}^{-1},{\mu }_s=10.0{\mathrm{mm}}^{-1},g=0.9)$ and $({\mu }_a=0.01{\mathrm{mm}}^{-1},{\mu }_s=0.3{\mathrm{mm}}^{-1},g=0.9)$ ] were chosen to investigate how OPs affect the fluence distribution. MC simulations for each set of OPs with three source radii (0.2, 1, $2\mathrm{mm}$ ) were conducted, and the fluences were compared with those calculated from the delta- ${\mathrm{P}}_1$ model. In Fig. 5 , relative difference in fluence is defined as:

Fig. 5

Relative difference in fluence versus radial distance as obtained from the delta- ${\mathrm{P}}_1$ model without the contribution of the second term in Eq. 23 compared with MC for different source sizes. The radii of the source spheres are 0.2, 1, and $2\mathrm{mm}$ , respectively, (a) and (b) were obtained with OPs (0.01,10,0.9) and (0.01,0.3,0.9), respectively.

In the following discussion, we classify the medium according to its value of reduced scattering coefficient, i.e., relatively low scattering in Table 3 and relatively high scattering in Table 4 . The reason for this classification is that it appears delta- ${\mathrm{P}}_1$ and ISM models show different fluence and phase behavior for high and low scattering media. Examining the low scattering medium, although not commonly encountered in biological tissues, allows us to thoroughly investigate various models.

Table 3

Optical properties: scattering is medium or low.

Table 4

Optical properties: scattering is medium or high.

3.3.

Forward Calculation Using the Delta- $P_{\mathit{1}}$ Model in the Frequency Domain

Forward calculation using the delta- ${\mathrm{P}}_1$ model was performed to describe the fluence and phase delay relative to the source. A collection of OPs with medium and low scattering as listed in Table 3 were used. Figure 6 plots the percent relative fluence and phase difference compared with those from MC for the medium with ${\mu }_a=0.005{\mathrm{mm}}^{-1}$ , ${\mu }_s^{\prime }=1{\mathrm{mm}}^{-1}$ . It demonstrates that SDA and the delta- ${\mathrm{P}}_1$ model both describe the phase and fluence well for the region between 10 and $50\mathrm{mm}$ from the source. Closer than $10\mathrm{mm}$ , the delta- ${\mathrm{P}}_1$ model still performs better than does SDA but gives errors up to 25% for phase and 15% for fluence. As the albedo decreases, similar results hold until the ${\mu }_s^{\prime }∕{\mu }_a$ ratio reaches around 20, as shown in Fig. 7 for the medium with ${\mu }_a=0.04{\mathrm{mm}}^{-1}$ and ${\mu }_s^{\prime }=0.80{\mathrm{mm}}^{-1}$ . The delta- ${\mathrm{P}}_1$ model gives less than 5% error in phase for the entire region. However, the delta- ${\mathrm{P}}_1$ model gives less than 10% error in fluence only for the region close to the source (for this albedo, less than $40\mathrm{mm}$ ). It appears that the discrepancy in fluence compared to MC simulation becomes larger as the radial distance increases and grows indefinitely. The SDA shares features of the delta- ${\mathrm{P}}_1$ model: increasing discrepancy in fluence with increasing radial distances and good prediction of phase, except in the region very close to the source. When the ${\mu }_s^{\prime }∕{\mu }_a$ ratio is decreased to close to or less than 1, shown in Fig. 8 for the medium ${\mu }_a=0.03{\mathrm{mm}}^{-1}$ , ${\mu }_s^{\prime }=0.05{\mathrm{mm}}^{-1}$ , the fluence estimated by the delta- ${\mathrm{P}}_1$ model is in error by less than 25% for radial distances less than about $25\mathrm{mm}$ . SDA’s description of fluence is in error by more than 25% at $5\mathrm{mm}$ and up to 85% just at the source. Hence the delta- ${\mathrm{P}}_1$ model is better than SDA in terms of describing fluence close to the source for media with low albedo. The description of phase by the delta- ${\mathrm{P}}_1$ model in this case is inaccurate with maximum deviation 350% close to the source, leveling off at 10% beyond $15\mathrm{mm}$ . We have plotted the absolute phase against the radial distance and discovered that the phase was negative for distance less than $4\mathrm{mm}$ . This is unphysical because it means that photons were detected prior to their being launched from the source. This is indicative of the limitations of the collimation approximation [Eq. 18]. From the given OPs, ${\mu }_t^{*}$ can be calculated to be $0.125{\mathrm{mm}}^{-1}$ . The collimation approximation requires $r\ge 2∕{\mu }_t^{*}=16\mathrm{mm}$ for the delta- ${\mathrm{P}}_1$ model to be valid. Phase description by SDA is in error by less than 10% for distances greater than $15\mathrm{mm}$ . Thus for low scattering media, the delta- ${\mathrm{P}}_1$ model gives a better description of fluence than SDA for low and high albedos, but the phase description is worse than that of SDA for low albedos, for all radial distances.

Fig. 6

Relative difference in fluence and phase calculated from SDA and delta- ${\mathrm{P}}_1$ models compared with Monte Carlo simulation $(g=0.9)$ as a function of position. ${\mu }_a=0.005{\mathrm{mm}}^{-1}$ , ${\mu }_s^{\prime }=1.00{\mathrm{mm}}^{-1}$ .

Fig. 7

Relative difference in fluence and phase calculated from SDA and delta- ${\mathrm{P}}_1$ models compared with Monte Carlo simulation $(g=0.9)$ as a function of position. ${\mu }_a=0.04{\mathrm{mm}}^{-1}$ , ${\mu }_s^{\prime }=0.80{\mathrm{mm}}^{-1}$ .

Fig. 8

Relative difference in fluence and phase calculated from SDA and delta- ${\mathrm{P}}_1$ models compared with Monte Carlo simulation $(g=0.9)$ as a function of position. ${\mu }_a=0.03{\mathrm{mm}}^{-1}$ , ${\mu }_s^{\prime }=0.05{\mathrm{mm}}^{-1}$ .

The preceding paragraph discussed how the delta- ${\mathrm{P}}_1$ model performed for media with medium and low scattering. In reality, we often encounter biological tissues with relatively high scattering. A collection of OPs with medium and high scattering is shown in Table 4. The delta- ${\mathrm{P}}_1$ model and SDA demonstrated similar trends, as shown in Fig. 7 (with ${\mu }_a=0.04{\mathrm{mm}}^{-1}$ and ${\mu }_s^{\prime }=0.80{\mathrm{mm}}^{-1}$ ), for $1<{\mu }_s^{\prime }∕{\mu }_a<20$ , except that phase description by the delta- ${\mathrm{P}}_1$ model shows improved precision compared to SDA close to the source, but a somewhat worse performance in the far region. For ${\mu }_s^{\prime }∕{\mu }_a>20$ , the delta- ${\mathrm{P}}_1$ model rapidly converges to SDA. For high scattering media, the delta- ${\mathrm{P}}_1$ model and SDA give similar descriptions of fluence for all radial distances. Phase calculated by the delta- ${\mathrm{P}}_1$ model is better than that by SDA close to the source, but worse in the far region.

3.4.

Forward Calculation Using the Isotropic Similarity Model in the Frequency Domain

First we used MC to investigate if the similarity relation [OPs $({\mu }_a,{\mu }_s,g)$ replaced by the equivalent OPs $({\mu }_a,{\mu }_s^{\prime },0)$ , where ${\mu }_s^{\prime }=(1-g){\mu }_s$ ] is a reasonable approximation for describing fluence and phase. We investigated how fluence and phase behaved if $g$ changed from 0.9 to 0 for a fixed ratio ${\mu }_s^{\prime }∕{\mu }_a=100$ , where ${\mu }_a=0.01{\mathrm{mm}}^{-1}$ and ${\mu }_s^{\prime }=1.00{\mathrm{mm}}^{-1}$ . Six chosen values of $g$ were 0.9, 0.7, 0.5, 0.3, 0.1, and 0.0. Since the fluence ranges over several orders of magnitude, we defined the relative fluence (phase), shown in Fig. 9 , as the ratio of the difference of absolute fluence (phase) from MC using $g=0.9$ and that from MC using $g=0.0$ (0.3, 0.7) to the absolute fluence (phase) from MC using $g=0.9$ . The maximum relative difference of fluence of 11% occurred at less than $0.5\mathrm{mm}$ from the source for $g=0$ . The error rapidly dropped to less than 2% at $2\mathrm{mm}$ . Beyond that distance, the relative difference in fluence is less than 1% for all values of $g$ . Similar behavior was observed for the relative phase, for which the maximum difference of 2% also occurred $2\mathrm{mm}$ from the source. Larger fluctuations for distances over $50\mathrm{mm}$ were due to the intrinsic noise of the MC simulations. We also examined how the relative difference in fluence and phase (percent difference in fluence obtained in the medium with $g=0.9$ and that obtained in the medium with $g=0.0$ ) behaved as ${\mu }_s^{\prime }∕{\mu }_a$ changes. Table 3 summarizes relevant optical parameters and the values of ${\mu }_s^{\prime }∕{\mu }_a$ , which spanned more than three orders of magnitude. For ${\mu }_s^{\prime }∕{\mu }_a>10$ , a typical graph that demonstrates the relative difference of fluence and phase is shown in Fig. 10 for ${\mu }_a=0.04{\mathrm{mm}}^{-1}$ , ${\mu }_s^{\prime }=0.80{\mathrm{mm}}^{-1}$ . The difference of phase has a maximum of 14% at the source and decreases to less than 1% beyond $5\mathrm{mm}$ . The difference in fluence is 10% at the source, rapidly decreases to 1% at about $2\mathrm{mm}$ , then gradually increases with a certain slope. However, for ${\mu }_s^{\prime }∕{\mu }_a>100$ , this slope is virtually zero. As ${\mu }_s^{\prime }∕{\mu }_a$ decreases, the slope slowly increases until ${\mu }_s^{\prime }∕{\mu }_a=10$ , where different behavior is observed (see next). It should be noted that the increase is so slow that for distances of practical interest $(<20\mathrm{mm})$ , the maximum error in fluence is only a few percent. It appears that the relative errors for phase are smaller than those for fluence for high scattering coefficients. Figure 11 plotted with ${\mu }_a=0.01{\mathrm{mm}}^{-1}$ and ${\mu }_s^{\prime }=0.03{\mathrm{mm}}^{-1}$ shows the typical trend for the relative differences of fluence and phase when $0.5<{\mu }_s^{\prime }∕{\mu }_a<10$ . The difference in phase reaches its maximum value at about $5\mathrm{mm}$ and gradually decreases away from the source. The difference of fluence reaches its maximum value at about $9\mathrm{mm}$ and also gradually decreases with distance. It is interesting to note that for this domain, the relative errors in phase exceed those in fluence, just the opposite to the conclusion drawn for media with high albedo.

Fig. 9

Relative difference in (a) fluence and (b) phase obtained from MC using $g=0.7$ (solid line), from MC using $g=0.3$ (dashed line), and from MC using $g=0.0$ (dash-dotted line), compared with that obtained from Monte Carlo simulation $(g=0.9)$ as a function of position. ${\mu }_a=0.01{\mathrm{mm}}^{-1}$ , ${\mu }_s^{\prime }=1.00{\mathrm{mm}}^{-1}$ .

Fig. 10

Relative difference in fluence and phase obtained from ISM model compared with that obtained from Monte Carlo simulation $(g=0.9)$ as a function of position. ${\mu }_a=0.04{\mathrm{mm}}^{-1}$ , ${\mu }_s^{\prime }=0.80{\mathrm{mm}}^{-1}$ .

Fig. 11

Relative difference in fluence (bottom curve) and phase (top curve) obtained from ISM model compared with that obtained from Monte Carlo simulation $(g=0.9)$ as a function of position. ${\mu }_a=0.01{\mathrm{mm}}^{-1}$ , ${\mu }_s^{\prime }=0.03{\mathrm{mm}}^{-1}$ .

Having explored the validity of the similarity relationships using MC, we then performed analytical calculations using Eq. 37 to obtain the amplitude and phase of the complex fluence for a number of OPs whose ${\mu }_s^{\prime }∕{\mu }_a$ ratios range from 1 to 100. It was found that the fluence amplitude calculated from the ISM model agreed with that obtained from MC $(g=0.0)$ for all these OPs across the distances investigated (0.1 to $60\mathrm{mm}$ ) with less than 0.5% error. However, the estimate of the phase using the analytical expressions was less satisfactory, agreeing well closer to the source but deviating from MC in the far region for relatively high ${\mu }_s^{\prime }$ . This supports the claim that the coefficient in the second-order term in Eq. 35 does play a role in estimating the phase. We also performed a numerical calculation of Eq. 33 using the double exponential formula. It was found that both analytical and numerical methods gave satisfactory descriptions of fluence and phase within $10\mathrm{mm}$ for all OPs listed in Tables 3, 4. The fluence and phase calculated by the analytical solution start to deviate from those of the numerical method beyond $10\mathrm{mm}$ when ${\mu }_s^{\prime }>0.3{\mathrm{mm}}^{-1}$ ( ${\mu }_a$ is usually less than ${\mu }_s^{\prime }$ ). To be cautious, numerical integration was invoked for this range of OPs.

3.5.

Inverse Calculation Using Four Models

Since we are interested in interstitial detection of light, detector locations are typically less than $1\mathrm{cm}$ from the source. In this work, two detector locations were chosen to be 4 and $9\mathrm{mm}$ for low absorption media as shown in Table 5 , and 3 and $6\mathrm{mm}$ for high absorption media as shown in Table 6 . Shorter distances were chosen for high absorption media because the signal becomes virtually undetectable in practice. For medium and low scattering media, recovery of reduced scattering and absorption coefficients is within 5% for ${\mu }_s^{\prime }∕{\mu }_a$ greater than 10 for all models (Table 5). When ${\mu }_s^{\prime }∕{\mu }_a$ is less than 10, the SDA and delta- ${\mathrm{P}}_1$ models break down completely, because SDA cannot give good estimates of fluence, and delta- ${\mathrm{P}}_1$ is unable to give correct phase values. The ISM model fares better, with the estimate of ${\mu }_a$ within 70% and ${\mu }_s^{\prime }$ within 64%. This suggests that even though the difference of fluence and phase is within 20%, the inverse calculation still gives fairly large errors in recovering OPs. Hence, inverse calculation is sensitive to the small variations of fluence and phase. For medium and high scattering media, all four models perform very well for ${\mu }_s^{\prime }∕{\mu }_a$ greater than 10. When the albedo is decreased, SDA gives a maximum error of about 25% for ${\mu }_s^{\prime }∕{\mu }_a$ as low as 1.11. The ISM model gives 17% error of the estimate of ${\mu }_s^{\prime }$ and 5.1% error of the estimate of ${\mu }_a$ . The delta- ${\mathrm{P}}_1$ model performs slightly better than diffusion theory but worse than the ISM. It is not surprising that the SMC model works very well for all types of media, since no approximation has been used in the forward calculation (small deviations are due to interpolation errors). It is generally thought that SDA does not perform well if ${\mu }_s^{\prime }∕{\mu }_a$ is not large, but from our results, the inverse calculation using SDA as a forward model can recover ${\mu }_s^{\prime }$ and ${\mu }_a$ both within 10% for ${\mu }_s^{\prime }∕{\mu }_a$ as low as 5. It can also be seen that SDA does not apply well for distances within several mean free paths, as shown in Table 5. Fortunately, most biological tissues have scattering coefficients greater than $0.5{\mathrm{mm}}^{-1}$ . Hence, diffusion theory should be adequate in most practical circumstances. It is interesting to note that the delta- ${\mathrm{P}}_1$ model behaves about the same as the diffusion theory in recovering OPs, in contrast to findings in other papers. The discrepancy may be due to the fact that: 1. the collimation approximation is not good enough for a point source; and 2. the measured quantity in our case is fluence instead of irradiance. For the latter, since the fluence is the integral of the radiance over $4\pi$ while the irradiance is the integral over some forward solid angle, it is likely that the delta- ${\mathrm{P}}_1$ model may provide a more accurate description of the radiance in the forward direction. Hence, although improvement over SDA has been reported for irradiance,^{7, 8, 23} this may not be true for fluence.

Table 5

Recovery of optical properties (medium and low scattering) for various μs′∕μa using four models. Two positions are 4 and 9mm from the source.

Table 6

Recovery of optical properties (medium and high scattering) for various μs′∕μa using four models. Two positions are 3 and 6mm from the source.