Journal of Biomedicine and Biotechnology

Volume 2008 (2008), Article ID 784354, 5 pages

http://dx.doi.org/10.1155/2008/784354

## An Approximate Numerical Technique for Characterizing Optical Pulse Propagation in Inhomogeneous Biological Tissue

Advanced Computing and Simulation Laboratory (AXL), Department of Electrical and Computer Systems Engineering, Monash University, Clayton, VIC 3800, Australia

Received 8 September 2007; Accepted 19 December 2007

Academic Editor: Daniel Howard

Copyright © 2008 Chintha C. Handapangoda and Malin Premaratne. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

An approximate numerical technique for modeling optical pulse propagation through weakly scattering biological tissue is developed by solving the photon transport equation in biological tissue that includes varying refractive index and varying scattering/absorption coefficients. The proposed technique involves first tracing the ray paths defined by the refractive index profile of the medium by solving the eikonal equation using a Runge-Kutta integration algorithm. The photon transport equation is solved only along these ray paths, minimizing the overall computational burden of the resulting algorithm. The main advantage of the current algorithm is that it enables to discretise the pulse propagation space adaptively by taking optical depth into account. Therefore, computational efficiency can be increased without compromising the accuracy of the algorithm.

#### 1. Introduction

Modeling of light propagation through biological tissues is important for many medical applications such as optical tomography for cancer detection [1] and noninvasive detection of diabetes mellitus [2]. Researchers have been working on modeling biological tissues over the last two decades [3, 4].

Light propagation through biological tissues can be modeled using the photon transport equation (PTE) [5, 6]. Several numerical models have been developed to solve the PTE over the recent years [7–10]. These models include techniques for solving the steady-state PTE [7], as well as the transient PTE for short pulse propagation [8–10]. Several different variations of PTE for inhomogeneous media have been proposed in the literature [6, 11–15]. However, most of these variations are not a result of fundamental errors or differences but due to different assumptions about the medium or wave field properties. For example, [6, 11, 13] look at spatially slowly varying isotropic refractive index profiles in their work. Interestingly, the approach given in [15] is formulated to accommodate geometric optics approximations but ignore the wavefront curvature in their derivation. Wavefront curvature in the context of slowly varying refractive index approximation is considered in [6]. Numerical considerations necessary to account for such slowly varying spatial refractive profiles are considered in [12, 14].

This paper presents a technique for modeling the light propagation through weakly scattering and absorbing media by solving the three-dimensional transient PTE numerically. In a weakly scattering medium, the ray paths can be approximated by the Eikonal equation [16]. First, a set of ray paths is calculated depending on the refractive index profile of the medium. There are several existing methods to do this [16–18]. In this paper, we have briefly described a ray tracing technique proposed by Sharma et al. [16]. The next step is to solve the RTE written in terms of the arc length [6] on each of these paths. The Laguerre Runge-Kutta-Fehlberg method [19] can be used for this purpose.

This paper is organized in four sections. Section 2 presents the formulation of the numerical technique. Section 3 contains some results obtained using this technique and a discussion on these results. Section 4 includes the conclusions.

#### 2. Formulation

The photon transport equation for a medium with a spatially varying isotropic refractive index, in standard spherical coordinates, is [6]

where is the position vector of a point on a path of a ray, is the arc length along a ray, , is the time variable, is the intensity, is the refractive index profile, is the speed of light in vacuum, and are the principal radii of curvatures of the geometrical wave-fronts, is the attenuation coefficient with , is the absorption coefficient and is the scattering coefficient, is the phase function, and represents sources inside the medium.

The path of rays in a medium with a spatially varying refractive index is given by the Eikonal equation [16–18]:

where , , , and are unit vectors along , , and axes, respectively, and . We use the ray model of light here, based on geometric optic techniques, which neglects the concept of wavelength (i.e., ) [20]. Thus, this approximation will be valid for modeling light propagation through inhomogeneous media in which the properties vary very slowly compared to the wavelength of the incident light. Also, the geometric optic techniques assume that the field behaves locally as a plane wave and that the intensity does not change rapidly [20].

The technique proposed in this paper is to first solve (2) to obtain a set of possible ray paths and then solve the PTE, (1), for each of these paths. The main advantage of the current algorithm is that it enables to discretise the pulse propagation space adaptively by taking optical depth into account. Therefore, computational efficiency can be increased without compromising the accuracy of the algorithm [21]. Several techniques had been introduced for solving the Eikonal equation by various research groups [16–18]. We adopt the technique proposed by Sharma et al. [16] because it uses Runge-Kutta integration, which will be used again later to solve photon transport equation in discrete ordinate method setting. In this method, a new variable is introduced such that . Then, (2) can be written as

The optical ray vector is defined as

where , , and . Equation (3) can be written in matrix format as

where Thus, (5) can be solved using the Runge-Kutta algorithm starting from a known point . That is, (5) with the initial condition will successively generate points along the path [16]. Therefore, if we select a set of starting points and initial directions, , we can obtain a set of ray paths by numerically integrating (5). For example, such ray tracing for Maxwell's fish-eye gives the paths as shown in Figure 1 [22].

Next step is to solve the PTE, (1), for a weakly scattering medium on each of the above sets of paths, by numerically integrating (5). While tracing the ray paths from the above algorithm, the PTE can be solved to find the intensity at each point on the path.

First, we use the following transformation which maps the PTE to a moving coordinate system on ray paths:

where is the speed of light inside the medium along the ray path. Using (7) in (1) results in

In this paper, we consider plane waves so that the second term in the left-hand side of (8) vanishes. That is,

The Laguerre Runge-Kutta-Fehlberg (LRKF)method [19] can be used to solve (9) for the intensity at selected points on each ray path. The LRKF method is briefly described below. Since we solve (9) on a known ray path at a known point , , and are constants at that point. First, we use Gaussian quadrature [23] to approximate the integral:

where is the th quadrature point and is the corresponding Gaussian weight [23]. Then, the time dependence is represented using a Laguerre expansion [24]:

where

Using (11) in (10) and taking moments, we get

where is the Laguerre coefficient of the source term . Since we have previously traced the ray paths and chosen a set of points and values are known, thus, (13) can be solved using the Runge-Kutta-Fehlberg algorithm [25].

#### 3. Results and Discussion

Figure 3 shows few ray paths for a medium with a refractive index profile given by

where . Figures 4 and 5 were obtained from the above algorithm. We have obtained these results for a pulse propagating through a single ray path. The Henyey-Greenstein phase function [27] was used for the simulation where is the asymmetry factor.

Figure 4 shows the variation of intensity with time at mm with varying . The graphs correspond to , , , and . Other parameters such as the scattering coefficient and the absorption coefficient were kept constant for all the three graphs. The condition corresponds to the isotropic scattering case while represents forward scattering. This is illustrated by the above four graphs.

Figure 5 shows the variation of the forward intensity at different locations, that is, corresponding to different values (in mm), with a constant asymmetry factor . It can be clearly seen from this figure that the intensity reduces with increasing distance due to scattering and absorption. Also, the pulse is shifted in time as shown.

#### 4. Conclusion

This paper introduced a novel approximate numerical model to solve the photon transport equation with inhomogeneous refractive index and inhomogeneous scattering and absorption cross-sections for weakly scattering biological tissue. The proposed technique consists of two steps: first, the Eikonal equations describing the geometric-optic ray paths are solved using an efficient Runge-Kutta routine to construct a set of possible photon migration paths through the inhomogeneous tissue sample. Thereafter, the transient photon transport equation, written in terms of the arc length, is solved along each ray path to construct the optical intensity variation as time progresses. The main advantage of the current algorithm is that it enables to discretise the pulse propagation space adaptively by taking optical depth into account. Therefore, computational efficiency can be increased without compromising the accuracy of the algorithm. Computational efficiency becomes a bottle-neck when large, realistic simulations of optical pulse propagation in tissue are required. Therefore, current proposed method is very much suited for extensive computations work in time-resolved photon-diffusion tomography work.

#### References

- S. B. Colak, M. B. van der Mark, G. W. 'T Hooft, J. H. Hoogenraad, E. S. van der Linden, and F. A. Kuijpers, “Clinical optical tomography and NIR spectroscopy for breast cancer detection,”
*IEEE Journal on Selected Topics in Quantum Electronics*, vol. 5, no. 4, pp. 1143–1158, 1999. View at Publisher · View at Google Scholar - B. D. Cameron, H. W. Gorde, B. Satheesan, and G. L. Cote, “The use of polarized laser light through the eye for noninvasive glucose monitoring,”
*Diabetes Technology & Therapeutics*, vol. 1, no. 2, pp. 135–143, 1999. View at Publisher · View at Google Scholar - F. Liu, K. M. Yoo, and R. R. Alfano, “Ultrafast laser-pulse transmission and imaging through biological tissues,”
*Applied Optics*, vol. 32, no. 4, pp. 554–558, 1993. View at Google Scholar - G. Yao and L. V. Wang, “Theoretical and experimental studies of ultrasound-modulated optical tomography in biological tissue,”
*Applied Optics*, vol. 39, no. 4, pp. 659–664, 2000. View at Publisher · View at Google Scholar - A. D. Klose, U. Netz, J. Beuthan, and A. H. Hielscher, “Optical tomography using the time-dependent equation of radiative transfer—part 1: forward model,”
*Journal of Quantitative Spectroscopy and Radiative Transfer*, vol. 72, no. 5, pp. 691–713, 2002. View at Publisher · View at Google Scholar - M. Premaratne, E. Premaratne, and A. J. Lowery, “The photon transport equation for turbid biological media with spatially varying isotropic refractive index,”
*Optics Express*, vol. 13, no. 2, pp. 389–399, 2005. View at Publisher · View at Google Scholar - K. Stamnes, S. Tsay, W. Wiscombe, and K. Jayaweera, “Numerically stable algorithm for discrete-ordinate-method for radiative transfer in multiple scattering and emitting layered media,”
*Applied Optics*, vol. 27, no. 12, pp. 2502–2509, 1988. View at Google Scholar - A. D. Kim and M. Moscoso, “Chebyshev spectral methods for radiative transfer,”
*SIAM Journal on Scientific Computing*, vol. 23, no. 6, pp. 2074–2094, 2002. View at Publisher · View at Google Scholar - J. V. P. de Oliveira, A. V. Cardona, and M. T. M. B. de Vilhena, “Solution of the one-dimensional time-dependent discrete ordinates problem in a slab by
the spectral and ${\text{LTS}}_{\text{N}}$ methods,”
*Annals of Nuclear Energy*, vol. 29, no. 1, pp. 13–20, 2002. View at Publisher · View at Google Scholar - M. Sakami, K. Mitra, and P.-F. Hsu, “Analysis of light pulse transport through two-dimensional scattering and absorbing media,”
*Journal of Quantitative Spectroscopy and Radiative Transfer*, vol. 73, no. 2–5, pp. 169–179, 2002. View at Publisher · View at Google Scholar - T. Khan and H. Jiang, “A new diffusion approximation to the radiative transfer equation for scattering media
with spatially varying refractive indices,”
*Journal of Optics A: Pure and Applied Optics*, vol. 5, no. 2, pp. 137–141, 2003. View at Publisher · View at Google Scholar - H. Dehghani, B. Brooksby, K. Vishwanath, B. W. Pogue, and K. D. Paulsen, “The effects of internal refractive index variation in near-infrared optical tomography:
a finite element modelling approach,”
*Physics in Medicine and Biology*, vol. 48, no. 16, pp. 2713–2727, 2003. View at Publisher · View at Google Scholar - G. Bal, “Radiative transfer equations with varying refractive index: a mathematical perspective,”
*Journal of the Optical Society of America A*, vol. 23, no. 7, pp. 1639–1644, 2006. View at Publisher · View at Google Scholar - M. L. Shendeleva and J. A. Molly, “Scaling property of the diffusion equation for light in a turbid medium with varying refractive index,”
*Journal of the Optical Society of America A*, vol. 24, no. 9, pp. 2902–2910, 2007. View at Publisher · View at Google Scholar - L. Martí-López, J. Bouza-Domínguez, R. A. Martínez-Celorio, and J. C. Hebden, “An investigation of the ability of modified radiative transfer equations to accommodate laws of geometrical optics,”
*Optics Communications*, vol. 266, no. 1, pp. 44–49, 2006. View at Publisher · View at Google Scholar - A. Sharma, D. V. Kumar, and A. K. Ghatak, “Tracing rays through graded-index media: a new method,”
*Applied Optics*, vol. 21, no. 6, pp. 984–987, 1982. View at Google Scholar - B. Richerzhagen, “Finite element ray tracing: a new method for ray tracing in gradient-index media,”
*Applied Optics*, vol. 35, no. 31, pp. 6186–6189, 1996. View at Google Scholar - H. A. Ferwerda, “Radiative transfer equation for scattering media with a spatially varying refractive index,”
*Journal of Optics A: Pure and Applied Optics*, vol. 1, no. 3, pp. L1–L2, 1999. View at Publisher · View at Google Scholar - C. C. Handapangoda, M. Premaratne, L. Yeo, and J. Friend, “Laguerre Runge-Kutta-Fehlberg method for simulating laser pulse propagation in biological tissue,”
*IEEE Journal of Selected Topics in Quantum Electronics*, vol. 14, no. 1, 2008. View at Publisher · View at Google Scholar - M. Born and E. Wolf,
*Principles of Optics*, Cambridge University Press, Cambridge, UK, 7th edition, 1999. - M. Premaratne, “Numerical simulation of nonuniformly time-sampled pulse propagation in nonlinear fiber,”
*Journal of Lightwave Technology*, vol. 23, no. 8, pp. 2434–2442, 2005. View at Publisher · View at Google Scholar - A. D. Greenwood and J.-M. Jin, “A field picture of wave propagation in inhomogeneous dielectric lenses,”
*IEEE Antennas and Propagation Magazine*, vol. 41, no. 5, pp. 9–18, 1999. View at Publisher · View at Google Scholar - W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, “Integration of functions,” in
*Numerical Recipes in C++*, pp. 152–157, Cambridge University Press, Cambridge, UK, 2003. View at Google Scholar - M. Abramomitz and I. A. Stegun, “Orthogonal polynomials,” in
*Handbook of Mathematical Functions with Formulas, Graphs and Mathematical Tables*, pp. 773–784, Dover, New York, NY, USA, 1965. View at Google Scholar - S. C. Chapra and R. P. Canale, “Runge-Kutta methods,” in
*Numerical Methods for Engineers*, pp. 719–720, McGraw-Hill, New York, NY, USA, 4th edition, 2002. View at Google Scholar - M. P. Mengüç and R. Viskanta, “Radiative transfer in three-dimensional rectangular enclosures containing inhomogeneous,
anisotropically scattering media,”
*Journal of Quantitative Spectroscopy and Radiative Transfer*, vol. 33, no. 6, pp. 533–549, 1985. View at Publisher · View at Google Scholar - G. E. Thomas and K. Stamnes,
*Radiative Transfer in the Atmosphere and Ocean*, Cambridge University Press, Cambridge, UK, 1999.