# An Approximate Numerical Technique for

Authors: Chintha C. Handapangoda; Malin Premaratne
Journal: Journal of Biomedicine and Biotechnology
 (2008)

---

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


---

## Body

## 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. Section2
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](1)n(r)c∂∂tI(r,Ω,t)+(1R1(s)+1R2(s))I(r,Ω,t)+n2(r)∂∂s(I(r,Ω,t)n2(r))+σt(r)I(r,Ω,t)=σs(r)∫4πP(Ω,Ω′)I(r,Ω′,t)dΩ′+F(r,Ω,t),where r is the position
vector of a point on a path of a ray, s is the arc
length along a ray, Ω=dr/ds, t is the time
variable, I(r,Ω,t) is the
intensity, n(r) is the
refractive index profile, c is the speed of
light in vacuum, R1(s) and R2(s) are the
principal radii of curvatures of the geometrical wave-fronts, σt is the
attenuation coefficient with σt=σa+σs, σa is the
absorption coefficient and σs is the
scattering coefficient, P(Ω,Ω′) is the phase
function, and F(r,Ω,t) 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]:(2)dds(n(r)drds)=∇n(r), where r=xi+yj+zk, i, j, and k are unit
vectors along x, y, and z axes,
respectively, and ∇=((∂/∂x)i+(∂/∂y)j+(∂/∂z)k). We use the ray model of light here, based on
geometric optic techniques, which neglects the concept of wavelength λ (i.e., λ→0) [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 q is introduced
such that dq=ds/n. Then, (2) can be written as(3)d2rdq2=n(r)∇n(r).The optical ray vector is defined as(4)Q=drdq=nξi+nηj+nμk, where ξ=sinθcosϕ, η=sinθ⁢sinϕ, and μ=cos⁡θ. Equation (3) can be written in matrix format as(5)d2Rdq2=D(R), where (6)R=(xyz),Q=(nξnηnμ),D=n(∂n∂x∂n∂y∂n∂z).Thus, (5) can be solved using the Runge-Kutta
algorithm starting from a known point (R0,Q0). That is, (5) with the initial condition (R0,Q0) will
successively generate points (R1,Q1),(R2,Q2),…,(Rn,Qn) along the path
[16]. Therefore, if we select a set of starting points and initial directions, Q0, 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:(7)τ=t−sv, where v is the speed of
light inside the medium along the ray path. Using (7) in (1) results in(8)n2(r)∂∂s(I(r,Ω,τ)n2(r))+(1R1(s)+1R2(s))I(r,Ω,τ)+σt(r)I(r,Ω,τ)=σs(r)∫4πP(Ω,Ω')I(r,Ω',τ)dΩ'+F(r,Ω,τ). In this paper,
we consider plane waves so that the second term in the left-hand side of (8)
vanishes. That is,(9)n2(r)∂∂s(I(r,Ω,τ)n2(r))+σt(r)I(r,Ω,τ)=σs(r)∫4πP(Ω,Ω')I(r,Ω',τ)dΩ'+F(r,Ω,τ). 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 n(r), σt(r),
and σs(r) are constants
at that point. First, we use Gaussian quadrature [23] to approximate the
integral:(10)n2(r)∂∂s(I(r,Ω,τ)n2(r))+σt(r)I(r,Ω,τ)=σs(r)∑j=1qwjP(Ω,Ωj)I(r,Ωj,τ)+F(r,Ω,τ), where Ωj is the jth quadrature
point and wj is the
corresponding Gaussian weight [23]. Then, the time dependence is represented
using a Laguerre expansion [24]:(11)I(r,Ω,τ)=∑k=0NBk(r,Ω)Lk(τ), where(12)∫0∞Ln(τ)Lm(τ)e−τdτ=δmn. Using (11) in
(10) and taking moments, we get(13)n2(r)∂∂s(Bn(r,Ω)n2(r))+σt(r)Bn(r,Ω)=σs(r)∑j=1qwjP(Ω,Ωj)Bn(r,Ωj)+Fn(r,Ω), where Fn(r,Ω) is the Laguerre
coefficient of the source term F(r,Ω,τ). Since we have previously traced the ray paths and
chosen a set of points r and Ω values are
known, thus, (13) can be solved using the Runge-Kutta-Fehlberg algorithm [25].


## 3. Results and Discussion

Figure3 shows
few ray paths for a medium with a refractive index profile given byFigure 3
Few ray paths for a medium with a refractive index
profile given by (14).(14)n=n0e−x2, where n0=2. 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 g is the
asymmetry factor.Figure 4
Intensity
az=1 mm for different asymmetry factor (g) values.Figure 5
Forward
intensity at different locations for isotropic (g=0) scattering.Figure4 shows the variation of intensity with time at z=1 mm with varying g. The graphs correspond to g=0.8, g=0.7, g=0.6, and g=0. Other parameters such as the scattering coefficient
and the absorption coefficient were kept constant for all the three graphs. The
condition g=0 corresponds to
the isotropic scattering case while g=0.8 represents
forward scattering. This is illustrated by the above four graphs.Figure5 shows the variation of the forward intensity
at different locations, that is, corresponding to different z values (in mm),
with a constant asymmetry factor g=0. 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.



---

*Source: 784354-2008-02-06.xml*