1. Introduction

Diffuse reflectance spectroscopy has been extensively investigated as a promising tool for the non-invasive study, analysis, and characterization of biological tissues. Diffuse reflectance provides information about the absorption and scattering properties of light by the tissues under study. Correlation of these optical properties with histological and biochemical characteristics provides a wealth of information regarding tissue composition and biochemistry. This information may be also useful for the study of tissue morphology and physiology as well as for the early detection and diagnosis of disease.

Biomedical applications of diffuse reflectance spectroscopy include a wide variety of organs and tissue types such as colon [1–2], esophagus [3], stomach [4], bladder [5–6], cervix [7–8], ovaries [9], breast [10–11], brain [12], liver [13], pancreas [14], heart [15], oral tissues [16], and skin [17–18]. Based on diffuse reflectance, a number of biological tissue chromophores may be quantitatively assessed such as hemoglobin in its oxygenated and deoxygenated forms [2, 19], melanin [17], bilirubin [20] and cytochromes [13, 21]. In addition, light scattering properties of biological tissues may be determined and correlated with morphological characteristics such as epithelial cell nuclear size and density [22], or average scatterer size and density [2]. Furthermore, a number of physiological processes and characteristics may be monitored such as microcirculation [4, 14, 19], hemodynamics [23], erythema [24], and metabolism [21]. Finally, pathological conditions may be diagnosed such as precancerous lesions in various organs [1–3, 7, 8, 12, 16].

A typical experimental setup for measuring diffuse reflectance from biological tissues includes a light source, a spectrograph, a detector, and a fiber optic probe for delivery and collection of light. The measured diffuse reflectance spectra contain entangled information regarding the scattering and absorption of light in tissue. Spectral analysis requires a model describing the diffuse reflectance in terms of these optical properties. In particular, an inverse application of a model is needed i.e. starting from the diffuse reflectance spectra as input, the model should be capable of providing the absorption and scattering properties of the tissue under study. Even though there are several ways to model diffuse reflectance from biological tissues, most of them involve approximations to the radiative transport equation (such as diffusion theory) which incorporate potential inaccuracies and deviations into the final model [25–26]. One other frequently used approach involves the use of Monte Carlo simulation of light propagation in tissue [27], which, however, provides a solution in numerical form only and often requires long computational times.

In this article, starting from basic physical principles underlying the diffuse reflectance process, we propose a new modeling technique appropriate for the analysis of diffuse reflectance spectra from biological tissues. The model is characterized by its simplicity and effectiveness, requires only a one-time calibration on a tissue physical model (phantom) and can be adapted to various optical probe geometries. To minimize mathematical complexity, tissue is modeled as a semi-infinite turbid medium. Validity of the model is tested on tissue phantoms and its performance is further demonstrated with the analysis of diffuse reflectance spectra from human skin in vivo.

2. Diffuse reflectance theory

2.1 Simple “exponential model” for total diffuse reflectance

A simple physical picture of the diffuse reflectance process is outlined in Fig. 1. Light is incident on the surface of tissue, which is assumed to be a homogeneous semi-infinite turbid medium characterized by absorption and reduced scattering coefficients, µ_a and µ′_s, respectively. The tissue is also characterized by a refractive index, n, which usually lies in the range 1.35–1.45 for soft biological tissues i.e. it is slightly higher than the refractive index of water [28–30]. Part of the incident light is absorbed in tissue while the non-absorbed part is subjected to multiple scattering and eventually emerges from the tissue surface as diffuse reflectance.

 

Fig. 1. Physical picture of the diffuse reflectance process: light is incident on tissue at point x_i and after several successive scattering events the unabsorbed fraction of the incident light emerges from the tissue surface as diffuse reflectance at point x_f . The smooth line describes the “classical path” or “most probable path” of the light according to the path integral formalism.

Download Full Size | PDF

With R_t denoting the total diffuse reflectance from a semi-infinite medium, it can be easily shown, based on scale invariance arguments, that R_t =R_t (µ′_s /µ_a). The use of an optical probe for measuring the diffuse reflectance, R_p , breaks the scale invariance and introduces explicit dependencies on the reduced scattering and absorption coefficients, e.g. R_p =R_p (µ′ _s /µ _a ,µ′ _s ,µ _a ). The latter is the quantity of interest in realistic situations i.e. with the use of an optical probe. Nevertheless, investigation of R_t provides useful insight into the physics of the diffuse reflectance problem and can serve as a basis for modeling R_p .

The simplest way to calculate R_t from a homogeneous semi-infinite turbid medium is by assuming that the incident and backscattered light are exponentially attenuated according to the sum of the absorption and reduced scattering coefficients, (µ_a+µ′ _s ). Even though this approximation essentially ignores the integral term in the radiative transport equation [31], it is well established that exponential solutions of this type constitute an acceptable approximation, especially in one-dimensional geometries [31–32].

 

Fig. 2. Geometry of the “exponential” diffuse reflectance model. The intensity of the light propagating in tissue is exponentially attenuated according to the sum of the reduced scattering and absorption coefficients. Diffuse reflectance is scattered back being proportional to the reduced scattering coefficient.

Download Full Size | PDF

Figure 2 illustrates the geometry of the problem. A semi-infinite turbid medium is considered with z=0 on the surface where light is incident with an intensity I ₀ and is exponentially attenuated as it propagates within the medium i.e. $I=I_0{e}^{-\left({\mu }_a+{\mu }_s^{\prime }\right)z}$. At a given depth z, light is scattered back from a thin layer of thickness dz. The backscattered light is assumed proportional to the reduced scattering coefficient µ′ _s and travels back to the surface while being exponentially attenuated by a factor exp(-(µ _a +µ′ _s )z). The total diffuse reflectance is then given by:

Equation (1) is normalized such that R_t =1 for µ_a=0 and exhibits the correct functional dependence i.e. it is a function of µ′ _s /µ _a , as required by the scale invariance of the problem. In addition, it exhibits the correct dependence in the limits of very high absorption (R_t →0) and very low absorption (R_t →1).

Even though this simple derivation succeeds on a qualitative basis, to reach an acceptable quantitative agreement with other more accurate models a correction must be introduced by means of rescaling the reduced scattering coefficient by an empirical factor k. If the substitution µ′ _s →µ′ _s /k is made Eq. (1) becomes

Equation (2) constitutes a reasonable approximation in the range of optical properties that is of interest for biological tissues i.e. µ′ _s /µ _a =1-1000. Thus, this “exponential” model is noteworthy for its effectiveness and simplicity.

Equation (2) is valid when there is a refractive index match between the semi-infinite turbid medium and the surrounding medium. Assuming that the surrounding medium is air, with refractive index approximately equal to unity, it is possible to derive a general expression relating the total diffuse reflectance under refractive index mismatch conditions to that under index match conditions. This can be implemented by appropriately modifying the parameter k, allowing for a dependence on the refractive index.

Refractive index mismatch conditions cause total internal reflection for incident angles greater than the critical angle, θ _c =sin^-1(1/n), with respect to the vertical on the surface. This is equivalent to an additional external source of light with intensity equal to that of the internally reflected light, which in turn generates additional diffuse reflectance, part of which is internally reflected to create a new source, and so on. Figure 3 illustrates this process. The total diffuse reflectance, R_t , for the refractive index mismatch case may be obtained in terms of the index match total diffuse reflectance, R_t , and the fraction, q, of the diffuse reflectance which is internally reflected [34],

It is possible to calculate q by assuming an isotropic specific intensity for the light incident on the inner boundary [31, 35]. This is expressed by the cosθ term in Eq. (4):

Because some light is reflected back even for angles smaller than the critical angle θ _c Eq. (4) constitutes only an approximation. Exact values for the parameter q have been calculated and tabulated taking into account all relevant effects unaccounted for in Eq. (3). A more accurate expression may be derived in terms of an empirical fit [36]:

However, Eq. (4) is in good general agreement with Eq. (5), typically within 5%, which is usually acceptable for modeling diffuse reflectance from biological tissues.

 

Fig. 3. Calculation of the diffuse reflectance under refractive index mismatch conditions with q being the fraction of the internally reflected light acting as an additional source. This source creates additional diffuse reflectance, part of which undergoes total internal reflection to generate an additional source, and so on. The sum over infinite steps of this process gives the total diffuse reflectance, R_tm , under refractive index mismatch conditions.

Download Full Size | PDF

Equation (2) can be now generalized for the refractive index mismatch case. An expression may be derived for the dependence of the parameter k on the tissue refractive index. Let k_m denote the generalized k parameter for the refractive index mismatch case. From Eq. (3) follows that

Equations (2), (3), and (6) yield

Equation (6) can be written as

As expected, the refractive index mismatch condition tends to decrease the total diffuse reflectance. Equation (8) also exhibits the correct behavior in the limiting case n→∞ where R_tm →0 (all light internally reflected). Thus, Eqs. (2) and (8) share the same functional form and refractive index mismatch conditions may be implemented in a straightforward way with a simple rescaling of the parameter k.

2.2 Other models

Equations (2), (6) and (8) provide simple expressions for the total diffuse reflectance from semi-infinite turbid media. Other simple models also exist and it is instructive to examine and compare these models in more detail.

First, we consider a photon migration model [33]. This model employs a probabilistic concept to describe absorption and scattering of light. Photons are assumed to come in “bunches” or “packets”, with a fixed initial weight (typically equal to unity), which is gradually reduced by subsequent absorption events. Furthermore, a set of functions, f_n , is defined describing the probability of a photon packet to escape as diffuse reflectance after n scattering events in a turbid medium with zero absorption. The total diffuse reflectance is then given by Eq. (9),

with the normalization condition ${\sum }_{n=1}^{\infty }$ f_n =1 i.e. all photons escape as diffuse reflectance when there is zero absorption. Hence, f_n depends on the scattering properties of the medium only and it is independent of absorption. By making the reasonable approximation that f_n has an exponential dependence on n, such that f_n =e ^-n/k/k with k a constant, and by inserting this in Eq. (9), we arrive at the following expression for the total diffuse reflectance:

For µ′_s≫µ_a, Eq. (10) reduces to

Eq. (11) is identical to Eq. (2), which was derived in a distinctly different manner. Just as in the previous case, the constant k is not explicitly evaluated here either; it is still a free parameter to be independently determined.

Next, we consider a path integral model [37]. Just like the photon migration model, the path-integral formulation treats light scattering as a probabilistic process. The concepts of photon weight and zero-absorption photon path are used again, but the zero-absorption path probability, f_n , is here explicitly evaluated using the Feynman path integral formalism [38]. According to this formalism, the concept of the classical path is introduced, which is defined as the most probable path for a photon entering at a point x_i and exiting at a point x_f (Fig. 1). The probability for this path to occur may be written as a path integral, as follows:

where the scattering phase function is assumed Gaussian,

x and z are the coordinates spanning the surface and the depth in the medium respectively, θ is the scattering angle, t is the time and T is the total travel time for the photon inside the medium. By evaluating the path integral in Eq. (12) an expression for the photon escape probability, f_n , may be found,

with g=<cosθ> being the anisotropy of scattering. Combination of Eqs. (9) and (14) yields an expression for the total diffuse reflectance,

where β=2.68 is a constant. Note the different functional form of Eq. (15) as compared to Eqs. (2) and (11). In addition, the scattering anisotropy, g, and scattering coefficient, µ _s , appear explicitly in this expression, instead of being combined into µ′ _s .

Finally, we turn to diffusion theory. The diffusion approximation to the radiative transport equation is a well studied technique for calculating the diffuse reflectance from turbid media [31, 35] and has extensive biomedical applications. An analytical solution for the 3-dimensional diffusion equation in the case of an infinitely thin light beam incident on a semi-infinite turbid medium may be found using the Green’s function method [25, 39]. However, this method is mathematically complex and the final expression requires numerical evaluation. One of the simplest implementations of diffusion theory employs the method of images [26] according to which the total diffuse reflectance, R_t , is given by

with µ²=3(µ′ _s +µ _a )µ _a , z₀=1/(µ′ _s +µ _a ), and A=2n ²-1.

2.3 Model comparison

We now proceed to evaluate the relative performance of the models presented so far in the range 0.1<µ′ _s /µ _a <1000, typical for biological tissues.

To facilitate comparison of the models we also calculate the diffuse reflectance using the Monte Carlo technique [27]. This technique traces photons as they undergo multiple scattering and absorption within biological tissue and numerically evaluates the portion of the incident light which is remitted as diffuse reflectance. The Monte Carlo (MC) method involves no major approximations and is potentially more accurate compared to the models described above. The publicly available MC code [27] we employed uses the Henyey-Greenstein phase function and the asymmetry parameter was set to g=0.9 which is typical for soft biological tissues.

Figure 4 shows R_t as a function of µ′ _s /µ _a . Note that for all models, R_t →1 for µ′ _s /µ _a →∞, as expected. Also, for µ′ _s /µ _a =0 all models correctly predict R_t =0. Assuming the MC model to be the gold standard, we proceed to compare the four models. Figure 5 shows such a comparison with the MC model. The relative performance of each model is presented, divided by the MC results. Approximately 25 Monte Carlo simulations were conducted to cover the range shown in Figs. 4 and 5. The trace shown is a spline interpolation between the MC data points which is very accurate since the curve is smooth. For the exponential model, k=13 was used because this value provides a reasonable agreement with the other models and facilitates further comparison. For the path integral model g=0.9 was used as in the case of the MC simulations.

 

Fig. 4. Total diffuse reflectance as a function of the µ′ _s /µ _a , as calculated by (a) the diffusion model (red dashed line), (b) the exponential model (green dotted line), (c) the path integral model (blue dashed-dotted line), and (d) the Monte-Carlo model (black solid line). All models appear to be in good general agreement.

Download Full Size | PDF

 

Fig. 5. Total diffuse reflectance as a function of the µ′ _s /µ _a , as predicted by the models presented in Fig. 4. The curves shown are divided by the Monte-Carlo curve, in order to illustrate differences between the models. (a) diffusion (red dashed line), (b) exponential (green dotted line), and (c) path integral (blue dashed-dotted line). All models are in good general agreement as long as absorption is significantly smaller than scattering.

Download Full Size | PDF

Despite the simple and approximate way in which it was derived, the exponential model performs well, with a maximum deviation not exceeding approximately 20%. The diffusion model agrees very well with the MC model for small absorption; this is where the diffusion approximation provides its highest accuracy [31]. As expected, the largest deviation for the diffusion model is observed for high absorption which is, though, rarely encountered in biological tissues for visible and NIR light. Finally, the path integral model exhibits superior agreement with MC (down to 0.1%) for low absorption, but it deviates significantly in the high absorption region.

2.4 Model development

Following the short review of different approaches to diffuse reflectance modeling techniques, we make the proposition that the exponential model can serve as the simplest base for the development of a diffuse reflectance model. Eq. (2), however, is not suitable for describing diffuse reflectance measured with a fiber optic probe. This is because an optical probe measures only part of the diffuse reflectance exiting the tissue surface and hence breaks the scale invariance inherent in Eq. (2). A straightforward way to remove scale invariance is by introducing an explicit dependence on µ′ _s in Eq. (2) as follows:

where k ₁ and k ₂ are parameters depending on the geometrical characteristics of the optical probe as well as on the refractive indices of tissue and the surrounding medium and can be determined experimentally using a tissue physical model.

To model skin spectra, the absorption and reduced scattering coefficients were related to hemoglobin concentration, c* _Hb , hemoglobin oxygen saturation, α, melanin concentration, c_mel , and effective scatterer size, d_s , as follows [2, 17]:

The molecular extinction coefficient spectra of oxyhemoglobin, deoxyhemoglobin, and melanin are ${\epsilon }_{Hb{O}_2}\left(\lambda \right)$, ε _Hb (λ), ε _mel (λ) respectively. Equation (19) describes a linear dependence of the reduced scattering coefficient on the wavelength. This is a reasonable approximation which is supported by Mie theory calculations for spherical scatterers with Gaussian distribution in size [17]. The maximum and minimum wavelengths are λ_max=900 nm and λ_min=450 nm, respectively, and d ₀=0.0625µm is a constant.

3. Methods

3.1 Instrumentation

The experimental setup employed is shown in Fig. 6. A CCD spectrophotometer (USB2000, Ocean Optics, Dunedin, FL, USA), capable of collecting spectra in the 450–900 nm range was employed for diffuse reflectance measurements. Light was delivered and collected by means of a fiber optic probe (R200-7, Ocean Optics, Dunedin, FL, USA) which consisted of six 200 µm core diameter optical fibers for delivery and a single 200 µm core diameter optical fiber for collection. The delivery fibers were placed around the central collection fiber such that the center to center distance of the optical fibers was approximately 200 µm. Illumination was provided by a tungsten-halogen light source (HL-2000-FHSA-LL, Ocean Optics, Dunedin, FL, USA). All spectra were referenced to a reflectance calibration standard (WS-1, Ocean Optics, Dunedin, FL, USA). Diffuse reflectance spectra were normalized by dividing by the diffuse reflectance spectrum of the calibration standard. Typical integration time for spectral collection was 80 ms and signal to noise ratio was typically greater than 100:1.

 

Fig. 6. Experimental setup for diffuse reflectance measurements, which includes a tungsten-halogen light source, a compact CCD spectrophotometer, a fiber optic probe and a computer for data acquisition.

Download Full Size | PDF

3.2 Tissue physical model

Tissue phantoms were constructed using 1.0 µm diameter polystyrene beads (Polysciences, Inc, USA) and hemoglobin from bovine blood (Sigma, H2500), which was predominately methemoglobin. The beads were suspended in deionized water and hemoglobin solutions were prepared in distilled water. The scattering properties of the beads were calculated using Mie theory [40] and the absorption spectra of the hemoglobin solutions were calculated using the molar extinction coefficient of methemoglobin [41]. Diffuse reflectance measurements were made by placing the probe in contact with the phantom surface.

3.3 In vivo data

In vivo data were collected on seven healthy adult volunteers with skin type III. Diffuse reflectance spectra were measured on several normal skin sites on the volar forearm and also on melanocytic nevi found on the same volunteers. Spectra were measured by placing the optical probe in soft contact with the skin.

4. Results

Determination of k ₁ and k ₂ was performed on tissue phantoms. Figure 7 shows the measured dependence of diffuse reflectance on µ′ _s and µ _a . Note the approximately linear dependence of R_p on µ′ _s , which is an observation consistent with the findings of other studies employing fiber optic probes [42–43], while the dependence of R_p on µ _a shows a monotonically decreasing behavior as µ _a increases. The data shown in Fig. 7 can be described well by Eq. (17). A best fit of the phantom data to Eq. (17) yielded k ₁=0.025 mm^-1 and k ₂=0.057.

To further test the performance of Eq. (17) in extracting the optical properties of the medium under study, we measured diffuse reflectance on phantoms with varying optical properties. The phantom data were fitted to Eq. (17) and Table 1 summarizes typical results of this part of the study by showing the precision of Eq. (17) in predicting the actual optical property values of the phantoms.

Because there is more than one pair of the absorption and reduced scattering coefficient values corresponding to a given value of R_p , the results presented in Table 1 were found by fitting Eq. (17) to the experimentally measured value of R_p in a stepwise manner, first by keeping the absorption coefficient fixed and fitting to find the reduced scattering coefficient and then by keeping the reduced scattering coefficient fixed and fitting to find the absorption coefficient. Precision values shown in Table 1 actually represent the average precision of this two-step fitting process. The overall picture of the Table 1 results indicates that the optical properties can be recovered with a precision generally better than 6%. Note that these precision figures are generally better than those expected by inspecting Fig. 4. This is due to the increased flexibility afforded by Eq. (17), which includes two parameters for describing optical probe reflectance as opposed to a single parameter in Eq. (2) for the description of total reflectance.

 

Fig. 7. Measurements on tissue phantoms: diffuse reflectance as function of µ′ _s for zero absorption (top) and as a function of µ _α (bottom) for three different values of the reduced scattering coefficient: µ′ _s =3.2 mm^-1 (circles), µ′ _s =1.8 mm^-1 (triangles) and µ′ _s =0.65 mm^-1 (squares). Solid red lines represent the best fit to Eq. (17).

Download Full Size | PDF

Following testing on tissue phantoms, the model was applied to the analysis of diffuse reflectance spectra measured in vivo on human skin. Figure 8 shows typical fits to skin spectra, using Eqs. (17), (18) and (19), demonstrating the successful application of the model. The spectra exhibit absorption features characteristic of hemoglobin at approximately 540 nm and 580 nm and the melanocytic nevi spectrum, in particular, is characterized by smoothly decreasing reflectance with decreasing wavelength, which is typical of melanin absorption. In both spectra, the approximately linear dependence of diffuse reflectance on the wavelength in the 700–900 nm range is due to the linear dependence of the reduced scattering coefficient on the wavelength, expressed by Eq. (19). In this wavelength range, absorption by skin chromophores is minimal and scattering effects dominate. Table 2 summarizes the extracted parameter values measured on all volunteers. These figures are in good general agreement with the results of previous studies reporting on the optical properties of human skin [17, 44].

 

Fig. 8. Typical diffuse reflectance spectra of human skin measured in vivo. Black lines represent experimental data while red lines correspond to model fits.

Download Full Size | PDF

Table 1. Precision (%) for recovery of phantom optical properties using Eq. (17)

View Table | View all tables in this article

Table 2. Tissue parameters obtained from skin spectra analysis

View Table | View all tables in this article

5. Discussion and conclusions

We have shown that, based on insight gained from inspecting the physical picture of the diffuse reflectance process and the various available models, it is possible to derive a simple analytical model for describing diffuse reflectance from biological tissues. To our knowledge, this is the simplest model as compared to other models available in the literature for describing diffuse reflectance measured with an optical fiber probe. The model performs remarkably well in the analysis of spectra from tissue phantoms and in vivo human skin and its performance is comparable to that of other more complicated and rigorous models e.g. models based on diffusion theory or MC simulations.

We believe that the one-time phantom calibration required by the model does not constitute a serious disadvantage given that even more analytically-complicated and theoretically rigorous models usually do require at least one verification-calibration procedure on tissue phantoms while their utility in the biomedical optics field may be hampered by their complexity. In summary, the modeling approach described in this article presents the following major advantages:

(a) It describes diffuse reflectance collected by an optical probe through a simple expression with a minimal set of parameters and requires only a one-time calibration on a tissue phantom.

(b) It can be easily adapted to common fiber optic probe geometries through the initial calibration process. Once the probe response is determined for a given application, no further calibration or adjustment to the model is required.

An additional potential advantage of the model is the very fact that it successfully describes diffuse reflectance spectra with a small delivery-collection separation which is generally considered a weak point of other models based on diffusion theory.

In modeling diffuse reflectance, we have assumed tissue to be a semi-infinite, homogeneous medium. Even though this is clearly a first approximation and may even appear as a gross oversimplification for biological tissue, its applicability and utility has been tested and demonstrated by a number of studies [2, 17, 33]. In favor of the one-layer approach adopted in this article is also the fact that light does not penetrate very deep in tissue due to the small delivery-collection separation distance and thus does not reach tissue layers such as muscle and fat.

Nevertheless, a diffuse reflectance model that takes into account the multi-layer morphology of biological tissue may be more appropriate for accurate data analysis in several applications. The model presented here has the capacity to serve as the basis for the development of such a multi-layer model. In such a case, the advantages of its mathematical simplicity and the straightforward nature of the calibration procedure become very attractive and promising. We are currently exploring the generalization of the model to multi-layer tissue geometries.

One other very interesting potential application of the model, going beyond the use of an optical fiber probe, is diffuse reflectance spectral imaging. In such an application, real time processing of a very large number of spectra may be required and the simplicity of the present algorithm can significantly assist in reducing computing time.

Finally, the model is not limited to steady state diffuse reflectance and may be also useful for time domain or frequency domain applications of diffuse reflectance as well. It may also be adapted to model the spatial distribution of diffuse reflectance on the tissue surface. This array of capabilities makes the model even more attractive and opens the way for an extended range of applications.