Retrieving the optical parameters of biological tissues using diffuse reflectance spectroscopy and Fourier series expansions. I. theory and application.

The determination of optical parameters of biological tissues is essential for the application of optical techniques in the diagnosis and treatment of diseases. Diffuse Reflection Spectroscopy is a widely used technique to analyze the optical characteristics of biological tissues. In this paper we show that by using diffuse reflectance spectra and a new mathematical model we can retrieve the optical parameters by applying an adjustment of the data with nonlinear least squares. In our model we represent the spectra using a Fourier series expansion finding mathematical relations between the polynomial coefficients and the optical parameters. In this first paper we use spectra generated by the Monte Carlo Multilayered Technique to simulate the propagation of photons in turbid media. Using these spectra we determine the behavior of Fourier series coefficients when varying the optical parameters of the medium under study. With this procedure we find mathematical relations between Fourier series coefficients and optical parameters. Finally, the results show that our method can retrieve the optical parameters of biological tissues with accuracy that is adequate for medical applications.

1. Introduction

The determination of the optical parameters, which are defined by the absorption coefficient (μ), the scattering coefficient (μ), the refraction index (n) and the anisotropy factor (g), is of vital importance for the characterization of biological tissues using optical methods. In particular, for Diffuse Reflectance Spectroscopy is this determination is crucial. Many authors have presented different methods for retrieving optical parameters with different levels of difficulty and precision.

Farrell et al [1] present in their work a function for the reflectance with radial distribution through the approximation of the diffusion theory. Using a methodology of adjustment of the curve the least squares, this mathematical model allows the determination of optical parameters. In recent years models to predict the distribution of reflectance in turbid medium [2, 3, 4, 5] with different border conditions and levels of mathematical difficulty have been developed but our model allows us to perform the fitting of the diffuse reflectance curve regardless of the physical model, making use of the border conditions to calculate the retrieve of optical parameters. Other models extract optical properties (scattering and absorption coefficients) of the medium using small source-detector separations, for which the diffusion approximation is not valid [6, 7].

In this present work, we propose a new method to retrieve the optical parameters of biological tissues using the diffuse reflection with radial resolution data along with a Fourier series with adjustment by nonlinear least squares. Furthermore, using spectra generated by the Monte Carlo Multilayered Technique we can vary the values of the optical parameters and establish relationships between the Fourier series coefficients and the optical parameters which allow us to retrieve the optical parameters from real spectra.

2. Theoretical foundations

Fourier Series Expansion: the basic idea of these series is that any periodical function (T) can be expressed as a trigonometric sum of sines and cosines where their frequencies are multiples of the fundamental frequency ω0, but at the same intervals in which the function is defined. In our particular case, we do not have a periodical function. However, it could be considered a periodical function by parts, being defined as a range of period T equivalent to the radial distance to be studied. The idea of this proposal is to represent the reflectance as a linear combination of sines and cosines, called an order n trigonometric polynomial [8], Eq (1), where the expansion coefficients are determined through the adjustment of data. where r is the radial distance from the incident point to the point of exit of the light this distance is due to the multiple point of absorption and scattering within the sample before exiting. ω0 is the natural frequency defined for the period T (ω0 = 2π/T) and a, b are the Fourier coefficients.

The best adjustment R(r) with experimental data for to establish the relation of the Fourier series coefficients a and b with the optical parameters a(μ, μ, n, g) and b(μ, μ, n, g) must be found.

3. Methodology

To simulate diffuse reflection, the Monte Carlo method [9] was used. To do this, the important input parameters are: the absorption coefficient, the scattering coefficient, the anisotropy factor, the refraction index and the sample thickness. One million photons were used. Fig 1 shows the interaction between the incident beam and a sample of thickness d.

Fig. 1

Physical model of the simulation using the Monte Carlo method.

Taking the skin’s optical parameters as input values [1] we proceeded to use the Monte Carlo Multilayered (MCML) computing algorithm from the Oregon Medical Laser Center, Oregon Health and Sciences University [9], a program proven effective by other authors [10, 11], to obtain the diffused reflection with radial distribution in turbid media.

With the help of the Matlab Version 7.9.0.525 adjusting tools, the expansion of the Fourier series expansion was defined order 8. Using the results of the simulations obtained by Monte Carlo method as experimental data, we studied the behavior of Fourier series coefficients with the variation of the optical parameters and proceeded to propose an analytical relation (between the Fourier coefficients and the optical parameters) to facilitate the future retrieval of the optical parameters.

4. Results

Simulation of radial diffuse reflection: to evaluate the change of the diffuse reflection radial distribution using optical parameters, we first vary the scattering coefficient keeping the other parameters constant. Second, we vary the absorption coefficient without varying the other parameters. Third, we vary the refraction index, keeping the other optical parameters constant. All simulations are carried out with, 1,000,000 photons and with a sample thickness d = 200cm.

The first case simulated was by varying the scattering coefficient. We used as input parameters: the absorption coefficient (μ = 0.1mm−1), the anisotropy factor (g = 0.8), the refraction index (n = 1.4), varying the scattering coefficients over a range (μ = 10 – 200mm−1), In Fig 2 (A), as you can see, the diffuse reflectance intensity rises due to the increase of the scattering centers.

Fig. 2

Diffuse Reflectance when varying (A) the scattering coefficient (B) the absorption coefficient and (C) the refraction index.

When varying the absorption coefficient over a range (μ = 0.1 to 1mm−1) for this simulation, the input parameters are the scattering coefficient (μ = 75mm−1), the anisotropy factor (g = 0.8) and the refraction index (n = 1.4), Note that as the absorption coefficient increases, the reflection curve decreases progressively (see Fig 2 (B)). This was expected due to the increase of the absorption centers. Finally, the third case is for the refraction index. The variation was made in the interval 1.4 to 1.6 with the follows optical parameters: the scattering coefficient (μ = 75mm−1), the anisotropy factor (g=0. 8) and absorption coefficients (μ = 0.1mm−1). Note that as that parameter increases, the reflection curve decreases (see Fig 2 (C)).

Adjustment of the Reflection Curve with Trigonometric Functions: Once the simulations of the diffuse reflection resolution radial have been obtained, we do the trigonometric fit with least squares using Fourier series expansion. Fig 3 shows the radial reflection curve using the trigonometric polynomial of order 8, with an ω0 = 2.137(1/cm) and the root mean square of 0.9997.

Fig. 3

Trigonometric fitting of the reflection with radial resolution.

In Fig 4 we observe the radial reflectance curve fitted using Fourier Series (FS) and Farrell model [1] with Diffuse Approximation (DA). We observed that FS presents a good adjustment DA, but outside diffusion approximation, the adjustment by FS continued being effective; therefore, fitting the curve by diffuse reflectance Fourier series is more versatile and can be applied under different conditions.

Fig. 4

Curves adjustments radial reflectance (A)[σ = 0.9998;σ = 0.9999],(B)[σ = 0.9998;σ = 0.8029], where σ is standard deviation.

In Fig 5, we show the results for the variations of the Fourier series expansion coefficients a. Note that when increasing the scattering coefficient, the coefficients a1, a6 and a8 decrease but coefficients a2, a5 and a7 increase (see Fig 6). Then, two types of fit were made (linear and cubic) with the help of the program MATLAB, Version 7.9.0.525 computing tool CurveFitting.

Fig. 5

Variations of the Fourier coefficients a2, a5 and a7 related to the scattering.

Fig. 6

Variations of the Fourier coefficients a1, a6 and a8 related to the scattering.

In Table 1, parameters obtained in the adjustment are shown (see Fig 7). Then, for simplicity’s sake, we select a linear relation between both coefficients.

Table 1

Curve Fitting made for the scattering coefficient

Type	Equation	Adjustment Parameter	Standard Deviation
Linear	f(x) = p1x + p2	p1 = 9.344x109; p2 = − 8.199x1010	0.9721
Cubic	f(x) = p1x3 + p2x2 + p3x + p4	p1 = −2.436x105, p2 = 8.785x107p3 = 4.891x108, p4 = 1.264x1011	0.9855

Fig. 7

Variation of the Fourier coefficients a1 in relation to the scattering coefficient.

In addition, the behaviors of the Fourier series coefficients b were studied. The results are shows in Fig 8. Note that coefficients b2, b4 and b6, increase when the coefficient of scattering increases.

Fig. 8

Variations of the Fourier series coefficients b2, b4 and b6 in relation to the scattering coefficient.

It is important to note that coefficient a0 has been discarded from the study because its only contribution in the series is the displacement of the curve in the abscissa axis as seen in Fig 9.

Fig. 9

Increase of reflection when varying the scattering coefficient.

It is also noted that the first 5 coefficients a and b of the Fourier series expansion have been analyzed. The other terms were eliminated for their limited contribution in the series, though the adjustment made was of order 8 obtaining 16 Fourier coefficients (8 for a and 8 for b). Fig 10 shows the first five Fourier series coefficients. The nominal values are inappreciable in relation to the first terms.

Fig. 10

Relation of the nominal values of the Fourier series expansion coefficients.

In the case of the absorption coefficient in the Fig 11, the behavior of a1 is shown. Note an increase in a1 where the absorption coefficient increases. In Table 2, the adjustments made to the curve obtained in Fig 11 are shown in order to find the analytical relation between the coefficient a1 and the absorption coefficient, obtaning an analytical relation between both coefficients of an order 5 polynomial.

Fig. 11

Variation of the Fourier coefficient a1 in relation to the absorption coefficient.

Table 2

Curve Fitting made for the absorption coefficient

Type	Equation	Adjustment Parameter	Standard Deviation
Cubic	f(x) = p1x3 + p2x2 + p3x + p4	p1 = −1.533x1011, p2 = 2.749x1011p3 = −1.786x1011, p4 = 6.554x1011	0.9855
Order 5 polynomial	f(x) = p1x5 + p2x4 + p3x3 + p4x2 + p5x + p6	p1 = −1.44x1013, p2 = −4.032x1013p3 = 4.053x1013, p4 = −1.75x1013p5 = 2.88x1012, p6 = 5.38x1011	0.9925

In Fig 12, the behavior of the rest of the coefficients a of the Fourier series is shown, obtaining the same polynomial as the previous case in addition, a study of the variations of the coefficients b of the Fourier series expansion was carried out.

Fig. 12

Variation of the Fourier coefficient (A) a and (B) b with respect to the absorption coefficient.

In Figs. 13 (A) and 13 (B), the behavior of the coefficient a is shown where one can observe the curve made by the method of least squares for a order 2 polynomial in the case of the refraction index (see Table 3).

Fig. 13

Variation of the Fourier coefficient (A) a and (B) b with respect to respect of the refraction index.

Table 3

Curve Fitting made for the refraction index

an	Equation	Adjustment Parameter	Standard Deviation
Quadratic	f(x) = p1x2 + p2x + p1	p1 = −3.70x1012, p2 = 1.50x1012p3 = 1.69x1013	0.9934
Quadratic	f(x) = p1x2 + p2x + p1	p1 = 1.41x1012, p2 = −5.74x1012p3 = 6.44x1012	0.9933
Quadratic	f(x) = p1x2 + p2x + p1	p1 = 1.10x1011, p2 = −4.74x1011p3 = 5.379x1011	0.9948

Finally, no relationship was found between the Fourier series coefficients and the anisotropy factor. Table 4 shows the values of Fourier coefficients by varying the anisotropy factor; the slight changes in the values due to variations in the coefficients a and b can be attributed to the stochastic nature of Monte Carlo method.

Table 4

Adjustments of curves of coefficient a1 in function of the refraction index

Fourier coefficient	0.10	0.20	0.30	0.40	0.50
a1	−3.09x1012	−3.08x1012	−3.15x1012	−3.10x1012	−3.08x1012
b1	−1.88x1012	−1.87x1012	−1.92x1012	−1.88x1012	−1.87x1012
a2	1.18x1012	1.17x1012	1.20x1012	1.18x1012	1.17x1012
b2	2.2x10128	2.27x1012	2.32x1012	2.28x1012	2.27x1012
a3	9.97x1010	9.86x1010	1.02x1011	1.00x1011	9.94x1010
b3	−1.43x1012	−1.42x1012	−1.45x1012	−1.43x1012	−1.42x1012
a4	−3.55x1011	−3.53x1011	−3.61x1011	−3.56x1011	−3.53x1011
b4	5.01x1011	4.99x1011	5.10x1011	5.01x1011	4.98x1011
a5	1.80x1011	1.80x1011	1.83x1011	1.80x1011	1.79x1011
b5	−7.77x1010	−7.76x1010	−7.90x1010	−7.78x1010	−7.72x1010

5. Conclusion

Trigonometric adjustment is an effective method for the parametrization of the diffused reflection curves with radial distribution. A clear relationships between the Fourier series coefficients and the optical parameters is shown in the results. In the case of the scattering coefficient was obtained a linear behavior of the coefficients a and b of the Fourier series expansion. For the absorption coefficient there is a mathematical relationship of an order 5 polynomial and for the refraction index, the relationship is of an order 2 polynomial. Finally, no relation was found between expansion coefficients for the anisotropy factor, at least not the first ten expansions.