Open Access
Add to BibSonomy
Abstract
Light propagation in biological tissues is frequently modeled by the Monte Carlo (MC) method, which requires processing of many photon packets to obtain adequate quality of the observed backscattered signal. The computation times further increase for detection schemes with small acceptance angles and hence small fraction of the collected backscattered photon packets. In this paper, we investigate the use of a virtually increased acceptance angle for efficient MC simulation of spatially resolved reflectance and estimation of optical properties by an inverse model. We devise a robust criterion for approximation of the maximum virtual acceptance angle and evaluate the proposed methodology for a wide range of tissue-like optical properties and various source configurations.
© 2017 Optical Society of America
1. Introduction
Light backscattered from a turbid sample carries an abundance of information on its structure and chemical composition. As such, the measured backscattered light can be exploited for noninvasive analysis and characterization of various tissue abnormalities [1–6]. For this purpose, the backscattered light is frequently captured at multiple source-detector separations (SDSs), which yields the so-called spatially resolved reflectance (SRR), R(r). Such measurements can be conveniently performed by optical fiber probes [7–13] or imaging systems [14–17]. The additional spatial information can be exploited to estimate the optical properties independently of the wavelength [18]. Ideally, the number of optical properties that can be estimated at each wavelength is limited by the number of utilized SDS. In general, imaging systems offer a higher spatial resolution and more SDSs than optical fiber probes, where the measurements are typically limited to about ten SDSs. Furthermore, the contactless nature of the measurements conducted by imaging systems overcomes the influence of the contact pressure induced by optical fiber probes [19], which can alter the optical properties of the probed sample.
The nature of the acquired backscattered light strongly depends on the SDS. For SDS longer than the transport mean free path 1/μtr = 1/(μa + μ′s) (diffusive regime), the backscattered light undergoes many scattering events and hence the shape of the scattering phase function p(cos θ) is averaged out. Consequently, the reflectance can be adequately described by the absorption μa coefficient and the reduced scattering μ′s = μs (1 − g1) coefficient, which depends on the scattering coefficient μs and the first Legendre moment g1 of the scattering phase function also termed as the anisotropy factor. In contrast, for SDS shorter or comparable to the transport mean free path (termed also as the subdiffusive regime), the backscattered light undergoes only a few scattering events and is consequently more sensitive to the shape of the scattering phase function. To accommodate some additional information on the scattering phase function, a parameter γ = (1 − g2)/(1 − g1), comprising the first g1 and second Legendre moments g2 of the scattering phase function, was introduced [13,20]. Furthermore, it has been shown that the additional parameter γ also improves the estimation of optical properties from SRR by an inverse model [18]. Recently, a similar parameter σ that incorporates higher order Legendre moments of the scattering phase function was proposed by Bodenschatz et al. [21]. However, in terms of optical properties estimated from SRR, γ and σ yield similar results [21].
Light propagation in turbid sample such as biological tissue can be adequately described by the radiative transport equation (RTE) [22]. For a limited number of experimental setups, the RTE can be solved analytically [23,24]. Analytical solutions are also readily available for samples with dominant scattering at large SDS where the light field can be considered nearly isotropic. In this particular case, the RTE can be solved within the diffusion approximation [25–28]. However, for complex experimental settings, analytical solutions of the RTE are not available, and the Monte Carlo (MC) method is most frequently used to solve the RTE numerically. Due to the stochastic nature of the MC method, a large number of launched photon packets is required to achieve an adequate signal-to-noise ratio of the simulated quantity such as the backscattered light or reflectance. Moreover, the computation times significantly increase for detection schemes with small acceptance angles (i.e., low numerical apertures), which are typically used in imaging-based experimental settings.
In this paper, we present a framework that can be used to significantly reduce the computation time of the MC simulations of SRR by virtually increasing the acceptance angle of the detection scheme. The use of a virtually increased numerical aperture was already disscused by Thueler et al. [5], however the validation was limited and did not involve the use of inverse models. Additionally, Bargo et al. [29, 30] pointed out the importance of using the correct numerical aperture and showed that the collection efficiency of the optical fiber significantly depends on the optical properties of a turbid sample. Furthermore, the difference between the numerical aperture of the experimental setup and the numerical aperture used in the MC simulations is still frequently neglected in the existing literature [17]. Therefore, we first study the influence of the virtually increased acceptance angle on the relative errors of the MC simulated SRR as a function of the sample optical properties. Subsequently, we devise a criterion for estimation of the maximum virtual acceptance angle of the detection scheme, which retains small relative errors of the simulated SRR. The proposed methodology is explained in Section 2, where a detailed description of the light propagation model and of the inverse models used in this study are provided. In Section 3, we evaluate the methodology by a forward and inverse MC models for a wide range of optical properties that can be found in biological tissues, and compare the computation times as a function of the virtual acceptance angle.
2. Materials and methods
2.1. Light propagation model
Monte Carlo method [31,32] is considered as the gold standard in the biomedical community since it offers accurate description of the light propagation for an arbitrary tissue geometry and experimental setup within the context of RTE [33]. The methodology is particularly useful when a complex structure of the tissue and/or experimental setting must be taken into account [9,34]. In this study, an experimentally validated MC model for multi-layered tissue, implemented in the OpenCL™ standard for parallel programming [35], was used. The outline of the algorithm follows the publicly avaliable MC models MCML [36] and CUDA-MCML [37] where discrete absorption weighting is applied at each scattering event.
All the MC simulations were conducted under the semi-infinite geometry (nm = 1.33). The MC simulations are frequently run with a fixed number of launched photon packets regardless of the number or total weight of the backscattered photon packets collected through the detection scheme. This can lead to orders of magnitude variations in the total weight of the collected photon packets across turbid samples with different optical properties. With the aim to attain comparable backscattered signal levels for simulated SRR across the full range of biologically relevant optical properties, the MC simulations were terminated when the total remaining weight of the collected photon packets (Wtot) exceeded a predefined value:
where wi is the remaining weight of the i-th collected backscattered photon packet. For absorbing media wi is less then 1 and the number of collected photon packets Ncollected exceeds Wtot. In this study, the value of Wtot was set to 106 and was realized through consecutive runs of MC simulations by a fixed number of launched photon packets (107).2.2. Inverse model
The spatially resolved reflectance R(r) for a given experimental setting and sample with known optical properties can be predicted by a light propagation model f [20]:
where μ′s, μa, p(cos θ) are the tissue optical properties and G holds all the remaining details of the experimental setting. If the inverse model f−1 can be expressed analytically, then the tissue optical properties μa, μ′s, p(cos θ) can be easily estimated from R(r) [20]: However, in practice, analytical solutions of the light propagation model and the corresponding inverse are rarely available. Consequently, empirical models are frequently used to solve the inverse problem. For this purpose, a computationally very efficient look-up table-based (LUT) model [38,39] is frequently employed. Briefly, the LUT entries are reflectance values or SRR presented as a function of the relevant optical properties. The optical properties of the probed sample are derived by searching for the best match between the measured SRR and the LUT entries. The LUT can be populated by experimental measurements [40–42] of turbid phantoms with well defined optical properties or by MC simulations [41, 43]. In this study, we used an extended form of a simulated LUT that was introduced by Hennessy et al. [38]. The two-dimensional LUT defined over μ′s and μa was extended by an additional dimension γ, which accounts for the phase function information. The LUT was populated by MC simulations using the Gegenbauer kernel (GK) phase function [44] with two parameters ggk and αgk. The GK phase function was chosen for the wide range of γ values that can be obtained. A recent study [18] has also shown that in terms of the parameter γ, other phase functions such as the HG, MHG and Mie can be well approximated by the GK phase function.The inverse problem (Eq. 3) can be solved by capturing the backscattered light at different wavelengths (spectrally resolved reflectance) or at different SDS (spatially resolved reflectance). In the former case, additional wavelength-dependent regularization of the absorption and reduced scattering coefficients are required. The μa is usually modeled by a sum of the absorption spectra weighted by the respective volume fractions of the corresponding absorbing constituents comprising the sample. On the other hand, μ′s is frequently described by a parametric form such as μ′s (λ) = a(λ/λ0)−b, where a is the reduced scattering coefficient at a normalization wavelength λ0 and b the scattering power-law exponent related to the scatterer size. Unfortunately, the sample constituents are not always known a priori. Therefore, the reflectance captured at three or more different SDS ideally provides sufficient information to solve the inverse problem independently of the wavelength [18]. The fine spatial resolution of the SRR captured by imaging systems should be well suited for this task.
The inverse model used in this study comprises two steps. First, the measured R(r) is compared to all the LUT entries by using the selected criterion function. The obtained global minimum is then refined by local optimization utilizing linear interpolation of the LUT entries. The criterion function CF was defined as the sum of squared differences of the logarithmic SRR:
where i runs from 1 to Nr through all the SDSs of the SRR, Rm(r) is the measured SRR of the sample and RLUT (r) is the LUT entry or value interpolated from the LUT. The logarithmic scale of the criterion function was selected based on the findings presented in a recent study [18]. The LUT used in this study was populated for a wide range of optical properties that can be found in biological tissues [13,20,45,46]. The absorption coefficient μa was varied from 0 cm−1 to 20 cm−1 in 21 equally spaced steps, the reduced scattering coefficient μ′s from 5 cm−1 to 70 cm−1 in 25 equally spaced steps, and the phase function parameter γ from 1.60 to 2.30 in 15 equally spaced steps. The anisotropy factor g1 was held constant at 0.8.2.3. Virtually increased acceptance angle
Since the acceptance angle of a real detection scheme is limited, only a fraction of the total number of the backscattered photon packets from the observed turbid sample can be collected. In general, a detection scheme collects photon packets that are backscattered within the nominal acceptance angle θ0 of the front lens (Fig. 1(a)).
Fig. 1 Example of a typical imaging system with a narrow nominal acceptance angle θ0 (a). The nominal acceptance angle θ0 can be virtually increased in the MC simulations to capture a larger fraction of the backscattered photon packets (b). If the quotient between the R(θ0, r), captured at the nominal acceptance angle, and the R(θv, r), captured at the virtually increased acceptance angle θv, is approximately constant and independent of r, the virtual detection scheme can be used to reduce the computation time without introducing additional errors (c).
In order to accurately simulate the SRR captured by an imaging system, the full details of the measurement setting have to be carefully considered in the MC simulations. Furthermore, while the MC simulated SRR is computed on an absolute scale (total weight of the detected photon packets with respect to the total weight of the launched photon packets), the measured SRR are usually normalized to a calibration target such as Spectralon ®. Consequently, a calibration factor that connects the measured and the MC simulated SRR needs to be determined. For this purpose, measurements of optical phantoms with well-defined optical properties are usually used [47,48]. However, the multiplicative nature of the calibration factor can be also exploited to simplify or modify the experimental geometry in the MC simulations. The basic idea behind this study is to virtually increase the acceptance angle (Fig. 1(b)) of the detection scheme up to a value that is still compliant with the multiplicative nature of the calibration factor. In this way, the fraction of the detected backscattered photon packets can be increased, and consequently, the required computation time significantly reduced. Given that the multiplicative factor connecting the MC simulated R(θ0, r) captured at the nominal acceptance angle θ0 and the MC simulated R(θv, r) captured at the virtually increased acceptance angle θv is independent of r (Fig. 1(c)), the additional multiplicative factor k can be included in the calibration factor without any influence on the accuracy of the simulated SRR. Hence, the MC simulated R(θv, r) at a virtually increased acceptance angle can be approximated by the product of the MC simulated R(θ0, r) at the nominal acceptance angle and a multiplicative factor k:
In the case of real measurements, k is implicitly included in the calibration factor. On the other hand, in MC studies, the value of k can be estimated as the spatial average of the quotient between R(θv, r) and R(θ0, r): where Nr is the number of all SDSs.However, increasing the acceptance angle will at some point make the multiplicative factor k dependent on r and hence break the presumption. Therefore, two metrics for estimation of the maximum virtual acceptance angle θmax were introduced and compared. Firstly, the absolute relative SRR error (ARE) arising from the virtually increased acceptance angle was defined as:
Since the SRR signal at longer SDS decreases by a few orders of magnitude the ARE metric can become sensitive to the simulation noise. Consequently, a second error metric was formulated in a way that is less sensitive to the noise at longer SDS and was defined as a relative root mean square error rRMSE:The performance of the ARE and rRMSE metrics were evaluated for a turbid sample with μa = 2.0 cm−1, μ′s = 45.6 cm−1 and γ = 1.65. The computed AREs are presented in Fig. 2(a). It can be observed that ARE is not fully independent of r. At small r (in the subdiffusive regime), the ARE increases with the acceptance angle θv. Consequently, the selection of θv at small r can have a significant impact on the simulated SRR that goes beyond the spatially independent multiplicative factor k. In contrast, at larger r (in the diffuse regime), the ARE of the simulated SRR shows no significant dependence on the θv, but becomes increasingly affected by the simulation noise. Since the influence of simulation noise on the ARE at larger r can easily exceed the ARE values observed in the subdiffusive regime, deriving θmax from ARE would lead to inconsistent results. Consequently, the second metric, rRMSE, was evaluated. The integral nature of the metric leads to smooth monotonically increasing values of rRMSE as a function of θv (Fig. 2(b)). The value of θmax can be easily determined by a predefined threshold Et (Fig. 2(b)). For the purpose of this study, the threshold Et was set to 1%.
Fig. 2 Reflectance error arising from the virtually increased acceptance angle θv quantified by the spatial dependence of the absolute relative SRR error ARE (a) and the relative root mean square SRR error rRMSE (b). Results are presented for a nominal acceptance angle of 3° and for a turbid sample with absorption coefficient of 2.0 cm−1, reduced scattering coefficient of 45.6 cm−1, and phase function parameter γ of 1.65. A predefined 1% threshold value for Et is marked by a horizontal black dashed line and the estimated maximum virtual acceptance angle θmax with a vertical red dashed line.
2.4. Experimental setup and synthetic data sets
The simulated R(r) comprised 500 equally spaced (5 μm) reflectance points that followed from the geometry of the imaging system introduced in [16] where the fine spatial resolution was achieved with a 60 mm macro lens. The high spatial resolution comes at a cost of a narrow nominal acceptance angle θ0 of 3°. However, to cover different experimental settings, three values of θ0 (3°, 5°, 10°) were evaluated in this study. In the existing literature, normally incident beams of uniform spatial intensity and diameter of 200 μm [5,9,49] are most frequently used as a light source, however, some applications based on smaller and larger beam diameters can also be found [14,50]. The spatial intensity of laser sources is usually approximated by a Gaussian distribution [21]. Following the variety of sources used in the literature, the simulations were prepared for three uniform beams with diameters of 50 μm, 200 μm and 500 μm and similarly for three Gaussian beams with 1/e2 beam diameters of 50 μm, 200 μm and 500 μm. The R(r) of different beams were derived from MC simulations by convolution [51]. To exclude the incident beam, they were cropped at the outer edge of the uniform beam or at 1/e2 of the Gaussian beam. Finally, the SRR profiles were also cut off at the distance of 1500 μm from the edge of the beam for all types of sources.
The influence of the virtually increased acceptance angle θv was evaluated by observing the relative error of the MC simulated SRR and optical properties estimated by the inverse model. To avoid potential systematic errors arising from the properties of the phase function and to isolate the influence of θv on SRR, the synthetic data sets and the LUT entries were calculated by using the same scattering phase function. In this way, the errors of the optical properties estimated by the inverse model could be entirely attributed to the change in the acceptance angle. However, to objectively evaluate the inverse models, the optical properties of the synthetic data sets were not located on the LUT grid. For this purpose, synthetic data sets were computed on a grid of eight equally spaced values of μ′s, μa and γ that spanned [10 cm−1, 60 cm−1], [0.5 cm−1, 20 cm−1] and [1.65, 2.20], respectively.
3. Results and discussion
3.1. Maximum virtual acceptance angle
The values of θmax were first estimated for a wide range of optical properties by using the rRMSE metric and a 1% threshold Et. Dependence of the θmax on the turbid sample optical properties for a uniform beam of diameter 50 μm and nominal acceptance angle of θ0 = 3° is depicted in Fig. 3. The maximum virtual acceptance angle θmax strongly depends on the properties of the turbid sample.
Fig. 3 Maximum virtual acceptance angle θmax as a function of the absorption μa and reduced scattering μ′s coefficients at γ = 1.75 (a) and γ = 2.15 (b). The minimum values of θmax observed across the full range of μa and μ′s as a function of γ (c) and the corresponding values of μa (d) and μ′s (e).
Figure 3(a) shows the estimated θmax values as a function of μa and μ′s at γ = 1.75. The highest θmax values are observed for μa near 10 cm−1 and μ′s near 40 cm−1 where the acceptance angle can be virtually increased from 3° to 35°. However, even for low μ′s below 10 cm−1, where the smallest values of θmax are observed, the acceptance angle can be virtually increased up to 10°. Similarly, Fig. 3(b) shows the map of θmax at γ = 2.15. In this case, the acceptance angle can be increased to 16° and the lowest θmax is observed at high values of μ′s above 65 cm−1.
The dependence of θmax on the phase function parameter γ is presented in Fig. 3(c) where the minimum values observed across the full range of μa and μ′s are shown. For γ ⩽ 1.9, the minimum values of θmax are observed at μ′s below 10 cm−1 and low μa, while for γ > 1.9, the minimum values of θmax are located at μ′s above 65 cm−1 and at high μa (see Figs. 3(d) and 3(e)). The observed shift in location of θmax also explains the sharp change in θmax at γ = 1.9 that can be seen in Fig. 3(c).
For a 50 μm uniform beam and θ0 = 3°, the value of θmax was approximated to 9.5°, which corresponds to the minimum value of θmax observed in Fig. 3(c). Following the same procedure, values of θmax were estimated for all the remaining sources and are summarized in Table 1. For θ0 = 3°, θmax of 10° can be used for all the studied sources. The value of θmax further increases with the value of the nominal acceptance angle.
Table 1. Maximum virtual acceptance angles θmax (°) for different sources and nominal acceptance angles θ0.
3.2. Verification by inverse model
In the previous section, the maximum virtual acceptance angle θmax was estimated by allowing a 1% absolute relative SRR error rRMSE. Even though the selected error margin might appear low, it does not necessarily guarantee adequate performance of the proposed methodology when used to estimate the optical properties by an inverse model. Therefore, the influence of the virtually increased acceptance angle on the inverse model was evaluated on synthetic data sets, described in Section 2.4. The SRR of synthetic data sets was simulated by using the nominal acceptance angle of θ0 = 3°. In contrast, ten different LUT-based inverse models were prepared by varying the virtual acceptance angle. First, a reference LUT was prepared by setting θv to the nominal acceptance angle of 3°. Subsequently, nine more LUTs were computed for virtual acceptance angles of 10°, 20°, 30°, 40°, 50°, 60°, 70°, 80° and 90°. The multiplicative factor at each virtual acceptance angle k(θv) was estimated as the spatial average over all the optical properties in the LUT:
Influence of the θv on the inverse model was quantified by the root-mean-square error (RMSE) and the relative root-mean-square error (rRMSE): where i runs over all the synthetic data sets, y denotes one of μa, μ′s or γ, yi represents the true value and ŷi the corresponding estimated value.First, we analyzed the performance of the inverse model based on the reference LUT that was prepared with the nominal acceptance angle (Fig. 4 top row). These results present a baseline for the RMSE and rRMSE, since the obtained errors can be attributed solely to the properties of the inverse model and are not affected by the change in the acceptance angle. The obtained rRMSE values for μa, μ′s and γ are under 1% and are comparable to the results from our recent study [18].
Fig. 4 Scatter plots of the true and estimated values of the absorption coefficient μa (first column) reduced scattering coefficient μ′s (second column) and γ (third column) for virtual acceptance angles θv of 3° (first row), 10° (second row) and 40° (third row) obtained for a 50 μm uniform beam and nominal acceptance angle θ0 = 3°. Ideal relation between the true and estimated values is denoted by a dashed gray line.
Results presented in the second row of Fig. 4 were obtained by an inverse model based on a LUT computed at θv = 10°. The estimated rRMSE for μa, μ′s and γ are slightly higher, but still significantly under the 1% margin that is frequently used to identify successful measurements [18]. These results are consistent with the initial hypothesis, that small nominal acceptance angles can be significantly increased in the MC simulations without adding error to the SRR. The results also confirm that the proposed error metric from Section 2.3 can be used to estimate the maximum virtual acceptance angle.
Finally, in the last row of Fig. 4, the value of θv was increased to 40°, which exceeds by four times the estimated maximum virtual acceptance angle. As expected, the proposed methodology breaks down and significant errors are introduced into the optical properties estimated by the inverse model. The rRMSE of the inverse model now grows above 1% margin for all the three optical properties, in particular for γ. These observations could be explained by the fact that larger values of θv lead to higher errors in the SRR at shorter SDS (subdiffusive regime), where γ is much more sensitive to the SRR information than the other two optical properties. Especially lower γ values are more scattered which is consistent with the results from Fig. 3(c) where lower θv was predicted for lower γ values.
The values of RMSE and rRMSE for all the ten LUT-based inverse models are summarized in Table 2. As expected, the RMSE and rRMSE are increasing monotonically with the value of θv.
Table 2. Root-mean-square error (RMSE) and the relative root-mean-square error (rRMSE) of the estimated absorption μa and reduced scattering μ′s coefficients and γ as a function of the virtual acceptance angle θv for a 50 μm uniform beam and nominal acceptance angle θ0 = 3°.
All results presented in Fig. 4 and Table 2 were obtained for a 50 μm uniform beam and a nominal acceptance angle of θ0 = 3°. The same experiments were conducted for all the remaining source configurations defined in Section 2.4. The obtained rRMSE are presented in Fig. 5 and the obtained RMSE in Fig. 6. In general, the results are not significantly affected by the shape and size of the beam.
Fig. 5 Relative root-mean-square error (rRMSE) of the estimated μa (a), μ′s (b) and γ (c) as a function of the virtual acceptance angle θv obtained for different sources. Dashed line marks 1% rRMSE.
Fig. 6 Root-mean-square error (RMSE) of the estimated μa (a), μ′s (b) and γ (c) as a function of the virtual acceptance angle θv obtained for different sources.
3.3. Time complexity of Monte Carlo simulations
Reducing the computation time required to populate the LUT by MC simulations, is the most prominent aspect of the proposed methodology. Since a larger acceptance angle of the detection scheme leads to a higher fraction of collected backscattered photon packets, fewer photon packets need to be processed to attain the same total weight of the backscattered photon packets Wtot. In this study, Wtot was fixed to 106 and all simulations were performed on the same GPU device (Nvidia Tesla K40).
The normalized distributions of the number of launched photon packets as a function of the virtual acceptance angle (3°, 10°, 40°) are shown in Fig. 7. The distributions are mostly governed by the range of optical properties and somewhat also by the stochastic nature of the MC method. In general, a larger virtual acceptance angle requires fewer launched photon packets in order to attain a predefined total remaining weight of the detected backscattered photon packets. For example, increasing the acceptance angle from the nominal value of 3° to its estimated maximum value of 10°, leads to more then tenfold reduction in the number of launched photon packets. Since the computation time of GPU-based MC simulation is approximately linearly related to the number of launched photons packets, the computation time should scale down by a similar factor. Table 3 summarizes the computation times required to prepare the 3D LUTs for different θv. Increasing the nominal acceptance angle θ0 = 3° to the estimated maximum value θmax = 10°, reduced the initial computation time from 18.8 days (452.2 h) to only about one day and a half (35.9 h), which is consistent with the observed drop in the number of launched photon packets. Increasing the acceptance angle beyond 10° results in further reduction of the computation time (Table 3), however, at the cost of increased error in the optical properties estimated by the inverse model. Figure 8 summarizes the required computation times and the maximum rRMSE of the estimated optical properties as a function of the virtual acceptance angle. If a higher than 1% rRMSE of the estimated optical properties is acceptable, the virtual acceptance angle can be further increased. Since the presented values are valid for a wide range of optical properties, additional gains are possible for reduced subsets of optical properties, which might be more realistic when conducting measurements on a certain biological tissue.
Fig. 7 Normalized distribution of the number of launched photon packets over the full range of optical properties as a function of the virtual acceptance angle θv and for a fixed total weight of the backscattered photons packets Wtot = 106. The median value of each distribution is denoted by a vertical dashed line.
Fig. 8 Required computation time and relative root-mean-square error (rRMSE) of the estimated optical properties as a function of the virtual acceptance angle θv.
Table 3. The computation times required to populate a LUT as a function of the virtual acceptance angle θv.
4. Conclusion
This study offers a simple yet effective approach to reduce the computation time of Monte Carlo (MC) simulations of spatially resolved reflectance. The proposed methodology assumes that the spatially resolved reflectance (SRR) can be, up to a certain acceptance angle of the detection scheme, sufficiently modeled by a multiplicative factor that is independent of the source detector separation (SDS). In this way, the nominal acceptance angle of the detection scheme can be virtually increased in the MC simulations, which results in a higher fraction of detected backscattered photon packets and consequently shorter computation time.
In this study, the errors introduced by using the virtually increased acceptance angle in the MC simulations were investigated extensively by analyzing the simulated SRR and optical properties estimated by an inverse model. We found that errors are much higher in the reflectance captured at shorter SDS, while the acceptance angle does not significantly affect the reflectance captured at longer SDS. We devised a criterion that allows robust estimation of the maximum virtual acceptance angle considering the acceptable level of error in the simulated SRR. The results showed that the proposed criterion allows estimation of the reduced scattering coefficient μ′s, absorption coefficient μa and phase function related parameter γ by an inverse model without introducing additional significant errors. For the given experimental setup and for a wide range of optical properties, the nominal acceptance angle of the detection scheme could be increased from 3° to 10°. Consequently, the relative root-mean-square errors of the optical properties estimated by the inverse model increased only slightly, by 0.2% for μa, by 0.0% for μ′s and by 0.3% for γ. On the other hand, the increased acceptance angle resulted in a thirteen fold shorter average computation time, reducing the time required to prepare a 3D lookup table of SRR from 18.8 days (452.2 h) to only about one and a half day (35.9 h).
Funding
The authors acknowledge the financial support from the Slovenian Research Agency for research core funding No. P2-0232 and projects No. J2-8173, J2-7118, J2-7211 and J2-5473.
Acknowledgments
The authors acknowledge the Slovenian Initiative for National Grid (SLING) for providing the access to the GPU-cluster.
Disclosures
The authors declare that there are no conflicts of interest related to this article.
References and links
1. N. Bodenschatz, S. Lam, A. Carraro, J. Korbelik, D. M. Miller, J. N. McAlpine, M. Lee, A. Kienle, and C. MacAulay, “Diffuse optical microscopy for quantification of depth-dependent epithelial backscattering in the cervix,” J. Biomed. Opt. 21, 066001 (2016). [CrossRef]
2. D. M. McClatchy, E. J. Rizzo, W. A. Wells, P. P. Cheney, J. C. Hwang, K. D. Paulsen, B. W. Pogue, and S. C. Kanick, “Wide-field quantitative imaging of tissue microstructure using sub-diffuse spatial frequency domain imaging,” Optica 3, 613 (2016). [CrossRef] [PubMed]
3. J. J. Bravo, K. D. Paulsen, D. W. Roberts, and S. C. Kanick, “Sub-diffuse optical biomarkers characterize localized microstructure and function of cortex and malignant tumor,” Opt. Lett. 41, 781 (2016). [CrossRef] [PubMed]
4. P. Usenik, M. Bürmen, A. Fidler, F. Pernuš, and B. Likar, “Automated Classification and Visualization of Healthy and Diseased Hard Dental Tissues by Near-Infrared Hyperspectral Imaging,” Appl. Spectrosc. 66, 1067–1074 (2012). [CrossRef]
5. P. Thueler, I. Charvet, F. Bevilacqua, M. St. Ghislain, G. Ory, P. Marquet, P. Meda, B. Vermeulen, and C. Depeursinge, “In vivo endoscopic tissue diagnostics based on spectroscopic absorption, scattering, and phase function properties,” J. Biomed. Opt. 8, 495–503 (2003). [CrossRef] [PubMed]
6. E. Salomatina, B. Jiang, J. Novak, and A. N. Yaroslavsky, “Optical properties of normal and cancerous human skin in the visible and near-infrared spectral range,” J. Biomed. Opt. 11, 064026 (2006). [CrossRef]
7. M. Sharma, R. Hennessy, M. K. Markey, and J. W. Tunnell, “Verification of a two-layer inverse Monte Carlo absorption model using multiple source-detector separation diffuse reflectance spectroscopy,” Biomed. Opt. Express 5, 40–53 (2013). [CrossRef]
8. T. Y. Tseng, C. Y. Chen, Y. S. Li, and K. B. Sung, “Quantification of the optical properties of two-layered turbid media by simultaneously analyzing the spectral and spatial information of steady-state diffuse reflectance spectroscopy,” Biomed. Opt. Express 2, 901 (2011). [CrossRef] [PubMed]
9. P. Naglič, F. Pernuš, B. Likar, and M. Bürmen, “Limitations of the commonly used simplified laterally uniform optical fiber probe-tissue interface in Monte Carlo simulations of diffuse reflectance,” Biomed. Opt. Express 6, 3973 (2015). [CrossRef]
10. S. C. Kanick, U. A. Gamm, M. Schouten, H. J. C. M. Sterenborg, D. J. Robinson, and A. Amelink, “Measurement of the reduced scattering coefficient of turbid media using single fiber reflectance spectroscopy: Fiber diameter and phase function dependence,” Biomed. Opt. Express 2, 1687–1702 (2011). [CrossRef] [PubMed]
11. M. Bregar, M. Bürmen, U. Aljančič, B. Cugmas, F. Pernuš, and B. Likar, “Contact pressure–aided spectroscopy,” J. Biomed. Opt. 19, 020501 (2014). [CrossRef]
12. U. Utzinger and R. R. Richards-Kortum, “Fiber optic probes for biomedical optical spectroscopy,” J. Biomed. Opt. 8, 121–147 (2003). [CrossRef] [PubMed]
13. F. Bevilacqua and C. Depeursinge, “Monte Carlo study of diffuse reflectance at source–detector separations close to one transport mean free path,” J. Opt. Soc. Am. A 16, 2935 (1999). [CrossRef]
14. F. Foschum, M. Jäger, and A. Kienle, “Fully automated spatially resolved reflectance spectrometer for the determination of the absorption and scattering in turbid media,” Rev. Sci. Instrum. 82, 103104 (2011). [CrossRef] [PubMed]
15. H. Cen and R. Lu, “Optimization of the hyperspectral imaging-based spatially-resolved system for measuring the optical properties of biological materials,” Opt. Express 18, 17412 (2010). [CrossRef] [PubMed]
16. M. Ivančič, P. Naglič, F. Pernuš, B. Likar, and M. Bürmen, “Extraction of optical properties from hyperspectral images by Monte Carlo light propagation model,” “Proc. SPIE 9706, 97061 (2016).
17. I. Nishidate, T. Ishizuka, A. Mustari, K. Yoshida, S. Kawauchi, S. Sato, and M. Sato, “Evaluation of Cerebral Hemodynamics and Tissue Morphology of In Vivo Rat Brain Using Spectral Diffuse Reflectance Imaging,” Appl. Spectrosc. 69, 03702816657569 (2016).
18. P. Naglič, F. Pernuš, B. Likar, and M. Bürmen, “Estimation of optical properties by spatially resolved reflectance spectroscopy in the subdiffusive regime,” J. Biomed. Opt. 21, 095003 (2016). [CrossRef]
19. B. Cugmas, M. Bregar, M. Bürmen, F. Pernuš, and B. Likar, “Impact of contact pressure–induced spectral changes on soft-tissue classification in diffuse reflectance spectroscopy: Problems and solutions,” J. Biomed. Opt. 19, 037002 (2014). [CrossRef]
20. K. W. Calabro and I. J. Bigio, “Influence of the phase function in generalized diffuse reflectance models: Review of current formalisms and novel observations,” J. Biomed. Opt. 19, 075005 (2014). [CrossRef]
21. N. Bodenschatz, P. Krauter, A. Liemert, and A. Kienle, “Quantifying phase function influence in subdiffusively backscattered light,” J. Biomed. Opt. 21, 035002 (2016). [CrossRef]
22. L. V. Wang and H. I. Wu, Biomedical Optics: Principles and Imaging, 1st ed. (Wiley-Interscience, Hoboken, New Jersey, USA, 2007).
23. A. Liemert and A. Kienle, “Exact and efficient solution of the radiative transport equation for the semi-infinite medium,” Sci. Rep. 3, 2018 (2013). [CrossRef] [PubMed]
24. E. Vitkin, V. Turzhitsky, L. Qiu, L. Guo, I. Itzkan, E. B. Hanlon, and L. T. Perelman, “Photon diffusion near the point-of-entry in anisotropically scattering turbid media,” Nat. Commun. 2, 587 (2011). [CrossRef] [PubMed]
25. T. J. Farrell, M. S. Patterson, and B. Wilson, “A diffusion theory model of spatially resolved, steady-state diffuse reflectance for the noninvasive determination of tissue optical properties in vivo,” Med. Phys. 19, 879–888 (1992). [CrossRef] [PubMed]
26. A. Kienle and M. S. Patterson, “Improved solutions of the steady-state and the time-resolved diffusion equations for reflectance from a semi-infinite turbid medium,” J. Opt. Soc. Am. A 14, 246 (1997). [CrossRef]
27. P. Naglič, L. Vidovič, M. Milanič, L. L. Randeberg, and B. Majaron, “Applicability of diffusion approximation in analysis of diffuse reflectance spectra from healthy human skin,” Proc. SPIE 9032, 90320N (2013). [CrossRef]
28. E. L. Hull and T. H. Foster, “Steady-state reflectance spectroscopy in the P3 approximation,” J. Opt. Soc. Am. A 18, 584 (2001). [CrossRef]
29. P. R. Bargo, S. A. Prahl, and S. L. Jacques, “Collection efficiency of a single optical fiber in turbid media,” Appl. Opt. 42, 3187–3197 (2003). [CrossRef] [PubMed]
30. P. R. Bargo, S. A. Prahl, and S. L. Jacques, “Optical properties effects upon the collection efficiency of optical fibers in different probe configurations,” IEEE J. Sel. Top. Quantum Electron. 9, 314–321 (2003). [CrossRef]
31. S. A. Prahl, M. Keijzer, S. L. Jacques, and A. J. Welch, “A Monte Carlo Model of Light Propagation in Tissue,” “SPIE Series Vol. 5, 102–111” (1989).
32. L. Wang and S. L. Jacques, “Monte Carlo modeling of light transport in multi-layered tissues in standard C,” The University of Texas, MD Anderson Cancer Center, Houston (1992).
33. C. Zhu and Q. Liu, “Review of Monte Carlo modeling of light transport in tissues,” J. Biomed. Opt. 18, 050902 (2013). [CrossRef]
34. B. Majaron, M. Milanič, and J. Premru, “Monte Carlo simulation of radiation transport in human skin with rigorous treatment of curved tissue boundaries,” J. Biomed. Opt. 20, 015002 (2015). [CrossRef] [PubMed]
35. J. E. Stone, D. Gohara, and G. Shi, “OpenCL: A Parallel Programming Standard for Heterogeneous Computing Systems,” IEEE Des. Test 12, 66–73 (2010).
36. L. Wang, S. L. Jacques, and L. Zheng, “MCML—Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Meth. Prog. Bio. 47, 131–146 (1995). [CrossRef]
37. E. Alerstam, T. Svensson, and S. Andersson-Engels, “Parallel computing with graphics processing units for high-speed Monte Carlo simulation of photon migration,” J. Biomed. Opt. 13, 060504 (2008). [CrossRef]
38. R. Hennessy, S. L. Lim, M. K. Markey, and J. W. Tunnell, “Monte Carlo lookup table-based inverse model for extracting optical properties from tissue-simulating phantoms using diffuse reflectance spectroscopy,” J. Biomed. Opt. 18, 037003 (2013). [CrossRef] [PubMed]
39. G. M. Palmer and N. Ramanujam, “Monte Carlo-based inverse model for calculating tissue optical properties Part I: Theory and validation on synthetic phantoms,” Appl. Opt. 45, 1062 (2006). [CrossRef] [PubMed]
40. B. S. Nichols, N. Rajaram, and J. W. Tunnell, “Performance of a lookup table-based approach for measuring tissue optical properties with diffuse optical spectroscopy,” J. Biomed. Opt. 17, 057001 (2012). [CrossRef] [PubMed]
41. N. Rajaram, T. H. Nguyen, and J. W. Tunnell, “Lookup table–based inverse model for determining optical properties of turbid media,” J. Biomed. Opt. 13, 050501 (2008). [CrossRef]
42. P. Naglič, M. Bregar, F. Pernuš, B. Likar, and M. Bürmen, “Accuracy of experimental data and Monte Carlo simulation lookup table-based inverse models for assessment of turbid media optical properties with diffuse reflectance spectroscopy,” Proc. SPIE 9333, 933310 (2015). [CrossRef]
43. I. Fredriksson, M. Larsson, and T. Strömberg, “Inverse Monte Carlo method in a multilayered tissue model for diffuse reflectance spectroscopy,” J. Biomed. Opt. 17, 047004 (2012). [CrossRef] [PubMed]
44. L. O. Reynolds and N. J. McCormick, “Approximate two-parameter phase function for light scattering,” J. Opt. Soc. Am. 70, 1206 (1980). [CrossRef]
45. A. N. Bashkatov, E. A. Genina, and V. V. Tuchin, “Optical properties of skin, subcutaneous, and muscle tissues: A review,” J. Innov. Opt. Heal. Sci. 04, 9–38 (2011). [CrossRef]
46. S. L. Jacques, “Optical properties of biological tissues: A review,” Phys. Med. Biol. 58, 5007 (2013). [CrossRef]
47. B. W. Pogue and M. S. Patterson, “Review of tissue simulating phantoms for optical spectroscopy, imaging and dosimetry,” J. Biomed. Opt. 11, 041102 (2006). [CrossRef] [PubMed]
48. R. Michels, F. Foschum, and A. Kienle, “Optical properties of fat emulsions,” Opt. Express 16, 5907 (2008). [CrossRef] [PubMed]
49. B. Yu, A. Shah, V. K. Nagarajan, and D. G. Ferris, “Diffuse reflectance spectroscopy of epithelial tissue with a smart fiber-optic probe,” Biomed. Opt. Express 5, 675–689 (2014). [CrossRef] [PubMed]
50. A. Eshein, W. Wu, T.-Q. Nguyen, A. J. Radosevich, and V. Backman, “A fiber optic probe to measure spatially resolved diffuse reflectance in the sub-diffusion regime for in-vivo use,” Proc. SPIE 9703, 970317 (2016). [CrossRef]
51. L. Wang, S. L. Jacques, and L. Zheng, “Conv—convolution for responses to a finite diameter photon beam incident on multi-layered tissues,” Comput. Meth. Prog. Bio. 54, 141–150 (1997). [CrossRef]
