**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.

## 1.

## Introduction

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
${\mu}_{a}$
, and the reduced scattering coefficient
${\mu}_{s}^{\prime}$
. 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 ${\mu}_{t}^{\prime}={\mu}_{a}+{\mu}_{s}^{\prime}$ , 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
${\mu}_{t}^{\prime}$
and
${\mu}_{\mathrm{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
$(\sim 0.5\phantom{\rule{0.3em}{0ex}}\mathrm{cm})$
. While these factors do not severely limit the characterization of
${\mu}_{\mathrm{eff}}$
, they can lead to significant difficulties in accurately determining
${\mu}_{t}^{\prime}$
. 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
${\mu}_{\mathrm{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 Patterson^{12} 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 study^{13} 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
$(\sim 1\phantom{\rule{0.3em}{0ex}}\mathrm{cm})$
. 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.5\phantom{\rule{0.3em}{0ex}}\mathrm{cm}$
.

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 (
$\sim 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
$\sim 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 Vitkin^{10} that employs an additional
${\mu}_{\mathrm{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 ${\mu}_{\mathrm{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 ${\mu}_{a}$ and ${\mu}_{s}^{\prime}$ are known exactly, and 2. the inversion is performed free of “experimental” errors in sensor positioning or detection angle. Both the unconstrained and ${\mu}_{\mathrm{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.

## 2.

## Methods and Theory for Radiance Optical Property Determination

## 2.1.

### Radiative Transfer Equation

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

## Eq. 1

$$L(\stackrel{\u20d7}{r},\widehat{s})+\widehat{s}\cdot \nabla L(\stackrel{\u20d7}{r},\widehat{s})+({\mu}_{s}+{\mu}_{a})L(\stackrel{\u20d7}{r},\widehat{s})={\mu}_{s}{\int}_{4\pi}p(\widehat{s},{\widehat{s}}^{\prime})L(\stackrel{\u20d7}{r},\widehat{s})d{\omega}^{\prime}+S(\stackrel{\u20d7}{r},\widehat{s}).$$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^{15} demonstrated that the Henyey-Greenstein (H-G) phase function^{16} 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{\u20d1}{r},\widehat{s})$
,
$p(\widehat{s},{\widehat{s}}^{\prime})$
, and
$S(\stackrel{\u20d7}{r},\widehat{s})$
in spherical harmonics to the
$N$
’th order.^{17} The
$N=1$
expansion expresses the radiance
$L(\stackrel{\u20d1}{r},\widehat{s})$
as the linear combination of the isotropic fluence
$\varphi \left(\stackrel{\u20d1}{r}\right)$
and a linearly anisotropic flux
$\stackrel{\u20d1}{j}\left(\stackrel{\u20d1}{r}\right)$
^{18}:

In spherical geometry, the resulting diffusion theory (DT) radiance expression for a point source is^{18}:

## Eq. 2

$${L}_{\mathrm{P}1}(\stackrel{\u20d7}{r},\theta )={L}_{\mathrm{DT}}(\stackrel{\u20d7}{r},\theta )=\frac{{P}_{o}}{4\pi}\frac{1}{4\pi rD}[1+3(D\u2215r+\sqrt{{\mu}_{a}D})\mathrm{cos}\phantom{\rule{0.2em}{0ex}}\theta ]\mathrm{exp}(-{\mu}_{\mathrm{eff}}r).$$Here, $D=1\u22153({\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\u22152}$ is the effective attenuation coefficient, and $\theta $ is the angle between the radiance direction $\widehat{s}$ and the outgoing radial normal vector $\stackrel{\u20d1}{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^{10} 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.^{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
$\delta $
-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 is^{19}:

## Eq. 3

$${L}_{P3}(r,\widehat{s})={P}_{o}\sum _{l=0}^{3}\frac{2l+1}{4\pi}[{C}^{\prime}{h}_{l}(-{\nu}^{-}){Q}_{l}(-{\nu}^{-}r)+{D}^{\prime}{h}_{l}(-{\nu}^{+}){Q}_{l}(-{\nu}^{+}r)]{P}_{l}(\widehat{s}\cdot \widehat{r}),$$^{10}Finally, the ${P}_{l}$ 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
${\mu}_{s}^{\prime}$
. However, Hull and Foster^{19} 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$
.^{19}

## 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.^{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

$$\varphi \left(r\right)=\frac{{P}_{o}}{4\pi Dr}\phantom{\rule{0.2em}{0ex}}\mathrm{exp}(-{\mu}_{\mathrm{eff}}r).$$^{21, 22}However, absolute fluence measurements are

*not*required to determine ${\mu}_{\mathrm{eff}}$ , since the slope is a relative measurement. Ideally, two measurements at different spatial positions are required to obtain ${\mu}_{\mathrm{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, ${\mu}_{a}$ and ${\mu}_{s}^{\prime}$ , of turbid media using the angular radiance distribution $L(\theta ,\stackrel{\u20d1}{r})$ at position $\stackrel{\u20d1}{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 $\mathrm{cos}\left(\theta \right)$ now yields a slope $K$ , given by:

## Eq. 5

$$K=\frac{1}{({\mu}_{s}^{\prime}+{\mu}_{a})}[1\u2215r+\sqrt{3{\mu}_{a}({\mu}_{s}^{\prime}+{\mu}_{a})}]=3\frac{{\mu}_{a}}{{\mu}_{\mathrm{eff}}^{2}}(\frac{1}{r}+{\mu}_{\mathrm{eff}}),$$^{10}

## 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^{10} 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.^{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 suggested^{10} 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.^{10}

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.^{23} 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.^{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
$\sigma =100\sqrt{C}\u2215C$
, 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{\u20d1}{r},\theta )$ were recorded for each spatial position $\stackrel{\u20d1}{r}$ and each detection angle $\theta =a\phantom{\rule{0.2em}{0ex}}\mathrm{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{\u20d1}{r},\theta )$ was then calculated according to the equation:

## Eq. 6

$$L(\stackrel{\u20d1}{r},\theta )=\frac{1}{{N}_{o}}\frac{N(\stackrel{\u20d1}{r},\theta )}{\left(4\pi {r}^{2}dr\right)\left[2\pi \phantom{\rule{0.2em}{0ex}}\mathrm{sin}\left(\theta \right)d\theta \right]}\frac{1}{{\mu}_{a}},$$^{18, 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).

## 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:

## Eq. 7

$${\chi}^{2}({\mu}_{a},{\mu}_{s}^{\prime})=\sum _{i=1}^{N\left(\theta \right)}{[L{(r,{\theta}^{i})}_{\mathrm{MC}}-{L}_{\mathrm{forward}}{(r,{\theta}^{i},{\mu}_{a},{\mu}_{s}^{\prime})}_{P3}]}^{2}.$$^{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 ${\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\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ and $0.001<{\mu}_{s}^{\prime}<100\phantom{\rule{0.3em}{0ex}}{\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\u22152}$ . Here, the known ${\mu}_{\mathrm{eff}}$ was calculated using the MC forward properties.

## 3.

## 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,
${\mu}_{s}^{\prime}$
was held constant at
$10\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$
while
${\mu}_{a}$
ranged from 0.001 to
$10\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$
, such that the transport albedo
${a}^{\prime}={\mu}_{s}^{\prime}\u2215{\mu}_{s}^{\prime}+{\mu}_{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
$1\phantom{\rule{0.3em}{0ex}}\mathrm{cm}$
, 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
${\mu}_{\mathrm{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.005\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$
. 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}

Beyond this lower limit and up to
${\mu}_{a}$
, values as high as
$\sim 3\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$
, fits at both sensor positions returned absorption coefficients to within 15% of the expected value (worst case
$\sim 25\%$
). Fits in the 0.05 to
$1.75\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$
absorption range were 10% more accurate with 0.5-cm sensor data than 1.0-cm data, whereas fits in the 1.75 to
$3\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$
range were 15% less accurate. This indicates that accurate
${\mu}_{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
$\sim 10\%$
(worst case 16%) for data at both sensor positions up to a
${\mu}_{a}$
of
$\sim 3\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$
.

For ${\mu}_{a}>3\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ , errors of greater than $\sim 30\%$ were obtained for both ${\mu}_{a}$ and ${\mu}_{s}^{\prime}$ . However, even at high ${\mu}_{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 ${\mu}_{\mathrm{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 ${\mu}_{a}$ and ${\mu}_{s}^{\prime}$ recovered to within $\sim 10\%$ between the ${\mu}_{a}$ range of 0.05 to $3\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ . However, while improvements in recovered optical properties were observed, two sensor fits still exhibited lower $\left(0.005\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}\right)$ and upper $(\sim 3\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1})$ 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 $\theta $ ) 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 ( ${\mu}_{\mathrm{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 ${\mu}_{\mathrm{eff}}$ . The improvements in optical property recovery are substantial. This is demonstrated in Figs. 2c and 2d, where it can be seen that both ${\mu}_{a}$ and ${\mu}_{s}^{\prime}$ are now recovered to within $\sim 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 ${\mu}_{\mathrm{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 ${\mu}_{\mathrm{eff}}$ constraint method is preferable compared to fitting two radiance slopes. Practically, ${\mu}_{\mathrm{eff}}$ can be determined via literature values or calculated from experimental measurements as described in Sec. 5.4.

## 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 ${\mu}_{a}$ and an underestimation of ${\mu}_{s}^{\prime}$ (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
${\mu}_{\mathrm{eff}}$
method in this section, and compare our results with the single sensor unconstrained approach.

## 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\phantom{\rule{0.3em}{0ex}}\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\phantom{\rule{0.3em}{0ex}}\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\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ and variable ${\mu}_{a}$ of 0.01, 0.1, $1\phantom{\rule{0.3em}{0ex}}{\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 $ :

## Eq. 8

$$\mathrm{\Delta}r\left(\theta \right)={r}_{p}(1-\mathrm{cos}\phantom{\rule{0.2em}{0ex}}\theta ).$$Probe sizes ranging up to $1\phantom{\rule{0.3em}{0ex}}\mathrm{mm}\phantom{\rule{0.2em}{0ex}}\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\phantom{\rule{0.3em}{0ex}}{\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\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ ) was used to fit the radiance data in Fig. 4. All forward generated data were normalized to the radiance value at $0\phantom{\rule{0.3em}{0ex}}\mathrm{deg}$ . The resulting error $\epsilon $ was then calculated as:

## Eq. 9

$$\epsilon =\frac{{\mu}_{\mathrm{fit}}-{\mu}_{\mathrm{true}}}{{\mu}_{\mathrm{true}}}\times 100\%,$$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\phantom{\rule{0.3em}{0ex}}\mathrm{cm}$ . The entire radiance dataset between 0 and $180\phantom{\rule{0.3em}{0ex}}\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.

The effects of this change are significant. For a typical probe radius of $\sim 0.3\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ $\left(600\phantom{\rule{0.3em}{0ex}}\mathrm{\mu}\mathrm{m}\phantom{\rule{0.3em}{0ex}}\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\phantom{\rule{0.3em}{0ex}}{\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\phantom{\rule{0.3em}{0ex}}{\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\phantom{\rule{0.3em}{0ex}}{\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\phantom{\rule{0.3em}{0ex}}{\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\u2215{\mu}_{s}^{\prime}[1\u2215r]$ . 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
^{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
$30\phantom{\rule{0.3em}{0ex}}\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, ^{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
${\mu}_{\mathrm{eff}}$
constrained fits. For a NA less than
$\sim 0.25$
(acceptance angle
$\sim 15\phantom{\rule{0.3em}{0ex}}\mathrm{deg}$
), errors in recovered
${\mu}_{a}$
were less than 20% and 2% for unconstrained and
${\mu}_{\mathrm{eff}}$
constrained fits, respectively.

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.^{10} 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.

## Experimental Evaluation in Phantoms

## 5.1.

### Experimental Setup

Figure 8 shows a schematic of the 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\phantom{\rule{0.3em}{0ex}}\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
^{13} 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\phantom{\rule{0.3em}{0ex}}\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\phantom{\rule{0.3em}{0ex}}\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\phantom{\rule{0.3em}{0ex}}\mathrm{deg}$ (facing away from the source). A complete set of 91 readings was acquired in approximately $90\phantom{\rule{0.3em}{0ex}}\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\phantom{\rule{0.3em}{0ex}}\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\phantom{\rule{0.3em}{0ex}}\mathrm{cm}$ from the source, while spatial radiance data (at a fixed collection angle $\theta $ of $90\phantom{\rule{0.3em}{0ex}}\mathrm{deg}$ ) were collected by translating the source at distances between 0.5 to $1.5\phantom{\rule{0.3em}{0ex}}\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\phantom{\rule{0.3em}{0ex}}{\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,^{26} batch differences can result in a wide range of
${\mu}_{s}^{\prime}$
.^{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 ${\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\phantom{\rule{0.3em}{0ex}}\mathrm{deg},r\left)\right]$ photovoltage $V$ between 0.5 and $1.5\phantom{\rule{0.3em}{0ex}}\mathrm{cm}$ in 0.05-cm increments, and determining the slope of a plot of $\mathrm{ln}(V\cdot r)$ versus $r$ . According to the PM technique, plotting ${\mu}_{\mathrm{eff}}^{2}$ as a function of ${\mu}_{a,dye}$ gives:

## Eq. 10

$${\mu}_{\mathrm{eff}}^{2}=3{\mu}_{a,\mathrm{dye}}{\mu}_{s,IL}^{\prime}+3{\mu}_{a,IL}{\mu}_{s,IL}^{\prime},$$^{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 ${\mu}_{s,IL}^{\prime}$ of $\sim 9.43\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ and a ${\mu}_{a,IL}$ of $\sim 0.0049\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ at $690\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ .

## 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.

Trends observed experimentally were similar to those seen in the MC analysis. With the exception of the lower limit of ${\mu}_{a}=0.01\phantom{\rule{0.3em}{0ex}}{\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\phantom{\rule{0.3em}{0ex}}{\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\phantom{\rule{0.3em}{0ex}}{\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\phantom{\rule{0.3em}{0ex}}\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.^{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
${\mu}_{\mathrm{eff}}$
was calculated from experimental measurements at *multiple* distances (0.5 to
$1.5\phantom{\rule{0.3em}{0ex}}\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.

## 6.

## 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
$1\phantom{\rule{0.3em}{0ex}}\mathrm{cm}$
, 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
${\mu}_{\mathrm{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
$\delta $
-P1 approximation^{20} by Hayakawa
^{29} and You, Hayakawa, and Venugopalan^{30} 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 Grobe^{31} 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.

## Acknowledgments

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.