## Abstract

Interstitial determination of the tissue optical properties is important in biomedicine, especially for interstitial laser therapies. Continuous wave (CW) radiance techniques which examine light from multiple directions have been proposed as minimally invasive methods for determining the optical properties under an interstitial probe arrangement. However, both the fitting algorithm based on the P3 approximation and the analytical method based on the diffusion approximation (DA), which are currently used recovery algorithms, cannot extract the optical properties of tissue with low transport albedos accurately from radiance measurements. In this paper, we proposed an incomplete P3 approximation for the radiance, the P3_{in} for short, which is the asymptotic part of the solution for the P3 approximation. The relative differences between the P3_{in} and the P3 were within 0.48% over a wide range of clinically relevant optical properties for measurements at source detector separations (SDS) from 5 mm to 10 mm and angles from $0\xb0$ to $160\xb0$. Based on the P3_{in}, we developed an analytical method for extracting the optical properties directly using simple expressions constructed from the radiance measurements at only two SDSs and four angles. The developed recovery algorithm was verified by simulated and experimental radiance data. The results show that both the absorption and reduced scattering coefficients were recovered accurately with relative errors within 5.28% and 3.86%, respectively, from the simulated data and with relative errors within 10.82% and 10.67%, respectively, from the experimental data over a wide range of albedos from 0.5 to 0.99. Since the developed P3_{in}-based radiance technique can obtain the optical properties rapidly from the measurements at only two SDSs and four angles, it is expected to be used for *in vivo* and *in situ* determination of the optical properties in online treatment planning during laser therapies.

© 2017 Optical Society of America

## 1. Introduction

Quantification of the tissue optical properties has long been an area of extensive research in biomedicine. Specifically, the *in vivo* determination of the optical properties can be used to calculate the fluence distribution of light for both pretreatment planning and online optimization during photodynamic therapy or thermal therapy [1–3]. Besides, the determination of the optical properties also can be applied for the early cancer detection of tissue according to the differences in the optical properties between cancerous and normal tissues, caused by the variations in extracellular and nuclear structures and chromophore compositions [4, 5].

The continuous wave (CW) techniques are commonly employed to determinate the optical properties due to relatively low cost and simple equipment. Numerous studies have been conducted for noninvasively obtaining the optical properties of tissue such as skin from the spatially resolved surface reflectance measurements [6, 7]. However, applications of such methods are limited by shallow light penetration and inappropriate for diagnosis of internal organs such as the prostate. In such cases, the minimally invasive measurements are desirable and there has been a good deal of work done on this topic. One part of groups employed the interstitial spectroscopy techniques by a single encapsulated probe for determining the optical properties. Bargo *et al.* implemented a spectrally resolved CW diffuse reflectance method where only two fibers (one source and one detector) spaced 2.5 mm apart are used for the determination of the optical properties interstitially [8]. However, such technique requires extensive calibration with a large library of optical phantoms [9]. Baran *et al.* developed a custom fiber optic probe which consisted of six side-firing fibers axially separated by 2 mm to interstitially obtain spatially and spectrally resolved reflectance data from turbid media [9]. The optical properties were extracted by a Monte Carlo (MC)-based fitting algorithm. However, the relatively small number of source detector separations (SDS) available by this probe places potential limits on the accuracy of the optical property determination [10]. Another part of groups inserted the source and detector fibers into the tissue at a prescribed spacing using guidance templates. Zhu *et al.* measured the ratio of light fluence rate to source power along a linear channel at a fixed distance of 5 mm and the diffusion approximation (DA) was used for fitting the measured data to extract the optical properties [11, 12]. However, an integrating sphere which is difficult to be implemented *in vivo* is required for calibrating an isotropic detector and a light source power [13].

Recently, a CW radiance technique measuring light from multiple directions has been proposed as another minimally invasive method for determining the optical properties under an interstitial probe arrangement, since such method can minimize the required number of SDSs [13–17]. Dickey *et al.* first demonstrated the capability to obtain the optical properties of liquid phantoms from radiance measurement data at a single SDS by using the P3 approximation for fitting [14]. However, their research was limited to a single optical property pair and a set of experimental parameters. Chin *et al.* had further investigated the determination of the optical properties using the P3 approximation for fitting from the MC simulated and experimentally measured radiance caused by an isotropic source [15, 16]. However, the shape-based fitting algorithm disclosed a significant breakdown in estimating the optical properties for the transport albedos, ${a}^{\prime}$ defined as ${a}^{\prime}={{\mu}^{\prime}}_{s}/({{\mu}^{\prime}}_{s}+{\mu}_{a})$, less than ~0.9 due to the lack of inherent uniqueness in fitting algorithms, where ${\mu}_{a}$ and ${{\mu}^{\prime}}_{s}$ are the absorption coefficient and the reduced scattering coefficient, respectively. Grabtchak *et al.* derived two analytical expressions based on the DA to extract the effective attenuation coefficient, ${\mu}_{eff}$, and the diffusion coefficient, *D*, directly from radiance measurements caused by an isotropic source [13, 17]. Although such an analytical method could avoid the problems of lack of inherent uniqueness in fitting algorithms, the optical properties can be recovered finely just under the condition of ${{\mu}^{\prime}}_{s}>>{\mu}_{a}$ for which the DA is valid.

For the tissue with high absorption or low scattering that may be encountered in the clinical environment [18, 19], a higher order approximation for radiance is required for recovering the optical properties analytically. The P3 approximation, the P3 for short, is naturally chosen as a light propagation model to recover the optical properties. The P3 has a solution consisting of two parts; one is the transient part which contributes significantly only for the SDSs of less than two transport mean free paths, $2{{l}^{\prime}}_{t}$, where ${{l}^{\prime}}_{t}=1/({{\mu}^{\prime}}_{s}+{\mu}_{a})$, and the other is the asymptotic part which dominates for the large SDSs [20]. However, analytical expressions for recovering the optical properties cannot be constructed based on the P3. Besides, small SDSs less than $2{{l}^{\prime}}_{t}$ are not used in interstitial radiance measurements, generally. Thus, in this paper, we have proposed an incomplete P3 approximation for radiance, the P3_{in} for short, which consists of only the asymptotic part of the solution for the P3. The proposed P3_{in} has a simpler expression than the complete solution for the P3, but provides a higher accuracy than the DA for modeling light propagation in the turbid media with high absorption or low scattering. Based on the P3_{in}, we have developed analytical formulae for extracting the optical properties directly by using simple analytical expressions constructed from CW radiance measurements at only two SDSs and four angles.

## 2. Methodology

#### 2.1 The P3_{in} approximation for radiance

When the PN method is applied to solve the radiative transfer equation (RTE), the angular quantities are expanded in spherical harmonics, and for the case of spherical symmetry the radiance is given by Eq. (1) in a spherical coordinate [21]:

*L*is the radiance, $r$ is the position, $\theta $ is the angle between the direction of propagation and the direction of the position $r$ from the origin, ${S}_{0}$ is the source power, ${\phi}_{l}(r)$ is the

*l*-th order Legendre moment, and ${P}_{l}$ is Legendre polynomial of

*l*-th degree. In the case of the P3, the Legendre moments are given by the following set of the ordinary differential equations for $l=0,\mathrm{...},3$ [22]:

*l-*th order absorption coefficients with

*g*and ${\mu}_{s}$ indicating the anisotropy factor and scattering coefficient, respectively. By setting ${\phi}_{4}(r)=0$ in Eq. (2), the solution of the radiance in the P3 for an isotropic point source in an infinite medium is given as [20]:

In the Eq. (3), the terms associated with ${h}_{l}({\nu}^{+}){Q}_{l}(-{\nu}^{+}r)$ represent the transient part and the terms associated with ${h}_{l}({\nu}^{-}){Q}_{l}(-{\nu}^{-}r)$ represent the asymptotic part. As ${Q}_{l}(-{\nu}^{+}r)$ $\text{(}l=0,\mathrm{...},3\text{)}$ decay rapidly to be zero with the increase in *r* due to the relatively large value of ${\nu}^{+}$, the transient part mainly contributes to the total radiance only for the SDSs less than $2{{l}^{\prime}}_{t}$. The value of ${\nu}^{-}$ is relatively small, and ${Q}_{l}(-{\nu}^{-}r)$ $\text{(}l=0,\mathrm{...},3\text{)}$ decay slowly, thus the asymptotic part dominates for large SDSs.

Conventionally, the P3 approximation has been used for recovering the optical properties analytically from the radiance measurements. For each $l=0,\mathrm{...},3$, the values of ${S}_{0}[{C}^{\prime}{h}_{l}({\nu}^{-}){Q}_{l}(-{\nu}^{-}r)+{D}^{\prime}{h}_{l}({\nu}^{+}){Q}_{l}(-{\nu}^{+}r)]$ in Eq. (3) can be obtained by measuring the radiances at different angles with a fixed *r* to determine the values of ${P}_{l}(cos\theta )$. However, the values of ${S}_{0}{C}^{\prime}{h}_{l}({\nu}^{-})$ and ${S}_{0}{D}^{\prime}{h}_{l}({\nu}^{+})$ cannot be separately obtained by measuring the radiances at different SDSs to determine the values of ${Q}_{l}(-{\nu}^{\pm}r)$ because ${\nu}^{\pm}$ appearing in ${Q}_{l}(-{\nu}^{\pm}r)$ are unknown. Thus ${h}_{l}({\nu}^{-})$ and ${h}_{l}({\nu}^{+})$ cannot be obtained analytically, and from Eq. (5) it is deduced that the optical properties cannot be recovered in the form of analytical expressions.

Analytical expressions for recovering the optical properties cannot be constructed based on the P3, but small SDSs less than $2{{l}^{\prime}}_{t}$ are usually not used for interstitial radiance measurements. Thus, in this paper, we propose an incomplete P3 approximation for the radiance, the P3_{in} for short, which is the asymptotic part of the solution for the P3. Hereafter, ${\nu}^{-}$ is replaced by $\nu $ for simplicity, and the radiance by the P3_{in} is given as:

#### 2.2 Accuracy of the P3_{in} approximation

In order to evaluate the differences between the P3_{in} and the P3, we examined the influence of the abandoned transient part of the solution in the P3 by the relative differences, $\text{\epsilon}=|L(r,\theta )-{L}_{in}(r,\theta )|/L(r,\theta )$, at different SDSs and angles for different optical properties. In the analysis, the optical properties of prostate with $0.01\text{\hspace{0.17em}}{\text{mm}}^{-1}\le {\mu}_{a}\le 1.0\text{\hspace{0.17em}}{\text{mm}}^{-1}$ and $0.5\text{\hspace{0.17em}}{\text{mm}}^{-1}\le {{\mu}^{\prime}}_{s}\le 4.5\text{\hspace{0.17em}}{\text{mm}}^{-1}$ were considered [18, 19, 23]. The SDSs were chosen from 5 mm to 10 mm which are typical SDSs during interstitial laser therapies [16, 23].

Figure 1 plots the maximum values of $\text{\epsilon}$, ${\text{\epsilon}}_{\text{max}}$, in the optical property range mentioned above at SDSs from 5 mm to 10 mm and angles from $0\xb0$ to $180\xb0$. The results show that the ${\text{\epsilon}}_{\text{max}}$ values are very small for angles from $0\xb0$ to $160\xb0$ and the ${\text{\epsilon}}_{\text{max}}$ values fluctuate strongly for angles from $160\xb0$ to $180\xb0$. It also can be seen that the ${\text{\epsilon}}_{\text{max}}$ values decrease with the increase in *r*. Beside, Fig. 1 also shows that the sharp drops appear at the angls of about $40\xb0$, $115\xb0$ and $145\xb0$, since the P3_{in} and the P3 have almost the same values at such angles. The ${\text{\epsilon}}_{\text{max}}$ values for angles from $0\xb0$ to $160\xb0$ are within 0.48%, 1.94 × 10^{-2%}, 6.48 × 10^{-4%}, 4.78 × 10^{-5%}, 5.67 × 10^{-6%} and 6.76 × 10^{-7%} for the SDSs of 5, 6, 7, 8, 9 and 10 mm, respectively. Whereas, the ${\text{\epsilon}}_{\text{max}}$ values for angles from $160\xb0$ to $180\xb0$ are within 17.11%, 2.46%, 0.18%, 2.12 × 10^{-3%}, 7.71 × 10^{-5%} and 2.20 × 10^{-6%} for the SDSs of 5, 6, 7, 8, 9 and 10 mm, respectively. In summary, the relative differences between the P3_{in} and the P3 are within 0.48% for angles from $0\xb0$ to $160\xb0$ and within 17.11% for angles from $160\xb0$ to $180\xb0$for the optical property range mentioned above at SDSs from 5 mm to 10 mm. It indicates that the P3_{in} is acceptable for an alternative to the P3 as a light propagation model to recover the optical properties for angles from $0\xb0$ to $160\xb0$.

#### 2.3 Recovery of the optical properties based on the P3_{in} approximation

Based on the P3_{in} expressed by Eq. (6), this study formulates an analytical method to recover the optical properties from the radiance measurements, ${L}_{m}(r,\theta )$, at several pairs of $(r,\theta )$. From Eq. (5), it is deduced that the optical properties of the medium can be determined analytically if $\nu $, ${h}_{1}(\nu )$ and ${h}_{2}(\nu )$ were obtained. As the first step for obtaining $\nu $, ${h}_{1}(\nu )$ and ${h}_{2}(\nu )$, the terms of ${S}_{0}{C}^{\prime}{h}_{l}(\nu ){Q}_{l}(-\nu r)$ $(l=0,1,2)$should be determined. As the terms of ${S}_{0}{C}^{\prime}{h}_{l}(\nu ){Q}_{l}(-\nu r)$ with four different numbers of *l* in Eq. (6) are simply added, they can be separated from the radiance measurement data of ${L}_{m}(r,{\theta}_{i})$ $(i=0,\mathrm{...},3)$ at four angles of ${\theta}_{0}$, ${\theta}_{1}$, ${\theta}_{2}$,${\theta}_{3}$with a fixed *r*, as shown in Fig. 2, by a set of linear equations given as:

As the second step for obtaining $\nu $, ${h}_{1}(\nu )$ and ${h}_{2}(\nu )$, the ratio of ${f}_{0}(r,{\theta}_{0},{\theta}_{1},{\theta}_{2},{\theta}_{3})$ to ${f}_{0}({r}^{\prime},{\theta}_{0},{\theta}_{1},{\theta}_{2},{\theta}_{3})$ is taken from the radiance measurements at two SDSs, *r* and ${r}^{\prime}$, as shown in Fig. 2, to delete ${S}_{0}$ for avoiding the measurements of the absolute radiances:

Now $\nu $, ${h}_{1}(\nu )$ and ${h}_{2}(\nu )$ are obtained, and the optical properties can be determined from Eq. (5) as the following:

The advantages of this recovery algorithm of the optical properties are; (i) data fitting, a mass of measurement data and the absolute radiance measurements are not required, and (ii) both ${\mu}_{a}$ and ${{\mu}^{\prime}}_{s}$ are determined directly from the radiance measurements at two SDSs and four angles using the analytical formulas from Eq. (8) to Eq. (14).

## 3. Results

#### 3.1 Simulation evaluation using the P25 approximation

The developed algorithm of recovering the optical properties based on the P3_{in} was evaluated using simulated radiance data calculated by the P25 approximation for the radiance. The P25 approximation solution of the RTE produced by an isotropic source placed inside an infinitely extended turbid medium was used as the light propagation model to generate the simulated data because the excellent agreement between the P25 approximation solution and the results of a MC method has been verified [21, 24]. In this analysis, the values of 1 mm^{−1} and 0.5 mm^{−1} were chosen for ${{\mu}^{\prime}}_{s}$ as typical and low scattering cases, respectively, while ${\mu}_{a}$ ranged from 0.01 mm^{−1} to 0.5 mm^{−1}, so that the albedos ranged from 0.5 to 0.99. The chosen optical property sets roughly covered the range of the albedos expected in interstitial therapies [18, 19, 23, 25]. The radiance measurement data were generated at four angles of $40\xb0$, $80\xb0$, $120\xb0$ and $160\xb0$. This set of angles avoided the radiances at angles larger than $160\xb0$ where the differences between the P3_{in} and the P3 are relative remarkable. The SDS of ${r}^{\prime}$ was fixed at 5 mm, while 6 mm and 10 mm were chosen for the SDSs of *r*. At the same time, the DA was also employed to extract the optical properties from the same simulated radiance data in order to compare with the P3_{in} [13]. The results are presented for the recovery algorithms based on both the DA and the P3_{in}.

We evaluated the performance of the developed algorithm based on the P3_{in} in the case of ${{\mu}^{\prime}}_{s}=1.0\text{\hspace{0.17em}}{\text{mm}}^{-1}$ first. In Figs. 3(a) and 3(b), the results from the algorithm based on the DA show that the accuracy of the recovered optical properties decreases dramatically with the increase in the true ${\mu}_{a}$ and that the accuracy is improved only a little with the increase in *r*. Just in the case of ${\mu}_{a}\le 0.1\text{\hspace{0.17em}}{\text{mm}}^{-1}$, *i.e.*, ${a}^{\prime}\ge 0.91$, both ${\mu}_{a}$ and ${{\mu}^{\prime}}_{s}$ are recovered better with the relative errors (REs) less than 14.75% where the REs are defined as $|{\mu}_{a\_r}-{\mu}_{a\_t}|/{\mu}_{a\_t}$, $|{{\mu}^{\prime}}_{s\_r}-{{\mu}^{\prime}}_{s\_t}|/{{\mu}^{\prime}}_{s\_t}$ for ${\mu}_{a}$ and ${{\mu}^{\prime}}_{s}$, respectively, with ${\mu}_{a\_r}$, ${\mu}_{a\_t}$, ${{\mu}^{\prime}}_{s\_r}$ and ${{\mu}^{\prime}}_{s\_t}$ indicating the recovered ${\mu}_{a}$, the true ${\mu}_{a}$, the recovered ${{\mu}^{\prime}}_{s}$ and the true ${{\mu}^{\prime}}_{s}$. The likely explanation is that the DA provides a poor description of the radiance in the case of a low albedo even at large SDSs [15]. In contrast to the DA, Figs. 3(a) and 3(b) demonstrate that the developed recovery algorithm based on the P3_{in} is able to recover both ${\mu}_{a}$ and ${{\mu}^{\prime}}_{s}$ accurately. A small improvement in the accuracy of the recovered optical properties has also been observed with the increase in *r* similarly to the DA. The REs of the recovered ${\mu}_{a}$ and ${{\mu}^{\prime}}_{s}$ based on the P3_{in} are less than 1.30% and 0.96%, respectively, for the albedos ranging from 0.67 to 0.99.

We also evaluated the performance of the developed algorithm for the case of ${{\mu}^{\prime}}_{s}=0.5\text{\hspace{0.17em}}{\text{mm}}^{-1}$ which may be encountered during interstitial laser therapies [18, 19]. Figures 4(a) and 4(b) show that the recovered ${\mu}_{a}$ and ${{\mu}^{\prime}}_{s}$ based on the P3_{in} are also very close to true values, whereas the recovered ${\mu}_{a}$ and ${{\mu}^{\prime}}_{s}$ based on the DA exhibit significant errors. The recovery algorithm based on the DA disclosed a significant breakdown in recovering ${\mu}_{a}$ with REs larger than 20.72%, although ${{\mu}^{\prime}}_{s}$ is recovered with the REs within 9.54% for the true ${\mu}_{a}$ less than 0.05 mm^{−1}. In contrast to the DA, the REs of the recovered ${\mu}_{a}$ and ${{\mu}^{\prime}}_{s}$ based on the P3_{in} are less than 5.28% and 3.86%, respectively, for the albedos ranging from 0.5 to 0.98. Thus the results show that the recovery algorithm based on the P3_{in} is able to recover the optical properties accurately in the cases of both typical and low scattering.

#### 3.2. Experimental evaluation by liquid phantoms

Figure 5 presents a schematic of the relevant experimental setup. An $850\text{\hspace{0.17em}}\text{\mu m}$ diameter isotropic spherical diffuser (Medlight, Switzerland) was coupled to a 675nm laser (Thorlabs, USA) with an optical attenuator (EXFO, Canada) to illuminate a liquid phantom interstitially. The source fiber was attached to a holder mounted on a two-dimensional precision translation stage that allowed changing the source fiber position relative to the radiance detector probe. The radiance detector probe was constructed similarly to that of Barajas *et al.* by attaching a $\text{500}\text{\mu m}$ right-angle prism to a cleaved plane cut $475\text{\hspace{0.17em}}\text{\mu m}$ fiber with a numerical aperture of 0.22 [26]. The detector fiber was mounted on a motorized rotation stage that allowed complete $360\xb0$ rotation. The detector fiber was connected to a photomultiplier tube (PMT) (Hamamatsu, Japan) for data acquisition.

Similarly to the conditions in the evaluation using the simulated radiances, the radiance measurements were made at four angles of $40\xb0$, $80\xb0$, $120\xb0$ and $160\xb0$. The SDS of ${r}^{\prime}$ was fixed at 5 mm, while 6 mm and 10 mm were chosen for the SDSs of *r*. A gate time of 0.5 s and the average measurements of 10 times were chosen for the PMT. All measurements were background subtracted and were conducted in the center region of a 10 × 10 × 10 cm^{3} black acrylic box which minimized reflections from the box boundary. Tissue simulating phantoms were prepared using a combination of 20% Intralipid and India ink for ${{\mu}^{\prime}}_{s}$ of 1 mm^{−1} and 0.5 mm^{−1} and ${\mu}_{a}$ ranging from 0.01 mm^{−1} to 0.5 mm^{−1}. The concentrations of ink were determined according to the absorbance measured by a spectrometer while the concentrations of Intralipid were determined according to the equations derived by van Staveren *et al* [27].

The developed recovery algorithm based on the P3_{in} was evaluated by the experimental radiance data in the case of typical scattering of ${{\mu}^{\prime}}_{s}=\text{\hspace{0.17em}}1.0\text{\hspace{0.17em}}{\text{mm}}^{-1}$ in Figs. 6(a) and 6(b). The trends obtained by the phantom measurements were similar to those shown in Figs. 3(a) and 3(b) by the corresponding simulation. The results show that the REs of the recovered optical properties based on the DA increase dramatically with the increase in the true ${\mu}_{a}$ and that only for ${\mu}_{a}\le 0.1\text{\hspace{0.17em}}{\text{mm}}^{-1}$, *i.e.*, ${a}^{\prime}\ge 0.91$, both ${\mu}_{a}$ and ${{\mu}^{\prime}}_{s}$ were recovered better with the REs within 16.12%. The results using the recovery algorithm based on the P3_{in} show that significant improvements in the accuracy of the recovered optical properties have been obtained when compared with those based on the DA. Although the REs of the recovered optical properties from the phantom experiments are slightly larger than those from the simulations due to the systematic errors in measurements, the REs of the recovered ${\mu}_{a}$ and ${{\mu}^{\prime}}_{s}$ based on the P3_{in} are less than 10.39% and 9.86%, respectively, for the albedos ranging from 0.67 to 0.99.

The developed recovery algorithm based on the P3_{in} was also evaluated by the experimental radiance data in the case of low scattering of ${{\mu}^{\prime}}_{s}=\text{\hspace{0.17em}}0.5\text{\hspace{0.17em}}{\text{mm}}^{-1}$ in Figs. 7(a) and 7(b). It is observed that the recovery algorithm based on the DA disclosed a breakdown in recovering ${\mu}_{a}$ with the REs larger than 17.07%. ${{\mu}^{\prime}}_{s}$ were recovered finely only for the true ${\mu}_{a}$ less than $\text{0}\text{.05}{\text{mm}}^{-1}$ with the REs within 8.15%. In contrast to the DA, the REs of the recovered ${\mu}_{a}$ and ${{\mu}^{\prime}}_{s}$ based on the P3_{in} are less than 10.82% and 10.67%, respectively, for the albedos ranging from 0.5 to 0.98.

In summary, the recovery algorithm based on the DA can recover the optical properties from simulated and experimental data with the REs within 16.12% only for ${a}^{\prime}\ge 0.91$, whereas the recovery algorithm based on the P3_{in} can recover the optical properties with REs within 10.82% over a wide range of albedos from 0.5 to 0.99.

## 4. Discussion

Although the recovery algorithm based the P3_{in} greatly improved the accuracy of the recovered optical properties from simulated and experimental radiances when compared with the DA, the P5_{in} approximation as a light propagation model to recover the optical properties was also explored to examine whether the accuracy of the recovered optical properties could be further improved. Similarly to the P3_{in}, the solution of the P5_{in} is the asymptotic part of the solution for the radiance by the P5 approximation. The recovery algorithm based on the P5_{in} was developed and the optical properties were extracted using the analytical expressions constructed from the radiance measurements at two SDSs and six angles.

The recovery algorithms based on the P3_{in} and the P5_{in} were evaluated by the simulated radiance data calculated by the P25 for the optical property range of $0.01\text{\hspace{0.17em}}{\text{mm}}^{-1}\le {\mu}_{a}\le 1.0\text{\hspace{0.17em}}{\text{mm}}^{-1}$ and $0.5\text{\hspace{0.17em}}{\text{mm}}^{-1}\le {{\mu}^{\prime}}_{s}\le 4.5\text{\hspace{0.17em}}{\text{mm}}^{-1}$. ${r}^{\prime}$ and *r* were chosen as 5 mm and 6 mm, respectively. Four angles of $40\xb0$, $80\xb0$, $120\xb0$ and $160\xb0$ for the P3_{in} and six angles of $40\xb0$, $70\xb0$, $90\xb0$, $120\xb0$, $140\xb0$ and $160\xb0$ for the P5_{in} were chosen.

Figures 8(a)-8(d) plot the REs of the recovered optical properties based on the P3_{in} and the P5_{in}, respectively. The REs of the recovered ${\mu}_{a}$ and ${{\mu}^{\prime}}_{s}$ based on the P3_{in} are less than 6.59% and 7.49%, respectively, whereas those based on the P5_{in} are less than 4.19% and 4.62%, respectively. Although more radiance measurements are required for optical property recovery based on the P5_{in}, minor improvement in the accuracy of the recovered optical properties by the P5_{in} was obtained when compared with the P3_{in}.

Here, we present a recovery algorithm based on the P3_{in} that is capable of recovering ${\mu}_{a}$ and ${{\mu}^{\prime}}_{s}$ with the REs within 10.82% and 10.67%, respectively, over a wide range of albedos from 0.5 to 0.99. The REs of the recovered optical properties based on the P3_{in} were comparable to the errors found in previous studies of interstitial optical property recovery. Zhu *et al.* demonstrated the recovery of ${\mu}_{a}$ and ${{\mu}^{\prime}}_{s}$ with standard (and maximum) deviations of 8% (15%) and 18% (32%), respectively [12]. The method developed by Chin *et al.* resulted in the recovery of the optical properties to within ~20% for the albedos larger than ~0.9 from relative radiance measurements [9, 16]. Baran *et al.* demonstrated that ${\mu}_{a}$ and ${{\mu}^{\prime}}_{s}$ were recovered with mean errors of 9% and 19%, respectively, for the albedos larger than 0.95 [9]. Grabtchak *et al.* showed that they were able to recover ${\mu}_{eff}$ and *D* within 2% and 40%, respectively, using a liquid phantom of Intralipid-1%, but they did not comment on the accuracy of the recovered ${\mu}_{a}$ and ${{\mu}^{\prime}}_{s}$ directly [13]. It indicates that the developed algorithm based on the P3_{in} in this paper can recover the optical properties of turbid media with albedos as low as 0.5 accurately when compared with the results found in previous studies.

The developed P3_{in}-based radiance technique has three advantages over other CW techniques of interstitial optical property determination [11, 12, 16]. Firstly, for the P3_{in}-based radiance technique, neither absolute measurement nor a *priori* information is required, making the technique easier to be implemented in interstitial probe arrangements for clinical applications. Secondly, the developed method that only single wavelength measurements are employed, and resultantly both characterization and treatment of target tissues can be conducted by one laser. Thirdly, the shape-based fitting algorithms require large number of measurement data, whereas the developed recovery algorithm just requires the radiance measurements at two SDSs and four angles.

In theory, CW radiance measurements at two SDSs and four angles are sufficient to determine the optical properties of the homogeneous medium using the developed recovery algorithm based on the P3_{in}. However, for the inhomogeneous tissue sample in the clinical environments, more radiance measurement data are required to obtain a reliable result of optical property recovery. Similar to the arrangement of PDT for prostate [23], the minimally invasive measurements for tissue sample in clinical application can be conducted by a source catheter and two detector catheters which can be inserted into the tissue at a prescribed spacing using guidance templates. In order to minimize the invasiveness during the interstitial measurements, radiance measurements are still carried out at two fixed SDSs. Instead, radiance measurements at multiple angles can be obtained by rotating the detector fibers. The detecting angles can be chosen from $0\xb0$ to $160\xb0$by the step of $5\xb0$. For optical property recovery, CW radiance measurements at two SDSs and at four angles are employed to extract a set of optical properties where the four angles are randomly selected from four angle ranges of $0\xb0~40\xb0$, $41\xb0~80\xb0$, $81\xb0~120\xb0$ and $121\xb0~160\xb0$, separately. The process of optical property recovery can be repeated many times for statistical analysis.

## 5. Conclusion

In this paper, we have proposed the incomplete P3 approximation, the P3_{in} for short, which is the asymptotic part of the solution for the radiance by the P3 approximation. The proposed P3_{in} provides the simpler formulation than the P3 for recovering the optical properties and the higher degree approximation than the DA for modeling light propagation in the medium with a low albedo. Based on the P3_{in}, we have developed an analytical formulation method for extracting the optical properties directly from the radiance measurements at only two SDSs and four angles. The developed algorithm can recover ${\mu}_{a}$ and ${{\mu}^{\prime}}_{s}$ of the turbid media with albedos as low as 0.5 accurately with the REs within 10.82% and 10.67%, respectively, which demonstrates the effectiveness of the developed recovery algorithm based on the P3_{in} for recovering the optical properties of the turbid medium with a high absorption or a low scattering.

The ultimate goal of our work is to interstitially determine the optical properties of internal organs such as the prostate. Patient invasiveness limits the clinically available number of sources and detectors employed during interstitial laser treatments. Hence, it is interesting to explore whether ${\mu}_{a}$ and ${{\mu}^{\prime}}_{s}$ can be extracted, separately from relative radiance measurements at a single small SDS in the future. Thus, a source fiber and a detector fiber can be encapsulated in a probe. By replacing the laser with a broadband white light source and replacing the PMT with a spectrometer, CW interstitial radiance spectroscopy measurements can be carried out using the probe to determine the chromophore concentrations and the reduced scattering coefficient spectrum of tissue.

## Funding

The authors acknowledge the funding supports from the National Natural Science Foundation of China (81401453, 81371602, 61475115, 61475116, 61575140, 81571723, 81671728); Tianjin Municipal Government of China (15JCZDJC31800, 16JCZDJC31200, 17JCZDJC32700,17JCQNJC12700); the 111 Project (B07014).

## Acknowledgments

The authors acknowledge the insightful suggestions and revisions given by Prof. Yukio Yamada from University of Electro-Communications, Japan.

## References and links

**1. **G. Shafirstein, D. Bellnier, E. Oakley, S. Hamilton, M. Potasek, K. Beeson, and E. Parilov, “Interstitial photodynamic therapy-a focused review,” Cancers (Basel) **9**(2), 12 (2017). [PubMed]

**2. **T. M. Baran, “Recovery of optical properties using interstitial cylindrical diffusers as source and detector fibers,” J. Biomed. Opt. **21**(7), 77001 (2016). [PubMed]

**3. **V. K. Nagarajan and B. Yu, “Monitoring of tissue optical properties during thermal coagulation of *ex vivo* tissues,” Lasers Surg. Med. **48**(7), 686–694 (2016). [PubMed]

**4. **R. Nachabé, D. J. Evers, B. H. Hendriks, G. W. Lucassen, M. van der Voort, J. Wesseling, and T. J. Ruers, “Effect of bile absorption coefficients on the estimation of liver tissue optical properties and related implications in discriminating healthy and tumorous samples,” Biomed. Opt. Express **2**(3), 600–614 (2011). [PubMed]

**5. **Y. Pu, W. Wang, M. Al-Rubaiee, S. K. Gayen, and M. Xu, “Determination of optical coefficients and fractal dimensional parameters of cancerous and normal prostate tissues,” Appl. Spectrosc. **66**(7), 828–834 (2012). [PubMed]

**6. **R. Watté, B. Aernouts, R. Van Beers, and W. Saeys, “Robust metamodel-based inverse estimation of bulk optical properties of turbid media from spatially resolved diffuse reflectance measurements,” Opt. Express **23**(21), 27880–27898 (2015). [PubMed]

**7. **S.-Y. Tzeng, J.-Y. Guo, C.-C. Yang, C.-K. Hsu, H.-J. Huang, S.-J. Chou, C.-H. Hwang, and S.-H. Tseng, “Portable handheld diffuse reflectance spectroscopy system for clinical evaluation of skin: a pilot study in psoriasis patients,” Biomed. Opt. Express **7**(2), 616–628 (2016). [PubMed]

**8. **P. R. Bargo, S. A. Prahl, T. T. Goodell, R. A. Sleven, G. Koval, G. Blair, and S. L. Jacques, “In vivo determination of optical properties of normal and tumor tissue with white light reflectance and an empirical light transport model during endoscopy,” J. Biomed. Opt. **10**(3), 034018 (2005). [PubMed]

**9. **T. M. Baran, M. C. Fenn, and T. H. Foster, “Determination of optical properties by interstitial white light spectroscopy using a custom fiber optic probe,” J. Biomed. Opt. **18**(10), 107007 (2013). [PubMed]

**10. **T. M. Baran, J. D. Wilson, S. Mitra, J. L. Yao, E. M. Messing, D. L. Waldman, and T. H. Foster, “Optical property measurements establish the feasibility of photodynamic therapy as a minimally invasive intervention for tumors of the kidney,” J. Biomed. Opt. **17**(9), 98002 (2012). [PubMed]

**11. **A. Dimofte, J. C. Finlay, and T. C. Zhu, “A method for determination of the absorption and scattering properties interstitially in turbid media,” Phys. Med. Biol. **50**(10), 2291–2311 (2005). [PubMed]

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

**13. **S. Grabtchak and W. M. Whelan, “Separation of absorption and scattering properties of turbid media using relative spectrally resolved cw radiance measurements,” Biomed. Opt. Express **3**(10), 2371–2380 (2012). [PubMed]

**14. **D. J. Dickey, R. B. Moore, D. C. Rayner, and J. Tulip, “Light dosimetry using the P3 approximation,” Phys. Med. Biol. **46**(9), 2359–2370 (2001). [PubMed]

**15. **L. C. L. Chin, W. M. Whelan, and I. A. Vitkin, “Information content of point radiance measurements in turbid media: implications for interstitial optical property quantification,” Appl. Opt. **45**(9), 2101–2114 (2006). [PubMed]

**16. **L. C. L. Chin, A. E. Worthington, W. M. Whelan, and I. A. Vitkin, “Determination of the optical properties of turbid media using relative interstitial radiance measurements: Monte Carlo study, experimental validation, and sensitivity analysis,” J. Biomed. Opt. **12**(6), 064027 (2007). [PubMed]

**17. **S. Grabtchak, L. G. Montgomery, and W. M. Whelan, “Feasibility of interstitial near-infrared radiance spectroscopy platform for *ex vivo* canine prostate studies: optical properties extraction, hemoglobin and water concentration, and gold nanoparticles detection,” J. Biomed. Opt. **19**(5), 057003 (2014). [PubMed]

**18. **D. Piao, K. E. Bartels, Z. Jiang, G. R. Holyoak, J. W. Ritchey, G. Xu, C. F. Bunting, and G. Slobodov, “Alternative transrectal prostate imaging: a diffuse optical tomography method,” IEEE. J. Sel. Top. Quant. **16**(4), 715–729 (2010).

**19. **T. Svensson, E. Alerstam, M. Einarsdóttír, K. Svanberg, and S. Andersson-Engels, “Towards accurate *in vivo* spectroscopy of the human prostate,” J. Biophotonics **1**(3), 200–203 (2008). [PubMed]

**20. **E. L. Hull and T. H. Foster, “Steady-state reflectance spectroscopy in the P3 approximation,” J. Opt. Soc. Am. A **18**(3), 584–599 (2001).

**21. **A. Liemert and A. Kienle, “Analytical Green’s function of the radiative transfer radiance for the infinite medium,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. **83**(3 Pt 2), 036605 (2011). [PubMed]

**22. **T. Khan and A. Thomas, “Comparison of PN or spherical harmonics approximation for scattering media with spatially varying and spatially constant refractive indices,” Opt. Commun. **255**(1), 130–166 (2005).

**23. **J. Li and T. C. Zhu, “Determination of *in vivo* light fluence distribution in a heterogeneous prostate during photodynamic therapy,” Phys. Med. Biol. **53**(8), 2103–2114 (2008). [PubMed]

**24. **S. Grabtchak, T. J. Palmer, F. Foschum, A. Liemert, A. Kienle, and W. M. Whelan, “Experimental spectro-angular mapping of light distribution in turbid media,” J. Biomed. Opt. **17**(6), 067007 (2012). [PubMed]

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

**26. **O. Barajas, Å. M. Ballangrud, G. G. Miller, R. B. Moore, and J. Tulip, “Monte Carlo modelling of angular radiance in tissue phantoms and human prostate: PDT light dosimetry,” Phys. Med. Biol. **42**(9), 1675–1687 (1997). [PubMed]

**27. **H. J. van Staveren, C. J. Moes, J. van Marie, S. A. Prahl, and M. J. van Gemert, “Light scattering in Intralipid-10% in the wavelength range of 400-1100 nm,” Appl. Opt. **30**(31), 4507–4514 (1991). [PubMed]