Translator Disclaimer
1 November 2007 Determination of the optical properties of turbid media using relative interstitial radiance measurements: Monte Carlo study, experimental validation, and sensitivity analysis
Author Affiliations +
Interstitial quantification of the optical properties of tissue is important in biomedicine for both treatment planning of minimally invasive laser therapies and optical spectroscopic characterization of tissues, for example, prostate cancer. In a previous study, we analyzed a method first demonstrated by Dickey et al., [Phys. Med. Biol. 46, 2359 (2001)] to utilize relative interstitial steady-state radiance measurements for recovering the optical properties of turbid media. The uniqueness of point radiance measurements were demonstrated in a forward sense, and strategies were suggested for improving performance under noisy experimental conditions. In this work, we test our previous conclusions by fitting the P3 approximation for radiance to Monte Carlo predictions and experimental data in tissue-simulating phantoms. Fits are performed at: 1. a single sensor position (0.5 or 1 cm), 2. two sensor positions (0.5 and 1 cm), and 3. a single sensor position (0.5 or 1 cm) with input knowledge of the sample's effective attenuation coefficient. The results demonstrate that single sensor radiance measurements can be used to retrieve optical properties to within ~20%, provided the transport albedo is greater than ~0.9. Furthermore, compared to the single sensor fits, employing radiance data at two sensor positions did not significantly improve the accuracy of recovered optical properties. However, with knowledge of the effective attenuation coefficient of the medium, optical properties can be retrieved experimentally to within ~10% for an albedo greater or equal to 0.5.



The determination of tissue optical properties has long been a subject of interest in the field of biophotonics. In particular, in-vivo measurements of tissue optical properties can be incorporated into treatment planning models for either pretreatment planning or online optimization of novel interstitial laser therapies such as photodynamic therapy (PDT) and laser interstitial thermal therapy (LITT). 1, 2, 3, 4 For therapies such as PDT and LITT, the properties of interest are the absorption coefficient μa , and the reduced scattering coefficient μs . Although surface reflectance measurements can be utilized for noninvasive characterization of tissues such as breast and skin, such methods are limited by light penetration when diagnosing internal organs such as the prostate. In such cases, minimally invasive fluence (i.e., integrated light intensity) point probes inserted directly into the organ of interest have commonly been utilized.5, 6

Light attenuation one or two mean free paths from the source is characterized by the transport attenuation coefficient μt=μa+μs , and is primarily a function of the reduced scattering coefficient, while attenuation far from the source is described by the effective attenuation coefficient,

is a nonunique combination of both absorption and scattering.7

Reflectance techniques provide a well-sampled dataset of radially resolved surface measurements. This allows for the assessment of both μt and μeff . As such, optical properties can be uniquely determined in reflectance geometry using relative spatially resolved reflectance measurements. 7, 8 However, unlike the noninvasive geometry of reflectance techniques, interstitial characterization of optical properties is complicated by two factors.

First, patient invasiveness limits the number of clinically practical sources and sensors that can be employed during minimally invasive treatments. Second, interstitial placement of sources and sensors are typically performed using guidance templates (similar to those employed for brachytherapy), which utilize a square grid with equidistant spaced insertion points. Since the sources and sensors are inserted parallel to each other, their closest radial separation is largely determined by the minimum grid spacing (0.5cm) . While these factors do not severely limit the characterization of μeff , they can lead to significant difficulties in accurately determining μt . Hence, a particular challenge of interstitial geometries is to devise strategies for accurately determining optical properties using measurements far from the source, where knowledge of μeff alone cannot be used to uniquely separate absorption and scattering.9, 10

Absolute (calibrated) and multispectral measurements have been developed to overcome the limitations of interstitial geometry.5, 11 These techniques typically require measurements at multiple source-detector positions to separate optical properties. However, the necessity for spatially resolved data results in a larger sampling volume and, overall, poorer spatial characterization of the tissues, which may be of clinical consequence when significant heterogeneity exists in the target organ.11 The development of techniques that minimize the number of required spatial measurements and provide local information on optical properties may, therefore, allow for improved spatial resolution and greater specificity in patient planning. One promising interstitial method, reported by Xu and Patterson12 employed frequency-domain measurements that utilized only relative measurements of intensity and phase shift at two positions (relative to a modulated point source) to accurately determine optical properties over a wide range of clinically relevant optical properties and sensor positions. In general, the greater the information content of the measurements, the easier it is to separate optical properties.

Recently, point radiance measurements have been proposed as a potential alternative to the previously described fluence-based methods. The concept, first described by Dickey, 13 measures radiance data at a single spatial location as a function of detection angle relative to a point source, and fits the data using the P3 approximation to the radiative transport equation for radiance to separate absorption and scattering properties. Here, the additional information provided by radiance measurements comes from resolving the fluence into separate directional components. Since radiance is collected by rotation instead of spatial translation, the technique provides a potential means to improve the information content of steady-state measurements, while simultaneously reducing invasiveness.

While the initial study13 indicated the significant potential of radiance-based optical property determination, a number of important issues remain to be investigated. For example, the experimental work was performed only for a single high-albedo optical property set and at a distal sensor position (1cm) . These conditions may not always be relevant to the clinical environment where albedos may be as low as 0.5 and sensors positioned as close as 0.5cm .

To partially address these issues, Chin, Whelan, and Vitkin 10 recently performed a forward analysis using the P3 approximation to establish conditions of uniqueness for point radiance information, with and without the presence of noise, over a broader transport albedo range ( 0.5 to 1). Their work demonstrated that under ideal noiseless conditions, the absorption and scattering properties of tissue were separable using relative radiance measurements at a single position. However, the same analysis demonstrated that the addition of noise significantly degrades quantification accuracy, particularly at low albedos (less than 0.85 ). Given the sensitivity to noise, the same authors presented possible strategies for improving optical property retrieval. Their analysis is yet to be verified by inversion of either synthetic (i.e., Monte Carlo) or experimental data.

In this study, we further explore the potential of radiance measurements for optical property characterization. Specifically, we evaluate the accuracy of the P3 approximation as a radiance inversion model for optical property determination over a wide albedo range. We further test the methodology initially proposed by Chin, Whelan, and Vitkin10 that employs an additional μeff constraint to the radiance fits that greatly improves the inversion performance compared to unconstrained fits.

This work is organized as follows. In Sec. 2, we review the theory and methods utilized in the study. The basic concepts of radiance uniqueness as pertaining to optical property retrieval are examined, and the fundamental basis of the μeff constraint approach is presented. In Sec. 3 we evaluate the limits of the P3 approximation as an inversion model for recovering optical properties using point radiance measurements, by performing fits to synthetically generated Monte Carlo (MC) data. The analysis allows a direct assessment of the P3 algorithm because: 1. the optical properties μa and μs are known exactly, and 2. the inversion is performed free of “experimental” errors in sensor positioning or detection angle. Both the unconstrained and μeff -constrained techniques are evaluated. After testing the radiance techniques under “ideal” conditions, in Sec. 4 we consider potential errors that may arise due to systematic differences from the ideal case. Using simulated data, we perform a sensitivity analysis of the radiance inversion algorithm to experimentally relevant finite probe size and numerical aperture. The results indicate the important experimental factors that must be accounted for during radiance fitting. In Sec. 5, we experimentally confirm the constrained and unconstrained radiance approach in tissue-simulating phantoms. We also compare the radiance approach to other techniques that have been employed for determining optical properties under interstitial geometry. Finally, we conclude with a brief discussion of the implications of this work.


Methods and Theory for Radiance Optical Property Determination


Radiative Transfer Equation

Light propagation in turbid media has been shown to be well described by the radiative transport equation (RTE)14:

Eq. 1

In Eq. 1 we introduce the radiance L(r,ŝ) , which is a scalar quantity representing the radiated power per unit area per unit solid angle ω at a position r due to photons traveling in direction ŝ . Integration of the radiance over all solid angles gives the fluence ϕ(r)=4πL(r,ŝ)dω , which describes the total light intensity converging on position r . Both the radiance and fluence arise from a source distribution S(r,ŝ) irradiating the turbid medium. The resulting light distribution is a function of the medium’s optical properties.

In the RTE, the optical properties of interest are the absorption coefficient μa and scattering coefficient μs , which are the probability of absorption and scattering per unit path length, respectively. In addition, the directionality of scattering must also be considered. This is described by the scattering phase function p(ŝ,ŝ) , which is the probability of a photon traveling in direction ŝ scattering into direction ŝ . While precise information on p(ŝ,ŝ) in tissues would be desirable, Jacques15 demonstrated that the Henyey-Greenstein (H-G) phase function16 provides a convenient analytic approximation of single-scattering angular distributions for most tissues. All expressions employed in this work utilized the H-G phase function.


P1 Approximation

Although exact analytical solutions can be obtained for Eq. 1, they exist only for idealized geometries that are too restrictive for most practical applications. Simplifying assumptions are necessary to obtain analytic solutions to the RTE.

The most commonly employed approach, the PN approximation, expands L(r,ŝ) , p(ŝ,ŝ) , and S(r,ŝ) in spherical harmonics to the N ’th order.17 The N=1 expansion expresses the radiance L(r,ŝ) as the linear combination of the isotropic fluence ϕ(r) and a linearly anisotropic flux j(r) 18:

The N=1 expansion is also known as the diffusion approximation to the RTE.

In spherical geometry, the resulting diffusion theory (DT) radiance expression for a point source is18:

Eq. 2

LP1(r,θ)=LDT(r,θ) =Po4π14πrD[1+3(Dr+μaD)cosθ]exp(μeffr).
An illustration of the described geometry is shown in Fig. 1 .

Fig. 1

Geometry for the radiance L , fluence ϕ , positive flux j+ , and negative flux j at position r relative to a point source in a scattering medium with optical properties μa and μs . Examples of the radiance are shown for detection angles of θ=0 and 180deg . The detection angle for the radiance traveling in direction ŝ is determined relative to the radial normal vector r .


Here, D=13(μs+μa) is the diffusion coefficient, Po is the source power, μeff=[3μa(μa+μs)]12 is the effective attenuation coefficient, and θ is the angle between the radiance direction ŝ and the outgoing radial normal vector r .

In Eq. 2 we have also introduced the reduced scattering coefficient μs=(1g)μs , where g=4πp(ŝ,ŝ)(ŝŝ)dω is also known as the anisotropy factor. The reduced scattering coefficient allows the scattering coefficient μs in an anisotropically scattering medium to be replaced by an equivalent isotropic scattering parameter μs . As such, the parameter μs can be thought of as the probability per unit path length that the direction of a photon will be completely randomized. In most practical applications, it is μs and not μs that is the quantity of interest.

The diffusion approximation provides a simple framework for analyzing point radiance measurements (see Secs. 2.4 through 2.5). However, for most practical interstitial situations, the N=3 expansion (i.e., P3 approximation) is needed to accurately describe the radiance10 in a forward sense.


P3 Approximation

In the current work, the P3 approximation was chosen as an inversion model due to its validity in describing radiance under anisotropic conditions typical of interstitial geometries.10 The choice of the P3 approximation was also based on the computational convenience of the higher order similarity relations developed by Hull and Foster.19 An additional consideration for choosing the P3 approximation is that inversion strategies based on diffusion theory are intuitively and directly applicable to the higher order P3 approximation. Other models such as the δ -P1 approximation,20 the modified spherical harmonics method,12 or even the P5 approximation,17 are also potential candidates for inversion applications.

The P3 radiance solution for a point source in spherical geometry is19:

Eq. 3

where ν and ν+ are analagous to the effective attenuation coefficient μeff and reduced transport coefficient μt=μa+μs , respectively, in the diffusion approximation.10 Finally, the Pl are Legendre polynomials of order l that describe the angular shape of each N ’th order mode of the radiance, while the terms in square brackets [ ] are the relative contribution of each mode. Equation 2 can be obtained from Eq. 3 if the latter is truncated to first order and the optical properties are set to the limits of the diffusion approximation.10

Note that unlike the diffusion approximation, scattering can no longer be described simply by μs . However, Hull and Foster19 have shown that a simple similarity relation can still be used to directly relate higher N ’th order expansions of the anisotropy factor to μs . This relationship is given by μs(N)=σ(N)μs with σ(N=0,1,2,3)=(0,1,1.85,2.6) . This simplification allows the scattering to be described by the more commonly employed μs while maintaining accurate predictions of the radiance for g as low as 0.7 .19


Optical Property Determination Using Relative Fluence or Radiance Measurements: an Analysis Using the P1 (Diffusion) Approximation

In this section, we focus on the P1 (diffusion) approximation to the RTE, and examine the ability of either fluence or radiance measurements to uniquely determine the absorption and scattering properties of turbid media.

The rationale for such an analysis is two-fold. First, while higher order approximations to the RTE (e.g., P3 approximation) provide a more accurate description of the fluence or radiance, in some cases, experimental uncertainty can limit the accuracy of the measurements to first order.10 This is particularly true of interstitial measurements where the spatial sampling and measurement numbers are sparse. By determining first-order conditions for unique optical property determination, we can devise inverse algorithms that are potentially less sensitive to experimental uncertainties. Second, first-order solutions provide simpler analytical expressions compared to higher order solutions that allow the physics of the problem to be clearly identified. Furthermore, because the higher order solutions to the RTE are an extension of the lower order solutions, strategies derived from the P1 approximation are also applicable to the P3 approximation.

We first examine the case of spatially resolved fluence measurements and then proceed with an analysis of angularly resolved radiance measurements. Starting with Eq. 2, integration of the radiance over all solid angles yields an expression for the fluence:

Eq. 4

Equation 4 can be rewritten: ln[ϕ(r)r]=μeffr+ln[Po(4πD)] , implying that a plot of ln[ϕ(r)r] versus r yields a straight line with the slope equal to the negative of the effective attenuation coefficient of the medium. Since the intercept contains the Po term, absolute fluence measurements are required to uniquely separate μa and μs . This often requires complex calibration methods to correctly relate the measured photovoltage from a photodetector to the absolute fluence.21, 22 However, absolute fluence measurements are not required to determine μeff , since the slope is a relative measurement. Ideally, two measurements at different spatial positions are required to obtain μeff . However, in practice, multiple measurements are preferable to improve signal to noise.

In the case of radiance, we are interested in characterizing the optical properties, μa and μs , of turbid media using the angular radiance distribution L(θ,r) at position r . Working within the P1 radiance formulation, we can derive an analogous radiance term to the μeff parameter used to describe fluence. To avoid the necessity of absolute radiance measurements, we normalize Eq. 2 to the radiance at an arbitrary detection angle θn . Plotting the radiance as a function of cos(θ) now yields a slope K , given by:

Eq. 5

and an intercept C that is independent of the optical properties and only a function of θn .10


Strategies for Improving Optical Property Determination Using Single or Multiple Sensors: an Analysis Using the P1 Approximation

While improving modeling accuracy may benefit radiance optical property retrieval, the gains from more accurate models may potentially be offset by uncertainties from common experimental errors. To compensate for these limitations, we have derived and investigated a novel strategy that utilizes both radiance and spatial attenuation information.

From Eq. 5 we observe that various combinations of μa and μs can result in the same K value. Chin, Whelan, and Vitkin10 have previously demonstrated using Eq. 5 that a family of similar solution sets (i.e., different combinations of μa and μs ) exist within the diffusion formulation that yield identical relative radiance distributions.

In their analysis, P1-derived similar solution sets to radiance distributions calculated using the P3 approximation resulted in similar but separated (unique) radiance curves. The authors concluded that the added information of the higher order modes potentially allows for the unique separation of μa and μs . However, in the presence of typical experimental noise, the same solution sets were virtually indistinguishable.10 As such, optical property recovery using relative radiance measurements at a single sensor position can suffer from significant uncertainty due to experimental noise.

Based on Eq. 5, it was also suggested10 that, using radiance measurements, improved recovery of μa and μs can be obtained either by: 1. fits to the radiance slopes K at two different sensors positions, or 2. knowledge of both K and μeff . Of the two methods, the combination of a measured point radiance curve at a known position and a known μeff was shown to be most robust in providing convergence to a unique set of optical properties.10

Either method can be implemented using a minimum of two radiance sensors at different spatial positions ( r1 and r2 ). Practically, large separations of the two sensor positions are preferred. This is because as r1r2 , both methods converge to the information content of a single sensor. The ideal separation distance will be dependent on the optical properties of the medium of interest. However, while the radiance slope changes linearly as a function of sensor positions [Eq. 5], the fluence changes exponentially [Eq. 4]. As such, μeff can be determined with greater sensitivity at shorter r1 and r2 separations. This implies that while relatively large sensor separations are required for method 1, method 2 is feasible over a larger range of source-sensor separations and overall is preferable for interstitial applications.

In contrast to the single-sensor strategies, a particularly appealing property of the μeff radiance strategy is that it is less dependent on the accuracy of the particular forward model employed. In addition, the strategy is based on measurable first-order radiance and attenuation parameters. This is an important distinction, as first-order measurements are likely to be: 1. robust to experimental errors, and 2. less dependent on higher order models that, in some cases, can be computationally intensive.23 Such characteristics are particularly important for online interstitial applications where experimental errors can be significant and rapid calculations are preferable.


Monte Carlo Simulations of Interstitial Point Radiance at a Single Wavelength

A previously described Monte Carlo (MC) simulation for a point source in an infinite medium was used to generate synthetic radiance data.24 The Monte Carlo method propagates individual photons based on the optical properties of the simulated medium to generate a statistical representation of the resulting light distribution and is commonly employed as a gold standard of the RTE solutions. However, the accuracy of the method can be compromised by random noise when an inadequate number of photons are sampled. In this work, the simulation was divided into spatial bins with spherical symmetry. At locations far from the source and for highly attenuating media, statistical noise limits the accuracy of the MC data. To provide a measure of the statistical noise, we defined the percent standard deviation as σ=100CC , with C being the total number of interactions at a given bin. For the radiance data used in our analysis, σ was typically less than 0.3%, with a worst case σ of 1% at the 0- and 180-deg detection angles. Depending on the optical properties investigated, as few as 105 and as many as 109 launched photons were employed for a given simulation.

An anisotropy factor g of 0.9 was fixed for all simulations, typical of tissue, employing the commonly used H-G scattering phase function. All simulations employed detection angles with 1.8-deg angular resolution and spatial bins with 0.1-cm spatial resolution. The number of absorbed photons N(r,θ) were recorded for each spatial position r and each detection angle θ=acos(μ̂x,y,zn̂) . Here, μ̂x,y,z is the vector describing the direction of photon propagation, while n̂ is the spherical normal vector relative to the source.

The radiance as a function of detection angle L(r,θ) was then calculated according to the equation:

Eq. 6

where No is the total number of simulated photon packets and the two bracketed terms in the denominator are the voxel volume and solid angle volume, respectively. In Eq. 6, the division by μa converts the absorbed power density to the radiance. To check the accuracy of the simulations, we compared absolute radiance curves with MC results reported by other groups18, 25 and found good quantitative agreement. To further reduce noise, all data was smoothed using a three-point median filter from the Image Processing Toolbox of Matlab (Mathworks, Natick, Massachusetts).


Data Fitting

The goal of the fitting algorithms is to optimize the vector of parameters u=u(μa,μs,r,θ) , such that the forward generated radiance Lforward(u) matches the MC generated radiance L(r,θi)MC . Here, r and θ are assumed known and are thus fixed during fitting. The optical properties μa and μs are left as free parameters. A Levenberg-Marquardt algorithm from the Matlab 6.12 function fmincon (Mathworks) was used to optimize the values of μa and μs that minimized the chi-squared function defined as:

Eq. 7

In all cases, the P3 approximation was chosen as a forward model, as it was shown to provide an accurate description of the radiance over a wide range of optical properties and sensor positions typical of interstitial applications.10 However, we stress that the use of the P3 approximation does not obviate the conclusions from the P1 analysis performed in Sec. 2.3, as the P3 approximation is an extension of the P1 model.

While initial guesses for μa and μs close to the true solution aid in expediting convergence, to avoid bias, all initial guesses for μa and μs were purposely chosen to be unphysical and set to negative values. The optical properties were constrained such that 0.0001<μa<20cm1 and 0.001<μs<100cm1 . Where it was assumed that knowledge of the effective attenuation coefficient was available, an additional constraint was applied such that μa and μs were restricted to obey the equation: μeff=[3μa(μs+μa)]12 . Here, the known μeff was calculated using the MC forward properties.


Monte Carlo Study of Optical Property Recovery Using the P3 Approximation

To test the performance of the P3 approximation described before (Sec. 2.3), we generated a series of synthetic radiance curves using the MC simulation described in Sec. 2.6, and fit the data using the inversion algorithm described in Sec. 2.7. In this analysis, μs was held constant at 10cm1 while μa ranged from 0.001 to 10cm1 , such that the transport albedo a=μsμs+μa ranged from 0.5 to 0.999. The chosen optical property sets roughly span the range of expected transport albedos of tissue that may be relevant in interstitial therapies.12 Results are presented for single sensor fits at 0.5 and 1cm , as well as two sensor fits at the same positions. These positions are typical sensor locations for interstitial laser treatments. In the two sensor fit, only the shape of the radiance curves are considered and μeff is not employed as a constraint.

Figures 2a and 2b plot the recovered absorption coefficient and reduced scattering coefficient as a function of the actual absorption coefficient employed to generate the MC data. Several interesting observations can be seen from Fig. 2. Starting with the single sensor fits, we observe that the apparent lower limit for accurate recovery of absorption is approximately 0.005cm1 . This is predicted by Eq. 5, whereby radiance measurements lack the needed information content to separate optical properties as conditions approach the diffusion approximation.10

Fig. 2

(a) and (c) Recovered values of μa and (b) and (d) μs obtained from LP3 fits to Monte Carlo datasets; expected properties (solid line). (a) and (b) show unconstrained single and two sensor fits. (c) and (d) show single sensor fits with μeff constraint.


Beyond this lower limit and up to μa , values as high as 3cm1 , fits at both sensor positions returned absorption coefficients to within 15% of the expected value (worst case 25% ). Fits in the 0.05 to 1.75cm1 absorption range were 10% more accurate with 0.5-cm sensor data than 1.0-cm data, whereas fits in the 1.75 to 3cm1 range were 15% less accurate. This indicates that accurate μa recovery is dependent not only on the optical properties of the medium but also on the location of the radiance sensor.10 Interestingly, while the absorption coefficient demonstrated significant variations in accuracy, the reduced scattering coefficient was typically retrieved to within 10% (worst case 16%) for data at both sensor positions up to a μa of 3cm1 .

For μa> 3cm1 , errors of greater than 30% were obtained for both μa and μs . However, even at high μa , only small differences were observed between the P3 and MC generated radiance curves (data not shown). This implies that for single sensor radiance data, very small errors in predicted radiance potentially result in very large errors in recovered optical properties during fitting. Greater stability is expected by employing a μeff constraint or extending the data set to two sensor locations.

Figures 2a and 2b show the recovered optical properties when the radiance slope at two different sensor positions is fit simultaneously. Improvements were obtained due to the additional data with both μa and μs recovered to within 10% between the μa range of 0.05 to 3cm1 . However, while improvements in recovered optical properties were observed, two sensor fits still exhibited lower (0.005cm1) and upper (3cm1) limits of accuracy similar to single sensor fits. This is not surprising, since both single and two sensor fits are based on the radiance slope (a function of θ ) and should suffer similar limitations in uniqueness and model inaccuracy. To further improve the stability of the radiance fits, we can incorporate information on the spatial change in the radiance ( μeff constraint). Such information is independent of the radiance slope and should therefore not suffer from the same limitations.

We repeated the prior analysis constraining the radiance fit at single sensor positions of 0.5- or 1-cm positions with a known μeff . The improvements in optical property recovery are substantial. This is demonstrated in Figs. 2c and 2d, where it can be seen that both μa and μs are now recovered to within 1.4% . Even more striking is that the lower and upper limits that existed for the radiance-slope-based technique are no longer apparent. The results indicate that radiance fits performed using a μeff constraint provide accurate recovery of optical properties over a wider range of sensor positions and reduced albedos typical of interstitial applications, and that employing a μeff constraint method is preferable compared to fitting two radiance slopes. Practically, μeff can be determined via literature values or calculated from experimental measurements as described in Sec. 5.4.


Effects of Experimental Deviations from Theoretical Behavior

The synthetic MC data represent an idealized environment where the fundamental limits of the P3 approximation as an inversion model could be tested. However, while fits to MC generated radiance curves resulted in accurate optical property retrieval, initial experimental validation in tissue-simulating phantoms using Eq. 3 demonstrated a systematic overestimation of μa and an underestimation of μs (data not shown). As such, we hypothesized that these errors resulted from potential discrepancies in the P3 model (from the idealized situation) that were encountered during experiments.

In this section, we briefly discuss the differences between the idealized and experimental cases, and analyze potential fitting errors that may arise if such discrepancies are unaccounted for. Furthermore, to focus on the effects of such experimental deviations, we generate radiance curves (with various systematic errors to reflect real experimental situations) using the P3 approximation (instead of MC), and employ the same P3 model to recover optical properties. Such an analysis ensures that the resulting errors in recovered optical properties are not due to the assumptions inherent in the forward model, but are due solely to any systematic error introduced. For simplicity, we focus on the μeff method in this section, and compare our results with the single sensor unconstrained approach.


Probe Dimensions

Experimentally, radiance information is acquired by rotation of a sensor at different collection angles about a single radial position. However, due to finite probe size, a radial shift in position occurs (relative to the 0-deg detection angle) that reaches a maximum discrepancy at 180deg . To investigate the potential errors due to probe dimensions, forward radiance distributions were generated at two sensor positions (0.5 and 1.0cm ) using the P3 approximation [Eq. 3] for three combinations of optical properties, typical of biological tissues. To simulate high and low transport albedos situations, a constant μs=10cm1 and variable μa of 0.01, 0.1, 1cm1 were selected. The effects of finite probe dimensions were simulated by assuming that the probe size results in a radial shift Δr that is a function of the probe diameter rp and detection angle θ :

Eq. 8

Figure 3 illustrates the geometry under consideration.

Fig. 3

Geometry for a radial shift Δr due to rotation of a finite-sized radiance sensor of dimensions rp . (a) At a detection angle of θ=0deg , the collection face of the sensor is located at position r . (b) At a detection angle of θ=180deg , the collection face of the sensor shifts to position r+Δr .


Probe sizes ranging up to 1mmdiam were evaluated, representative of the range of right commercially available angle-etched fibers (i.e., Poly-micro Technologies, Phoenix, Arizona).

Figure 4 shows the resulting radiance distribution when Eq. 8 is incorporated into the analytical P3 calculations. Significant differences in radiance are observed between the idealized and finite size cases, where a progressively steeper slope is predicted as the probe diameter increases. As an example, we focus on the μa=0.01cm1 dataset. It is clear that the radiance curve resulting from a 0-mm probe diameter (top dashed line) decreases much more slowly as a function of detection angle compared to the curve resulting from a 1-mm probe diameter (bottom dashed line). The same trend is observed for the other two absorption values. To examine the errors that may result from neglecting probe size on the recovered optical properties, the idealized P3 expression (with rp=0mm ) was used to fit the radiance data in Fig. 4. All forward generated data were normalized to the radiance value at 0deg . The resulting error ε was then calculated as:

Eq. 9

where μfit and μtrue are the fitted and true absorption or scattering properties, respectively.

Fig. 4

Change in radiance distribution at r=1cm due to finite probe dimensions for three different optical property sets. For each family of curves, the top dashed line is the idealized case (0-mm probe diameter), while the bottom dashed line was generated using a probe diameter of 1mm . Ten different curves are shown for each optical property set with probe radii between 0 and 0.055cm .


Figure 5 shows the relative error in the recovered absorption coefficient and reduced scattering coefficient, with and without μeff constraints, as a function of increasing probe size for a sensor positioned at 1cm . The entire radiance dataset between 0 and 180deg was employed in the fits. Some general observations can be made from Fig. 5. First, errors in absorption are slightly greater than scattering. Second, finite probe size results in an overestimation of μa and an underestimation of μs . This result is easily explained by Fig. 4, where the effect of the finite probe size effectively steepens the radiance curve. This results in a radiance curve that is representative of an artificially smaller transport albedo.

Fig. 5

Error in recovered (a) and (c) μa and (b) and (d) μs if finite probe size is unaccounted for during fitting. Fits are for radiance data at 1cm (a) and (b) without and (c) and (d) with a μeff constraint. A μs of 10cm1 was employed for all cases. Similar trends were observed at 0.5cm .


The effects of this change are significant. For a typical probe radius of 0.3mm (600μmdiam) , errors of greater than 150% can result for the unconstrained case, compared to a 20% error for the constrained case. For μa<0.1cm1 , μa is highly sensitive to probe radius, whereas μs errors are under 10% . Interestingly, for μa=1cm1 , the opposite trend is observed with μa errors decreasing to below 25% and μs errors increasing to 40% .

It should be noted that for the μa=1cm1 case, a μeff constraint does not significantly minimize the error due to finite probe size for either μs or μa . We suspect that for this particular optical property set, sufficient information in the radiance curve exists such that the optical properties can be recovered with similar accuracy as when a μeff constraint is employed.

In the lower absorption cases ( μa=0.01 , 0.1cm1 ), adding a μeff constraint significantly improved recovery of μa but provided little improvement in recovered μs . An examination of Eq. 5 provides some insight into these trends. From Eq. 5, we see that as μa0 , K1μs[1r] . This means that while μs can be directly recovered from K , μa can attain virtually any value and is highly sensitive to small errors. In this case, a μeff constraint provides significant improvements in the convergence of μa .

Therefore, if the improvement in recovered μa is significantly greater than the error due to finite probe size, an overall reduction in μa error is observed. By contrast, a similar improvement in μs would not be observed, since a μeff constraint provides little benefit in μs recovery when the absorption is low. Here, finite probe size is the dominant source of error.

The results demonstrate that probe dimensions must be accounted for during inversion for accurate optical property quantification. As such, in all experimental studies to follow, Eq. 8 was incorporated into the fitting algorithm.


Numerical aperture

Experimentally, radiance information is integrated over the numerical aperture (NA) of the sensor. Dickey 13 assumed that the radiance is well represented by the central angle of the radiance probe, and demonstrated good agreement with experimental data. However, a finite NA may produce a loss of angular resolution, and hence cannot be ignored if the radiance changes rapidly with detection angle. The effects of increasing NA were simulated by integrating Eq. 3 over a typical range of fiber NAs.

Figure 6 shows the resulting change in the shape of radiance distribution due to various NA ranging from 0 (idealized) to 0.5 (corresponding to acceptance angles between 0 and 30deg ). As shown in Fig. 6, the effect of finite NA is a slight flattening of the radiance curve. However, as verified by the previous experimental work of Dickey, 13 the overall effects of finite NA are generally negligible. This observation is also demonstrated in Fig. 7 , which shows the relative quantification error for the absorption coefficient and reduced scattering coefficient as a function of increasing NA for unconstrained and μeff constrained fits. For a NA less than 0.25 (acceptance angle 15deg ), errors in recovered μa were less than 20% and 2% for unconstrained and μeff constrained fits, respectively.

Fig. 6

Change in radiance distribution at 1cm due to different sized numerical apertures for three different optical property sets. For each optical property set, the top dashed line is the idealized case (NA=0deg) , while the bottom dashed line was generated using an NA of 0.5. Six different curves are shown for each optical property set with NA ranging between 0 and 0.5.


Fig. 7

Error in recovered (a) and (c) μa and (b) and (d) μs if numerical aperture is ignored during fitting. Fits are for radiance data at 1cm (a) and (b) without and (c) and (d) with a μeff constraint.


While the optical properties are both overestimated for the unconstrained case, an underestimation in μa and an overestimation in μs was observed for the constrained method. For the unconstrained case, the overestimation of both μa and μs are likely explained by optical similarity. This is because similar solution sets for radiance tend to maintain the same albedo.10 This is not the case for μeff , where a decrease in μa must be balanced by an increase in μs to maintain the same attenuation.

Note that for simplicity, we assume the index of refraction of the surrounding medium equal to 1, although the results can be scaled appropriately to the index of refraction of the medium of interest. Indeed, given average tissue properties with index of refraction 1.4 , a general lowering of the acceptance angle is expected to occur, because of tissue/probe refractive index similarity, leading to even further reductions in error. Based on these modest effects of NA, we have ignored them in subsequent analyses.


Experimental Evaluation in Phantoms


Experimental Setup

Figure 8 shows a schematic of the experimental setup.

Fig. 8

Experimental setup.


A 650-μm-diam isotropic source (Resonance Optics) coupled to a 690-nm laser was used to illuminate the phantom interstitially at 100- to 300-mW power levels. The source fiber was also attached to a precision translation stage (Thorlabs) that allowed for positioning of 0.5 to 1.5cm from the radiance probe with 0.01-cm resolution. The radiance probe was constructed according to the method of Dickey 13 by attaching a 700-μm right-angle prism (Melles Griot) to a cleaved plane cut fiber. The acceptance angle of the radiance sensor was measured to be 10deg (NA0.17) in air.

The collected radiance data was read using a photomultiplier tube (PMT) and secured by a fiber chuck into the center of a calibrated rotating optical stage (Melles Griot). The rotating optical stage was controlled by Labview software to rotate with angular increments of 2deg . At each angular position, an average of 300 photovoltage acquisitions was read into a PC and stored into an output file containing a complete rotation between 0 (facing the source) and 180deg (facing away from the source). A complete set of 91 readings was acquired in approximately 90s . Sensors were placed within a hollow bore needle and inserted into an acrylic template for positioning within the phantom. To avoid precession during rotation, a thin jacket with the exact dimensions of the needle was wrapped around the probe.

All measurements were made at the center of a 10×11×11-cm acrylic box. This allowed for the radiance probe to be positioned 4 to 5cm from the box boundaries, thereby mimicking an infinite medium geometry. The box was painted black to minimize reflections from the box interface and to limit external light contamination. Point (angular) radiance data were measured at 0.5 and 1cm from the source, while spatial radiance data (at a fixed collection angle θ of 90deg ) were collected by translating the source at distances between 0.5 to 1.5cm with 0.05-cm step size from the sensors. All experiments were performed in a darkened room. As such, the detected background signal was negligible in all cases (i.e., no background subtraction was required).


Phantom Material and Characterization

Turbid tissue simulating liquid phantoms were prepared using a combination of 20% Intralipid (Kabi-Pharmacia, Germany) for scattering, and added quantities of Naphthol Green dye (Sigma Aldrigde) for absorption. In our experiments, the Intralipid concentration was kept approximately constant, and small aliquots of dye with negligible volume (compared to the overall phantom volume) were added, such that μa ranged from 0 to 5cm1 . Note that unlike the MC simulation, we did not perform measurements at lower transport albedos (0.5) due to limitations in detected signal.

The absorption coefficient of the Naphthol Green stock solution was measured using a spectrophotometer. However, while estimates for the scattering properties of Intralipid exist in the literature,26 batch differences can result in a wide range of μs .6 Therefore, to establish the baseline scattering properties of the Intralipid solution, we utilized the poisoned moderator (PM) technique.18, 27

With the PM technique, quantities of dye with known absorption (measured using a spectrophotometer) are added to the phantom and the resulting μeff measured for each quantity of added absorption μa,dye . In our experiments, μeff was obtained by performing spatial scans of the radiance [L(θ=90deg,r)] photovoltage V between 0.5 and 1.5cm in 0.05-cm increments, and determining the slope of a plot of ln(Vr) versus r . According to the PM technique, plotting μeff2 as a function of μa,dye gives:

Eq. 10

where the slope is a function of the reduced scattering coefficient μs,IL of the Intralipid solution, and the intercept is a combination of the Intralipid scattering μs,IL and absorption μa,IL . The range of added μa,dye utilized for the PM calculations were between 0 and 0.06cm1 , such that the estimated ratio of μa,dyeμs,IL0.1 . Martelli 18 previously demonstrated that by using this quoted ratio, the optical properties of the target medium can be determined to within 2% using the PM technique. Linear regression of Eq. 10 was performed using the Excel spreadsheet package. For a 5% dilution of 20% Intralipid, the PM technique yielded a μs,IL of 9.43cm1 and a μa,IL of 0.0049cm1 at 690nm .



Figures 9a through 9d demonstrate the recovered values of μa and μs , respectively. In Figs. 9a and 9c (μa) , the dashed line represents the expected absorption measured using a spectrophotometer in the presence of no Intralipid. All data are background absorption subtracted using the μa determined from the PM technique. In Figs. 9b and 9d (μs) , the dashed line represents the expected μs determined using the PM method. Note that finite probe size was accounted for in the fits by incorporating Eq. 8.

Fig. 9

Recovered values of (a) and (c) μa and (b) and (d) μs obtained from LP3 fits to experimental measurements in tissue simulating Intralipid phantoms. Expected properties (solid line); (a) and (b) unconstrained single and two sensor fits; (c) and (d) single sensor fits with μeff constraint.


Trends observed experimentally were similar to those seen in the MC analysis. With the exception of the lower limit of μa=0.01cm1 , unconstrained single sensor fits recovered both optical properties to within 10% at the 0.5-cm position and to within 20% at the 1-cm position. Above 1cm1 , fits for both sensor positions displayed a significant breakdown in recovered μs , although the 1-cm sensor recovered μa accurately to within 30% , up to 4cm1 . The two sensor fitting algorithm (no μeff constraint) only showed minor improvements in accuracy compared to the single sensor fits. However, fits performed using a μeff constraint typically recovered both optical properties to within 5% (worst case 12%) over the entire optical property range investigated.

In summary, in the experimental setting, fitting algorithms based solely on the shape of the radiance curve (single or two sensor) can recover optical properties to within 20% for a reduced albedo range of 0.9 to 0.999 at sensor positions between 0.5 and 1cm . Furthermore, with the additional information regarding the spatial distribution of light (i.e., μeff constraint), radiance-based fitting algorithms are able to experimentally determine optical properties to within 12% over a larger albedo range of 0.67 to 0.999.


Comparison with Other Interstitial Techniques

It is helpful to compare the μeff radiance technique to other potential interstitial methods to gain some insight into its practical clinical implementation. Current fluence-based techniques require absolute calibration, source modulation, or multiple wavelengths to accurately extract optical properties. The proposed μeff radiance strategy requires only relative steady-state information and is, in terms of source and detection equipment, technologically simpler than absolute or frequency-domain methods. The fact that only single wavelength steady-state measurements are employed is a particular advantage since the same laser may be employed to both characterize and treat the targeted tissue. Although the radiance method requires mechanical rotation, fluence techniques generally require mechanical translation and as such offer similar technical difficulties and practical (clinical) limitations. In addition, rotational requirements may be further minimized by bundling multiple radiance sensors.

How can the μeff radiance technique be implemented practically? If the μeff of the targeted tissue is known a priori (i.e., via literature or clinical database), the approach is straightforward. In such cases, the accuracy of the recovered optical properties will be limited by the uncertainty of the μeff estimate.

Alternatively, μeff may be determined directly. This approach was utilized during experimental verification of the radiance technique in Secs. 5.1, 5.2, 5.3. It should be noted that most fluence-based approaches require experimental determination of μeff . Hence, the experimental evidence presented in this work indicates that our radiance method is, at the very least, comparable to such methods.

Practically, the assessment of μeff will require measurements of the radiance at a minimum of two sensor positions. In this case, the μeff constraint may be applied individually to constrain the radiance fits at each sensor position (for local characterization of optical properties) or to simultaneously constrain fitting at both sensor positions (for global determination of optical properties). Such an approach allows for an assessment of local variations in optical properties that may provide more patient-specific planning of interstitial treatments. In cases where one of the optical properties is expected to remain constant, the constant property (scattering or absorption) may be employed as a further constraint to allow for accurate differential determination of the changing property at individual sensor locations. As an example, such a case may be relevant to PDT applications where changes in the chromophore concentration occur both globally and locally, while the scattering is expected to remain virtually constant.

A two-sensor μeff radiance technique appears to be comparable to the two-sensor frequency-domain methodology of Xu and Patterson.12 Similar to the work of Xu and Patterson,12 the proposed technique is accurate over a wide range of properties. However, note that while the accuracy of our phantom studies was typically two times more accurate than the technique reported by Xu and Patterson,12 our μeff was calculated from experimental measurements at multiple distances (0.5 to 1.5cm ) and not from two positions. It is likely that somewhat larger errors will result when only two sensors are used to determine μeff .

Despite the current clinical limitations of single sensor radiance measurements, it is important to stress that we were able to experimentally recover both optical properties to within 20% over a very reasonable range of conditions typical of interstitial therapies. These promising results indicate the significant potential of the point radiance methodology. In particular, even in cases when errors in optical property recovery were greater than 20%, reasonable fits were often obtained with the experimental measurements. This suggests not a serious breakdown in the P3 methodology, but more a lack of inherent uniqueness for fitting algorithms based solely on single wavelength point steady-state radiance measurements. This is supported by the analysis of Eq. 5 where, to the first order, multiple pairs of absorption and scattering parameters can lead to essentially the same radiance, LP3(θ,r) curve. For the case of steady-state radiance measurements at a single wavelength, the additional information of a μeff constraint greatly improves inversion performance. However, if point radiance measurements are performed spectrally or are combined with other modalities such as frequency domain, similar improvements in uniqueness and experimental robustness may yet result.


Discussion and Conclusions

We investigate the role of relative radiance measurements for quantifying the optical properties of turbid media in interstitial geometry. Using Monte Carlo simulations and phantom experiments we show that at sensor positions 0.5 to 1cm , the P3 approximation is adequate for retrieving optical properties with albedo ranges of 0.9 to 0.99 and 0.6 to 0.999 using unconstrained and μeff constrained fits to radiance measurements, respectively. In addition, we perform a systematic analysis of potential modeling discrepancies that arise between idealized and experimental radiance measurements. Our analysis indicates that while errors due to a finite numerical aperture are small, significant errors may result if the finite dimensions of a radiance probe are not accounted for when fitting radiance data. A simple analytical expression has been introduced that properly accounts for errors related to probe dimensions. This work extends the initial study of Dickey, 13 who demonstrated optical property recovery for a much smaller subset of experimental conditions (single optical property set and sensor position).

Of fundamental interest is the source of the uniqueness of relative radiance information. Time-domain methods rely on differences in optical path length and fit a histogram of the photons arriving at a fixed source-detector distance to extract optical properties.28 Because radiance information provides similar quantification capabilities in the steady-state regime, it is possible that the radiance slope K is actually a representation of differences in average path length traveled between the measured detection angles. Extension of the P3 approximation or Monte Carlo simulations to include both time and radiance information may help clarify this issue. In addition, recent reports using the δ -P1 approximation20 by Hayakawa 29 and You, Hayakawa, and Venugopalan30 have demonstrated that the near- to far-field transition in phase and amplitude of spatially resolved fluence scans may be suggestive of the absolute value of the anisotropy factor g . Similarly, Menon, Su, and Grobe31 recently reported on a spatially resolved fluence method that employs a directional source to recover g . Hence, an interesting avenue of investigation may be to explore whether interstitial radiance data contain sufficient sensitivity to accurately recover the anisotropy factor.

It is envisioned that the described radiance technique may prove ideal for online monitoring and treatment planning of laser interstitial thermal therapies (LITT). In our laboratory, we are currently developing an optical toolkit for online treatment planning and monitoring of thermal damage. In particular, it has been demonstrated that real-time coagulation-induced changes in radiance information are particularly sensitive to coagulation events compared to fluence sensors.24 Using a two sensor radiance strategy, the native optical properties of the targeted medium could be characterized to optimize input power and heating times. The same radiance sensors could be employed to monitor the evolution of the coagulation boundary during heating. Finally, coagulated optical properties may be evaluated post-treatment for the planning of future treatments. While the current work focuses on point source geometry, extension of the P3 model to cylindrical geometries will be necessary for proper online clinical implementation of radiance optical property determination. Furthermore, since most interstitial sensors are inserted inside the target tissue via catheters, additional investigation will be needed to address the potential perturbation effects of the additional interfaces.


Financial support for this work was provided by the National Cancer Institute of Canada (with funds from the Canadian Cancer Society) and the Natural Sciences and Engineering Research Council of Canada. The authors would like to thank Drs. Brian Wilson and Dwayne Dickey for helpful discussions and insight.



S. Puccini, N. K. Bar, M. Bublat, T. Kahn, and H. Busse, “Simulations of thermal tissue coagulation and their value for the planning and monitoring of laser-induced interstitial thermotherapy (LITT),” Magn. Reson. Med., 49 351 –362 (2003). 0740-3194 Google Scholar


M. D. Altshuler, T. C. Zhu, J. Li, and S. M. Hahn, “Optimized interstitial PDT prostate treatment planning with the Cimmino feasibility algorithm,” Med. Phys., 32 3524 –3536 (2005). 0094-2405 Google Scholar


R. A. Weersink, A. Bogaards, M. Gertner, S. R. Davidson, K. Zhang, G. Netchev, J. Trachtenberg, and B. C. Wilson, “Techniques for delivery and monitoring of TOOKAD (WSTO9)-mediated photodynamic therapy of the prostate: clinical experience and practicalities,” J. Photochem. Photobiol., B, 79 211 –212 (2005). 1011-1344 Google Scholar


C. Holmer, K. S. Lehmann, J. Risk, A. Roggan, C. T. Germer, C. Reissfelder, C. Isbert, H. J. Buhr, and J. P. Ritz, “Colorectal tumors and hepatic metastases differ in their optical properties—relevance for dosimetry in laser-induced interstitial thermal therapy,” (2006) Google Scholar


L. Lilge, N. Pomerleau-Dalcourt, A. Douplik, S. H. Selman, R. W. Keck, M. Szkudlarek, M. Pestka, and J. Jankun, “Transperineal in vivo fluence-rate dosimetry in the canine prostate during SnET2-mediated PDT,” Phys. Med. Biol., 49 3209 –3225 (2004). 0031-9155 Google Scholar


A. Dimofte, J. C. Finlay, and T. C. Zhu, “A method for determination of the absorption and scattering properties interstitially in turbid media,” Phys. Med. Biol., 50 2291 –2311 (2005). 0031-9155 Google Scholar


T. J. Farrell, M. S. Patterson, and B. C. Wilson, “A diffusion theory model of spatially resolved, steady-state diffuse reflectance for the noninvasive determination of tissue optical properties in vivo,” Med. Phys., 19 879 –888 (1992). 0094-2405 Google Scholar


E. L. Hull, M. G. Nichols, and T. H. Foster, “Quantitative broadband near-infrared spectroscopy of tissue-simulating phantoms containing erythrocytes,” Phys. Med. Biol., 43 3381 –3404 (1998). 0031-9155 Google Scholar


D. R. Wyman, M. S. Patterson, and B. C. Wilson, “Similarity relations for the interaction parameters in radiation transport,” Appl. Opt., 28 5243 –5249 (1989). 0003-6935 Google Scholar


L. C. Chin, W. M. Whelan, and I. A. Vitkin, “Information content of point radiance measurements in turbid media: implications for interstitial optical property quantification,” Appl. Opt., 45 101 –114 (2006). 0003-6935 Google Scholar


T. C. Zhu, J. C. Finlay, and S. M. Hahn, “Determination of the distribution of light, optical properties, drug concentration, and tissue oxygenation in-vivo in human prostate during motexafin lutetium-mediated photodynamic therapy,” J. Photochem. Photobiol., B, 79 231 –241 (2005). 1011-1344 Google Scholar


H. Xu and M. S. Patterson, “Determination of the optical properties of tissue-simulating phantoms from interstitial frequency domain measurements of relative fluence and phase difference,” Opt. Express, 14 6485 –6501 (2006). 1094-4087 Google Scholar


D. J. Dickey, R. B. Moore, D. C. Rayner, and J. Tulip, “Light dosimetry using the P3 approximation,” Phys. Med. Biol., 46 2359 –2370 (2001). 0031-9155 Google Scholar


S. Chandrasekhar, Radiative Transfer, Univ. Clarendon Press, Oxford, UK (1950). Google Scholar


S. L. Jacques, “Origins of optical properties in the UVA, visible, and NIR regions,” OSA Trends Opt. Photonics Ser., 2 364 –371 (1996). 1094-5695 Google Scholar


L. G. Henyey and J. L. Greenstein, “Diffuse radiation in the galaxy,” Astrophys., 93 70 –83 (1941). 0571-7256 Google Scholar


G. W. Faris, “P(N) approximation for frequency-domain measurements in scattering media,” Appl. Opt., 44 2058 –2071 (2005). 0003-6935 Google Scholar


F. Martelli, M. Bassani, L. Alianelli, L. Zangheri, and L. Zaccanti, “Accuracy of the diffusion equation to describe photon migration through an infinite medium: numerical and experimental investigation,” Phys. Med. Biol., 45 1359 –1373 (2000). 0031-9155 Google Scholar


E. Hull and T. H. Foster, “Steady-state reflectance spectroscopy in the P3 approximation,” J. Opt. Soc. Am. A, 18 584 –599 (2001). 0740-3232 Google Scholar


V. Venugopalan, J. S. You, and B. J. Tromberg, “Radiative transport in the diffusion approximation: an extension for highly absorbing media and small source-detector separations,” Phys. Rev. E, 58 2395 –2407 (1998). 1063-651X Google Scholar


J. P. Marijnissen and W. M. Star, “Calibration of isotropic light dosimetry probes based on scattering bulbs in clear media,” Phys. Med. Biol., 41 1191 –1208 (1996). 0031-9155 Google Scholar


J. P. Marijnissen and W. M. Star, “Performance of isotropic light dosimetry probes based on scattering bulbs in turbid media,” Phys. Med. Biol., 47 2049 –2058 (2002). 0031-9155 Google Scholar


K. Rinzema, L. H. P. Murrer, and W. M. Star, “Direct experimental verification of light transport in an optical phantom,” J. Opt. Soc. Am. A, 15 2078 (1998). 0740-3232 Google Scholar


L. C. L. Chin, B. C. Wilson, W. M. Whelan, and I. A. Vitkin, “Radiance-based monitoring of the coagulation boundary during laser interstitial thermal therapy,” Opt. Lett., 29 959 –961 (2004). 0146-9592 Google Scholar


O. Bajaras, A. M. Ballangrud, G. G. Miller, R. B. Moore, and J. Tulip, “Monte Carlo modeling of angular radiance in tissue phantoms and human prostate: PDT light dosimetry,” Phys. Med. Biol., 42 1675 –1687 (1997). 0031-9155 Google Scholar


H. J. Van Staveren, C. J. M. Moes, J. Van Marle, S. A. Prahl, and M. J. C. van Gemert, “Light scattering in Intralipid-10% in the wavelength range of 4001100nm,” Appl. Opt., 30 4507 –4514 (1997). 0003-6935 Google Scholar


M. N. Iizuka, M. D. Sherar, and I. A. Vitkin, “Optical phantom materials for near infrared laser photocoagulation studies,” Lasers Surg. Med., 25 159 –169 (1999).<159::AID-LSM10>3.0.CO;2-V 0196-8092 Google Scholar


A. Liebert, H. Wabnitz, D. Grosenick, M. Moller, R. Macdonald, and H. Rinneberg, “Evaluation of optical properties of highly scattering media by moments of distributions of time of flight of photons,” Appl. Opt., 42 5795 –5792 (2003). 0003-6935 Google Scholar


C. K. Hayakawa, B. Y. Hill, J. S. Joon, F. Bevilacqua, J. Spanier, and V. Venugopalan, “Use of the δ-P1 approximation for recovery of optical absorption, scattering, and asymmetry coefficients in turbid media,” Appl. Opt., 43 4677 –4684 (2004). 0003-6935 Google Scholar


J. S. You, C. K. Hayakawa, and V. Venugopalan, “Frequency domain photon migration in the delta- P1 approximation: analysis of ballistic, transport, and diffuse regimes,” Phys. Rev. E, 72 021903 (2005). 1063-651X Google Scholar


S. Menon, Q. Su, and R. Grobe, “Determination of g and μ using multiply scattered light in turbid media,” Phys. Rev. Lett., 94 153904 (2005). 0031-9007 Google Scholar
©(2007) Society of Photo-Optical Instrumentation Engineers (SPIE)
Lee C. L. Chin, Arthur E. Worthington, William M. Whelan, and I. Alex Vitkin "Determination of the optical properties of turbid media using relative interstitial radiance measurements: Monte Carlo study, experimental validation, and sensitivity analysis," Journal of Biomedical Optics 12(6), 064027 (1 November 2007).
Published: 1 November 2007

Back to Top