2.1.

Monte Carlo Simulations

The light propagation was modeled by an MC stochastic approach derived from the code of Wang et al.²¹ and the CUDA-based implementation of Alerstam et al.²³ In our study, two optical fiber probe geometries that included five SDS were introduced.

The first geometry (SG) was based on a simple laterally uniform probe–tissue interface, which takes into account only the mismatch between the refractive indices of tissue ($n_m=1.33$) and optical fibers ($n_{\mathrm{fib}}=1.45$). Each simulation included ${10}^8$ photon packets, which were launched uniformly over the fiber opening and the solid angle defined by the numerical aperture ($\mathrm{NA}=0.22$). Additionally, the detection scheme was modified to account for the numerical aperture of the detector fibers. Due to the uniformity of the probe–tissue interface, the detection scheme around the source fiber was divided into $5\text{-}\mu \mathrm{m}$-thick concentric annular rings, from which the number of detected photon packets through a particular detection fiber was estimated in the postprocessing step. In this way, a better signal-to-noise ratio was obtained compared to the exact geometry of the individual detection fibers.

The second geometry (RG) was based on a realistic probe–tissue interface that additionally considers (1) the reflections from the stainless steel probe tip and (2) the refractive index mismatch between the black epoxy resin ($n_{\text{epoxy}}=1.6$) and the tissue (Fig. 1). Stainless steel reflectivity was set to 57% (measured value). Its effect on the SRR was already investigated and found significant.²⁴ The photon packets (${10}^9$) were launched following the scheme of the first (SG) geometry. In contrast, the detection scheme was constrained to individual detection fibers located at five different SDS.

Fig. 1

Linear array optical fiber probe used to acquire the SRR profiles in the subdiffusive regime. The materials and source-detector configuration of the probe are presented in different colors.

The MC simulations were conducted using a semi-infinite medium geometry. To speed up the simulations, all photon packets that drifted more than 0.8 cm laterally and axially from the source center were terminated. The two termination criteria were validated and proved consistent with the semi-infinite medium geometry.

2.2.

Scattering Phase Functions and Their Implementation in Monte Carlo Simulations

Since the SRR in the subdiffusive regime significantly depends on the phase function, several phase functions that have been proposed for biological tissues were considered. In addition to the common Henyey–Greenstein (HG) phase function $p_{\mathrm{HG}}$, we used the modified Henyey–Greenstein (MHG) phase function, which also accounts for the angular dependence of the Rayleigh scattering:¹²

Existing MC light propagation simulations primarily use the HG phase function to sample the scattering angles when a photon packet undergoes a scattering event. In this aspect, the HG phase function is convenient, since it allows analytical expression of the cumulative probability distribution (CPD) from which the scattering angles can be sampled. Analytical expression for the CPD can also be obtained for the GK phase function.²⁸ In contrast, the CPD of MHG and Mie phase functions can only be derived numerically, and thus, a lookup table-based approach was used for sampling the scattering angle.²⁹

2.3.

Inverse Monte Carlo Model

The IMC model based on the lookup table approach introduced by Palmer et al. and Hennessy et al.^30,31 was used to estimate the optical properties from the acquired SRR profiles $\mathbf{R}=R({\mathrm{SDS}}_1),R({\mathrm{SDS}}_2),\dots ,R({\mathrm{SDS}}_5)$, where $R({\mathrm{SDS}}_i)$ is the reflectance acquired at the $i$’th SDS. In both studies, the lookup tables were introduced by only considering ${\mu }_{\mathrm{a}}$ and ${\mu }_{\mathrm{s}}^{\prime }$. A similar basic IMC (b-IMC) methodology was adopted in the initial part of this study. In order to comply with biological tissue characteristics³², the lookup table values for each SDS ranged from 0 to $25{\mathrm{cm}}^{-1}$ for ${\mu }_{\mathrm{a}}$ and from 5 to $70{\mathrm{cm}}^{-1}$ for ${\mu }_{\mathrm{s}}^{\prime }$ using 30 uniform steps. The b-IMC model was based on the HG phase function with the anisotropy factor $g$ set to 0.8. A constant value of $g$ was selected in accordance with the first similarity relation for ${\mu }_{\mathrm{s}}^{\prime }$, which states that regardless of the individual values of $g$ and scattering coefficient ${\mu }_{\mathrm{s}}$, the same values of ${\mu }_{\mathrm{s}}^{\prime }$ result in the same reflectance.^33–35

It has been shown, however, that higher similarity relations are important when the SDS of the optical probe becomes comparable to the transport mean free path length.¹² By following the lookup table methodology, the b-IMC model was extended (e-IMC) by an additional similarity parameter $\gamma$ that was proposed by Bevilacqua and Depeursinge.¹² The e-IMC model lookup tables included a third dimension for $\gamma$, which spanned from 1.6 to 2.3 in 20 uniform steps, taking into account the biological variations.^16,22 The e-IMC model was based on the GK phase function that allowed the full span of $\gamma$ values at $g=0.8$. On the other hand, the MHG and HG phase functions are less useful, since they only support $\gamma$ values of up to 2, and HG phase function does not allow independent selection of $g$ and $\gamma$.

The simplified geometry (SG) was utilized by the IMC models when applied to the synthetic $\mathbf{R}$. In contrast, the realistic geometry (RG) was used, when the IMC models were applied to the measured $\mathbf{R}$.

As the number of independent measurements used in this study is relatively low (five different SDS), formulation of a new CF is essential for robust estimation of optical properties. In this study, two CFs were compared

2.4.

Synthetic Datasets of Spatially Resolved Reflectance Profiles

To objectively evaluate the performance of the IMC models, synthetic datasets $\mathcal{R}$ of SRR profiles $\mathbf{R}$ were computed according to the parameters in Table 1. Additionally, the synthetic datasets included computed $\mathbf{R}$ of turbid phantoms, later used for validation (see Sec. 2.6, Table 2). Each phantom comprised 150 synthetic $\mathbf{R}$ for wavelengths from 450 to 800 nm.

Table 1

Optical properties used for synthetic datasets R of SRR profiles R.

Table 2

Microsphere number density n, reduced scattering coefficient μs′ and absorption coefficient μa at 630 nm for each phantom P comprising 0.51- and 0.99-μm diameter polystyrene microspheres and, for the absorbing phantoms, also green molecular dye.

It should be noted that an arbitrary combination of $g$ and $\gamma$ cannot be obtained for all the used phase functions. For HG phase function, the relation between $g$ and $\gamma$ (i.e., $\gamma =g+1$) results in only six valid combinations of $g$ and $\gamma$ (Table 1, second line). In the case of MHG phase function, a mathematical constrain (i.e., $\gamma <1+g$) leads to only 14 combinations (Table 1, third line). Lastly, combinations of $g$ and $\gamma$, (0.75, 2.25) and (0.95, 1.65), are not valid for the GK phase function. Thus, only 18 valid combinations of $g$ and $\gamma$ can be obtained (Table 1, fifth line). The datasets ${\mathcal{R}}_{\mathrm{HG}\text{-}\mathrm{I}}$ and ${\mathcal{R}}_{\mathrm{GK}\text{-}\mathrm{I}}$ were introduced for the performance verification of the b-IMC and e-IMC inverse models, respectively. For these datasets, the error of estimated parameters should be close to zero, given a proper formulation of the CF is employed.

2.5.

Experimental Setup for Reflectance Measurements

SRR profiles in the subdiffusive regime were acquired by a custom-made linear array optical fiber probe (FiberTech Optica Inc., Ontario, Canada) with five SDSs (220, 440, 660, 880, and $1100\mu \mathrm{m}$) and an outer diameter of 6.0 mm (Fig. 1). The diameter of the optical fiber cores and the numerical aperture of the fibers was $200\mu \mathrm{m}$ and 0.22, respectively. The optical fiber probe was placed in a holder, which enabled repeatable repositioning of the probe within the cylindrical beaker containing the liquid turbid phantom. The cylindrical beaker with a diameter of 19.5 mm and a height of 30 mm was internally coated with a black matte paint to reduce reflections from the container walls and eliminate stray light pollution. Additionally, the optical fiber probe was held at a position far away from the beaker walls to satisfy the semi-infinite medium assumption for the turbid phantoms. A broadband halogen light source (AvaLight-Hal LS, Avantes, The Netherlands) was coupled to the source fiber. The collected light from each of the five detector fibers was transmitted to a custom-made multiplexer, which comprised a motorized linear stage that aligned the selected detector fiber with a fiber leading to the spectrometer (AvaSpec-2048TEC-FT, Avantes). The recorded signal from each detector fiber was corrected for the sensor dark current and normalized by the Spectralon white signal (Avantes). The SRR profiles were acquired in the wavelength range from 450 to 800 nm.

2.6.

Turbid Phantoms

Turbid phantoms with known optical properties are commonly used to objectively evaluate experimental setups and computational methodologies in biomedical optics.³⁷ In this study, we prepared 24 water-based turbid phantoms comprising a mixture of absorbers and scatterers. A green molecular dye found in fountain pen inks (Live Line Green) was used as the absorber. The molecular dye was thoroughly tested for its stability through several weeks and proved to be stable even when exposed to direct sunlight.³⁸ The scattering component of the dye was negligible due to its molecular nature. A 1-mm cuvette and a cuvette holder (CVH100/M, Thorlabs) were used to measure the absorption coefficient ${\mu }_{\mathrm{a}}$ of the molecular dye (without the scatterer) by observing the attenuation of a collimated light beam through the cuvette. The absorption coefficient was computed according to the Beer–Lambert law. Polystyrene microspheres (diameter $0.51\pm 0.01$ and $0.99\pm 0.03\mu \mathrm{m}$, Polysciences Inc.) were used as the scatterer. The respective ${\mu }_{\mathrm{s}}$ and phase functions were calculated according to the Mie theory and the narrow size distribution of the microspheres (nearly monodisperse) provided by the manufacturer. The wavelength dependence of the polystyrene refractive index was taken from Nikolov et al.³⁹

The prepared turbid phantoms were divided into two groups. The first group of phantoms ${\mathrm{P}}_1$ to ${\mathrm{P}}_{12}$ comprised only scatterers (phantoms ${\mathrm{P}}_1$ to ${\mathrm{P}}_6$ comprised $0.51$ and phantoms ${\mathrm{P}}_7$ to ${\mathrm{P}}_{12}$ comprised $0.99\text{-}\mu \mathrm{m}$ microspheres, Table 2). The second group of phantoms ${\mathrm{P}}_{13}$ to ${\mathrm{P}}_{24}$ comprised scatterers and green molecular dye (phantoms ${\mathrm{P}}_{13}$ to ${\mathrm{P}}_{18}$ comprised 0.51 and phantoms ${\mathrm{P}}_{19}$ to ${\mathrm{P}}_{24}$ comprised $0.99\text{-}\mu \mathrm{m}$ microspheres, Table 2). The values of ${\mu }_{\mathrm{a}}$ for the second group of phantoms ranged between 8.2 and $10.3{\mathrm{cm}}^{-1}$ at 630 nm (absorption peak). The microsphere number density $n$ and corresponding ${\mu }_{\mathrm{a}}$ and ${\mu }_{\mathrm{s}}^{\prime }$ for all the prepared turbid phantoms are listed in Table 2.

3.1.

Cost Function Evaluation

The performance and robustness of the proposed CFs were evaluated on the synthetic dataset ${\mathcal{R}}_{\mathrm{HG}\text{-}\mathrm{I}}$ by deploying the b-IMC model. The b-IMC model was selected for its dependence on only two parameters (${\mu }_{\mathrm{a}}$ and ${\mu }_{\mathrm{s}}^{\prime }$) and, hence, easier visualization of the CF (Fig. 2). The presented example ($\mathbf{R}$ from the dataset ${\mathcal{R}}_{\mathrm{HG}\text{-}\mathrm{I}}$ with the true values of ${\mu }_{\mathrm{a}}=12.5{\mathrm{cm}}^{-1}$ and ${\mu }_{\mathrm{s}}^{\prime }=18.9{\mathrm{cm}}^{-1}$) clearly shows that the CF significantly affects the estimation of optical properties by the IMC models. The ${\mathrm{CF}}_2$ [Eq. (4)] exhibits improved localization of the minimum with respect to the ${\mathrm{CF}}_1$. Moreover, the relative root-mean-square (RMS) error obtained by the b-IMC model using the ${\mathrm{CF}}_1$ for the entire dataset ${\mathcal{R}}_{\mathrm{HG}\text{-}\mathrm{I}}$ is 2.7% and 2.8% for ${\mu }_{\mathrm{a}}$ and ${\mu }_{\mathrm{s}}^{\prime }$, respectively. In contrast, the relative RMS error obtained by the ${\mathrm{CF}}_2$ was one order of magnitude lower, i.e., 0.32% and 0.33% for ${\mu }_{\mathrm{a}}$ and ${\mu }_{\mathrm{s}}^{\prime }$, respectively. For the sake of completeness, the performance of the two CFs was also compared by the e-IMC model. As with the b-IMC model, the relative RMS errors obtained by the e-IMC model based on the ${\mathrm{CF}}_1$ (5.6%, 4.4%, and 4.7% for ${\mu }_{\mathrm{a}}$, ${\mu }_{\mathrm{s}}^{\prime }$, and $\gamma$, respectively) were one order of magnitude higher than the errors obtained by the e-IMC model based on the ${\mathrm{CF}}_2$ (0.60%, 0.58%, and 0.76% for ${\mu }_{\mathrm{a}}$, ${\mu }_{\mathrm{s}}^{\prime }$, and $\gamma$, respectively). The obtained results clearly show the superiority of the ${\mathrm{CF}}_2$ over ${\mathrm{CF}}_1$. Consequently, ${\mathrm{CF}}_2$ was used in all subsequent calculations involving the b-IMC and e-IMC models.

Fig. 2

CF dependence on the values of optical properties in terms of absorption ${\mu }_{\mathrm{a}}$ and reduced scattering ${\mu }_{\mathrm{s}}^{\prime }$ coefficients for (a) ${\mathrm{CF}}_1$ from Eq. (3) and for (b) ${\mathrm{CF}}_2$ from Eq. (4).

3.2.

Influence of the Phase Function on the Subdiffusive Reflectance and on the Inverse Monte Carlo Models

Subdiffusive reflectance is known to significantly depend on the phase function. As explained in Sec. 1, the similarity parameter $\gamma$ has been suggested to account for some of the subdiffusive reflectance variability that arises from the phase functions. In other words, for a constant ${\mu }_{\mathrm{a}}$ and ${\mu }_{\mathrm{s}}^{\prime }$, the phase function-dependent changes in the subdiffusive reflectance should be reduced when the additional similarity parameter $\gamma$ is kept constant as well.

In order to study the beneficial effect of $\gamma$, we introduced relative variability maps of the subdiffusive reflectance as a function of ${\mu }_{\mathrm{a}}$ and ${\mu }_{\mathrm{s}}^{\prime }$ at a particular SDS. The relative variability maps were derived as a ratio between the standard deviation and the mean of the reflectance maps across different phase functions. Figure 3 shows the relative variability maps at three SDS of 220, 660, and $1100\mu \mathrm{m}$. In the top row (constant $g$), the HG, MHG, and GK phase functions were used with $g$ set to 0.8 and $\gamma$ varied from 1.7 to 1.9. In the bottom row (constant $\gamma$), the HG, MHG, and GK phase functions were used with $\gamma$ set to 1.7 and $g$ varied from 0.7 to 0.9. In comparison to $g$, $\gamma$ significantly reduces the variability of the reflectance for ${\mu }_{\mathrm{s}}^{\prime }$ exceeding $\sim 15{\mathrm{cm}}^{-1}$ when different phase functions are used. Moreover, although $g$ is frequently described in the literature as an important observable of the angular scattering distribution in turbid media, the subdiffusive reflectance is more dependent on the changes of $\gamma$ than $g$. Similar conclusions were recently made by Calabro and Bigio¹⁵ and Bodenschatz et al.⁴⁰ The reason for the increased sensitivity of the subdiffusive reflectance to the phase function observed at the low reduced scattering coefficients can be attributed to the longer photon mean free path lengths. In this case, the photons that contribute to the reflectance signal exhibit only a few scattering events, with one of those events likely to be a large-angle deflection. Therefore, reflectance is highly dependent on the large-angle section of the phase function, which is only approximately accounted for by $\gamma$. The effect can be observed in the reflectance variability map in Fig. 3 (bottom row) computed for phase functions that despite having a common $\gamma$ still exhibit significant differences in the large-angle section. To account for this shortcoming, Bodenschatz et al.⁴⁰ have proposed a parameter $\sigma$, which reduced the variability of reflectance for a constant ${\mu }_{\mathrm{a}}$ and ${\mu }_{\mathrm{s}}^{\prime }$ and different phase functions in the spatial frequency domain. However, in the same study, $\sigma$ was not found to significantly reduce the variability in the spatially resolved domain in comparison to $\gamma$. We believe that the performance of the e-IMC model for ${\mu }_{\mathrm{s}}^{\prime }$ below $\sim 15{\mathrm{cm}}^{-1}$ could be significantly improved by introducing higher order similarity parameters beyond $\gamma$,⁴¹ e.g., $\delta =(1-g_3)/(1-g_1)$,⁴² where $g_3$ is the third Legendre moment of the phase function. Keeping the values of $\gamma$ and $\delta$ constant should in principle significantly reduce the relative variability of subdiffusive reflectance computed for different phase functions at any given value of ${\mu }_{\mathrm{a}}$ and ${\mu }_{\mathrm{s}}^{\prime }$.

Fig. 3

Relative variability maps of the subdiffusive reflectance at (a) $220\mu \mathrm{m}$, (b) $660\mu \mathrm{m}$, and (c) $1100\mu \mathrm{m}$ SDS for different phase functions, as a function of the absorption ${\mu }_{\mathrm{a}}$ and reduced scattering ${\mu }_{\mathrm{s}}^{\prime }$ coefficients. In the top row (constant $g$), HG, MHG, and GK phase functions were used with $g=0.8$ and $\gamma$ varied from 1.7 to 1.9. In the bottom row (constant $\gamma$), HG, MHG, and GK phase functions were used with $\gamma =1.7$ and $g$ varied from 0.7 to 0.9.

Another interesting observation can be made at the smallest SDS, where a region of significantly decreased relative reflectance variability occurs. In this region, the reflectance seems unaffected by the choice of the phase function as long as ${\mu }_{\mathrm{a}}$ and ${\mu }_{\mathrm{s}}^{\prime }$ are kept constant. This observation, named as the “isobestic point,” was also made by Calabro and Bigio,¹⁵ and it is believed to occur at ${\mu }_{\mathrm{s}}^{\prime }$ $\mathrm{SDS}\approx 0.7$. In this study, the SDS of $220\mu \mathrm{m}$ evaluates to an isobestic point of $32{\mathrm{cm}}^{-1}$, which is consistent with the results in Fig. 3.

To evaluate the effect of the phase function on the estimation of optical properties, the synthetic dataset ${\mathcal{R}}_{\mathrm{GK}\text{-}\mathrm{I}}$ was analyzed by the b-IMC and e-IMC models. The ${\mathcal{R}}_{\mathrm{GK}\text{-}\mathrm{I}}$ dataset is based on the GK phase function with the value of $g$ set to 0.8 and $\gamma$ varied from 1.65 to 2.25, which is similar to the conditions used for producing the results in the top row of Fig. 3, where $g$ was constant. Figure 4 (top row) shows the relative errors obtained for the estimated ${\mu }_{\mathrm{a}}$ and ${\mu }_{\mathrm{s}}^{\prime }$ by the b-IMC model from the dataset ${\mathcal{R}}_{\mathrm{GK}\text{-}\mathrm{I}}$. The values of ${\mu }_{\mathrm{a}}$ and ${\mu }_{\mathrm{s}}^{\prime }$ used in ${\mathcal{R}}_{\mathrm{GK}\text{-}\mathrm{I}}$ are color-coded and can be deduced from the corresponding colorbars. It can be observed that the accuracy of the estimated optical properties decreases significantly due to the variations in $\gamma$. The relative errors are the lowest for $\gamma =1.8$, where the GK phase function simplifies to HG phase function. For all the other values of $\gamma$, the performance is significantly degraded. In contrast, a significant improvement in the accuracy of the optical properties estimated by the e-IMC model from the same dataset ${\mathcal{R}}_{\mathrm{GK}\text{-}\mathrm{I}}$ can be observed in Fig. 4 (bottom row). The relative errors of the estimated optical properties obtained by the e-IMC model are mostly within 2% for all of the selected values of $\gamma$.

Fig. 4

Relative errors of (a) the absorption ${\mu }_{\mathrm{a}}$ and (b) reduced scattering ${\mu }_{\mathrm{s}}^{\prime }$ coefficients estimated by the b-IMC model (top row) and e-IMC model (bottom row) for the synthetic reflectance dataset ${\mathcal{R}}_{\mathrm{GK}\text{-}\mathrm{I}}$. The true values of ${\mu }_{\mathrm{a}}$ and ${\mu }_{\mathrm{s}}^{\prime }$ used in ${\mathcal{R}}_{\mathrm{GK}\text{-}\mathrm{I}}$ are color coded and can be deduced from the colorbars.

The results of b-IMC model presented in Fig. 4 clearly suggest that using only the first similarity relation will lead to large relative errors of the estimated optical properties. For example, although the dataset ${\mathcal{R}}_{\mathrm{GK}\text{-}\mathrm{I}}$ is simulated using a phase function with a constant $g$, the changes in the reflectance due to changes in phase function arising from $\gamma$ cannot be sufficiently accounted for by the first similarity relation. This observation could explain one of the possible sources of errors in studies that utilized b-IMC-like models for optical fiber probes with similar SDS as used in this study.^31,43,44 We believe that the listed studies could be improved by using additional parameters such as $\gamma$ in addition to ${\mu }_{\mathrm{a}}$ and ${\mu }_{\mathrm{s}}^{\prime }$ in the IMC models to account for some of the phase function variability that influences the reflectance at small SDS.

3.3.

Performance of the Basic and Extended Inverse Monte Carlo Models on Synthetic Datasets of Spatially Resolved Reflectance Profiles

To further examine the effect of $\gamma$ on the estimated optical properties, the b-IMC and e-IMC models were extensively evaluated and compared on synthetic datasets of $\mathbf{R}$ (Sec. 2.5, Table 1). The variability of phase functions captured by the synthetic datasets should account for the variability observed in biological tissues. Consequently, similar errors can be expected for measured $\mathbf{R}$ of biological tissues.

Figure 5 shows the values of ${\mu }_{\mathrm{a}}$ and ${\mu }_{\mathrm{s}}^{\prime }$ estimated by the b-IMC and e-IMC models with respect to the true values used to compute the synthetic datasets. The performance of the e-IMC model is evidently superior to the b-IMC model, since the values obtained by the b-IMC model are significantly more dispersed along the ideal relationship line between the true and estimated parameter values. Moreover, the e-IMC model offers an additional estimate of the similarity parameter $\gamma$ [Fig. 5(e)]. Quantitative results obtained for each synthetic dataset in terms of the two IMC models are summarized in Table 3. As expected, both IMC models performed with relative RMS errors well below 1% when tested on the datasets ${\mathcal{R}}_{\mathrm{HG}\text{-}\mathrm{I}}$ and ${\mathcal{R}}_{\mathrm{GK}\text{-}\mathrm{I}}$, respectively. The two datasets were obtained by the same phase functions that were used to create the lookup tables. The performance of the b-IMC model clearly degrades for the dataset ${\mathcal{R}}_{\mathrm{HG}\text{-}\mathrm{II}}$ with $g$ varying between 0.7 and 0.95. In contrast, the e-IMC model does not seem to be significantly affected by the variation in $g$, which is in agreement with the results in Fig. 3 (bottom row). The e-IMC model performs similarly on the ${\mathcal{R}}_{\mathrm{GK}\text{-}\mathrm{II}}$ dataset. However, the RMS errors of the IMC models further increase when applied to the MHG and Mie phase function-based datasets (${\mathcal{R}}_{\mathrm{MHG}}$, pure scattering ${\mathrm{P}}_1$ to ${\mathrm{P}}_{12}$, and scattering and absorbing ${\mathrm{P}}_{13}$ to ${\mathrm{P}}_{24}$ datasets). Once more, the b-IMC model is much more affected by the changes in the phase function and exhibits twice the relative RMS error of the e-IMC model. The quantitative advantage of the e-IMC model over the b-IMC model can be observed in Table 3 that summarizes the results across all the datasets (all datasets).

Fig. 5

The values of (a) ${\mu }_{\mathrm{a}}$ and (b) ${\mu }_{\mathrm{s}}^{\prime }$ estimated by the b-IMC model and the values of (c) ${\mu }_{\mathrm{a}}$, (d) ${\mu }_{\mathrm{s}}^{\prime }$, and (e) $\gamma$ estimated by the e-IMC model with respect to the corresponding true values that were used to compute the synthetic datasets of SRR profiles. The ideal relationship is denoted by a black solid line. The synthetic datasets are presented with different colors and markers (see the legend).

Table 3

The RMS errors of the estimated optical properties obtained by the b-IMC and e-IMC models on synthetic datasets of SRR profiles R.

Despite the improved result, the performance of e-IMC model on the dataset of scattering and absorbing phantoms ${\mathrm{P}}_{13}$ to ${\mathrm{P}}_{24}$ is still significantly degraded. Previous results (bottom row in Fig. 3) point out that $\gamma$ does not fully eliminate the variability of reflectance for low reduced scattering coefficients (especially below $\sim 15{\mathrm{cm}}^{-1}$), when different phase functions are used. Since Fig. 3 only includes HG, MHG, and GK phase functions, we also constructed similar relative subdiffusive reflectance variability maps that depend only on the GK and Mie phase functions. Specifically, Fig. 6 shows the relative variability map for three different phase functions with $\gamma$ set to 2.1. The first-phase function was GK with $g=0.8$, which was also used to create the e-IMC lookup table. The second-phase function was Mie for $0.51\text{-}\mu \mathrm{m}$ polystyrene spheres at 705.8 nm, which yields $g=0.79$. Finally, the third-phase function was Mie for $0.99\text{-}\mu \mathrm{m}$ polystyrene spheres at 499.5 nm, which yields $g=0.93$. Figure 6 highly resembles the bottom row of Fig. 3, thus again suggesting that $\gamma$ cannot fully account for the variability of the phase function for low ${\mu }_{\mathrm{s}}^{\prime }$. Moreover, since the GK phase function is used in the e-IMC model, an SRR profile simulated by a different phase function with the same value of $\gamma$ could be different due to the phase function properties that cannot be accounted for by the first two Legendre moments encapsulated in $\gamma$. These differences between the e-IMC model and the SRR profiles then propagate into errors of the estimated ${\mu }_{\mathrm{a}}$, ${\mu }_{\mathrm{s}}^{\prime }$, and $\gamma$. The observation is further confirmed by removing the phantoms ${\mathrm{P}}_1$, ${\mathrm{P}}_7$, ${\mathrm{P}}_{13}$, ${\mathrm{P}}_{19}$, that exhibited ${\mu }_{\mathrm{s}}^{\prime }$ below $15{\mathrm{cm}}^{-1}$ over the majority of the wavelength range, from the summarized RMS errors ($\mathrm{All}-\{{\mathrm{P}}_1,{\mathrm{P}}_7,{\mathrm{P}}_{13},{\mathrm{P}}_{19}\}$; Table 3). As a result, the summarized RMS errors drop significantly and are now $0.12{\mathrm{cm}}^{-1}$ (3.2%), $0.56{\mathrm{cm}}^{-1}$ (2.0%), and 0.049 (2.3%) for ${\mu }_{\mathrm{a}}$, ${\mu }_{\mathrm{s}}^{\prime }$, and $\gamma$, respectively.

Fig. 6

The relative variability maps of the subdiffusive reflectance at (a) $220\mu \mathrm{m}$, (b) $660\mu \mathrm{m}$, and (c) $1100\mu \mathrm{m}$ SDS for different phase functions at a constant $\gamma =2.1$. The employed phase functions included the GK phase function with $g$ set to 0.8 (the value employed by the e-IMC lookup table) and Mie phase functions for 0.51- and $0.99\text{-}\mu \mathrm{m}$ polystyrene microspheres at 705.8 and 499.5 nm, respectively.

3.4.

Performance of the Extended Inverse Monte Carlo Model on Measured Spatially Resolved Reflectance Profiles

Finally, we evaluated the e-IMC model on measured $\mathbf{R}$ acquired from turbid phantoms (Sec. 2.6). To compare the measured reflectance to the MC simulated reflectance at a particular SDS, a calibration procedure was required. Briefly, the MC simulated reflectance is normalized against the initial number of launched photon packets, while the measured reflectance is normalized against a reflectance standard. Consequently, to account for the reflectance properties of the standard, the reflectance of the turbid phantoms with known optical properties is measured and simulated. The calibration factor is introduced as the ratio between the MC simulated and measured reflectance. Ideally, the calibration factor for a given wavelength and SDS should remain the same across all the turbid phantoms. In this study, the purely scattering phantoms ${\mathrm{P}}_1$ to ${\mathrm{P}}_{12}$ were used for computing the calibration factors. The variations of the calibration factors at all spectral bands and SDS were within 2%, and mostly arising from the variations in the measurements.

The performance of the e-IMC model was evaluated on the remaining turbid phantoms ${\mathrm{P}}_{13}$ to ${\mathrm{P}}_{24}$ with nonzero ${\mu }_{\mathrm{a}}$. The estimated values of the optical properties with respect to the corresponding true values are shown in Fig. 7 and are similar to those in Figs. 5(c)–5(e). Moreover, the RMS errors of ${\mu }_{\mathrm{a}}$, ${\mu }_{\mathrm{s}}^{\prime }$, and $\gamma$ obtained for the turbid phantoms are $0.38{\mathrm{cm}}^{-1}$ (15.6%), $0.71{\mathrm{cm}}^{-1}$ (5.8%), and 0.10 (4.8%), respectively, and tightly follow the results obtained for the corresponding synthetic datasets (Table 4). In addition, by removing the turbid phantoms ${\mathrm{P}}_{13}$ and ${\mathrm{P}}_{19}$ (${\mu }_{\mathrm{s}}^{\prime }$ below $15{\mathrm{cm}}^{-1}$) from the summarized RMS errors, nearly a twofold improvement in terms of relative RMS error is attained. Consequently, the RMS errors stay within the 10% margin, however, only for a subregion of optical properties where ${\mu }_{\mathrm{s}}^{\prime }$ exceeds $\sim 15{\mathrm{cm}}^{-1}$. The obtained results in this section are consistent with the results shown in Figs. 3 and 6, where $\gamma$ accounts for the phase functions variability only for ${\mu }_{\mathrm{s}}^{\prime }$ values exceeding $\sim 15{\mathrm{cm}}^{-1}$.

Fig. 7

The values of (a) ${\mu }_{\mathrm{a}}$, (b) ${\mu }_{\mathrm{s}}^{\prime }$, and (c) $\gamma$ estimated by the e-IMC model from the measured $\mathbf{R}$ of turbid phantoms with respect to the corresponding true values. The ideal relationship is denoted by a black solid line. Different phantoms are presented with different colors (see the legend).

Table 4

The RMS errors of the estimated optical properties obtained by the e-IMC models on synthetic (left) and measured (right) SRR profiles of turbid phantoms.

The related studies from Hennessy et al.³¹ and Rajaram et al.⁴³ estimated ${\mu }_{\mathrm{a}}$ and ${\mu }_{\mathrm{s}}^{\prime }$ with relative RMS errors of 0.74% and 1.74%, and 11.6% and 5.9%, respectively. Both studies used optical fiber probes with small SDS between 250 and $300\mu \mathrm{m}$. According to these results, the proposed e-IMC model gives slightly higher relative errors for the turbid phantoms ${\mathrm{P}}_{13}$ to ${\mathrm{P}}_{24}$. The results by Hennessy et al. are especially surprising since the employed IMC model was based on the HG phase function at a constant $g=0.85$ (similar to our b-IMC model), while the validation was performed on turbid phantoms containing polystyrene microspheres that follow the Mie phase function. We have shown that the reflectance at small SDS and constant $g$ can vary significantly, if different phase functions are used (Fig. 3). Moreover, we have observed a significant degradation of the b-IMC model on synthetic datasets (Table 3). One of the possible causes for a better performance of the IMC models in the related studies is the use of spectral constrains in terms of tissue chromophore concentrations and spectral dependence of the ${\mu }_{\mathrm{s}}^{\prime }$. Spectral constrains were required because only one SDS was utilized, thus the full reflectance spectrum was used to estimate only a few parameters. This, however, significantly reduced the number of optimized parameters in comparison to the number of available measurements. It is important to note that unlike the IMC models presented in the related studies, the e-IMC model does not require spectral constrains. In this way, the e-IMC model requires no prior knowledge of the chromophores and scatterers in the turbid medium to estimate ${\mu }_{\mathrm{a}}$, ${\mu }_{\mathrm{s}}^{\prime }$, and $\gamma$. As such, the e-IMC can be utilized only at particular wavelengths. In addition, spectral constrains can be easily introduced into the e-IMC model. We believe that this should in principle further reduce the relative error of the estimated quantities since there would be significantly more measurements and less parameters to optimize. For human tissues, this would most commonly include modeling ${\mu }_{\mathrm{a}}$ by the absorption spectra of chromophores, such as oxygenated and deoxygenated hemoglobin, melanin, carotenoids, while ${\mu }_{\mathrm{s}}^{\prime }$ can be modeled as ${\mu }_{\mathrm{s}}^{\prime }(\lambda )=A{(\lambda /{\lambda }_0)}^{-B}$.³² The spectral dependence of $\gamma$ in tissue is close to constant across the visible spectrum and could thus be adequately modeled by a linear or quadratic function.^16,45

Finally, the turbid phantoms used in this study to experimentally test the e-IMC model were based on polystyrene microspheres, which follow the Mie phase function. Polystyrene microspheres are very practical to use since they are available in standardized nearly monodisperse suspensions. In this way, the scattering coefficient and the phase function can easily be calculated. However, monodisperse solutions do not exactly mimic the tissue optical properties in terms of the shape of the phase function. As a result, the accuracy of the estimated optical properties for tissues can somewhat differ from the accuracy obtained for such turbid phantoms. Since biological tissues exhibit phase functions similar to GK or MHG, more accurate results can be expected by the e-IMC model (similar to the synthetic dataset ${\mathcal{R}}_{\mathrm{GK}\text{-}\mathrm{II}}$ or ${\mathcal{R}}_{\mathrm{MHG}}$ in Table 3).