1. Introduction

Many studies in biomedical optics have employed diffuse reflectance spectroscopy for the non-invasive characterization and analysis of biological tissues. The technique is simple and straightforward in its application and can be implemented with relatively inexpensive instrumentation. The measured spectra contain information about the absorption and scattering properties of tissues under study, which can potentially be used to obtain biochemical, morphological, and histochemical information. However, no exact analytical solution exists for the problem of light propagation in turbid media such as biological tissues. Early attempts to model diffuse reflectance from biological tissues were characterized by cumbersome analytical expressions based on the diffusion approximation [1,2]. Recently, several studies have introduced relatively simple models, which are designed for specific probe geometries [3,4], with some of them being empirical or semi-empirical in nature while others are based on a relatively simple solution of the diffusion approximation utilizing the method of images [5]. On the other hand, numerical methods such as Monte Carlo (MC) simulations [6] are accurate and they are not restricted to any particular light delivery and collection geometry but require long computation times and they are difficult to apply in the inverse mode since they provide no analytical solution. Many researchers have tried to address these shortcomings by coming up with various ways to accelerate the execution time of MC simulations [7–9]. More recently, significant acceleration using graphics processing units (GPU) has also been possible [10–12]. Overall, these improvements yield an advantage of several orders of magnitude faster execution, making MC simulations a very appealing technique. At least one study has recently been published utilizing a MC method in the forward mode in lieu of an analytical solution [13].

In this work, MC simulations are used to solve the problem of diffuse reflectance from a homogeneous semi-infinite turbid medium for arbitrary detector diameter geometry. To overcome the lack of an analytical solution using the MC method, a simple semi-empirical solution is developed based on the results of the MC simulations, which can be used to describe the MC results quite well. In addition, several interesting properties of diffuse reflectance are discussed and pointed out in terms of the semi-empirical modeling approach and the particular detector diameter. The overall modeling approach and results may be useful for the further modeling of reflectance from more complicated tissue geometries, including multi-layer morphological structures.

2. Methods

2.1. Monte Carlo simulations

Simulations of light propagation in turbid media such as biological tissues were performed using the publicly available MCML code [6]. A one-layer semi-infinite tissue model was assumed (implemented as a 1000 mm thick tissue layer) and incident light was modeled as an infinitely thin beam perpendicular to tissue. Reflectance was scored up to a distance of 10 mm away from the incident beam with a step size of 0.5 mm so as to simulate a wide range of delivery/collection separations corresponding to detectors with various diameters. Reflectance at various detector radii was calculated by integrating the total reflectance from r = 0 up to the actual radius r of the detector. Total reflectance (corresponding to r = ∞) was also scored. Reflectance was scored at various photon exit angles (with respect to the perpendicular on the tissue surface) corresponding to detectors with different numerical apertures. All possible exit angles from 0 to 90 degrees were covered in 6 different bins, each one 15 degrees wide. The refractive index of tissue was assumed equal to 1.4 and the refractive index of the medium above was assumed equal to 1.0. For simplicity, the refractive index of the detector was also assumed equal to 1.0 and was otherwise assumed to be an ideal detector (i.e. with no absorption, reflection, etc.). The scattering anisotropy coefficient of tissue was assume to be g = 0.9 (Heaney-Greenstein scattering phase function [6]). Table 1 summarizes the distinct values of the absorption and reduced scattering coefficient used in the simulations. These cover the entire typical relevant range for soft biological tissues. There were 5 × 13 = 65 runs of the MCML code in total, each one for each pair of absorption and reduced scattering coefficient values. Every MCML run was performed with 10 million incident photons. All reflectance results are reported as percentage (%) of the total number of incident photons. Figure 1 illustrates the delivery/collection geometry employed in the simulations.

Table 1. Optical parameter values used in the MC simulations

View Table | View all tables in this article

 

Fig. 1 Illumination/collection geometry employed in the simulations.

Download Full Size | PDF

2.2. Data analysis

Data analysis and fitting was performed using the TableCurve2D software package (Systat) which utilized the Levenberg-Marquardt minimization method. TableCurve2D enabled the identification of the simplest and more appropriate functions to describe the data by fitting the data to a library of thousands of trial functions and categorizing the results according to the goodness of fit. The exact functional form of Eqs. (3)–(8) was determined in this way.

3. Results

3.1. General model

Diffuse reflectance calculated using the MCML code for various detector diameters is shown in Fig. 2 as a function of the reduced scattering coefficient. Note that for smaller detector diameters, reflectance exhibits an approximately linear dependence on the reduced scattering coefficient and this holds for all values of the absorption coefficient. The solid lines in Fig. 2 represent a fit to the MC data using the model described by Eq. (1),

 

Fig. 2 Reflectance as a function of the reduced scattering coefficient (MC data) for various detector radii, r, and various values of the absorption coefficient. (a) ${\mu }_a$ = 0.001 mm⁻¹, (b) ${\mu }_a$ = 0.01 mm⁻¹, (c) ${\mu }_a$ = 0.1 mm⁻¹, and (d) ${\mu }_a$ = 1 mm⁻¹. Solid lines represent a fit to the model described by Eq. (1).

Download Full Size | PDF

In Eq. (1), R is the reflectance while ${\mu }_a$and ${{\mu }^{\prime }}_s$ are the absorption and reduced scattering coefficients, respectively. This is a simple expression that can be derived from basic physics [14] for the description of total reflectance (i.e. infinite detector radius). Here, the parameters k₁ and k₂ depend on the absorption coefficient and detector radius such that Eq. (1) can be used to describe the reflectance collected by an optical detector with arbitrary radius. These parameters can be determined by fitting the MC data to Eq. (1) (Fig. 2). Figure 3 shows the resulting values of k₁ and k₂ as a function of detector radius for different values of the absorption coefficient. Overall, the dependence of both parameters on the detector radius and the absorption coefficient is quite smooth and once partially determined, both k₁ and k₂ can be fully determined using spline interpolation and a lookup table technique [15] to evaluate the reflectance collected using an optical detector of arbitrary radius over the entire range of optical parameter values relevant for biological tissue.

 

Fig. 3 k₁ and k₂ parameter values obtained from the fitting results shown in Fig. 2, using Eq. (1). The solid lines represent spline interpolation between the data points.

Download Full Size | PDF

A more straightforward modeling approach can be identified by focusing on the dependence of reflectance on absorption. Figure 4 shows the reflectance for detectors with various diameters as a function of the absorption coefficient with the solid lines representing MC data fits to Eq. (2),

 

Fig. 4 Reflectance as function of the absorption coefficient for detectors with various radii. The solid lines represent fits to the MC data using Eq. (2) (${{\mu }^{\prime }}_s$ = 1.8 mm⁻¹ for the data shown here).

Download Full Size | PDF

As can be seen, Eq. (2) describes the MC data very well. All data shown in Fig. 4 are for ${{\mu }^{\prime }}_s$ = 1.8 mm⁻¹ with different values of the reduced scattering coefficient producing similar results. This confirms that the reflectance when using an optical detector with finite diameter, exhibits a simple inverse dependence on the absorption coefficient. We have also previously discovered this dependence experimentally using a small-diameter optical probe [14]. This observation simplifies modeling of reflectance considerably. What is needed in order to completely determine a general reflectance model based on Eq. (2) is to establish the exact dependence of the two parameters k₁ and k₂ on the reduced scattering coefficient and the detector radius r, i.e., $k_1({{\mu }^{\prime }}_s,r)$and$k_2({{\mu }^{\prime }}_s,r)$.

Figure 5 shows the values of these last two parameters as a function of the detector radius. The solid lines in Fig. 5 represent fits to Eqs. (3) and (4) for k₁ and k₂, respectively:

 

Fig. 5 Dependence of the k₁ and k₂ parameters of Eq. (2) on the detector radius. The solid lines represent fits to Eqs. (3) and (4) for k₁ and k₂, respectively.

Download Full Size | PDF

Equations (3) and (4) include four parameters to be determined:$a_1$,$b_1$,$a_2$and $b_2$. The last step in fully developing the model is to determine the dependence of these four parameters on the reduced scattering coefficient. It turns out that the dependence is again simple enough such that concise expressions can be found describing these parameters (Fig. 6 and Eqs. (5)–(8)):

 

Fig. 6 Dependence of the parameters$a_1$,$b_1$,$a_2$and $b_2$ of Eqs. (3) and (4) on the reduced scattering coefficient. The solid lines represent fits to Eqs. (5)–(8).

Download Full Size | PDF

Thus, an analytical expression for the reflectance using a detector with arbitrary radius can be fully determined in terms of Eq. (2) and eight empirically determined parameters which are summarized in Table 2 . This expression combines the simplicity of a concise analytical expression with the accuracy of the MC method. It is important to point out that due to the semi-empirical nature of the model it is only valid within the optical parameter ranges specified in Table 1: ${\mu }_a$ = 0.001–10.0 mm⁻¹ and ${{\mu }^{\prime }}_s$ = 0.6–3.0 mm⁻¹.

Table 2. Model parameter values

View Table | View all tables in this article

3.2. Small diameter detector model

As already mentioned (Fig. 2 above) the dependence of reflectance on the reduced scattering coefficient is approximately linear for small detector diameters. This effect has been previously observed experimentally both by us and other researchers [3,4,14,16] and simplifies greatly modeling of diffuse reflectance when using small diameter detectors. To model reflectance from such a small diameter detector, an even simpler analytical expression can be found described by Eq. (9):

Figure 7 shows that Eq. (9) can describe the reflectance when using a detector with 0.5 mm radius remarkably well. In contrast to the general model developed in the previous section, the parameters k₁ and k₂ have fixed constant values and they are independent of the reduced scattering coefficient (k₁ = 6.1 and k₂ = 1.55 in this case). In addition, Fig. 7(a) demonstrates again here the inverse dependence of reflectance on the absorption coefficient, which was also pointed out in the previous section (Fig. 4). It is interesting to note that we have experimentally observed previously that Eq. (9) describes very well the reflectance measured with a standard six-around-one 200 micron optical fiber probe [14].

 

Fig. 7 Diffuse reflectance scored using an optical detector with a radius of 0.5 mm showing the dependence of reflectance on (a) the absorption coefficient (reflectance shown is normalized to unity at the smallest value of the absorption coefficient), and (b) the reduced scattering coefficient. Solid lines represent fits to Eq. (9).

Download Full Size | PDF

Another interesting observation can be pointed out here regarding the inverse dependence of reflectance on the absorption coefficient. It can be seen from Figs. 4 and 7(a) that this dependence (expressed by Eqs. (1), (2), and (9)) describes the reflectance remarkably well for detector radii smaller than 5 mm but starts to deviate from the MC data for larger detector diameters. In this large detector diameter case, and also for a detector with infinite diameter (total reflectance collected) the MC data can be better described by Eq. (10)

Figure 8 shows an example of the excellent description that Eq. (10) provides for the total reflectance (r = ∞) compared to Eq. (2), which does not perform that well for large detector diameters and total reflectance. Equation (10) is actually very similar to that obtained analytically for the total diffuse reflectance using the path integral method (see review in [14] and references therein). The parameter values used to obtain the fits shown in Fig. 8 were k₁ = 82.66 and k₂ = 15.65 in the case of Eq. (2) and k₁ = 94.96 and k₂ = 3.31 in the case of Eq. (10).

 

Fig. 8 Reflectance for a detector with infinite radius and best fits using Eqs. (2) and (10). Equation (10) gives a better fit compared to Eq. (2), and this is also true for large detector diameters (${{\mu }^{\prime }}_s$ = 1.8 mm⁻¹ for the data shown here).

Download Full Size | PDF

3.3. Angular dependence effects

All analysis and modeling work presented in the previous sections is based on the reflectance integrated over all exit angles. Some optical detection schemes such as those employing optical fibers typically collect diffuse reflectance within a limited exit angle range depending on the actual numerical aperture of the optical fibers used in probe construction as well as the tissue refractive index. Thus, maximum collection angles may lie typically in the 10–20 degrees range, measured from the perpendicular to the tissue surface, while other detectors such as integrating spheres are theoretically capable of collecting all light exiting up to angles of 90 degrees. Most of the published diffuse reflectance models neglect exit angle effects [3–5] and in this section, we set out to evaluate such effects using our MC results.

Figure 9 shows the differences between a fiber optical probe detector and an integrating sphere detector. The vertical axis “efficiency” shown in Fig. 9 is the ratio of light collected within the 0–15 degree exit angle range versus the full 0–90 degree range. As can be seen, this efficiency is approximately constant for small values of the absorption coefficient but starts dropping when absorption becomes comparable to the reduced scattering coefficient. But even in the high absorption case the spectral distortions introduced are not higher than about 15% and they are more prominent for smaller values of the reduced scattering coefficient, as expected.

 

Fig. 9 Effects of detector numerical aperture as a function of the absorption coefficient for three different values of the reduced scattering coefficient: (a)${{\mu }^{\prime }}_s$ = 0.6 mm⁻¹, (b)${{\mu }^{\prime }}_s$ = 1.8 mm⁻¹, and (c) ${{\mu }^{\prime }}_s$ = 3.0 mm⁻¹. The efficiency parameter shown in the vertical axis represents the percent ratio of reflectance scored in the 0–15 deg exit angle range over the total reflectance scored over the entire 0–90 deg exit angle range. Solid lines represent spline interpolation between data points.

Download Full Size | PDF

4. Discussion

The problem of diffuse reflectance from a semi-infinite turbid medium has been addressed by many researchers in the field of biomedical optics. However, to the knowledge of the authors, there is no detailed treatment available based on MC calculations for a detector scheme with arbitrary diameter. The modeling approach presented here is characterized by its accuracy (an attribute of the MC method), its generality, and its simplicity. In addition, it goes beyond the restrictions of the MC method by providing a simple semi-empirical analytical expression for the reflectance.

The modeling geometry employed in this work is potentially useful for modeling realistic complex optical probe geometries. A similar approach has been successfully employed in the past [17] to model reflectance collected by fiber optic probe (using an “effective” maximum radius that best models the probe). However, that work was based on diffusion theory [5]. Compared to that work, the present work presents significant advantages mainly due to the fact that it is not limited by the inaccuracies and applicability limitations of the diffusion approximation.

It is interesting to compare the present work with another MC-based modeling technique [18]. This scaling MC method offers significantly faster computational times without sacrificing the modeling accuracy associated with standard MC methods. The advantage of our approach compared to such a method is the use of an analytical expression, albeit semi-empirical in nature, while retaining the accuracy associated with the MC method. The disadvantages of our approach lie in the restrictions associated with the specific illumination/collection geometry and the semi-infinite tissue approximation. However, these disadvantages do not seem to impose serious limitations because the modeling approach presented here can be potentially extended to more general delivery/collection geometries as well as multilayer tissue geometries by following a similar methodology in model development. We hope that we will be able to report on such generalizations of the present modeling approach in the future.

The specific delivery/collection geometry of the present model (i.e. an infinitely thin delivery beam and idealized detector with fixed radius) may at first appear far from realistic delivery/collection schemes. Nevertheless, we must point out that it may be quite useful in many practical situations. As already mentioned above, a similar approach has been successfully applied using a fiber-optic probe [17] and other similar approaches have also been reported [14,19,20]. In addition, the present approach may be useful in modeling integrating sphere detectors and/or delivery/collection geometries with free space illumination and a camera for reflectance collection.

Regarding the similar previous published models, it is worth mentioning those that employ a modeling approach similar to that described by Eq. (9) ([14,21–25] by us and [19,20] by other researchers) to describe reflectance collected using optical fiber probes. In particular, the simple inverse dependence of reflectance on the absorption coefficient (illustrated in Figs. 4 and 7(a)) is worth mentioning because it greatly simplifies modeling and, to the authors’ knowledge, it has not been previously reported by other researchers. This, in combination with the simple linear dependence of reflectance on the reduced scattering coefficient (which has been reported and confirmed by other researchers previously) makes for the very simple expression described by Eq. (9). Equations (1), (2), and (9) may indeed seem very simple and, in fact, they may even look too simplistic to some researchers. Nevertheless, they appear to be quite promising and effective and we would like to urge researchers in the field of biomedical optics to investigate their utility further. At least two other publications [19,20] have already reported on the successful application of a reflectance model similar to that described by Eq. (9).

Equation (1) relies on a look-up table technique which has already been suggested by other researchers as a convenient modeling tool in tissue reflectance [15]. In addition, it is possible to substitute the lookup table technique with an empirical fit of the data as is illustrated for Eq. (2). The problem in such a case is to identify the appropriate analytical expression that best describes the data. We have shown here that it is possible to develop such a semi-empirical model based on Eq. (2) using this approach, as we have also done previously for a two-layer reflectance model [25]. In both cases, the complexity of the model is significantly lower and the ease of model application and accuracy are potentially greater compared to other available reflectance models such as those based on diffusion theory [5,17].

Investigation of the angular dependence of diffuse reflectance revealed that the finite numerical aperture of a fiber optic probe does not introduce any appreciable spectral distortions as long as absorption is lower that scattering. Even in the limit of high absorption (which is sometimes present in biological tissues) the error introduced is no larger than about 15%. These results are consistent with the practice of many researchers who neglect exit angle effects in their modeling by using models that are directly based on the total reflectance at all exit angles [3–5]. When accurate measurements are required though, and absorption is high, the angular dependence of reflectance must be taken into account.

Another important point is that the same modeling approach described here using Eqs. (1), (2), and (9) can be applied to experimental data on tissue phantoms. In this way, a more accurate, reliable, and realistic calibration of the model can be achieved by either tuning the actual model parameter values or even coming up with new model analytical expressions such as those described by Eqs. (1), (2), (9), and (10). We have already employed this modeling approach using tissue phantoms in the development of both one and two-layer tissue reflectance models [14,21,25] and we hope that this approach will serve as the basis for more accurate and robust modeling of tissue reflectance in the future.

5. Conclusions

This work presented a simple semi-empirical tissue reflectance model useful for an optical detector with arbitrary diameter. The development of the model is based on MC simulations results and thus it is not restricted by the limitations of the diffusion approximation. In addition, this modeling approach utilizes two simple characteristics of the reflectance: (1) an inverse dependence on the absorption coefficient, and (2) a linear dependence on the reduced scattering coefficient (for small diameter detectors). We expect that the modeling approach presented here will facilitate further tissue analysis and investigations using reflectance spectroscopy and that it will be useful to researchers in the field of biomedical optics.