1.

Introduction

The in vivo determination of tissue optical properties (the absorption coefficient ${\mu }_a$ and the reduced scattering coefficient ${\mu }_s^{\prime }$) and the study of light propagation in the biological tissues are very important in a variety of biomedical fields because of their potential applications in photodynamic therapy,^1,2 glucose monitoring,^3,4 diffuse optical tomography,^5,6 and cancer detection.^7,8

Among the optical measurement methods, steady-state or continuous wave spatially resolved diffuse reflectance measurement has been extensively adopted for quantitative tissue characterization and disease diagnosis.^9–12 Estimations of optical properties based on spatially resolved diffuse reflectance measurements, which is also called the inverse problem, can be cataloged into two types according to the dependence of the knowledge of the light propagation in biological tissues. One type of the inverse problem, which needs a priori knowledge including partial least square (PLS),¹³ neural networks (NN)^14,15 support vector machines (SVM),¹⁶ and look-up table (LUT) approach.^17,18 However, PLS does not provide accurate predictions, as it is unable to account for the intrinsic nonlinear relationships in light transport problems.¹⁹ NN requires a number of accurate correlations between optical properties and spatially resolved diffuse reflectance measurements to develop training sets that can initially “teach the neurons”.¹⁵ LUT is based solely on experimental measurements of calibration standards to provide high accuracy in optical property determination even at high absorption levels using a nonlinear optimization routine.^17,18 The new SVM-based regression algorithm for determining the absorption and reduced scattering coefficients can be implemented rapidly and provide more accurate predictions than LUT.¹⁶ Though the methods of NN, LUT, and SVM can obtain good prediction accuracy, they all need high quality calibration sets, i.e., requiring a lot of premodeling work.

Another type of inverse problem is based on the light transportation model called the forward model. The most commonly adopted method within this form of the inverse problem is the nonlinear fitting methods based on either an analytical approximation of the radiative transfer equation (RTE) such as diffusion approximation (DA)²⁰ or the Monte Carlo model.²¹ Thus, the accuracy of the nonlinear fitting methods in determining optical properties is highly influenced by the veracity of the analytical and numerical models. DA, as a low-order spherical harmonics to RTE (PN approximation), is widely used in the biomedical field. However, DA has several limitations and is unsuitable for application in tissues where absorption is high in comparison with scattering and the SDSs are small, e.g., the application of noninvasive blood-glucose monitoring using light in the region of 1000 to 2000 nm.^22,23 To overcome these limitations of DA, Hull and Foster²⁴ derived the Green’s function of the steady-state RTE in the P3 approximation and demonstrated that the P3 approximation models the radiance in highly absorbing media or in the region close to the source with an accuracy superior to that of DA. To improve the DA model and to simplify the P3 approximation for highly absorbing media, a modified diffusion approximation model called the hybrid diffusion-P3 approximation (HDPA) was further investigated by Hull and Foster.²⁴ Tian et al.²⁵ used the HDPA model to derivate the optical properties.

Based on the previously mentioned studies, the derivation of optical properties based on DA, P3, and HPDA have been investigated. However, the inverse calculation of optical properties for large reduced albedo $a^{\prime }[a^{\prime }={\mu }_s^{\prime }/({\mu }_s^{\prime }+{\mu }_a)]$ range and small SDSs is still a difficult problem. The aim of this paper was to propose a joint derivation method to accurately determine the ${\mu }_a$ and ${\mu }_s^{\prime }$ at large $a^{\prime }$ range from 0.5 to 0.99 and small SDSs using DA, P3, and HPDA. To apply the joint derivation method, we first used the best forward model of ${\mu }_a$ to derive ${\mu }_a$ and ${\mu }_s^{\prime }$ simultaneously, and then used the derived ${\mu }_a$ as the best result of ${\mu }_a$. Secondly, we used the best forward model of ${\mu }_s^{\prime }$ to derive ${\mu }_a$ and ${\mu }_s^{\prime }$ simultaneously, and then used the derived ${\mu }_s^{\prime }$ as the best result of ${\mu }_s^{\prime }$. In Sec. 2, the aforementioned forward models of both one-point and two-point sources, and the inverse calculation method are introduced. In Secs. 3 and 4, the applicable regions of the forward models are investigated and a joint derivation method is proposed. The corresponding errors of the joint derivation method in deriving optical properties are subsequently given for the five ranges of $a^{\prime }$, which can be decided based on the steady-state spatially resolved diffuse reflectance measurement.

2.

Forward Models and the Inverse Calculation of Optical Properties

2.1.

Forward Models

The migration of photons in biological tissues can be described by the RTE, while disregarding any wave effect:

Eq. (1)

$$\widehat{s}·\nabla L(r,\widehat{s})=-({\mu }_a+{\mu }_s)L(r,\widehat{s})+{\mu }_s{\int }_{4\pi }p(\widehat{s}·{\widehat{s}}^{\prime })L(r,\widehat{s})\mathrm{d}{\widehat{s}}^{\prime }+S(r,\widehat{s}),$$

where $L(r,\widehat{s})$ is the radiance at position $r$, and in the direction $\widehat{s}$. $S(r,\widehat{s})$ describes the distribution of the source. $p(\widehat{s}·{\widehat{s}}^{\prime })$ is the scattering phase function for photons anisotropically scattering from direction ${\widehat{s}}^{\prime }$ to direction $\widehat{s}$ using the Henyey–Greenstein phase function.²⁶ ${\mu }_a$ and ${\mu }_s$ are the absorption coefficient and the scattering coefficient, respectively.

The spherical harmonics equation, which is one of the derivation methods of RTE, is commonly used to provide an approximation for the RTE. DA is one of the most widely used models called the P1 approximation and the Green’s function solution of DA, at the situation of an isotropic steady-state source an infinite medium, is given by:

Eq. (2)

$${\mathrm{\varphi }}_{\mathrm{DA}}(r)=\frac{1}{4\pi D}\frac{\exp (-{\mu }_{\mathrm{eff}}r)}{r},$$

where $\mathrm{\varphi }={\int }_{4\pi }L(r,\widehat{s})\mathrm{d}\widehat{s}$ is the fluence at position $r$, ${\mu }_{\mathrm{eff}}={[3{\mu }_a({\mu }_a+{\mu }_s^{\prime })]}^{1/2}$, and $D={[3({\mu }_a+{\mu }_s^{\prime })]}^{-1}$ are the effective attenuation coefficient and the diffusion coefficient, respectively. ${\mu }_s^{\prime }={\mu }_s(1-g)$ is the reduced scattering coefficient with $g$ being the anisotropic factor.

P3 approximation is the three-order approximation of the RTE using the spherical harmonics method. HDPA, the Green’s function of which had a good agreement with P3, can be obtained using the Green’s function of DA, but replacing ${\mu }_{\mathrm{eff}}$ and $D$ with $v^{-}$ and $D_{\mathrm{asym}}$, respectively.²⁴ $v^{-}$ and $D_{\mathrm{asym}}$ are used as the coefficients of P3 approximation. The Green’s function of HDPA, corresponding to Eq. (2) is described as follows:

Eq. (3)

$${\mathrm{\varphi }}_{\mathrm{HDPA}}(r)=\frac{1}{4\pi D_{\mathrm{asym}}}\frac{\exp (-v^{-}r)}{r},$$

where

Eq. (4)

$$\begin{gathered} v^{-}=\frac{1}{\sqrt{18}}{(\beta -\sqrt{{\beta }^2-{\gamma }_a})}^{\frac{1}{2}} \\ \beta \equiv 27{\mu }_a({\mu }_a+{\mu }_s^{\prime })+28{\mu }_a({\mu }_a+{\mu }_s^{\prime }\delta )+35{\mu }_a({\mu }_a+{\mu }_s^{\prime }\gamma ){\mu }_a({\mu }_a+{\mu }_s^{\prime }\delta ) \\ {\gamma }_a\equiv 3780({\mu }_a+{\mu }_s^{\prime })({\mu }_a+{\mu }_s^{\prime }\gamma )({\mu }_a+{\mu }_s^{\prime }\delta ) \\ \gamma =(1-g_2)/(1-g_1),\delta =(1-g_3)/(1-g_1) \\ {D}_{\mathrm{asym}}={\mu }_a/{(v^{-})}^2. \end{gathered}$$

In Eq. (4), $g_1$, $g_2$, and $g_3$ are the first-moment, second-moment, and third-moment of the phase function, respectively. For the Henyey–Greenstein phase function, $\gamma =1+g$ and $\delta =1+g+g^2$.

For a semi-infinite medium, the derivation of DA or HDPA involves characterizing the source terms and satisfying appropriate boundary conditions. To represent the pencil beam in terms of simpler source distributions, isotropic point sources have been developed to describe the distribution of sources. The distributions of point sources can be approximated to the dipole moment (one-point source) or to the dipole and quadruple moments (two-point sources) as illustrated in Fig. 1(a) and 1(b), respectively.^24,25

Fig. 1

Approximation for point sources and the extrapolated boundary condition (a) one-point source and (b) two-point sources.

To obtain the solutions of the radiance emitted from a semi-infinite medium, appropriate boundary conditions must be prescribed at the interface between the surrounding media and the biological tissue. The extrapolated boundary condition (EBC) is one of the boundary conditions used for a semi-infinite scattering medium,^27,28 shown in Fig. 1. To satisfy the extrapolated boundary condition, a negative image source is introduced above the tissue surface with the same magnitude to the corresponding real sources at $z^{\prime }=-(2z_b+z_0)$ for one-point source or $z_1=-(2z_b+z_{01})$ and $z_2=-(2z_b+z_{02})$ for two-point sources, respectively, shown in Fig. 1. Thus, for the semi-infinite medium in EBC, the source terms and the sources $q_1(z)$ (one-point source case) and $q_2(z)$ (two-point sources case), are the combination of the real and image sources, which can be expressed as:²⁵

Eq. (5)

$$\begin{gathered} q_1(z)=a^{\prime }[\delta (\rho ,z-z_0)-\delta (\rho ,z+2z_b+z_0)] \\ {q}_2(z)=\frac{1}{2}{a}^{\prime }[\delta (\rho ,z-z_{01})-\delta (\rho ,z+2z_b+z_{01})] \\ +\frac{1}{2}{a}^{\prime }[\delta (\rho ,z-z_{02})-\delta (\rho ,z+2z_b+z_{02})], \end{gathered}$$

where $z_0=1/{\mu }_t^{\prime }$, $z_{01}=2/{\mu }_t^{\prime }$, and $z_{02}=0$. $z_b$ is the distance between the real and extrapolated boundaries, which corresponds to 2AD and $2{\mathrm{AD}}_{\mathrm{asym}}$ for DA and HDPA, respectively. $A=(1+R_{\mathrm{eff}})/(1-R_{\mathrm{eff}})$, and $R_{\mathrm{eff}}$ is the effective reflection coefficient.

The radiance fluence calculated from the DA and the HDPA with EBC can be expressed as:²⁹

Eq. (6)

$$\begin{gathered} {\varphi }_{\mathrm{DA}-\mathrm{EBC}}(\rho ,z)={\mathrm{\varphi }}_{\mathrm{DA}}(\rho ,z)\otimes q(z) \\ {\varphi }_{\mathrm{HDPA}-\mathrm{EBC}}(\rho ,z)={\mathrm{\varphi }}_{\mathrm{HDPA}}(\rho ,z)\otimes q(z), \end{gathered}$$

where $\rho$ is the distance from the source to the detector.

The reflectance can then be given by:²⁷

Eq. (7)

$$\begin{gathered} R_{\mathrm{DA}-\mathrm{EBC}}(\rho )={\int }_{2\pi }\mathrm{d}\mathrm{\Omega }[1-R_{\mathrm{fres}}(\theta )]\frac{1}{4\pi }[{\varphi }_{\mathrm{DA}}(\rho ,z=0) \\ +3D\frac{\partial {\varphi }_{\mathrm{DA}}(\rho ,z=0)}{\partial z}\cos \theta ]\cos \theta \\ {R}_{\mathrm{HDPA}-\mathrm{EBC}}(\rho )={\int }_{2\pi }d\mathrm{\Omega }[1-R_{\mathrm{fres}}(\theta )]\frac{1}{4\pi }[{\varphi }_{\mathrm{HDPA}}(\rho ,z=0) \\ +3D_{\mathrm{asym}}\frac{\partial {\varphi }_{\mathrm{HDPA}}(\rho ,z=0)}{\partial z}\cos \theta ]\cos \theta , \end{gathered}$$

where $R_{\mathrm{fres}}(\theta )$ is the Fresnel reflection coefficient for a photon with an incident angle $\theta$ relative to the normal to the boundary.

3.

Comparison of Six Theoretical Models for Determining Optical Properties

In this section, we evaluate the results from the same inverse algorithm as described in Sec. 2.2, but used six forward models (DA, HDPA, and P3 in both one and two-point sources). The range of $a^{\prime }$ is from 0.5 to 0.99, and the anisotropic factor $g$ and the refractive index $n$ are assumed to be 0.9 and 1.4, respectively. The range of SDSs is between 0.4 and 8 mm.

In the MC simulation, a pencil photon beam is assumed to normally irradiate upon the semi-infinite turbid medium. The Henyey–Greenstein phase function is applied to calculate the scattering angle. The spatial resolution of the steady-state MC simulations is 0.1 mm, and distances up to 30 mm were scored. The number of the incidence photons was ${10}^8$.

3.1.

Best Forward Model for Different Reduced Albedo Ranges

In this section, the “experimental” spatially resolved diffuse reflectance was calculated from MC with the “real” ${\mu }_a$ being 0.01, 0.1, 0.5, and $1{\mathrm{mm}}^{-1}$, and ${\mu }_s$ being $10{\mathrm{mm}}^{-1}$, and $g$ being 0.9. For the spatially resolved measurements, multiple SDSs are required. As other researchers investigated, the aforementioned forward models had different applicable SDSs.²⁴ Thus, the evaluation of the forward models should appoint the SDS range. Here we define the smallest SDS in a measurement as ${\rho }_{\text{start}}$.

Figures 2 and 3 show the spatially resolved diffuse reflectance and compare the diffuse reflectance between MC simulation and the six forward models at different absorption coefficients (${\mu }_a=0.1{\mathrm{mm}}^{-1}$ and ${\mu }_a=0.5{\mathrm{mm}}^{-1}$).

Fig. 2

The value (a) and the variety (b) of diffuse reflectance $R$ at the detection range of 0.4∶0.1∶8 mm from the source. The “real” optical properties are ${\mu }_a=0.1{\mathrm{mm}}^{-1}$, ${\mu }_s^{\prime }=1{\mathrm{mm}}^{-1}$, $g=0.9$, $n=1.4$.

Fig. 3

The value (a) and the variety (b) of diffuse reflectance $R$ at the detection range of 0.4∶0.1∶8 mm from the source. The “real” optical properties are ${\mu }_a=0.5{\mathrm{mm}}^{-1}$, ${\mu }_s^{\prime }=1{\mathrm{mm}}^{-1}$, $g=0.9$, $n=1.4$.

Figures 4 and 5 show the derived ${\mu }_a$ and ${\mu }_s^{\prime }$ at different ${\rho }_{\text{start}}$, which varies from 0.4 to 4 mm at an interval of 0.1 mm. In the figures, SDSs are from ${\rho }_{\text{start}}$ to 8 mm at the spatial resolution of 0.1 mm. Figures 4 and 5 indicate that the derivation accuracy of ${\mu }_a$ and ${\mu }_s^{\prime }$ depends on the forward models used. When $a^{\prime }$ is 0.99 (${\mu }_a=0.1{\mathrm{mm}}^{-1}$, ${\mu }_s^{\prime }=1{\mathrm{mm}}^{-1}$), as shown in Fig. 4(a), the best models to derive ${\mu }_a$ are DA and HDPA of two-point sources. The derivation error is less than 10% at the ${\rho }_{\text{start}}$ range from 0.6 to 4 mm, and is even less than 1% at the ${\rho }_{\text{start}}$ range from 1.6 to 3.5 mm. However, in this $a^{\prime }$ range, the best models to derive ${\mu }_s^{\prime }$ are DA and HDPA of one-point source. The derivation error is less than 7% at the ${\rho }_{\text{start}}$ range from 0.4 to 4 mm. When $a^{\prime }$ equals 0.91 (${\mu }_a=0.1{\mathrm{mm}}^{-1}$, ${\mu }_s^{\prime }=1{\mathrm{mm}}^{-1}$), as shown in Fig. 4(b), the best model to derive ${\mu }_a$ and ${\mu }_s^{\prime }$ is HDPA of one-point source. The derivation errors of ${\mu }_a$ and ${\mu }_s^{\prime }$ are less than 10% at the ${\rho }_{\text{start}}$ range from 0.7 to 4 mm and the ${\rho }_{\text{start}}$ range from 1.7 to 4 mm, respectively.

Fig. 4

Absorption and reduced scattering coefficients derived from nonlinear regressions of DA, HDPA, and P3 for the steady-state spatially resolved reflectance with the extrapolated-boundary condition at the detection range of ${\rho }_{\text{start}}:0.1:8\mathrm{mm}$ from the source. The “real” ${\mu }_a$ and ${\mu }_s^{\prime }$ are (a) ${\mu }_a=0.01{\mathrm{mm}}^{-1}$, ${\mu }_s^{\prime }=1{\mathrm{mm}}^{-1}$ and (b) ${\mu }_a=0.1{\mathrm{mm}}^{-1}$, ${\mu }_s^{\prime }=1{\mathrm{mm}}^{-1}$.

Fig. 5

Absorption and reduced scattering coefficients derived from nonlinear regressions of DA, HDPA, and P3 for the steady-state spatially resolved reflectance with the extrapolated-boundary condition at the detection range of ${\rho }_{\text{start}}:0.1:8\mathrm{mm}$ from the source. The “real” ${\mu }_a$ and ${\mu }_s^{\prime }$ are (a) ${\mu }_a=0.5{\mathrm{mm}}^{-1}$, ${\mu }_s^{\prime }=1{\mathrm{mm}}^{-1}$ and (b) ${\mu }_a=1{\mathrm{mm}}^{-1}$, ${\mu }_s^{\prime }=1{\mathrm{mm}}^{-1}$.

Figure 5(a) and 5(b) shows similar results to Fig. 4 but with decreased values of $a^{\prime }$ at 0.67 (${\mu }_a=0.5{\mathrm{mm}}^{-1}$, ${\mu }_s^{\prime }=1{\mathrm{mm}}^{-1}$) and 0.5 (${\mu }_a=1{\mathrm{mm}}^{-1}$, ${\mu }_s^{\prime }=1{\mathrm{mm}}^{-1}$), respectively. From Fig. 5(a), where ${\mu }_a$ is equal to $0.5{\mathrm{mm}}^{-1}$ and with other optical properties kept constant from Fig. 4, it can be found that HDPA of two-point sources is more effective than other models in deriving ${\mu }_a$ and ${\mu }_s^{\prime }$. The relative error in deriving ${\mu }_a$ is between 0 and 5% at the ${\rho }_{\text{start}}$ range from 0.4 to 2.6 mm. The derivation error about ${\mu }_s^{\prime }$ is between 0 and 5% at the ${\rho }_{\text{start}}$ range from 0.5 to 2.6 mm. From Fig. 5(b), where ${\mu }_a=1{\mathrm{mm}}^{-1}$, it can be found that the one source DA model for deriving the absorption coefficient can obtain good results. The derivation error from this model is less than 6% at the ${\rho }_{\text{start}}$ range from 0.7 to 2.4 mm. For the derivation of ${\mu }_s^{\prime }$, the HDPA model of two sources is better than other models. However, the smallest derivation error is still near 20% at the ${\rho }_{\text{start}}$ range from 0.4 to 1.2 mm.

According to this analysis, these theoretical models have different applicability for different $a^{\prime }$ ranges. For a more detailed study of the applicability of these models, 19 sets of optical properties were selected, which all have the same ${\mu }_s^{\prime }$ ($1{\mathrm{mm}}^{-1}$) but different ${\mu }_a$. For ${\mu }_a$ equal to 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, and $1{\mathrm{mm}}^{-1}$ respectively, we can obtain the reduced albedos ranging from 0.5 to 0.99. By comparing the derivation error, the best model to derive ${\mu }_a$ and ${\mu }_s^{\prime }$ at different ranges of $a^{\prime }$ using the six theoretical models is summarized in Table 1. First, this means that the best forward model of ${\mu }_a$ is used to derive ${\mu }_a$ and ${\mu }_s^{\prime }$ simultaneously and then the derivation result of ${\mu }_a$ is used to be the best result of ${\mu }_a$; secondly, the best forward model of ${\mu }_s^{\prime }$ is used to derive ${\mu }_a$ and ${\mu }_s^{\prime }$ simultaneously, and then the derivation result of ${\mu }_s^{\prime }$ is used to be the best result of ${\mu }_s^{\prime }$.

Table 1

The best model to derive the absorption coefficient and reduced scattering coefficient at different ranges of the reduced albedo.

	0.96<=a′<0.99	0.91<=a′<0.96	0.71<=a′<0.91	0.62<=a′<0.71	0.5<=a′<0.62
${\mu }_a$	HDPA, two sources	HDPA, one source	DA, one source	HDPA, two sources	DA, one source
${\mu }_s^{\prime }$	HDPA, one source	HDPA, one source	DA, two sources	HDPA, two sources	HDPA, two sources

3.2.

Best ${\rho }_{\text{start}}$ Range for Different Best Forward Model

As the discussion in Sec. 3.1 reveals, the forward models also have different applicable ${\rho }_{\text{start}}$ ranges that will yield the small derivation error. Thus in this section, the derivation error for the best forward model under different ${\rho }_{\text{start}}$ is investigated.

For the best models listed in Table 1, Fig. 6 shows the derivation error of the aforementioned 19 sets of optical properties at different ${\rho }_{\text{start}}$ range from ${\rho }_{\text{start}}$ to 8 mm with the spatial resolution of 0.1 mm. In Fig. 6, the 19 sets of “real” ${\mu }_a$ and ${\mu }_s^{\prime }$ are shown in the abscissa axis. The hollow column represents the ${\rho }_{\text{start}}$ range of ${\mu }_a$ and the solid column represents the ${\rho }_{\text{start}}$ range of ${\mu }_s^{\prime }$. The corresponding derivation errors of the 19 sets of optical properties are indicated on the top of each column. From Fig. 6, it can be found that the best ${\rho }_{\text{start}}$ range, at which the derivation error is small for the best forward model has no obvious regularity. From Figs. 4 and 5, or other 19 sets of optical properties, we observed that all of the curves varied smoothly during the best ${\rho }_{\text{start}}$ range. Thus, an effective way to find the best ${\rho }_{\text{start}}$ range for a best forward model is to select the smooth range of the curve in Fig. 4 or 5, at which ${\rho }_{\text{start}}$ varying from 0.4 to 4 mm.

Fig. 6

The ${\rho }_{\text{start}}$ (the ordinate) and the corresponding derivation error range (the upper abscissa) for the 19 sets of optical properties (the lower abscissa) using the best model listed in Table 1. For all the 19sets of optical properties, the real reduced scattering coefficient is $1{\mathrm{mm}}^{-1}$. The hollow column represents the ${\rho }_{\text{start}}$ range of ${\mu }_a$ and the solid column represents the ${\rho }_{\text{start}}$ range of ${\mu }_s^{\prime }$.

3.3.

Comparison Between DA One Source and Best Forward Model

Because the DA one source model is widely used, it is necessary to compare the results between the DA one source model and the best model in Table 1 under their respective best ${\rho }_{\text{start}}$ ranges. Figures 7 and 8 give the comparison results in view of the derivation of ${\mu }_a$ and ${\mu }_s^{\prime }$, respectively. The hollow column represents the derivation error range using the DA one source model under its best ${\rho }_{\text{start}}$ range and the solid column represents those using the corresponding best models shown in Table 1, which are also indicated on top of the figures. From the figures, we find that the derivation error from the best model is less than the DA one source model. Thus, we conclude that it is helpful to use the corresponding best model to derive the optical properties at different reduced albedo ranges, especially at small SDSs.

Fig. 7

Comparison of the errors for the absorption coefficient derivation calculated from the DA one source model and the best model, respectively. The “real” ${\mu }_a$ is showing in the abscissa axis, and ${\mu }_s^{\prime }=1{\mathrm{mm}}^{-1}$. The hollow column represents the results using the DA one source model and the solid column represents the results using the best model shown in Table 1.

Fig. 8

Comparison of the errors for the reduced scattering coefficient derivation calculated from the DA one source model and the best model, respectively. The “real” ${\mu }_a$ is showing in the abscissa axis, and ${\mu }_s^{\prime }=1{\mathrm{mm}}^{-1}$. The hollow column represents the results using the DA one source model and the solid column represents the results using the best model showing in Table 1.

4.

Joint Derivation Method for Determining Optical Properties in Large Reduced Albedo Range

In addition to the case of ${\mu }_s^{\prime }=1{\mathrm{mm}}^{-1}$, we also studied other cases of the best forward model at ${\mu }_s^{\prime }=1.3{\mathrm{mm}}^{-1}$, ${\mu }_s^{\prime }=1.6{\mathrm{mm}}^{-1}$, and ${\mu }_s^{\prime }=2{\mathrm{mm}}^{-1}$ ($0.5\le a^{\prime }\le 0.99$). In the case of ${\mu }_s^{\prime }=2{\mathrm{mm}}^{-1}$, results show that the best forward model to determine the optical properties at different ranges of the reduced albedo is the same as that of ${\mu }_s^{\prime }=1{\mathrm{mm}}^{-1}$, except in the region of $0.95\le a^{\prime }\le 0.98$. When $0.95\le a^{\prime }\le 0.98$ and ${\mu }_s^{\prime }=2{\mathrm{mm}}^{-1}$, results indicated that the best forward model to derive the absorption coefficient is the P3 approximation of one-point source. Figure 9 shows one of the results of $a^{\prime }=0.975$, ${\mu }_a=0.05{\mathrm{mm}}^{-1}$, and ${\mu }_s^{\prime }=2{\mathrm{mm}}^{-1}$. Specific results from the other cases are not given; however, it is important to note that in the case of $0.95\le a^{\prime }\le 0.98$, if we adopt the best forward model in Table 1, which is the HDPA of two sources, the increasing derivation error of the absorption coefficient is still less than 10%. The situation of ${\mu }_s^{\prime }=1.3{\mathrm{mm}}^{-1}$ and ${\mu }_s^{\prime }=1.6{\mathrm{mm}}^{-1}$ is similar with the case of ${\mu }_s^{\prime }=2{\mathrm{mm}}^{-1}$. All results suggest that the best forward models listed in Table 1 are universally applicable.

Fig. 9

Absorption and reduced scattering coefficients derived from nonlinear regressions of DA, HDPA, and P3 for the steady-state spatially resolved reflectance with the extrapolated-boundary condition at the detection range of ${\rho }_{\text{start}}:0.1:8\mathrm{mm}$ from the source. The “real” ${\mu }_a$ and ${\mu }_s^{\prime }$ are ${\mu }_a=0.05{\mathrm{mm}}^{-1}$, and ${\mu }_s^{\prime }=2{\mathrm{mm}}^{-1}$.

Because different $a^{\prime }$ ranges have different corresponding best forward models, it is difficult to generate a universal derivation method with only one forward model. Based on these observations, a joint derivation method is proposed. The joint derivation method uses the best forward models in varying ranges of $a^{\prime }$ to derive the ${\mu }_a$ and ${\mu }_s^{\prime }$ simultaneous, and then uses the derivation result of their own best forward model as the best result of ${\mu }_a$ and ${\mu }_s^{\prime }$, respectively, based on the determination of the $a^{\prime }$ range and the best forward models being chosen. The scheme to determine and apply the joint derivation method is shown in Fig. 10. First, the range of $a^{\prime }$ is determined from the slope $k$ and the intercept $d_{x=1}$ of the ${\mathrm{Log}}_{10}(R_{\mathrm{meas}})$ with the ($R_{\mathrm{meas}}$) being the diffuse reflectance either measured from experiment or calculated from simulation. Then the corresponding best forward model, as shown in Sec. 3.1, is determined. Next, an initial set of the input parameters (${\mu }_a$ and ${\mu }_s^{\prime }$) is used in the best forward model and then the diffuse reflectance calculated from the forward model ($R_{\text{model}}$) is generated. The sum of square errors between the $R_{\text{model}}$ and $R_{\mathrm{meas}}$ is computed. The input parameters are then iteratively updated until the sum of square errors is minimized. The trust-region-reflective algorithm provided by MATLAB optimization toolbox is used as the optimization algorithm.

Fig. 10

The scheme of the joint derivation method.

It can be observed from Fig. 10 that the utilization of the joint derivation method is limited by the knowledge of the range of $a^{\prime }$. Because the optical properties are often unknown for a sample of tissue, it is necessary to investigate the possibility of the division of $a^{\prime }$ range through the diffuse reflection measurement. For $a^{\prime }$ range from 0.5 to 0.99 when ${\mu }_s^{\prime }$ equals $1{\mathrm{mm}}^{-1}$, the diffuse reflectance $R_{\mathrm{MC}}$ in different SDSs was calculated from MC. Figure 11 shows the first-order linear regression of ${\mathrm{Log}}_{10}(R_{\mathrm{MC}})$, that has a good linearity at each reduced albedo to SDS. The slope $k$ and the intercept $d_{x=1}$ of each line decreases as the reduced albedo decreases, i.e., different lines have different $k$ and $d_{x=1}$ values. Figure 11 indicates that the $k$ and $d_{x=1}$ may be the values that can be used initially to determine the range of $a^{\prime }$. To achieve this, the $k$ and $d_{x=1}$ at different $a^{\prime }$ ranges should not overlap at the same time.

Fig. 11

First-order linear regression of ${\mathrm{Log}}_{10}(R_{MC})$ in the case of $0.5\le a^{\prime }\le 0.99$ and ${\mu }_s^{\prime }=1{\mathrm{mm}}^{-1}$.

To validate the prediction of $a^{\prime }$ with Fig. 11, the medium with different ${\mu }_s^{\prime }$ (from 1 to $2{\mathrm{mm}}^{-1}$ at the interval of $0.1{\mathrm{mm}}^{-1}$) were investigated. Figure 12 shows the calculated $k$ and $d_{x=1}$ when $a^{\prime }$ equals 0.99, 0.96, 0.91, 0.71, 0.62, and 0.5. It can be imaged that the $k$ and $d_{x=1}$ of any $a^{\prime }$ range listed in Table 1 should be located between adjacent lines in Fig. 12. The value of $k$ and $d_{x=1}$ can be obtained from the measured or simulated diffuse reflectance ($R_{\mathrm{meas}}$) and generally the range of $a^{\prime }$ can be determined accordingly. However, sometimes, the overlapping may happen, which means that the two $a^{\prime }$ regions possessing the same value of $k$ may also be having the same value of $d_{x=1}$. To determine the range of $a^{\prime }$ with the $k$ and $d_{x=1}$, it is important to find out if the overlapping occurred. Through the analysis of the data in Fig. 12, it can be found that the simultaneous overlapping of $k$ and $d_{x=1}$ for $1{\mathrm{mm}}^{-1}\le {\mu }_{\mathrm{s}}^{\prime }\le 2{\mathrm{mm}}^{-1}$ happens only between the regions of $0.71\le a^{\prime }\le 0.91$ and $0.62\le a^{\prime }\le 0.71$, and between the regions of $0.62\le a^{\prime }\le 0.71$ and $0.5\le a^{\prime }\le 0.62$. Table 2 shows the overlapping region with this situation. We can also find that the overlapping ranges are adjacent. Thus, one can just choose one of the overlapping $a^{\prime }$ ranges as the desired one. Because the increase of the derivation error caused by the misusage of the adjacent forward models is less than 10%, the increase of the derivation error caused by this optional choice of $a^{\prime }$ range is also less than 10%.

Fig. 12

The slope $k$ and intercept $d_{x=1}$ of ${\mathrm{Log}}_{10}(R_{MC})$ in the case of $1{\mathrm{mm}}^{-1}\le {\mu }_s^{\prime }\le 2{\mathrm{mm}}^{-1}$.

Table 2

The values of k and dx=1 at different ranges of the reduced albedo and different reduced scattering coefficients.

Based on this discussion, a joint derivation method was proposed. To validate and apply the joint derivation method, the diffuse reflectance $R_{\mathrm{MC}}$ was generated by Monte Carlo simulation with the known true optical properties, and then comparisons were given between the best forward models determined by the proposed method and by the optical property derivation accuracy.

Firstly, the best models are determined with the method proposed in this paper. Figure 13 shows the first-order linear regression of ${\mathrm{Log}}_{10}(R_{\mathrm{MC}})$, where the “real” ${\mu }_a$ of sample 1 (recombination line) and sample 2 (dotted line) are 0.1 and $0.5{\mathrm{mm}}^{-1}$, respectively. The other optical properties of the two samples are ${\mu }_s^{\prime }=1.5{\mathrm{mm}}^{-1}$, $g=0.9$, and $n=1.4$. From Fig. 13, we obtain that the $k$ and $d_{x=1}$ of sample 1 are $-0.4616$ and 0.1697, respectively, whereas those of sample 2 are $-0.8463$ and $-0.2605$, respectively. From Fig. 12 and Table 3, it can be seen that the $k$ of sample 1 spans three regions of $a^{\prime }$, which are $0.71\le a^{\prime }\le 0.91$, $0.91\le a^{\prime }\le 0.96$, and $0.96\le a^{\prime }\le 0.99$. The $d_{x=1}$ of sample 1 spans two regions of a', which are $0.91\le a^{\prime }\le 0.96$, and $0.96\le a^{\prime }\le 0.99$. Then the desired $a^{\prime }$ range is the sharing region that the $k$ and $d_{x=1}$ share the same region of a' for the same ${\mu }_s^{\prime }$. From Table 3, it can be determined that the sharing region of sample 1 is $0.91\le a^{\prime }\le 0.96$, whereas that of sample 2 is $0.71\le a^{\prime }\le 0.91$. According to these regions of $a^{\prime }$, we can find from Table 1 that the best forward models for ${\mu }_a$ and ${\mu }_s^{\prime }$ of sample 1 are both the HDPA of one source, and the best forward models for ${\mu }_a$ and ${\mu }_s^{\prime }$ of sample 2 are the DA model of one source and the DA model of two sources, respectively.

Fig. 13

First-order linear regression of ${\mathrm{Log}}_{10}(R_{MC})$. The “real” ${\mu }_a$ of the sample 1 (recombination line) is $0.1{\mathrm{mm}}^{-1}$, and that of the sample 2 (dotted line) is $0.5{\mathrm{mm}}^{-1}$. The other optical properties of the two samples are ${\mu }_s^{\prime }=1.5{\mathrm{mm}}^{-1}$, $g=0.9$, and $n=1.4$.

Table 3

The a′ range determined from Fig. 12 for the k and dx=1 of sample 1 and sample 2.

Secondly, to validate the determined best models, the best models of sample 1 and sample 2 are directly decided by comparing the derivation accuracy of ${\mu }_a$ and ${\mu }_s^{\prime }$. Figure 14 shows the derived ${\mu }_a$ and ${\mu }_s^{\prime }$ of sample 1 and sample 2 with the method described in Sec. 2.2. We find that the best models determined using the $a^{\prime }$ range are the same as the results indicated in Fig. 14.

Fig. 14

The derived absorption and reduced scattering coefficients for (a) sample 1 and (b) sample 2 with the method described in Sec. 2.2. The “real” ${\mu }_a$ and ${\mu }_s^{\prime }$ of sample 1 are ${\mu }_a=0.1{\mathrm{mm}}^{-1}$, ${\mu }_s^{\prime }=1.5{\mathrm{mm}}^{-1}$ and these of sample 2 are ${\mu }_a=0.5{\mathrm{mm}}^{-1}$, ${\mu }_s^{\prime }=1.5{\mathrm{mm}}^{-1}$.

5.

Discussion and Conclusion

The accurate determination of the optical properties of tissues is very important in a variety of diagnostic and therapeutic procedures. In this paper, we proposed a joint derivation method for determining the optical properties accurately. Through the study of some cases of different ${\mu }_s^{\prime }$, we obtained a general conclusion and determined the best theoretical derivation model for different ranges of $a^{\prime }$. The best ${\rho }_{\text{start}}$ range and the corresponding derivation error range to derive ${\mu }_a$ and ${\mu }_s^{\prime }$ at small SDSs (between 0.4 and 8 mm) and large $a^{\prime }$ range (between 0.5 and 0.99) are also discussed. From the results, it can be concluded that the effect of optimization using the best forward model is significant.

However, in the actual process of measuring the optical properties, the value of $a^{\prime }$ cannot be determined because optical properties are unknown. Thus, how to choose the range of $a^{\prime }$ is an important issue. Figure 12 shows a method using the $k$ and $d_{x=1}$ of each ${\mathrm{Log}}_{10}(R_{\mathrm{meas}})$ to determine the range of $a^{\prime }$. Although there exists the situation of the overlapping regions of $k$ and $d_{x=1}$ at the same time for different $a^{\prime }$ ranges, the increase of the derivation error using the falsely best forward model is less than 10% in the case of $1{\mathrm{mm}}^{-1}\le {\mu }_s^{\prime }\le 2{\mathrm{mm}}^{-1}$. Thus, a complete joint derivation method which is an effective method to improve the derivation accuracy is shown in Fig. 10.

In addition, the trust-region-reflective algorithm of nonlinear least-squares (nonlinear data-fitting) methods is mainly applied for inversion calculations in this paper. Other algorithms, such as the Levenberg-Marquardt algorithm and the Gauss-Newton algorithm, are also widely used. We compared these three methods for inversion calculations and found that the derivation results from these three methods are similar and the conclusion in Table 1 was unchanged. Besides, the whole study focused on small source-detector separations because the diffuse reflectance attenuates quickly as the SDS increased for small reduced albedo. The condition of different largest SDSs (being 6, 7, 8, 9 mm, etc.) and different spatial resolutions (being 0.1, 0.2, and 0.3 mm) are discussed to validate the joint derivation method. All the results show that the result of Table 1 is a general or universal conclusion. In the future, the joint derivation method will be applied to the phantom experiment, and the results will be presented.