1. Introduction

Through the years diffuse reflectance (DR) spectroscopy has proved itself as a useful method for non-invasive study of biological tissues. The reflectance spectra can be efficiently collected with optical fiber probes [1] or imaging systems [2] and sufficiently described through the framework of radiative transport equation (RTE) [3]. For particular experimental settings [4,5] or large source-detector separations with dominant scattering [6] the RTE can be solved analytically. However, for complex experimental settings, the RTE needs to be solved numerically, which is commonly acomplished by the Monte Carlo (MC) method. Due to the stochastic nature of the MC method, a large number of photon packets has to be processed to obtain adequate quality of the simulated quantity. Since such simulations can quickly become time-consuming, the structure and geometry of the tissue and optical fiber probes is frequently described with some simplifications. The optical fiber probe tip model (PTM) is often approximated (Fig. 1) by a laterally uniform boundary with mismatched refractive indices [7–10], which allows radially symmetric detection schemes that lead to shorter simulation times. This approach can be valid in some cases, but in general, it could lead to significant errors in the MC simulated reflectance spectra [1,11,12]. To overcome these limitations, additional parameters of the fiber probe tip, such as the reflectivity of the stainless steel and refractive index of the black epoxy fill has to be taken into account. Unfortunately, using the full details of the probe tip in the MC simulations breaks down the radial symmetry of the detection scheme allowing collection of the photon packets only through the individual optical fibers. Consequently, the simulation time required to obtain a good quality spatially resolved reflectance (SRR) can quickly increase by more than two orders of magnitude.

 

Fig. 1. Simplified PTM (left) considers only mismatch between the refractive indices of the sample and optical fibers, while realistic PTM (right) takes into account all the details of the optical fiber probe tip.

Download Full Size | PDF

In this paper, we propose a deep artificial neural network-based (ANN) reflectance regression model (RRM) that maps reflectance, simulated with a simplified PTM, to reflectance, valid for realistic PTM. In Section 2, we present the light propagation model, reflectance datasets, inverse models for estimation of optical properties, the deep ANN-based RRM, experimental setup and description of turbid phantoms. Section 3 summarizes the performance evaluation of RRM, validation through the use of inverse models for estimation of optical properties and experimental validation by liquid turbid phantoms.

2. Materials and methods

2.1 Workflow

The basic idea of the methodology presented in this paper is summarized in Fig. 2. The SRR obtained by fast MC simulations utilizing a simplified PTM, which inaccurately describes the optical fiber probe tip and hence introduces errors into the simulated SRR, is rapidly transformed through the RRM to obtain accurate SRR. The accurate SRR can be also directly computed with MC simulations that utilize a realistic PTM, however such MC simulations are extremely slow. Detailed implementation of the RRM is discussed in Section 2.4 and the accuracy of the RRM is presented in Fig. 5. For additional numerical and experimental validation, the regressed reflectance was used to prepare inverse models for estimation of optical properties. Implementation details of the inverse models are presented in Section 2.5 and accuracy of the estimated optical properties is discussed in Section 3. A detailed workflow of the methodology that includes the RRM and inverse models is presented in Fig. 3.

 

Fig. 2. Basic workflow of the proposed methodology. SRR that can be efficiently simulated with the simplified PTM but consequently subject to significant simulation errors is rapidly regressed to accurate SRR that could be directly obtained by much slower Monte Carlo simulations utilizing the realistic PTM. The regressed reflectance is used to substantially speed-up (up to 480 fold) the preparation of inverse models for estimation of optical properties from measured and calibrated SRR.

Download Full Size | PDF

 

Fig. 3. Illustration of the framework that includes training of the ANN-based RRMs that rapidly map SRR from the simplified PTM to the realistic PTM and training of the inverse models for estimation of optical properties from the SRR valid for the simplified or realistic PTMs. The $\textrm {R}_{\textrm {S}}$ and $\textrm {R}_{\textrm {R}}$ stand for reflectance at one source-detector separation, valid for the simplified and realistic PTMs, respectively. The $\textrm {SRR}_{\textrm {S}}$ and $\textrm {SRR}_{\textrm {R}}$ stand for the spatially resolved reflectance valid for the simplified and realistic PTM, respectively. The $\textrm {OP}$ stands for optical property. (left) The SRRs simulated with the simplified PTM are passed through the inverse models for estimation of optical properties from the SRR valid for the simplified PTM. (right) The SRRs simulated with the realistic PTM are passed through the inverse models for estimation of optical properties from the SRR valid for the realistic PTM. Ideally, given that the RRM does not introduce additional errors into the regressed SRR, the estimation errors of the two inverse models should be similar and hence the inverse model of the simplified PTM serves as a performance baseline.

Download Full Size | PDF

2.2 Light propagation model

In this study, the light propagation in sample was modeled by an experimentally validated OpenCL^TM implementation of the MC method for a layered medium that follows the existing open source implementations [13,14]. All simulations of reflectance were prepared for an optical fiber probe with an outer diameter of 6.1 mm and six linearly placed optical fibers with a diameter of 200 µm and core-to-core spacing of 220 µm. Due to small source-detector separation (SDS), we took into account the so-called subdiffusive quantifiers of the scattering phase function (SPF), namely $\gamma$ [15] and $\delta$ [16], as it has been shown, that they are important for accurate description of reflectance close to the source. Parameter $\gamma$ is defined as a function of the first $g_1$ and second $g_2$ Legendre moments of SPF $\gamma = (1-g_2)/(1-g_1)$ and $\delta$ is defined as a function of the first $g_1$ and third $g_3$ Legendre moments of SPF $\delta = (1-g_3)/(1-g_1)$. The first Legendre moment $g_1$ is commonly termed as the anisotropy factor $g$. Bodenschatz et al. demonstrated that higher-order Legendre moments carry information of large-angle scattering [17]. Two different approaches were used in the MC simulations to formulate the PTM that followed the work of Naglič et al. [18] (Fig. 1).

The first, simplified PTM, assumed a laterally uniform interface that took into account only the mismatch between the refractive indices of the sample ($n_{samp}$=1.33) and optical fibers ($n_{fib}$=1.452) [9]. Each MC simulation included at least $10^{8}$ photon packets that were launched uniformly with respect to the numerical aperture (NA=0.22) and diameter (200 µm) of the optical fiber core. To attain an acceptable and consistent signal-to-noise ratio of the SRR collected through the optical fibers, the MC simulations were terminated when the total weight of the photon packets, collected through the optical fiber with the largest SDS exceeded 1000. Due to the lateral uniformity of the PTM, the backscattered light was collected through a radially symmetric detection scheme (annual rings) centered at the position of the source optical fiber. The SRR collected through the 5 remaining fibers was derived from the output of the radially symmetric detection scheme.

The second, realistic PTM, accounted for the reflections from the polished stainless steel (0.5), the refractive index of the optical fibers ($n_{fib}$ = 1.452) [9] and the refractive index of the black epoxy resin ($n_{epoxy}$ = 1.6) (EpoTek 353ND Black, Epoxy Technology Inc., MA, USA), which were found to significantly affect the simulated reflectance [11,19]. Reflectivity of the stainless steel was measured by observing the attenuation of a reflected laser diode beam (670 nm) (Newport Optical Power Meter 1936-C and Thermopile Sensor 818P-001-12) and the reflectivity of the black epoxy resin was calculated by Fresnel’s equations. Once the photon packet hits the polished stainless steel or black epoxy surface of the probe tip, it reflects with a probability that equals reflectivity, otherwise it is considered completely absorbed Since the realistic PTM is not laterally uniform, the radial symmetry of the detection scheme breaks down and hence the reflectance was collected only through the individual optical fibers. To attain a signal-to-noise ratio of the simulated SRR similar to the simplified PTM, the minimum number of launched photon packets was increased by the ratio between the surface area of the most distant optical fiber and the corresponding annual ring of the radially symmetric detection scheme used for the simplified PTM. Consequently, at least $40 \times 10^{8}$ photon packets were launched and the same MC termination criterion was used as for the simplified PTM.

2.3 Datasets

In this subsection we define the SRR datasets that were used to prepare and evaluate the RRM that maps the SRR valid for the simplified PTM to the SRR valid for the realistic PTM and the inverse models for estimation of optical properties. An overview of the datasets is provided in Table 1 and Fig. 3 shows the use of individual datasets within the proposed framework. Datasets denoted with letter T were utilized to train the ANN-based RMMs and inverse models, datasets denoted with letter V were utilized to control (validate) the training process and datasets denoted with letter E were used to test (independently evaluate) the trained ANNs. The range of optical properties was selected based on the values commonly reported for biological tissues [1,20–23]. However, the selected range can be customized to accommodate sample and application-specific requirements.

The points along the absorption coefficient ($\mu _a$) and reduced scattering coefficient ($\mu ^{\prime }_{s}$) axes are equally spaced and form a regular grid in the $\mu _a - \mu ^{\prime }_{s}$ plane. However, the points in the plane formed by the $\gamma$ and $\delta$ subdiffuse SPF quantifiers are not located on a regular grid but are instead uniformly distributed in a way that includes the points that are laying on the boundary of the $\gamma -\delta$ domain. We accomplish this by first uniformly distributing sample points along the $\gamma$ axis within the selected range of $\gamma$ values. At each $\gamma$ sample point we compute intersections with the domain boundary of the GK [24] scattering phase function in the $\gamma -\delta$ plane. The number of uniformly distributed sample points along the $\delta$ axis between the two intersections is selected in a way that keeps the distance between the sample points under a predefined value $d\delta$. Consequently, the total number of points in the $\gamma -\delta$ plane is specified for each dataset. The distributions of sample points for all of the optical properties are presented in Fig. 4. All MC simulations of SRR were based on the Gegenbauer kernel (GK) [24] SPF. In the following, we provide a short description of the three groups of datasets summarized in Table 1.

 

Fig. 4. Visualization of optical properties defined by the datasets summarized in Table 1. Red line shows the valid domain of $\gamma$ and $\delta$ subdiffusive quantifiers for the GK scattering phase function.

Download Full Size | PDF

Table 1. Summary of datasets that includes the total number of SRR samples ($\textrm {N}_{\textrm {SRR}}$), range of the optical properties ($\mu _{a}$, $\mu ^{\prime }_{s}$, $\gamma$ and $\delta$), source of the dataset (MC simulations or ANN-based reflectance regression model (RRM)) and probe tip model (PTM) complexity (simplified (S) or realistic (R)).

View Table | View all tables in this article

2.4 Reflectance regression model (RRM)

The RRMs that map SRR from the simplified PTM to the realistic PTM were based on deep ANNs. One ANN model was trained for each SDS. The five inputs of the ANN-based RRMs were formed by the reflectance value at the selected SDS and the corresponding absorption coefficient $\mu _{a}$, reduce scattering coefficient $\mu ^{\prime }_{s}$ and SPF quantifiers $\gamma$ and $\delta$. The output of each ANN-based RRM was defined as the relative difference between the reflectance values at the selected SDS that were simulated with the simplified and the realistic PTM. The regression ANNs were trained using the $\textrm {T}_{\textrm {RRM},\textrm {S}}$ dataset as the input and the relative difference between $\textrm {T}_{\textrm {RRM},\textrm {S}}$ and $\textrm {T}_{\textrm {RRM},\textrm {R}}$ as the output. A small three-layer deep ANN topology was used with 8, 4 and 1 nodes in the individual layers. A softplus activation function was applied to the first two layers and a hard-sigmoid activation function to the output layer. The training datasets $\textrm {T}_{\textrm {RRM},\textrm {S}}$ and $\textrm {T}_{\textrm {RRM},\textrm {R}}$ were gradually reduced in size as summarized in Table 2 and shown in Fig. 4 to obtain an estimate for the smallest dataset that is required to train sufficiently accurate RRMs.

The training of the ANNs was based on the Adam optimizer [25] using mean squared error cost function and fixed learning rate of $10^{-3}$. Due to the stochastic nature of the training process, 10 RRMs were trained for each SDS and dataset size. Subsequently, the model with the lowest regression error on the $\textrm {E}_{\textrm {RRM}, \textrm {S}}$ dataset was selected.

Table 2. Summary of the decimated $\textrm {T}_{\textrm {RRM},\textrm {S}}$ and $\textrm {T}_{\textrm {RRM},\textrm {R}}$ datasets used to train the ANN-based RRMs that map SRR from the simplified PTM to the realistic PTM.

View Table | View all tables in this article

2.5 Inverse models

Dedicated ANN-based inverse models for estimation of optical properties from SRR were prepared for the simplified PTM and for the realistic PTM. The inverse models were prepared following our methodology presented in [26]. The inverse models for the simplified and realistic PTM were trained using the $\textrm {T}_{\textrm {IM}, \textrm {S}}$ and $\textrm {T}_{\textrm {IM},\textrm {R}}$ datasets, respectively. The inputs of each inverse model were formed by the 5 SRR values collected through the optical fibers and one output represented the estimated optical property. A five-layer deep ANN topology was used with 30, 20, 10, 5 and 1 nodes in the individual layers. A linear activation function was used in the first layer, a hyperbolic tangent activation function was used in the hidden layers and a hard-sigmoid activation function was used in the output layer. The training of the ANNs was based on the Adam optimizer [25] using mean squared error cost function and adaptive learning rate that exponentially decayed with the number of training epochs $10^{-3} \cdot e^{-epoch/1500}$. The training process was limited to 10000 epochs and was stopped early if the value of the cost function on a small validation dataset, $\textrm {V}_{\textrm {IM}, \textrm {S}}$ for the simplified PTM and $\textrm {V}_{\textrm {IM}, \textrm {R}}$ for the realistic PTM, did not sufficiently improve for 500 consecutive training epochs. Due to the stochastic nature of the training process, 5 inverse models were trained for each optical property and PTM. Subsequently, the inverse model with the lowest estimation error was selected. For notation simplicity, the inverse models prepared for the simplified PTM will be referred to as $\textrm {IM}_{\textrm {S}}$ and the inverse models prepared for the realistic PTM will be referred to as $\textrm {IM}_{\textrm {R,n}}$, where n stands for the number of SRRs used to train the underlying RRM that produced the correspondent $\textrm {T}_{\textrm {IM},\textrm {R}}$ dataset (e.g. $\textrm {IM}_{\textrm {R,1215}}$)

2.6 Turbid phantoms

We prepared five liquid turbid phantoms comprising a mixture of absorbers and scatterers following the procedure described in [18]. A green molecular dye found in fountain pen inks (Live Line Green) was used as the absorber. A 10 mm cuvette (Hellma Optik GmbH, Jena, Germany) and a cuvette holder (CUV-UV/VIS, Avantes, Apeldoorn, The Netherlands) were used to measure the absorption coefficient $\mu _{a}$ of diluted molecular dye by applying Beer-Lambert law to the attenuation of a collimated light beam, transmitted through the cuvette. Aqueous suspensions of $0.525 \pm 0.01$ µm and $2.56 \pm 0.04$ µm polysterene particles (microParticles GmbH, Berlin, Germany) were used as scatterers. The scattering coefficient and SPF quantifiers $\gamma$ and $\delta$ of the liquid suspensions were computed using the Mie theory and the wavelength dependence of the polystyrene refractive index was taken from [27]. The MC simulations of turbid phantoms required numerical sampling of the scattering angles, which was implemented as detailed in [28]. The optical properties of the five turbid phantoms are summarized in Table 3.

Table 3. The range of optical properties ($\mu _a$, $\mu ^{\prime }_{s}$, $\gamma$ and $\delta$) spanned by the five liquid turbid phantoms (TP) that were prepared to experimentally validate the inverse models and the average diameter $d$ and standard deviation $\sigma$ of polystyrene particles in the individual turbid phantoms.

View Table | View all tables in this article

2.7 Experimental setup

Experimental validation was based on liquid turbid phantoms with known optical properties. The turbid phantoms were poured into a cylindrical beaker with a diameter of 19.5 mm and a height of 30 mm, coated with a black matte paint to reduce reflections from walls and stray light pollution. A custom-made linear array optical fiber probe (FiberTech Optica Inc., Ontario, Canada) with a 6.1 mm stainless steel housing and six 200 µm optical fibers (220, 440, 660, 880, and 1100 µm SDS) was dipped about 5 mm deep into the turbid phantoms. The turbid phantoms were illuminated by the first optical fiber that was coupled to a broadband halogen light source (AvaLight-Hal LS, Avantes, Apeldoorn, The Netherlands). The backscattered light collected through the remaining five optical fibers was dispersed onto a camera (acA1920-150um, Basler, Ahrensburg, Germany) by a fiber-coupled imaging spectrograph (ImSpector V10E, 400-1000 nm, Specim, Spectral Imaging Ltd., Oulu, Finland). The measured spectra were corrected for the dark response of the sensor array and normalized by the dark corrected reflectance of a standard diffuse tile (WS-2, Avantes, Apeldoorn, The Netherlands). A more detailed description of the measurement setup can be found in [29].

3. Results and discussion

The training process of RRMs was repeated 10 times for each dataset size summarized in Table 2 and visualized in Fig. 4. The best RRM was selected based on the mean absolute relative error between the regressed $\textrm {E}_{\textrm {RRM}, \textrm {S}}$ and simulated $\textrm {E}_{\textrm {RRM}, \textrm {R}}$ SRR datasets. The relative errors of the RRM as a function of the number of SRR samples in the training set is summarized in Fig. 5. The SRR regression errors begin to significantly increase if the training dataset is reduced to less than 225 SRR samples. The error maps in Fig. 6 provide visualization of the reflectance errors as a function of the optical properties. The relative errors shown in the maps were computed with respect to the reflectance values obtained by RRMs that were trained by the full size of the training dataset with 4335 SRR samples. The calculated errors are unbiased and start to significantly increase only for the two smallest training datasets with merely 135 and 36 SRR samples. To further validate the performance of the RRMs that map SRR from the simplified to the realistic PTM, we trained seven sets of inverse models for estimation of optical properties, namely $\textrm {IM}_{\textrm {S}}$, $\textrm {IM}_{\textrm {R,4335}}$, $\textrm {IM}_{\textrm {R,1215}}$, $\textrm {IM}_{\textrm {R,375}}$, $\textrm {IM}_{\textrm {R,225}}$, $\textrm {IM}_{\textrm {R,135}}$, $\textrm {IM}_{\textrm {R,36}}$. The errors of optical properties $\mu _a$, $\mu ^{\prime }_{s}$ and $\gamma$ estimated by the inverse models were independently evaluated using the $\textrm {E}_{\textrm {IM}, \textrm {S}}$ dataset for $\textrm {IM}_{\textrm {S}}$ and $\textrm {E}_{\textrm {IM}, \textrm {R}}$ dataset for $\textrm {IM}_{\textrm {R,4335}}$ through $\textrm {IM}_{\textrm {R,36}}$ (Fig. 7). The evaluation results confirm, that a small training set of merely 225 SRR samples is enough to prepare sufficiently accurate RRMs and from there inverse models ($\textrm {IM}_{\textrm {R}}$) for estimation of optical properties.

 

Fig. 5. Distribution of the relative reflectance errors for the realistic PTM at 5 SDS obtained by the RRMs as a function of the training dataset size. The total number of SRR samples in the training datasets varies from 4335 (left) to 36 (right).

Download Full Size | PDF

 

Fig. 6. Absolute relative error maps of reflectance that was regressed from SRR valid for the simplified PTM to SRR valid for the realistic PTM as a function of the number of SRR samples that were used to train the ANN-based RRMs. The relative errors are provided for three SDS, namely (first row) 220 µm, (second row) 660 µm and (third row) 1100 µm.

Download Full Size | PDF

 

Fig. 7. Root mean square errors (RMSE) and relative root mean square errors (RRMSE) of the absorption coefficient $\mu _a$, reduced scattering coefficient $\mu ^{\prime }_{s}$ and SPF quantifier $\gamma$, estimated by inverse models for the realistic PTM as a function of the size (4335, 1215, 375, 225, 135 and 36 SRR samples) of the dataset that was used to train the underlying RRM. Green line shows the performance of the inverse model for the simplified PTM that does not depend on the properties of the RRM and thus serves as a baseline.

Download Full Size | PDF

To highlight the computational efficiency of the proposed regression framework, we measured the MC simulation times required to prepare SRRs of comparable quality at 16 points representing the corners of the $\textrm {T}_{\textrm {IM}, \textrm {S}}$ dataset using the simplified and realistic PTMs. The MC simulations for the simplified PTM were terminated when the total weight of the photon packets collected through the annual ring representing the optical fiber with the largest SDS exceeded 1000. Likewise, the MC simulations for the realistic PTM were terminated when the total weight of the photon packets collected through the optical fiber with the largest SDS exceeded 40000. The two termination criteria produced reflectance of comparable quality on the most distant optical fiber. The MC simulations on Nvidia^TM Geforce GTX 1070 GPU took 242 s for the simplified PTM and 9533 s for the realistic PTM. The ratio between the two computational times $\approx 40$ tightly follows the ratio between the surface areas of the optical fiber and of the corresponding annual ring. The time required by the RRM to map the SRR from the simplified to the realistic PTM was found negligible. Since an accurate RRM that maps SRR from the simplified PTM to the realistic PTM can be derived from a small training dataset of merely 225 SRR samples, the total time required to prepare $\textrm {IM}_{\textrm {R,225}}$ is further reduced by a factor of $108500 {/} 225\approx 480$ (given that a training dataset for a simplified PTM is readily available).

From the timing data (Nvidia GTX 1070) we can estimate that the preparation of the datasets required to deploy the RRM in our study takes 19 days for $\textrm {T}_{\textrm {IM},\textrm {S}}$, 1.5 days for $\textrm {T}_{\textrm {RRM},\textrm {R}}$ and 1 hour for $\textrm {T}_{\textrm {RRM},\textrm {S}}$ (assuming that the datasets used to train the RRM contain 225 SRR samples). In contrast, without using the proposed RRM, the simulations of $\textrm {T}_{\textrm {IM},\textrm {R}}$ dataset, required to prepare the inverse models for realistic PTM, would take about 748 days.

To demonstrate the additional noise reduction attained with the proposed framework, we compared the spread of 100 reflectances simulated with the realistic PTM to the spread of 100 reflectances simulated with the simplified PTM and regressed to reflectances valid for the realistic PTM. The simulations were performed with $10^{7}$ photon packets for a selection of optical properties that exhibits the highest simulation noise: $\mu _a$ = 12 cm⁻¹, $\mu ^{\prime }_{s}$ = 5 cm⁻¹, $\gamma = 2.12$ and $\delta = 3.0$. The results presented in Fig. 8 show that the spread of reflectances valid for the realistic PTM is much smaller for the regressed reflectances that benefit from the radially symmetric detection scheme used in the MC simulations for simplified PTM. Note that the average values of the simulated and regressed reflectances match almost perfectly.

 

Fig. 8. Spread of 100 reflectance points simulated with the realistic PTM (red) and spread of 100 reflectance points simulated with the simplified PTM and regressed to the realistic PTM (blue). The reflectance points were simulated for optical properties ($\mu _a$ = 12 cm⁻¹, $\mu ^{\prime }_{s}$ = 5 cm⁻¹, $\gamma = 2.12$ and $\delta = 3.0$) that were found to produce the noisiest reflectance.

Download Full Size | PDF

Finally, we validated the proposed framework experimentally. For this purpose, two sets of measurements were collected from nonabsorbing and absorbing liquid turbid phantoms with well-defined optical properties summarized in Table 3. The nonabsorbing liquid phantoms were prepared using 2.56 µm polystyrene spheres and the absorbing phantoms were prepared using 0.525 µm polystyrene spheres and a green molecular absorber. The calibration factors that map the MC simulated SRR to the measured SRR were computed by dividing the MC simulated SRRs of nonabsorbing liquid phantoms and the corresponding measured SRRs. Ideally, the values of calibration factors at a given SDS should be the same for all the turbid phantoms. The values of calculated calibration factors for the 5 SDSs are shown in Fig. 9. The top row shows the calibration factors calculated from reflectance of optical phantoms simulated with the simplified (left) and realistic PTM (right) at each SDS. The bottom row shows the corresponding coefficients of variation (CV) for the calibration factors. One can clearly see that the calibration factors computed from the reflectance simulated with the simplified PTM do not match within a particular SDS and are biased in comparison to the calibration factors computed from the reflectance simulated with the realistic PTM. From there it follows that the CV of the calibration factors computed from the reflectance simulated with the simplified PTM is nearly 10%, while the CV of the calibration factors computed from the reflectance simulated with the realistic PTM is merely 1%. Further analysis and discussion on the sources of the observed mismatch and bias of the calibration factors computed from the reflectance simulated with the simplified PTM can be found in [11]. Consequently, the calibration factors at each SDS were smoothed by fitting a second-order polynomial to the calibration factors as a function of wavelength. The calibrated measured reflectances of the absorbing liquid phantoms $\textrm {Mie}_{\textrm {a},1}$ and $\textrm {Mie}_{\textrm {a},2}$ were passed through the inverse models $\textrm {IM}_{\textrm {R,225}}$. The SRR of the absorbing turbid phantoms for the realistic PTM were also simulated using the expected values of optical properties and the Mie or GK SPFs. The parameters of GK SPF $a_{GK}$ and $g_{GK}$ were calculated by matching the values of $\gamma$ and $\delta$ of the Mie SPF that holds for the liquid phantoms. Note, that some pairs of $\gamma$ and $\delta$ values cannot be attained by the GK SPF and were thus omitted from the SRR simulations. The optical properties estimated from the simulated SRR served as a baseline for the optical properties estimated from the measured SRR. The estimated values of optical properties are shown in Fig. 10 and the estimation errors are summarized in Table 4.The obtained RMSE and RRMSE errors of the estimated optical properties are similar for the simulated and the measured SRR, which indicates that the inverse models are unbiased. A slightly higher error of the estimated optical properties can be observed for the $\textrm {Mie}_{\textrm {a},2}$ turbid phantom (measured and simulated SRR), in particular at longer wavelengths, where the values of $\mu ^{\prime }_{s}$ drop bellow 9 cm⁻¹. The slightly increased errors could be explained by the influence of the SPF on the reflectance that is stronger at low $\mu ^{\prime }_{s}$ values, where the GK SPF that was used to simulate the training datasets of the inverse models cannot sufficiently account for the details of the Mie SPF that holds for the turbid phantoms.

 

Fig. 9. (top) Reflectance calibration factors for five SDS calculated from reflectance simulated with the simplified and realistic PTM for three non-absorbing liquid turbid phantoms $\textrm {Mie}_{\textrm {n},1}$ , $\textrm {Mie}_{\textrm {n},2}$ and $\textrm {Mie}_{\textrm {n},3}$. (bottom) Coefficient of variation (CV) for the calibration factors at each SDS for the simplified and realistic PTM.

Download Full Size | PDF

 

Fig. 10. Optical properties of the $\textrm {Mie}_{\textrm {a},1}$ and $\textrm {Mie}_{\textrm {a},2}$ turbid phantoms, predicted by $\textrm {IM}_{\textrm {R,225}}$. Red lines - optical properties estimated from SRR simulated with the Mie SPF using the realistic PTM, blue lines - optical properties estimated from the measured SRR, black lines - true/expected values of the optical properties.

Download Full Size | PDF

Table 4. Root mean square (RMSE) and relative root mean square (RRMSE) errors of optical properties estimated for the $\textrm {Mie}_{\textrm {a},1}$ (top rows) and $\textrm {Mie}_{\textrm {a},2}$ (bottom rows) absorbing turbid phantom from the SRR simulated with the GK and Mie SPFs using the realistic PTM, and from measured SRR.

View Table | View all tables in this article

Accurate MC simulation of reflectance for optical fiber probes are challenging and require careful consideration of all the details in the simulation model, in particular if used for estimation of optical properties by inverse models [11]. At small SDS that are frequently used with optical fiber probes, the inverse models need to take into account additional scattering phase function quantifiers such as $\gamma$ and $\delta$, which are related to second and third Legendre moments of the SPF, respectively. The additional parameters significantly increase the number of MC simulations required to prepare inverse models for estimation of optical properties. Furthermore, if a realistic description of the PTM is also used in the MC simulations, the symmetry of the detection scheme breaks down and the number of photon packets required for good quality SRRs further increases by about two orders of magnitude. Consequently, the total simulation time required to prepare datasets becomes excessively long and impractical even when using the latest generation GPUs. To overcome these limitations, many existing studies have used a simplified PTM [9,10,30–32], eventhough it has been shown [1,11,12] that such simplifications lead to large simulation errors of SRR. The proposed framework overcomes these limitations, and allows simulation of SRR using a realistic PTM at the speed and quality that matches MC simulations with a simplified PTM. The additional computational steps are based on efficient and accurate deep ANN-based RRMs, with negligible computational overhead. The proposed framework can be particularly beneficial for applications that utilize a linear layout of optical fibers in the probe, since the effect of the PTM on the reflectance [11] has been shown to increase with the SDS.

We believe that the presented framework can be applied to different measurement settings and in general take the efficiency and quality of MC simulations to a new level.

4. Conclusions

In this paper we present an efficient framework for MC simulation of SRR, that accounts for all the details of the optical fiber probe tip. For this purpose, we used deep neural networks that map the SRR, simulated with the simplified PTM to SRR, simulated with the realistic PTM. The results of the RRM were extensively validated through the use of inverse models for estimation of optical properties from simulated SRR and measured SRR of turbid phantoms. The results show that 225 SRR samples are sufficient to derive an accurate deep ANN-based RRM, which is about 480-fold less than the number of SRRs (108500) required to train a typical inverse model for estimation of optical properties. The remaining SRR samples can be efficiently simulated with the simplified PTM and rapidly regressed to SRR valid for the realistic PTM. Using the proposed methodology, it is possible to prepare high-quality inverse models for realistic PTM in approximately the same time as for simplified PTM.