Virtually increased acceptance angle for efficient estimation of spatially resolved reflectance in the subdiffusive regime: a Monte Carlo study

: 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 eﬃcient MC simulation of spatially resolved reﬂectance 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 conﬁgurations.


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][2][3][4][5][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][8][9][10][11][12][13] or imaging systems [14][15][16][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 − g 1 ) coefficient, which depends on the scattering coefficient µ s and the first Legendre moment g 1 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 − g 2 )/(1 − g 1 ), comprising the first g 1 and second Legendre moments g 2 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.

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 (n m = 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 (W tot ) exceeded a predefined value: where w i is the remaining weight of the i-th collected backscattered photon packet. For absorbing media w i is less then 1 and the number of collected photon packets N collected exceeds W tot . In this study, the value of W tot was set to 10 6 and was realized through consecutive runs of MC simulations by a fixed number of launched photon packets (10 7 ).

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  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 g gk 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 N r through all the SDSs of the SRR, R m (r) is the measured SRR of the sample and R LUT (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 g 1 was held constant at 0.8.

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 N r 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 E t (Fig.  2(b)). For the purpose of this study, the threshold E t was set to 1%.

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/e 2 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/e 2 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.

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 E t . 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. 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.

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 γ, y i 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].
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 .
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.

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 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  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 W tot . In this study, W tot was fixed to 10 6 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.    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.

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 thirteenfold 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.