## Abstract

Most analytical methods for describing light propagation in turbid medium exhibit low effectiveness in the near-field of a collimated source. Motivated by the Charge Simulation Method in electromagnetic theory as well as the established discrete source based modeling, we herein report on an improved explicit model for a semi-infinite geometry, referred to as “Virtual Source” (VS) diffuse approximation (DA), to fit for low-albedo medium and short source-detector separation. In this model, the collimated light in the standard DA is analogously approximated as multiple isotropic point sources (VS) distributed along the incident direction. For performance enhancement, a fitting procedure between the calculated and realistic reflectances is adopted in the near-field to optimize the VS parameters (intensities and locations). To be practically applicable, an explicit 2VS-DA model is established based on close-form derivations of the VS parameters for the typical ranges of the optical parameters. This parameterized scheme is proved to inherit the mathematical simplicity of the DA approximation while considerably extending its validity in modeling the near-field photon migration in low-albedo medium. The superiority of the proposed VS-DA method to the established ones is demonstrated in comparison with Monte-Carlo simulations over wide ranges of the source-detector separation and the medium optical properties.

© 2015 Optical Society of America

## 1. Introduction

In many diffuse optical spectroscopy (DOS) applications in medicine, such as endoscopic diagnosis and mesoscopic imaging, it has been desirable to establish a mathematically tractable model for near-field photon migration in turbid medium, especially with low albedo. The currently used diffusion approximation (DA) [1], *i.e.*, the P_{1} approximation to the Boltzmann radiative transport equation (RTE), is simple-to-implement but valid only in medium with high transport albedo: _{${\alpha}^{\prime}\left(={{\mu}^{\prime}}_{s}/{{\mu}^{\prime}}_{t}\right)>0.98$}, and for a source-detector separation (SDS) _{$\rho $} larger than about one transport mean free pathlength (mfp) _{${l}_{t}(\text{=}1/{{\mu}^{\prime}}_{t})$} [1], where _{${{\mu}^{\prime}}_{s}=(1-g){\mu}_{s}$} and _{${{\mu}^{\prime}}_{t}={\mu}_{a}+{{\mu}^{\prime}}_{s}$} are the reduced scattering and the total interaction coefficients, respectively, with _{$g$}, _{${\mu}_{a}$} and _{${\mu}_{s}$} being the anisotropy factor, absorption and scattering coefficients, respectively. To cope with the disadvantage, many improved schemes have been developed to more accurately describe the photon migration behavior in the near-field of low-albedo medium, while circumventing the computational complexity of the full P_{N} and discrete-coordinates (S_{N}) approximations [2]. The simplified spherical harmonics (SP_{N}, _{$\text{N}\ge 3$}) model approximates the RTE by coupled diffusion-like equations based on three-dimensional generalizations of the one-dimensional spherical harmonics (P_{N}) equations, and its semi-analytical solution is derived with universal geometry applicability [3]. In addition to the methods that take the high-order scattering effects into account, some others are reported based on correcting the DA to incorporate the scattering directionality: The Delta-Eddington (δ-P_{1}) approximation adds a Dirac-δ term to both the phase function and radiance expansions to explicitly decompose the light field into ballistic and diffuse components in the P_{1} approximation [4]; The hybrid DA and P_{3} (Hybrid-DA-P_{3}) modifying the DA by replacing the diffusion coefficient _{$D\left(=1/(3{{\mu}^{\prime}}_{\text{t}})\right)$} with the asymptotic diffusion coefficient _{${D}_{\text{asym}}$} of the P_{3} [5]; The phase-function-corrected DA (PFC-DA) takes advantage of the δ-isotropic phase function to derive a two-part expression for the diffuse reflectance, but restricted to the boundary condition of refractive-index-match [6,7]. Nevertheless, all these modified models are still less effective as applied in the near-field measurement, provided that the single isotropic point source model is assumed. To compensate the DA solution under an anisotropic point source, a semi-sphere source model was proposed with its radius and strength empirically derived depending on the optical properties of the medium [8].

For the measurement with a collimated source, in addition to the methods that modify the fidelity of the photon migration, those that are based on improved models of the collimated source have been also explored, among which the most conventional one is to take an isotropic source of line distribution along the incident direction of the collimated beam with exponentially decaying intensity: _{$I(z)={{\mu}^{\prime}}_{s}\mathrm{exp}\left(-{{\mu}^{\prime}}_{t}z\right)$}, *i.e.*, the line-source model [9]. To further remedy the scattering effect that might distort the exponentially decaying distribution, a hybrid numerical-analytical method has been surveyed with an isotropic line source of parameterized Gaussian distribution. This modified model was proved to be closer to Monte-Carlo (MC) simulations after 20 scattering interactions but presents notable errors in the near-field due to insufficient scattering of the incident light [10]. To be more physically realistic in the near field, one advanced but complex way is to model the collimated source as an isotropic continuous one with its intensity distribution relating to the phase function, such as the Henyey–Greenstein, and the reduced radiance component, which determine how the collimated light is scattered into the diffusion process [11]. In contrast to the continuous-distribution source models, it is also reasonable to employ discrete isotropic point sources to achieve both the improved description of the collimated source and mathematical conciseness of the applications. For example, Eason *et al*. proposed a multipole-expansion method that applies the Taylor's theorem to extend the integral expression of the DA solution to the exponentially decaying line source model into a series approximation to multiple dipoles (*i.e.*, isotropic point sources) [12]. Another simple source model, termed as the double-source, adopts only two isotropic point sources locating along the incident direction with the same dipole and quadrupole moments to an origin at the surrounding-tissue interface, and has been further combined with the Hybrid-DA-P_{3} approximation [5]. The results with this model remain significantly erroneous for _{$\rho /{l}_{t}<0.5$} and _{${\alpha}^{\prime}=0.99$}. As an alternative to the explicit source approximations, an implicit modeling can be conducted by embedding the realistic collimated source into the boundary condition. This scheme in general makes it difficult to express the solution to resultant equation in a closed form [13].

To enhance the photon migration modeling in the near-field (_{$0<\rho /{l}_{t}\le 5$}) of a collimated source as well as in low-albedo medium, a “Virtual Source (VS)” model is herein proposed in combination with the DA approximation, for the semi-infinite geometry. This model, abbreviated as the VS-DA, essentially extends the established discrete point source methods, by similarly modeling the collimated source as multiple isotropic point sources (VSs) and optimizing their locations and intensities for an effective fitting between the model calculation and the realistic data (simulated by the MC) over the near-field. Benefiting from the combined VS parameterization and P_{1} approximation strategy, the VS-DA model eventually offers an enhanced modeling of the near-field photon-migration and meanwhile retains the mathematical conciseness of the DA with an isotropic point source. To facilitate application of the method, an empirical strategy for explicitly deducing the VS parameters is specifically addressed for the commonly used range of the optical properties. The superior performance of the VS-DA to the established methods, *e.g.*, the δ-P_{1}, PFC-DA and Hybrid-DA-P_{3} *etc.*, is validated in comparison with the MC data. Also, it is demonstrated that further improved modeling of the near-field photon-migration can be potentially achieved by incorporating the VS methodology into the high-order approximation models.

## 2. Theory

#### 2.1 VS approximation to a collimated source

For a semi-infinite medium impinged by a collimated point source perpendicular to the boundary, *i.e.*, along the depth (z-direction), as shown in Fig. 1, the diffuse fluence-rate _{$\Phi (r)$} and the radiant flux _{$J(r)$} in the steady-state δ-P_{1} approximation can be expressed as [4]

_{${\mu}_{s}^{*}={\mu}_{s}(1-{g}^{2})$},

_{$Q(r,\widehat{z})$}describes the source with space position

_{$r=(x,y,z)$}, and

_{$\widehat{z}$}is the inward normal to the boundary. It is readily seen that, Eq. (1) would reduce to the standard DA equation:

_{$D{\nabla}^{2}\Phi (r)-{\mu}_{a}\Phi (r)=-{Q}^{*}(r,\widehat{z})$}, as its source item is modified into a form of

_{${Q}^{*}(r,\widehat{z})={\mu}_{s}^{*}Q(r,\widehat{z})$}. In addition, with an isotropic point source employed, the radiant flux given by Eq. (2) reduces to the Fick’s Law used with the standard DA:

_{$J(r)=-D\nabla \Phi (r)$}. Intuitively, the above formulations imply that, the physical essentials behind the δ-P

_{1}approximation is to effectively express directional nature of the collimated beam through the modification of the source item in the DA model. This strategy lends itself to a computationally-efficient remedy to the DA for low-albedo medium and also motivates the introduction of the VS for improved modeling of the near-field photon migration.

The VS conceptually originates from the Charge Simulation Method (CSM) in electromagnetic theory [14, 15]. The theorem claims that the electric field formed by a actual charge source of either continuous distribution or other complex type is equivalent to that formed together by a number of discrete point charges placed outside the region of interest, with their quantities and positions determined by satisfying the boundary conditions at a selected number of boundary sites [14, 15], which finally leads to a computational simplicity of the model due to availability of the analytical solution for the point source stimulation. Extending this principle to optical field, the realistic collimated light source directed perpendicularly to the medium surface can be analogously approximated to a group of discrete visual sources of point shape (_{${S}_{1}$} to _{${S}_{N}$}) along the z-direction, as shown in Fig. 1. Although, in a broad sense, there are several types of sources for choice, the isotropic point-shaped sources are herein adopted as the VS for the mathematical simplicity, as with the CSM scenarios. Besides the source substitution, the VS scheme also follows the uniqueness theorem as applied in the CSM, assuming that if the boundary conditions uniquely hold for the diffuse reflection generated by the equivalent point sources at the finite selected boundary sites, it must be continuously satisfied on the range covered by the selected sites even the whole boundary [14, 15]. Since the VS parameters including the number, intensities and positions of the constituent point sources, as well as a selected boundary locations can be optimized, the VS model has in principle the potential to obtain more accurate diffuse reflectance than the double- and line-source ones, in both the near- and far-fields, provided that an enough number of the VSs are adopted.

#### 2.2 VS*-*DA model

*-*

We calculate the boundary diffuse reflectance $R\left(\rho \right)$ under the partial-current boundary condition (PCBC) since it is closer to the actual application as compared to the extrapolated boundary condition [16]

*n*

_{t}being the relative refractive index of the tissue and assumed to be a constant; $\Phi $ is the fluence-rate that, for the DA with a single isotropic point source, is given as [16]

According to the uniqueness criteria, the intensities and positions of the VSs can be determined by satisfying the boundary conditions at multiple detection positions. Upon these values, the diffuse reflectance on the boundary (and the fluence-rate within the medium) can be then obtained by the superposition theorem

_{${I}_{n}$}is the intensity of the

*n*-th VS, noted as

_{${S}_{n}$}, and

_{${R}_{n}\left(\rho \right)$}is the diffuse reflectance at a SDS

_{$\rho $}generated by

_{${S}_{n}$}, and calculated by Eq. (3). The selection of VS number

*N*should be a tradeoff between the computation complicity and accuracy requirement of the modeling. To calculate the depth

_{${z}_{n}$}and intensity

_{${I}_{n}$}of the VSs, the following VS equation can be adoptedwhere ${R}_{\text{REF}}\left(\rho \right)$ is the reference data given by an advanced models suitable for describing the near-field light propagation, such as the MC simulation or the RTE. Considering ${R}_{\text{VS}}$ might not strictly equal to ${R}_{\text{MC}}$ over the whole boundary, the symbol “

_{$\to $}” is adopted instead of “=” to imply an application of the model fitting strategy. In addition, the location and number of “detecting” sites also affect the accuracy of the VS-DA model, a large “detecting” number

*M*is preferred. In this sense, the solution to the VS equation can be obtained in principle by a least-squares optimization procedure under appropriate constraints on the VS parameters

*M*detecting locations on the boundary with the SDSs:

_{${\rho}_{1},{\rho}_{2},\dots ,{\rho}_{M}$}, respectively;

_{$I={\left[{I}_{1},{I}_{2},\cdots ,{I}_{N}\right]}^{\text{T}}$}and $z={\left[{z}_{1},{z}_{2},\cdots ,{z}_{N}\right]}^{\text{T}}$ numerate the intensities and depths of the

*N*VSs, respectively; ${I}_{\text{total}}=1-{\left(1-{n}_{\text{t}}\right)}^{2}/{\left(1+{n}_{\text{t}}\right)}^{2}$ is the penetrating intensity of the incident collimated beam normal to the surface; $\Lambda =\text{diag}\left(\left[{R}_{\text{REF}}^{2}\left({\rho}_{1}\right),\cdots ,{R}_{\text{REF}}^{2}\left({\rho}_{M}\right)\right]\right)$.

To solve the nonlinear optimization problem in Eq. (7), a Newton-Raphson iteration scheme is used for its computational simplicity

_{${J}_{I}$}and

_{${J}_{z}$}are the Jacobian matrices with regard to the intensities and depths of the VSs, respectively. The resultant linear equations can be universally solved using an algebraic reconstruction technique. The constraints in Eq. (7) are achieved by simply employing a substitution of

_{${I}_{1}={I}_{\text{total}}-{\displaystyle {\sum}_{n=2}^{N}{I}_{n}}$}in Eq. (8) and forcing

_{${z}_{n}\ge 0$}after each iteration. The Lin’s Concordance Correlation Coefficient (

*r*) is employed as a similarity criterion to stop the iteration procedure as well as to evaluate the model effectiveness, which strikes a balance between a correlation measure insensitive to location differences and a measure of location discrepancy [17]

_{c}_{${R}_{\text{REF(Model)}}$}, respectively. By the definition,

*r*, is scaled to between −1 and 1, with “1” denoting the best fitting and “-1” the worst. We can thus terminate the iterative fitting process in Eq. (8), as

_{c}*r*stops increasing between the two successive iterations.

#### 2.3 Determination of VS parameters in close-form

For the VS-DA applications, it is convenient to calculate the VS parameters, *i.e.*, the intensities and locations of the VSs, in an explicit way, given the optical properties of the medium. First of all, the practically required number of the VSs, *N*, needs to be determined. In theory, the accuracy of the VS-DA approximation improves with increase in the VS number. However, this normally leads to degraded uncertainty of the solution, which has to be compensated by either more reference data or some priors. On the other hand, as illustrated by banana-shaped areas of different grayscales in Fig. 1, it is readily understood that, the VSs at increasing depths (except for those that approach the surface) exert decreasing effects on the detection at a certain SDS, due to the approximately exponential decays of both the penetrating beam intensity along the depth and the detected diffuse reflectance with the SDS. Thus, modeling with reasonable accuracy is probably attained with a limited number of the VSs that are optimized for over an appropriate range of _{$\rho $}. Without loss of generality, we start a trial from the DA with two VSs, referred to as the 2VS-DA model. Secondly, with an ultimate aim of establishing a VS-DA based DOS technique, solutions to the VS equation are found for the typical optical properties of biological tissue in the visible or near-infrared wavelength range, *i.e.*, _{$0.01\le {\mu}_{a}\le 0.5{\text{mm}}^{\text{-1}}$}, _{$5\le {\mu}_{s}\le 50{\text{mm}}^{\text{-1}}$}, and _{$0.7\le g\le 0.9$} [18,19], to further construct explicit expressions of the VS depths and intensities as functions of the optical properties, as described below. The MCML program was used to execute the Monte Carlo simulation to calculate *R*_{REF} [20], where the spatial resolution is set to _{${l}_{t}/10$} and the actual numbers of photons launched varies from 10^{7} to 10^{8} in terms of the medium albedo, respectively. *M* ( = 40) detecting locations are used spanning the near-field range with an equal spacing.

Suppose that the intensity of the 2nd VS is *k* times that of the 1st one, *i.e.*,_{${I}_{2}=k{I}_{1}$}, then we have _{$I=\left[{I}_{\text{total}}/\left(1+k\right),k{I}_{\text{total}}/\left(1+k\right)\right]$}. The above relationship transforms the _{${I}_{1}$}- and/or _{${I}_{2}$}-regarding fitting procedures to the only _{$k$}-regarding one. With the above substitution, the whole parameter set, *i.e.*, *k*, *z*_{1} and *z*_{2}, are calculated by the fitting procedure for the above typical optical properties, and listed in Table 1 for four different _{${\mu}_{a}$} and _{${{\mu}^{\prime}}_{s}$} combinations. With the observation of _{${z}_{1}\approx 0$} from the table, the total fitting parameters in the 2VS-DA model would finally reduce to *z*_{2} and *k* by fixing *z*_{1} to a small value. In practice, it is reasonable to set *z*_{1} to the mean of the four fitted values, *i.e.*, *z*_{1} = 8.29 _{$\times $}10^{−3} *l _{t}*.

To discern the relationship of the fitted parameter set of the 2VS-DA model, {*k*, *z*_{2}}, with the medium optical properties, Fig. 2 presents the two quantities as functions of ${\mu}_{s}$, ${\mu}_{a}$ and *g*, respectively, with the two others fixed. It is seen that *z*_{2} (or *k*) is approximately in inverse and direct proportion to _{${\mu}_{s}$} and _{${\mu}_{a}$}, respectively, and dependent on *g* in a nearly quadratic nonlinear way. Based on these observations, both *z*_{2} and *k* are assumed to take the same explicit expression of

For the determination of the above coefficient sets ${A}_{1-7}=\left\{{A}_{1},{A}_{2},\cdots ,{A}_{7}\right\}$, a fitting procedure between the solutions to the VS equation, Eq. (7), and the explicit function, Eq. (10), is adopted, where 6 combinations of the absorption and scattering coefficients and 3 anisotropy factors selected uniformly from the typical ranges of the absorption and scattering coefficients, *i.e.*, _{$0.01{\text{mm}}^{\text{-1}}\le {\mu}_{a}\le 0.5{\text{mm}}^{\text{-1}}$}, _{$5{\text{mm}}^{\text{-1}}\le {\mu}_{s}\le 50{\text{mm}}^{\text{-1}}$}, and *g* = 0.7, 0.8, and 0.9 are used, respectively, as the calibration points. ${A}_{1-7}=\left\{{A}_{1},{A}_{2},\cdots ,{A}_{7}\right\}$ for the *z*_{2}- and *k*- explicit expressions are obtained under the 18 combinations of the selected optical properties, respectively

For assessment of the above 2VS-DA explicit expressions, the fitting errors, ${\text{\epsilon}}_{\xi}$_{$\left(=\left({\xi}_{\text{im}}-{\xi}_{\text{ex}}\right)/{\xi}_{\text{ex}}\right)$}, of $\xi $$\left(\in \left\{{z}_{2},k\right\}\right)$ are calculated over the typical ranges of the absorption and scattering coefficients with *g* = 0.9, to compare the discrepancy between the implicit 2VS-DA solution ${\xi}_{\text{im}}$ and the explicit calculation ${\xi}_{\text{ex}}$, as illustrated in Fig. 3. The mean of squared residuals, _{$\overline{\epsilon}=\sqrt{{\displaystyle {\sum}_{i=1}^{{N}_{c}}{\left({\text{\epsilon}}_{\xi}^{i}\right)}^{2}}/{N}_{c}}$} (where _{${N}_{c}$} is the sample number), is also employed to evaluate the overall performance of the fitting procedure. It is seen from the figure, that _{${\text{\epsilon}}_{k}$} and _{${\text{\epsilon}}_{{z}_{\text{2}}}$} are within ± 10% and ± 1.5%, respectively, fairly conforming the effectiveness of the explicit solutions. Since the above errors are calculated under a constraints of ${\alpha}^{\prime}\ge 0.75$, the results also suggests an effective albedo range of _{$0.75\le {\alpha}^{\prime}\le 0.998$} for the 2VS-DA application, as will be discussed in Sec. 3.2. Finally, the sum of squared residual $\overline{\epsilon}$ is calculated to be 0.63% and 4.73% for *z*_{2} and *k*, respectively.

With the fitted VS parameters back to Eq. (5), the steady-state diffuse reflectance in the 2VS-DA model with PCBC is given as

## 3. Results

We first compare the proposed 2VS-DA with the three established models, *i.e.*, Hybrid-DA-P_{3} with PCBC, PFC-DA with extrapolated boundary condition (EBC) and δ-P_{1} with PCBC. Then, the 2VS-DA models are assessed for different albedo, anisotropic factor *g*, and detector numerical aperture (NA). The same MC configurations and the detecting locations are used as those in the above subsection. The relative refractive index of the tissue *n*_{t} = 1.4 are assumed for all the following results.

#### 3.1 Comparison of 2VS-DA with Hybrid-DA-P_{3}, PFC-DA, and δ-P_{1}

As shown in Fig. 4, the diffuse reflectances, $R\left(\rho \right)$, in the near-field are calculated from Eq. (12), as well as from the Hybrid-DA-P_{3}, PFC-DA, and δ-P_{1}, and compared to that from MC simulation ${R}_{\text{MC}}\left(\rho \right)$, for a wide range of the optical properties: _{$0.01\le {\mu}_{a}\le 0.3{\text{mm}}^{\text{-1}}$} and _{$0.769\le {\alpha}^{\prime}\le 0.99$}, and a fixed anisotropy factor of _{$g=0.9$}. For a fair comparison, the isotropic line-source model is employed in the three established models (Hybrid-DA-P_{3}, PFC-DA, and δ-P_{1}) for more accurate description of the collimated source. The deviations of the four models reflectances from ${R}_{\text{MC}}$ are also calculated based on a relative error: ${\text{\epsilon}}_{R}\left(\rho \right)=\left|R\left(\rho \right)-{R}_{\text{MC}}\left(\rho \right)\right|/{R}_{\text{MC}}\left(\rho \right)$. To be concise, the results are presented for only three different combinations of _{${\mu}_{a}$} and _{${\alpha}^{\prime}$}, representing the normal (_{${\mu}_{a}=0.01{\text{mm}}^{\text{-1}},{\alpha}^{\prime}=0.99$}), the high absorption (_{${\mu}_{a}=0.3{\text{mm}}^{\text{-1}},{\alpha}^{\prime}=0.9$}), as well as the high absorption and low-albedo (_{${\mu}_{a}=0.3{\text{mm}}^{\text{-1}},{\alpha}^{\prime}=0.769$}) situations, respectively. It is noteworthy that the minimum of the abscissa axis “0” actually denotes a tiny SDS of _{$\rho /{l}_{t}=0.01$} in the Fig. 4 (also all the following figures). To intuitively quantify the agreement of the four models with the MC reference, the Concordance Correlation Coefficients _{${r}_{c}$} is listed in Table 2.

For the high cases of _{${\alpha}^{\prime}\ge 0.9$}, Fig. 4 (a) and (b) show that, the relative error _{${\text{\epsilon}}_{R}$} provided by the 2VS-DA is less than 6% over the full range of _{${\rho /l}_{t}$}, and that the correlation coefficients _{${r}_{c}$} in Table 2 correspondingly answers to the errors, with those of the 2VS-DA proximal to “1” for all the combinations. It is also shown that, although the performances of the Hybrid-DA-P_{3} and δ-P_{1} approximations are slightly degraded compared to the 2VS-DA in the field of _{$1<{\rho /l}_{t}\le 5$}, they are notably erroneous in the nearly null SDS field of _{${\rho /l}_{t}\le 1$}. For the low albedo case of _{${\alpha}^{\prime}=0.769$}, as illustrated in Fig. 4 (c), the maximum error of the 2VS-DA approximation is above 10% in _{$\rho /{l}_{t}\le 1$} and *r _{c}* correspondingly reduces to below 0.984, indicating the less applicability of the 2VS-DA in the low-albedo case. On the other hand, the Hybrid-DA-P

_{3}approximation well accommodates the low-albedo medium for

_{$1<{\rho /l}_{t}\le 5$}, with

*r*approaching “1”, but again exhibits evident error for the nearly null SDS field. In summary, comparing to the established models with the exponentially decaying line-source, the improved near-field performance of the 2VS-DA is theoretically interpreted as the modified source description to allow both the VS intensities and positions to be optimized for the full boundary fitting. The similar effect can also be found in the Hybrid-DA-P

_{c}_{3}approximation for the double-source model, although the performance is inferior to the line source case due to the inefficient characterization of the source scattering [5]. As to the poor performance of the PFC-DA model, it is most probably attributed to the non-matched boundary condition, where a relative refractive index of 1.4 is realistically used.

For fair comparison, we additionally calculate in Fig. 5 the deviations of the PFC-DA and 2VS-DA to MC under the matched (*n*_{t} = 1) boundary condition for the same three optical property combinations as the above. For the high albedo cases of _{${\alpha}^{\prime}\ge 0.9$}, the PFC-DA model appealingly presents the excellent performance with _{${\text{\epsilon}}_{R}\le 10\%$} near the point of source entry (_{${\rho /l}_{t}\le 0.2$}), but significantly increased outside this range, especially for the low albedo case of _{${\alpha}^{\prime}=0.769$}. In contrast, the 2VS-DA exhibits the same performance as the PFC-DA near the point of source entry for all the three albedos, and is significantly superior to the PFC-DA over the whole near-field range.

To further validate the VS-DA superiority, the fluence-rates, $\Phi \left(\rho ,z\right)$, in the near-field are calculated from the 2VS-DA, Hybrid-DA-P3, PFC-DA and δ-P1, for the same three combinations of the optical properties under the mismatched boundary condition, and compared to the MC counterpart${\Phi}_{\text{MC}}$ in terms of the deviation:${\text{\epsilon}}_{\Phi}\left(\rho ,z\right)=\left|\Phi \left(\rho ,z\right)-{\Phi}_{\text{MC}}\left(\rho ,z\right)\right|/{\Phi}_{\text{MC}}\left(\rho ,z\right)$, respectively, as shown in Fig. 6. In agreement with the reflectance comparison, the 2VS-DA similarly outperforms the four established models in describing the near-field fluence-rate, with the improvement depending on the albedo. The inferior performance of the model for the low-albedo case, as shown in Fig. 4(c), might be interpreted as the result of the degraded match of the boundary conditions.

#### 3.2 Accommodation of 2VS-DA to different albedos

It is observed in Fig. 4(c) and Table 2, that all the models show a degraded near-field performance with the increased _{${\text{\epsilon}}_{R}$} for the high absorption (_{${\mu}_{a}=0.3{\text{mm}}^{-1}$}) and/or the low albedo (_{${\alpha}^{\prime}=0.769$}) case. Although the 2VS-DA is still applicable for the nearly null SDS field (_{${\rho /l}_{t}\le 1$}), the correlation coefficient *r _{c}* of the 2VS-DA drops to 0.984. This observation implies that, the directionality and exponential-decay nature of the beam transmission as described in the Bouguer’s law are restored with the scattering vanishing, and it would be increasingly difficult for the VSs to effectively model a collimated source [4]. To further address the applicable range of

_{${\alpha}^{\prime}$}for the 2VS-DA, Fig. 7 illustrates the relative errors of the model, for three different albedos of

_{${\alpha}^{\prime}=$}0.769, 0.85, and 0.9, but with a fixed absorption coefficient of

_{${\mu}_{a}=0.3{\text{mm}}^{-1}$}. All the errors for the three albedos are estimated below 13% over the whole near-field range, and overall go down with the albedo increasing. The appearance of “ripple wave” in the Fig. 7 is ascribed to the uncertainty in the MC data, since, with the increment of

_{${\mu}_{s}$}and

_{$\rho $}, the 2VS-DA keep generating smooth predictions while the shot-noise in the MC data ${R}_{\text{MC}}$ rapidly increasing. From the conclusion that the DA can be in principle made applicable for the low-albedo medium by modifying the source term, it is reasonable to deduce that the VS-DA approximation would be optimized to accommodate the low-albedo medium with improved performance.

#### 3.3 Performance improvement with increased VS number

To improve the applicability of the VS-DA to the near-field and lower-albedo medium, a natural way is to increase the VS number, which delivers more freedom to be fitted for improved capture of the near-field photon scattering patterns. For example, with incorporating three VSs in to the DA, briefly the 3VS-DA, the similar optimization procedure as for the 2VS-DA can be used, where the model parameters are solved as _{$I=\left[{I}_{\text{total}}/\left(1+{k}_{1}+{k}_{2}\right),{k}_{1}{I}_{\text{total}}/\left(1+{k}_{1}+{k}_{2}\right),{k}_{2}{I}_{\text{total}}/\left(1+{k}_{1}+{k}_{2}\right)\right]$} with (_{${k}_{1-2}=\left\{2.1071,\text{3}\text{.1532}\right\}$}, _{${z}_{1-3}=\left\{0.0\text{mm},1.0129\text{mm},2.002\text{mm}\right\}$}) for the low-albedo scenario:_{${\mu}_{a}=0.3{\text{mm}}^{-1}$} and _{${\alpha}^{\prime}=0.769$}, and (_{${k}_{1-2}=\left\{2.6078,\text{0}\text{.9998}\right\}$}, _{${z}_{1-3}=\left\{0.0\text{mm},1.348\text{mm},1.0778\text{mm}\right\}$}) for the normal tissue: ${\mu}_{a}=0.01{\text{mm}}^{-1}$ and _{${\alpha}^{\prime}=0.99$}. The relative errors of the 2VS-DA and the 3VS-DA approximations for the two optical property sets are compared, respectively, in Fig. 8. As expected, a performance improvement is evident for the 3VS-DA, especially in the near null SDS range, where the relative errors of the 3VS-DA fall below 3% for both the two cases, contrasting to the peak error of 13% for the 2VS-DA for the first set and 5% for the second. Nevertheless, a compromise should be made in practice between fitting accuracy and the VS number, depending on the task details.

#### 3.4 Sensitivity to anisotropy factor

Although, as with the standard DA, the anisotropy factor is eliminated in the explicit expression of the VS-DA through the fitting procedure for determine the VS parameters, the model is in nature *g*-dependent. Therefore it is worthy of assessing the sensitivity of the 2VS-DA to the anisotropic factor. For the implementation, the relative deviation _{${\widehat{\sigma}}_{R}={\sigma}_{R}/\overline{R}$} (where $\overline{R}$ and ${\sigma}_{R}$ is the mean and standard deviation of the diffuse reflection, respectively) is calculated as a function of _{$\rho /{l}_{t}$} over the typical *g*-range of 0.7-0.9 with an increment of 0.05, for two sets of the optical properties representing the high-absorption/low-albedo and normal tissues: (_{${\mu}_{a}=0.3{\text{mm}}^{-1}$},_{${\alpha}^{\prime}=0.769$}), and (_{${\mu}_{a}=0.01{\text{mm}}^{-1}$},_{${\alpha}^{\prime}=0.99$}), respectively, and compared among the 2VS-DA, δ-P_{1} and MC, as shown in Fig. 9. As expected, the relative deviations for both the scenarios show fair consistency between the 2VS-DA and MC, which implicitly validates the model applicability in the near-field. In contrast, the δ-P_{1} is notably inconsistent with the MC as a result of the modeling discrepancy in the near-field, as illustrated in Fig. 4.

#### 3.5 Influence of detector NAs

Since the radiance exhibits the varying angular distribution with the SDS due to insufficient scattering of the incident collimated light, the measured reflectance can be considerably affected by the detector NA in the near-field, as shown in Fig. 10, where the MC reflectances ${R}_{\text{MC}}^{\left(\text{NA}\right)}\left(\rho \right)$ are generated with the four conventional NAs of 0.05, 0.12, 0.22 and 0.37, and compared with the 2VS-DA data by calculating the relative error _{${\text{\epsilon}}_{R}\left(\rho \right)$} (_{$=\left|{R}_{\text{MC}}^{\left(\text{NA}\right)}\left(\rho \right)-{R}_{\text{VS}}\left(\rho \right)\right|/{R}_{\text{MC}}^{\left(\text{NA}\right)}\left(\rho \right)$}), for the low-albedo (_{${\mu}_{a}=0.3{\text{mm}}^{-1}$}, _{${\alpha}^{\prime}=0.769$}) and normal (_{${\mu}_{a}=0.01{\text{mm}}^{-1}$},_{${\alpha}^{\prime}=0.99$}) tissues, respectively. A great discrepancy between the MC and 2VS-DA data and notable dependence on the SDS are observed from the error profiles. To more effectively remedy the SDS-dependent NA effect, the VS-DA model in Eq. (5) is modified to the following form of linear relation

_{${\beta}_{n}$}is the scaling factor of the

*n*-th VS to be determined. This strategy of respectively scaling each VS retains the mathematical simplicity while releasing a freedom to appropriately count the NA-induced variation of the coupling efficiency with the SDS. Based on the VS-DA model that is essentially obtained as the solution for

_{$\text{NA}=\propto $},

_{${\beta}_{n}$}can be found by simply solving the following linear least-squares optimization

Table 3 lists the fitting results of _{$\beta $} for the two optical property sets. It is shown in Figs. 9(c) and (d) that, after the compensation for the NA effect, the relative errors _{${{\epsilon}^{\prime}}_{R}\left(\rho \right)$} (_{$=\left|{R}_{\text{MC}}^{\left(\text{NA}\right)}\left(\rho \right)-{R}_{\text{VS}}^{\left(\text{NA}\right)}\left(\rho \right)\right|/{R}_{MC}^{\left(\text{NA}\right)}\left(\rho \right)$}) all reduce to 27% and 25% for the two optical properties, respectively. Again, the notable “ripple wave” of _{${\text{\epsilon}}_{R}\left(\rho \right)$} for detector NA = 0.05 is explained as result of the MC statistical insufficiency with less received photons. In practice, to construct an explicit expression of the VS-DA model, the closed form for calculating _{${\beta}_{n}$} can be derived using a fitting procedure over the typical range of the optical properties (_{${\mu}_{a}$}, _{${\alpha}^{\prime}$}, *g*), analogous to the determination of the VS parameters (*I _{n}*,

*z*) as described in Sec. 2.3.

_{n}## 4. Discussion and conclusion

For both the normal and high-absorption/low-albedo medium, all the above assessments demonstrate that the VS-DA is superior to the established models in characterizing the photon transport nature in the near-field of a collimated source, and retains to be computationally efficient. In principle, this enhancement can be ascribed to the modified accuracy of the proposed model in delineating the incomplete scattering effect of a collimated source as well as the profound directionality of the photon migration in the near-field, which are comprehensively achieved by distributing multiple isotropic point sources along the depth and allowing both their intensities and positions to be optimized for the full near-field boundary fitting between the model prediction and the MC simulation.

As a result, the VS-DA model could be fairly applied to endoscopic and mesoscopic optical spectroscopy and imaging, such as laminar optical tomography or endoscopic diffuse optical detection where rapid fitting procedures or inverse problems are required using the data from short SDS detection [21]. For applications in the near-field with the typical optical properties, it is seen from Fig. 4 that a reasonably small model error of _{${\text{\epsilon}}_{R}<13\%$} can be achieved with the 2VS-DA. Nevertheless, as demonstrated in Fig. 8(a), the 3VS-DA would provide a better choice for applications in the near-field with the low albedo medium, though the computation cost is increased to a certain extent and more reference data might be required for the well-posedness of the VS equation.

In the VS-DA modeling, to guarantee the reliability of the model, the MC reference data should be generated with the signal-to-noise ratio (SNR) being more than 40dB over the whole near-field range. This means that the minimum photon number received at the *M* detecting locations exceeds 10000. Since the simulation is conducted for the near-field, the required SNR can be achieved without special computational burdens. In addition to improving the SNR, the reliability might be also enhanced by densely arranging the detecting locations in the nearly null SDS range for better capture of the rapid change in the reflectance. Concerning the inheritance of the VS-DA to the conventional DA, 0.5*l _{t}* is recommended as the lower limit of the applicable SDS, below which the inaccuracy of fluence-rate is particularly evident as seen in Fig. 6(a).

To apply the proposed method to the measurement of the optical properties or biomedical optical imaging, where a prior information about the background optical properties is unknown, a proper VS-DA model can be established beforehand by firstly solving the VS equation to optimize the VS parameters (intensity and position) over a predicted range of the optical properties, and then fitting these parameters to Eq. (10) to obtain their explicit expressions, as discussed in Sec. 2.3.

To take advantage of the high order approximations in more accurate modeling of the photon transport in the near field, the VS model can in principle be incorporated into the Hybrid-DA-P_{3} and δ-P_{1}, briefly the 2VS-Hybrid-DA-P_{3} and 2VS-δ-P_{1}, respectively, for two the sources. As an example, the relative errors in the near-field reflectances with the 2VS-based high order approximations are contrasted to that with the 2VS-DA in Fig. 11, for the three optical property sets as used in Sec. 3.1. The results clearly validate the further enhanced performance of the 2VS-based high-order approximations in modeling the near field photon transport.

In these scenarios, although the diffuse photon density in Eq. (4) should be calculated by the Hybrid-DA-P_{3} or δ-P_{1} and Eq. (4) be modified according to the corresponding boundary condition, the remaining steps to determinate the VS parameters are the same as that described in Sec. 2.2.

## Acknowledgments

The authors acknowledge the funding supports from the National Natural Science Foundation of China (81271618, 81371602, 61475115, 61475116, 81401453), the Research Fund for the Doctoral Program of Higher Education of China (20120032110056), and the Tianjin Municipal Government of China (13JCZDJC28000, 12JCQNJC09400).

## References and links

**1. **A. Ishimaru, *Wave Propagation and Scattering in Random Media* (Academic, 1978).

**2. **A. D. Klose, V. Ntziachristos, and A. H. Hielscher, “The inverse source problem based on the radiative transfer equation in optical molecular imaging,” J. Comput. Phys. **202**(1), 323–345 (2005). [CrossRef]

**3. **L.-M. Zhang, J. Li, X. Yi, H.-J. Zhao, and F. Gao, “Analytical solutions to the simplified spherical harmonics equations using eigen decompositions,” Opt. Lett. **38**(24), 5462–5465 (2013). [CrossRef] [PubMed]

**4. **W. M. Star, J. P. Marijnissen, and M. J. C. van Gemert, “Light dosimetry in optical phantoms and in tissues: I. Multiple flux and transport theory,” Phys. Med. Biol. **33**(4), 437–454 (1988). [CrossRef] [PubMed]

**5. **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). [CrossRef]

**6. **E. Vitkin, V. Turzhitsky, L. Qiu, L. Guo, I. Itzkan, E. B. Hanlon, and L. T. Perelman, “Photon diffusion near the point-of-entry in anisotropically scattering turbid media,” Nat. Commun. **2**, 587–605 (2011). [CrossRef] [PubMed]

**7. **R. J. Zemp, “Phase-function corrected diffusion model for diffuse reflectance of a pencil beam obliquely incident on a semi-infinite turbid medium,” J. Biomed. Opt. **18**(6), 067005 (2013). [CrossRef] [PubMed]

**8. **S. Fantini, M. A. Franceschini, and E. Gratton, “Effective source term in the diffusion equation for photon transport in turbid media,” Appl. Opt. **36**(1), 156–163 (1997). [CrossRef] [PubMed]

**9. **T. J. Farrell, M. S. Patterson, and B. Wilson, “A diffusion theory model of spatially resolved, steady-state diffuse reflectance for the noninvasive determination of tissue optical properties in vivo,” Med. Phys. **19**(4), 879–888 (1992). [CrossRef] [PubMed]

**10. **T. Spott and L. O. Svaasand, “Collimated light sources in the diffusion approximation,” Appl. Opt. **39**(34), 6453–6465 (2000). [CrossRef] [PubMed]

**11. **J. B. Domínguez and Y. Bérubé-Lauzière, “Light propagation from fluorescent probes in biological tissues by coupled time-dependent parabolic simplified spherical harmonics equations,” Biomed. Opt. Express **2**(4), 817–837 (2011). [PubMed]

**12. **G. Eason, A. R. Veitch, R. M. Nisbet, and F. W. Turnbull, “The theory of the back-scattering of light by blood,” J. Phys. D **11**(10), 1463–1479 (1978). [CrossRef]

**13. **M. Schweiger, S. R. Arridge, M. Hiraoka, and D. T. Delpy, “The finite element method for the propagation of light in scattering media: Boundary and source conditions,” Med. Phys. **22**(11), 1779–1792 (1995). [CrossRef] [PubMed]

**14. **N. H. Malik, “A review of the charge simulation method and its applications,” IEEE Trans. Electr. Insul. **24**(1), 3–20 (1989). [CrossRef]

**15. **K. Amano, “A charge simulation method for the numerical conformal mapping of interior, exterior and doubly-connected domains,” J. Comput. Appl. Math. **53**(3), 353–370 (1994). [CrossRef]

**16. **R. C. Haskell, L. O. Svaasand, T. T. Tsay, T. C. Feng, M. S. McAdams, and B. J. Tromberg, “Boundary conditions for the diffusion equation in radiative transfer,” J. Opt. Soc. Am. A **11**(10), 2727–2741 (1994). [CrossRef] [PubMed]

**17. **L. I.-K. Lin, “A concordance correlation coefficient to evaluate reproducibility,” Biometrics **45**(1), 255–268 (1989). [CrossRef] [PubMed]

**18. **S. L. Jacques, “Optical properties of biological tissues: a review,” Phys. Med. Biol. **58**(11), R37–R61 (2013). [CrossRef] [PubMed]

**19. **T. Lister, P. A. Wright, and P. H. Chappell, “Optical properties of human skin,” J. Biomed. Opt. **17**(9), 0909011 (2012). [CrossRef] [PubMed]

**20. **L. Wang, S. L. Jacques, and L. Zheng, “MCML--Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Methods Programs Biomed. **47**(2), 131–146 (1995). [CrossRef] [PubMed]

**21. **E. M. Hillman, D. A. Boas, A. M. Dale, and A. K. Dunn, “Laminar optical tomography: demonstration of millimeter-scale depth-resolved imaging in turbid media,” Opt. Lett. **29**(14), 1650–1652 (2004). [CrossRef] [PubMed]