2.1.

Radiative Transfer Equation

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

In the RTE, the optical properties of interest are the absorption coefficient ${\mu }_a$ and scattering coefficient ${\mu }_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(\widehat{s},{\widehat{s}}^{\prime })$ , which is the probability of a photon traveling in direction ${\widehat{s}}^{\prime }$ scattering into direction $\widehat{s}$ . While precise information on $p(\widehat{s},{\widehat{s}}^{\prime })$ in tissues would be desirable, Jacques¹⁵ demonstrated that the Henyey-Greenstein (H-G) phase function¹⁶ 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.

2.2.

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 ${\mathrm{P}}_N$ approximation, expands $L(\stackrel{⃑}{r},\widehat{s})$ , $p(\widehat{s},{\widehat{s}}^{\prime })$ , and $S(\vec{r},\widehat{s})$ in spherical harmonics to the $N$ ’th order.¹⁷ The $N=1$ expansion expresses the radiance $L(\stackrel{⃑}{r},\widehat{s})$ as the linear combination of the isotropic fluence $\varphi \left(\stackrel{⃑}{r}\right)$ and a linearly anisotropic flux $\stackrel{⃑}{j}\left(\stackrel{⃑}{r}\right)$ ¹⁸:

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

Fig. 1

Geometry for the radiance $L$ , fluence $\varphi$ , positive flux $j^{+}$ , and negative flux $j^{-}$ at position $r$ relative to a point source in a scattering medium with optical properties ${\mu }_a$ and ${\mu }_s^{\prime }$ . Examples of the radiance are shown for detection angles of $\theta =0$ and $180\mathrm{deg}$ . The detection angle for the radiance traveling in direction $\widehat{s}$ is determined relative to the radial normal vector $\stackrel{⃑}{r}$ .

Here, $D=1∕3({\mu }_s^{\prime }+{\mu }_a)$ is the diffusion coefficient, $P_o$ is the source power, ${\mu }_{\mathrm{eff}}=\left[3{\mu }_a\right({\mu }_a+{\mu }_s^{\prime }){]}^{1∕2}$ is the effective attenuation coefficient, and $\theta$ is the angle between the radiance direction $\widehat{s}$ and the outgoing radial normal vector $\stackrel{⃑}{r}$ .

In Eq. 2 we have also introduced the reduced scattering coefficient ${\mu }_s^{\prime }=(1-g){\mu }_s$ , where $g={\int }_{4\pi }p(\widehat{s},{\widehat{s}}^{\prime })(\widehat{s}\cdot {\widehat{s}}^{\prime })d\omega$ is also known as the anisotropy factor. The reduced scattering coefficient allows the scattering coefficient ${\mu }_s$ in an anisotropically scattering medium to be replaced by an equivalent isotropic scattering parameter ${\mu }_s^{\prime }$ . As such, the parameter ${\mu }_s^{\prime }$ 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 ${\mu }_s^{\prime }$ and not ${\mu }_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 radiance¹⁰ in a forward sense.

2.3.

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.¹⁰ The choice of the P3 approximation was also based on the computational convenience of the higher order similarity relations developed by Hull and Foster.¹⁹ 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 $\delta$ -P1 approximation,²⁰ the modified spherical harmonics method,¹² or even the P5 approximation,¹⁷ are also potential candidates for inversion applications.

The P3 radiance solution for a point source in spherical geometry is¹⁹:

Note that unlike the diffusion approximation, scattering can no longer be described simply by ${\mu }_s^{\prime }$ . However, Hull and Foster¹⁹ have shown that a simple similarity relation can still be used to directly relate higher $N$ ’th order expansions of the anisotropy factor to ${\mu }_s^{\prime }$ . This relationship is given by ${{\mu }_s^{\prime }}^{\left(N\right)}=\sigma \left(N\right){\mu }_s^{\prime }$ with $\sigma (N=0,1,2,3)=(0,1,1.85,2.6)$ . This simplification allows the scattering to be described by the more commonly employed ${\mu }_s^{\prime }$ while maintaining accurate predictions of the radiance for $g$ as low as $\sim 0.7$ .¹⁹

2.4.

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.¹⁰ 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:

In the case of radiance, we are interested in characterizing the optical properties, ${\mu }_a$ and ${\mu }_s^{\prime }$ , of turbid media using the angular radiance distribution $L(\theta ,\stackrel{⃑}{r})$ at position $\stackrel{⃑}{r}$ . Working within the P1 radiance formulation, we can derive an analogous radiance term to the ${\mu }_{\mathrm{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 ${\theta }_n$ . Plotting the radiance as a function of $\cos \left(\theta \right)$ now yields a slope $K$ , given by:

2.5.

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 ${\mu }_a$ and ${\mu }_s^{\prime }$ can result in the same $K$ value. Chin, Whelan, and Vitkin¹⁰ have previously demonstrated using Eq. 5 that a family of similar solution sets (i.e., different combinations of ${\mu }_a$ and ${\mu }_s^{\prime }$ ) 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 ${\mu }_a$ and ${\mu }_s^{\prime }$ . However, in the presence of typical experimental noise, the same solution sets were virtually indistinguishable.¹⁰ 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 suggested¹⁰ that, using radiance measurements, improved recovery of ${\mu }_a$ and ${\mu }_s^{\prime }$ can be obtained either by: 1. fits to the radiance slopes $K$ at two different sensors positions, or 2. knowledge of both $K$ and ${\mu }_{\mathrm{eff}}$ . Of the two methods, the combination of a measured point radiance curve at a known position and a known ${\mu }_{\mathrm{eff}}$ was shown to be most robust in providing convergence to a unique set of optical properties.¹⁰

Either method can be implemented using a minimum of two radiance sensors at different spatial positions ( $r_1$ and $r_2$ ). Practically, large separations of the two sensor positions are preferred. This is because as $r_1\to r_2$ , 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, ${\mu }_{\mathrm{eff}}$ can be determined with greater sensitivity at shorter $r_1$ and $r_2$ 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 ${\mu }_{\mathrm{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.²³ Such characteristics are particularly important for online interstitial applications where experimental errors can be significant and rapid calculations are preferable.

2.6.

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.²⁴ 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 $\sigma =100\sqrt{C}∕C$ , with $C$ being the total number of interactions at a given bin. For the radiance data used in our analysis, $\sigma$ was typically less than 0.3%, with a worst case $\sigma$ of $\sim 1\%$ at the 0- and 180-deg detection angles. Depending on the optical properties investigated, as few as ${10}^5$ and as many as ${10}^9$ 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(\stackrel{⃑}{r},\theta )$ were recorded for each spatial position $\stackrel{⃑}{r}$ and each detection angle $\theta =a\cos ({\widehat{\mu }}_{x,y,z}\cdot \widehat{n})$ . Here, ${\widehat{\mu }}_{x,y,z}$ is the vector describing the direction of photon propagation, while $\widehat{n}$ is the spherical normal vector relative to the source.

The radiance as a function of detection angle $L(\stackrel{⃑}{r},\theta )$ was then calculated according to the equation:

2.7.

Data Fitting

The goal of the fitting algorithms is to optimize the vector of parameters $u=u({\mu }_a,{\mu }_s^{\prime },r,\theta )$ , such that the forward generated radiance $L_{forward}\left(u\right)$ matches the MC generated radiance $L(r,{\theta }^i{)}_{MC}$ . Here, $r$ and $\theta$ are assumed known and are thus fixed during fitting. The optical properties ${\mu }_a$ and ${\mu }_s^{\prime }$ are left as free parameters. A Levenberg-Marquardt algorithm from the Matlab 6.12 function fmincon (Mathworks) was used to optimize the values of ${\mu }_a$ and ${\mu }_s^{\prime }$ that minimized the chi-squared function defined as:

While initial guesses for ${\mu }_a$ and ${\mu }_s^{\prime }$ close to the true solution aid in expediting convergence, to avoid bias, all initial guesses for ${\mu }_a$ and ${\mu }_s^{\prime }$ were purposely chosen to be unphysical and set to negative values. The optical properties were constrained such that $0.0001<{\mu }_a<20{\mathrm{cm}}^{-1}$ and $0.001<{\mu }_s^{\prime }<100{\mathrm{cm}}^{-1}$ . Where it was assumed that knowledge of the effective attenuation coefficient was available, an additional constraint was applied such that ${\mu }_a$ and ${\mu }_s^{\prime }$ were restricted to obey the equation: ${\mu }_{\mathrm{eff}}=\left[3{\mu }_a\right({\mu }_s^{\prime }+{\mu }_a){]}^{1∕2}$ . Here, the known ${\mu }_{\mathrm{eff}}$ was calculated using the MC forward properties.

4.1.

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 $180\mathrm{deg}$ . To investigate the potential errors due to probe dimensions, forward radiance distributions were generated at two sensor positions (0.5 and $1.0\mathrm{cm}$ ) 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 ${\mu }_s^{\prime }=10{\mathrm{cm}}^{-1}$ and variable ${\mu }_a$ of 0.01, 0.1, $1{\mathrm{cm}}^{-1}$ were selected. The effects of finite probe dimensions were simulated by assuming that the probe size results in a radial shift $\mathrm{\Delta }r$ that is a function of the probe diameter $r_p$ and detection angle $\theta$ :

Fig. 3

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

Probe sizes ranging up to $1\mathrm{mm}\mathrm{diam}$ 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 ${\mu }_a=0.01{\mathrm{cm}}^{-1}$ 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 $r_p=0\mathrm{mm}$ ) was used to fit the radiance data in Fig. 4. All forward generated data were normalized to the radiance value at $0\mathrm{deg}$ . The resulting error $\epsilon$ was then calculated as:

Fig. 4

Change in radiance distribution at $r=1\mathrm{cm}$ 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 $1\mathrm{mm}$ . Ten different curves are shown for each optical property set with probe radii between 0 and $0.055\mathrm{cm}$ .

Figure 5 shows the relative error in the recovered absorption coefficient and reduced scattering coefficient, with and without ${\mu }_{\mathrm{eff}}$ constraints, as a function of increasing probe size for a sensor positioned at $1\mathrm{cm}$ . The entire radiance dataset between 0 and $180\mathrm{deg}$ 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 ${\mu }_a$ and an underestimation of ${\mu }_s^{\prime }$ . 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) ${\mu }_a$ and (b) and (d) ${\mu }_s^{\prime }$ if finite probe size is unaccounted for during fitting. Fits are for radiance data at $1\mathrm{cm}$ (a) and (b) without and (c) and (d) with a ${\mu }_{\mathrm{eff}}$ constraint. A ${\mu }_s^{\prime }$ of $10{\mathrm{cm}}^{-1}$ was employed for all cases. Similar trends were observed at $0.5\mathrm{cm}$ .

The effects of this change are significant. For a typical probe radius of $\sim 0.3\mathrm{mm}$ $\left(600\mathrm{\mu }\mathrm{m}\mathrm{diam}\right)$ , errors of greater than 150% can result for the unconstrained case, compared to a $\sim 20\%$ error for the constrained case. For ${\mu }_a<0.1{\mathrm{cm}}^{-1}$ , ${\mu }_a$ is highly sensitive to probe radius, whereas ${\mu }_s^{\prime }$ errors are under $\sim 10\%$ . Interestingly, for ${\mu }_a=1{\mathrm{cm}}^{-1}$ , the opposite trend is observed with ${\mu }_a$ errors decreasing to below $\sim 25\%$ and ${\mu }_s^{\prime }$ errors increasing to $\sim 40\%$ .

It should be noted that for the ${\mu }_a=1{\mathrm{cm}}^{-1}$ case, a ${\mu }_{\mathrm{eff}}$ constraint does not significantly minimize the error due to finite probe size for either ${\mu }_s^{\prime }$ or ${\mu }_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 ${\mu }_{\mathrm{eff}}$ constraint is employed.

In the lower absorption cases ( ${\mu }_a=0.01$ , $0.1{\mathrm{cm}}^{-1}$ ), adding a ${\mu }_{\mathrm{eff}}$ constraint significantly improved recovery of ${\mu }_a$ but provided little improvement in recovered ${\mu }_s^{\prime }$ . An examination of Eq. 5 provides some insight into these trends. From Eq. 5, we see that as ${\mu }_a\to 0$ , $K\sim 1∕{\mu }_s^{\prime }[1∕r]$ . This means that while ${\mu }_s^{\prime }$ can be directly recovered from $K$ , ${\mu }_a$ can attain virtually any value and is highly sensitive to small errors. In this case, a ${\mu }_{\mathrm{eff}}$ constraint provides significant improvements in the convergence of ${\mu }_a$ .

Therefore, if the improvement in recovered ${\mu }_a$ is significantly greater than the error due to finite probe size, an overall reduction in ${\mu }_a$ error is observed. By contrast, a similar improvement in ${\mu }_s^{\prime }$ would not be observed, since a ${\mu }_{\mathrm{eff}}$ constraint provides little benefit in ${\mu }_s^{\prime }$ 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.

4.2.

Numerical aperture

Experimentally, radiance information is integrated over the numerical aperture (NA) of the sensor. Dickey ¹³ 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 $30\mathrm{deg}$ ). 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, ¹³ 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 ${\mu }_{\mathrm{eff}}$ constrained fits. For a NA less than $\sim 0.25$ (acceptance angle $\sim 15\mathrm{deg}$ ), errors in recovered ${\mu }_a$ were less than 20% and 2% for unconstrained and ${\mu }_{\mathrm{eff}}$ constrained fits, respectively.

Fig. 6

Change in radiance distribution at $1\mathrm{cm}$ 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 $(\mathrm{NA}=0\mathrm{deg})$ , 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) ${\mu }_a$ and (b) and (d) ${\mu }_s^{\prime }$ if numerical aperture is ignored during fitting. Fits are for radiance data at $1\mathrm{cm}$ (a) and (b) without and (c) and (d) with a ${\mu }_{\mathrm{eff}}$ constraint.

While the optical properties are both overestimated for the unconstrained case, an underestimation in ${\mu }_a$ and an overestimation in ${\mu }_s^{\prime }$ was observed for the constrained method. For the unconstrained case, the overestimation of both ${\mu }_a$ and ${\mu }_s^{\prime }$ are likely explained by optical similarity. This is because similar solution sets for radiance tend to maintain the same albedo.¹⁰ This is not the case for ${\mu }_{\mathrm{eff}}$ , where a decrease in ${\mu }_a$ must be balanced by an increase in ${\mu }_s^{\prime }$ 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 $\sim 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.

5.1.

Experimental Setup

Figure 8 shows a schematic of the experimental setup.

Fig. 8

Experimental setup.

A $650\text{-}\mathrm{\mu }\mathrm{m}\text{-}\mathrm{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.5\mathrm{cm}$ from the radiance probe with $\sim 0.01\text{-}\mathrm{cm}$ resolution. The radiance probe was constructed according to the method of Dickey ¹³ by attaching a $700\text{-}\mathrm{\mu }\mathrm{m}$ right-angle prism (Melles Griot) to a cleaved plane cut fiber. The acceptance angle of the radiance sensor was measured to be $\sim 10\mathrm{deg}$ $(\mathrm{NA}\sim 0.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 $2\mathrm{deg}$ . At each angular position, an average of $\sim$ 300 photovoltage acquisitions was read into a PC and stored into an output file containing a complete rotation between 0 (facing the source) and $180\mathrm{deg}$ (facing away from the source). A complete set of 91 readings was acquired in approximately $90\mathrm{s}$ . 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\times 11\times 11\text{-}\mathrm{cm}$ acrylic box. This allowed for the radiance probe to be positioned $\sim 4$ to $5\mathrm{cm}$ 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 $1\mathrm{cm}$ from the source, while spatial radiance data (at a fixed collection angle $\theta$ of $90\mathrm{deg}$ ) were collected by translating the source at distances between 0.5 to $1.5\mathrm{cm}$ 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).

5.2.

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 ${\mu }_a$ ranged from $\sim 0$ to $5{\mathrm{cm}}^{-1}$ . Note that unlike the MC simulation, we did not perform measurements at lower transport albedos $(\sim 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,²⁶ batch differences can result in a wide range of ${\mu }_s^{\prime }$ .⁶ 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 ${\mu }_{\mathrm{eff}}$ measured for each quantity of added absorption ${\mu }_{a,dye}$ . In our experiments, ${\mu }_{\mathrm{eff}}$ was obtained by performing spatial scans of the radiance $\left[L\right(\theta =90\mathrm{deg},r\left)\right]$ photovoltage $V$ between 0.5 and $1.5\mathrm{cm}$ in 0.05-cm increments, and determining the slope of a plot of $\ln (V\cdot r)$ versus $r$ . According to the PM technique, plotting ${\mu }_{\mathrm{eff}}^2$ as a function of ${\mu }_{a,dye}$ gives:

5.3.

Results

Figures 9a through 9d demonstrate the recovered values of ${\mu }_a$ and ${\mu }_s^{\prime }$ , respectively. In Figs. 9a and 9c $\left({\mu }_a\right)$ , 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 ${\mu }_a$ determined from the PM technique. In Figs. 9b and 9d $\left({\mu }_s^{\prime }\right)$ , the dashed line represents the expected ${\mu }_s^{\prime }$ 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) ${\mu }_a$ and (b) and (d) ${\mu }_s^{\prime }$ obtained from $L_{P3}$ 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 ${\mu }_{\mathrm{eff}}$ constraint.

Trends observed experimentally were similar to those seen in the MC analysis. With the exception of the lower limit of ${\mu }_a=0.01{\mathrm{cm}}^{-1}$ , unconstrained single sensor fits recovered both optical properties to within $\sim 10\%$ at the 0.5-cm position and to within $\sim 20\%$ at the 1-cm position. Above $\sim 1{\mathrm{cm}}^{-1}$ , fits for both sensor positions displayed a significant breakdown in recovered ${\mu }_s^{\prime }$ , although the 1-cm sensor recovered ${\mu }_a$ accurately to within $\sim 30\%$ , up to $\sim 4{\mathrm{cm}}^{-1}$ . The two sensor fitting algorithm (no ${\mu }_{\mathrm{eff}}$ constraint) only showed minor improvements in accuracy compared to the single sensor fits. However, fits performed using a ${\mu }_{\mathrm{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 $\sim 20\%$ for a reduced albedo range of 0.9 to 0.999 at sensor positions between 0.5 and $1\mathrm{cm}$ . Furthermore, with the additional information regarding the spatial distribution of light (i.e., ${\mu }_{\mathrm{eff}}$ constraint), radiance-based fitting algorithms are able to experimentally determine optical properties to within $\sim 12\%$ over a larger albedo range of 0.67 to 0.999.

5.4.

Comparison with Other Interstitial Techniques

It is helpful to compare the ${\mu }_{\mathrm{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 ${\mu }_{\mathrm{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 ${\mu }_{\mathrm{eff}}$ radiance technique be implemented practically? If the ${\mu }_{\mathrm{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 ${\mu }_{\mathrm{eff}}$ estimate.

Alternatively, ${\mu }_{\mathrm{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 ${\mu }_{\mathrm{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 ${\mu }_{\mathrm{eff}}$ will require measurements of the radiance at a minimum of two sensor positions. In this case, the ${\mu }_{\mathrm{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 ${\mu }_{\mathrm{eff}}$ radiance technique appears to be comparable to the two-sensor frequency-domain methodology of Xu and Patterson.¹² Similar to the work of Xu and Patterson,¹² 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,¹² our ${\mu }_{\mathrm{eff}}$ was calculated from experimental measurements at multiple distances (0.5 to $1.5\mathrm{cm}$ ) and not from two positions. It is likely that somewhat larger errors will result when only two sensors are used to determine ${\mu }_{\mathrm{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 $\sim 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, $L_{P3}(\theta ,r)$ curve. For the case of steady-state radiance measurements at a single wavelength, the additional information of a ${\mu }_{\mathrm{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.