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₁ 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₁) 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₁ approximation [4]; The hybrid DA and P₃ (Hybrid-DA-P₃) 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₃ [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\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₃ 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₁ 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₁, PFC-DA and Hybrid-DA-P₃ 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₁ approximation can be expressed as [4]

 

Fig. 1 Schematic of the VS approximation in semi-infinite geometry, where the multiple VSs, _{$S_1$} to _{$S_N$}, are adopted at depths _{$z_1$} to _{$z_N$}to replace the realistic collimated source. The banana-shaped areas with shoaling grayscales indicate the most probable paths of the photon transmission from the VSs to the detector _{$D_m$}(m = 1…M) but weaken effects of the VSs at increasing depths on the detection.

Download Full Size | PDF

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]

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

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

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⁷ to 10⁸ 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₁ and z₂, 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₂ and k by fixing z₁ to a small value. In practice, it is reasonable to set z₁ to the mean of the four fitted values, i.e., z₁ = 8.29 _$\times$10⁻³ l_t.

Table 1. The VS parameters in the 2VS-DA model fitted for different combinations of _{${\mu }_a$} and _{${{\mu }^{\prime }}_s$}

View Table | View all tables in this article

To discern the relationship of the fitted parameter set of the 2VS-DA model, {k, z₂}, 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₂ (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₂ and k are assumed to take the same explicit expression of

 

Fig. 2 The fitted parameters (z₂, k) of the 2VS-DA as functions of the optical properties: (a) _{$5\le {\mu }_s\le 50{\text{mm}}^{-1}$} with fixed _{${\mu }_a=0.01{\text{mm}}^{-1}$} and _$g=0.9$; (b) _{$0.01\le {\mu }_a\le 0.4{\text{mm}}^{-1}$} with fixed _{${\mu }_s=10{\text{mm}}^{-1}$} and _$g=0.9$; (c) _{$0.7\le g\le 0.9$} with fixed _{${\mu }_a=0.01{\text{mm}}^{-1}$} and _{${\mu }_s=10{\text{mm}}^{-1}$}.

Download Full Size | PDF

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₂- 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{ε}}_{\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{ε}}_{\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{ε}}_k$} and _{${\text{ε}}_{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₂ and k, respectively.

 

Fig. 3 The fitting errors of (a) k and (b) z₂ in the 2VS-DA, calculated for _{$0.01\le {\mu }_a\le 0.5{\text{mm}}^{-1}$}, _{$5\le {\mu }_s\le 50{\text{mm}}^{-1}$}, and _$g=0.9$

Download Full Size | PDF

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₃ with PCBC, PFC-DA with extrapolated boundary condition (EBC) and δ-P₁ 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₃, PFC-DA, and δ-P₁

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₃, PFC-DA, and δ-P₁, 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₃, PFC-DA, and δ-P₁) 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{ε}}_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.

 

Fig. 4 Reflectances _{$R\left(\rho \right)$} calculated with the 2VS-DA, Hybrid-DA-P₃, PFC-DA, and δ-P₁ in comparison with the MC (upper), as well as their relative errors _{${\text{ε}}_R$} (lower), for three optical property sets: (a) _{${\mu }_a=0.01{\text{mm}}^{-1}$} and _{${\alpha }^{\prime }=0.99({{\mu }^{\prime }}_s=1{\text{mm}}^{-1})$}; (b)_{${\mu }_a=0.3{\text{mm}}^{-1}$} and _{${\alpha }^{\prime }=0.9({{\mu }^{\prime }}_s=2.7{\text{mm}}^{-1})$}; (c) _{${\mu }_a=0.3{\text{mm}}^{-1}$} and _{${\alpha }^{\prime }=0.769({{\mu }^{\prime }}_s=1{\text{mm}}^{-1})$}. The inset figures illustrate _{$R\left(\rho \right)$} in the nearly null SDS field of _{$\rho /l_t<1$}.

Download Full Size | PDF

Table 2. Concordance Correlation Coefficients of the four models with the MC at _{$0<\rho /l_t\le 1$} and _{$1<\rho /l_t\le 5$}, and for the three different combinations of ${\mu }_a$ and _{${\alpha }^{\prime }$}, respectively.

View Table | View all tables in this article

For the high cases of _{${\alpha }^{\prime }\ge 0.9$}, Fig. 4 (a) and (b) show that, the relative error _{${\text{ε}}_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₃ and δ-P₁ 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₃ approximation well accommodates the low-albedo medium for _{$1<{\rho /l}_t\le 5$}, with r_c 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₃ 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{ε}}_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.

 

Fig. 5 Relative error _{${\text{ε}}_R$} of the 2VS-DA and PFC-DA for three optical property sets: (a) _{${\mu }_a=0.01{\text{mm}}^{-1}$} and _{${\alpha }^{\prime }=0.99$}; (b)_{${\mu }_a=0.3{\text{mm}}^{-1}$} and _{${\alpha }^{\prime }=0.9$}; (c) _{${\mu }_a=0.3{\text{mm}}^{-1}$} and _{${\alpha }^{\prime }=0.769$} for refractive-index-matched boundary.

Download Full Size | PDF

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{ε}}_{\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.

 

Fig. 6 Relative errors _{${\text{ε}}_{\Phi }$} of (a) the 2VS-DA, (b) Hybrid-DA-P₃, (c) PFC-DA, and (d) δ-P₁, for three optical property sets: (Top) _{${\mu }_a=0.01{\text{mm}}^{-1}$} and _{${\alpha }^{\prime }=0.99$}; (Middle)_{${\mu }_a=0.3{\text{mm}}^{-1}$} and _{${\alpha }^{\prime }=0.9$}; (Bottom) _{${\mu }_a=0.3{\text{mm}}^{-1}$} and _{${\alpha }^{\prime }=0.769$}, respectively.

Download Full Size | PDF

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{ε}}_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.

 

Fig. 7 Relative error _{${\text{ε}}_R$} of the 2VS-DA for _{${\alpha }^{\prime }=$}0.769, 0.85, and 0.9, respectively, with the absorption coefficient fixed at _{${\mu }_a=0.3{\text{mm}}^{-1}$}.

Download Full Size | PDF

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.

 

Fig. 8 Relative errors _{${\text{ε}}_R$} of the 2VS-DA and 3VS-DA for two optical property sets: (a) _{${\mu }_a=0.3{\text{mm}}^{-1}$} and _{${\alpha }^{\prime }=0.769$}; (b) _{${\mu }_a=0.01{\text{mm}}^{-1}$} and _{${\alpha }^{\prime }=0.99$}, respectively.

Download Full Size | PDF

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₁ 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₁ is notably inconsistent with the MC as a result of the modeling discrepancy in the near-field, as illustrated in Fig. 4.

 

Fig. 9 Comparison between the relative deviations _{${\widehat{\sigma }}_R$} of the 2VS-DA and δ-P₁ for assessing performance dependency on the anisotropy factor. The calculation is conducted for a range of anisotropy factors of _{$0.7\le g\le 0.9$}, and two optical property sets of (a) _{${\mu }_a=0.3{\text{mm}}^{-1}$} and _{${\alpha }^{\prime }=0.769$}; (b) ${\mu }_a=0.01{\text{mm}}^{-1}$ and _{${\alpha }^{\prime }=0.99$}, respectively.

Download Full Size | PDF

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{ε}}_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

 

Fig. 10 Relative errors _{${\text{ε}}_R\left(\rho \right)/{{\epsilon }^{\prime }}_R\left(\rho \right)$} of the 2VS-DA without (upper) and with (lower) compensation for detector NAs of 0.05, 0.12, 0.22, and 0.37, and for two combinations of the optical properties: (a) ${\mu }_a=0.3{\text{mm}}^{-1}$ and _{${\alpha }^{\prime }=0.769$}; (b) ${\mu }_a=0.01{\text{mm}}^{-1}$ and ${\alpha }^{\prime }=0.99$, respectively.

Download Full Size | PDF

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{ε}}_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_n) as described in Sec. 2.3.

Table 3. Fitting results of _{${\beta }_{1,2}$} with four typical NAs for two optical property sets.

View Table | View all tables in this article

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{ε}}_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.5l_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₃ and δ-P₁, briefly the 2VS-Hybrid-DA-P₃ and 2VS-δ-P₁, 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.

 

Fig. 11 Relative errors_{${\text{ε}}_R$} of the 2VS-DA, 2VS-Hybrid-DA-P₃ and 2VS-δ-P₁, for two optical property sets: (a)_{${\mu }_a=0.01{\text{mm}}^{-1}$} and _{${\alpha }^{\prime }=0.99$}; (b)_{${\mu }_a=0.3{\text{mm}}^{-1}$} and _{${\alpha }^{\prime }=0.9$}; (c) _{${\mu }_a=0.3{\text{mm}}^{-1}$} and _{${\alpha }^{\prime }=0.769$}, respectively, with the anisotropic factor fixed at g = 0.9.

Download Full Size | PDF

In these scenarios, although the diffuse photon density in Eq. (4) should be calculated by the Hybrid-DA-P₃ or δ-P₁ 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.