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

Figure 1: A set of possible ray paths in a Maxwell's fish-eye.

Figure 2: Light pulse incident on a tissue specimen.

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

Figure 3: Few ray paths for a medium with a refractive index profile given by (14).

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: Intensity a  mm for different asymmetry factor () values.

Figure 5: Forward intensity at different locations for isotropic () scattering.

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.