## Abstract

We present the first algorithm for solving the equation of radiative transfer (ERT) in the frequency domain (FD) on three-dimensional block-structured Cartesian grids (BSG). This algorithm allows for accurate modeling of light propagation in media of arbitrary shape with air-tissue refractive index mismatch at the boundary at increased speed compared to currently available structured grid algorithms. To accurately model arbitrarily shaped geometries the algorithm generates BSGs that are finely discretized only near physical boundaries and therefore less dense than fine grids. We discretize the FD-ERT using a combination of the upwind-step method and the discrete ordinates (*S _{N}*) approximation. The source iteration technique is used to obtain the solution. We implement a first order interpolation scheme when traversing between coarse and fine grid regions. Effects of geometry and optical parameters on algorithm performance are evaluated using numerical phantoms (circular, cylindrical, and arbitrary shape) and varying the absorption and scattering coefficients, modulation frequency, and refractive index. The solution on a 3-level BSG is obtained up to 4.2 times faster than the solution on a single fine grid, with minimal increase in numerical error (less than 5%).

© 2010 OSA

## 1. Introduction

There has long been an interest in modeling how light propagates in biological tissues [1]. The emergence of novel diagnostic and therapeutic methods that rely on lasers in the 1980’s and 1990’s has led to a need for accurately modeling light-tissue interaction. Various models have been developed that make use of ever-increasing computational power and new numerical methods. Fundamentally two different approaches have been pursued: (a) Monte Carlo (MC) modeling [2,3] and (b) numerical solutions to the equation of radiative transfer (ERT) [4–11] or its approximation, the diffusion equation (DE) [4,12–15]. MC methods employ statistical techniques to propagate a large number of photons through tissue. While highly accurate results can be achieved, MC methods are computationally very intensive. Finite-differences, finite-element, and finite–volume methods in which the ERT or DE is discretized and the resulting system of equations is solved numerically are, in general, less computationally expensive. However, the complexity of these codes can still lead to various computational challenges. Furthermore, it has been established that the ERT is more accurate and preferred over the DE when modeling light propagation through highly absorbing media and media with void like regions [7].

A major consideration when implementing non-MC methods is the form of the grid used in these calculations. The ERT is typically solved on either a structured or an unstructured grid. Structured and unstructured grids differ from one another by the method in which Euclidean space is discretized. Structured grids discretize Euclidean space into a set of hexahedron elements (quadrilaterals in two-dimensions), while typical unstructured grids discretize Euclidean space into a collection of tetrahedral elements (or triangles, quadrilaterals, etc., in two-dimensions). The connectivity between nodes in a structured grid is implicitly known; nodes always connect to form a cuboid of constant size and orientation. This *a priori* information is central to the finite differences method, where any given grid point is assumed to have neighboring grid points that also connect to form cuboids. Conversely, the connectivity of unstructured grid elements must be explicitly provided because unstructured grid nodes connect to form elements that vary in size, shape, and orientation [16]. Thus, algorithms that solve the ERT on unstructured grids are more complex than algorithms on structured grids because node connectivity information must be explicitly provided and processed. Generating unstructured grids can in its-self be an arduous task, often requiring third party applications, while generating structured grids is a relatively simple task [16]. Numerical algorithms on both unstructured and structured grids have been developed for solving the ERT. For example, the unstructured finite-element method has been used by Arridge et al. [13], Salah et al. [11], and Rasmussen et al. [17]. Unstructured finite-volume methods have been employed by Kim and Hielscher [9], Ren et al. [8], and Gu et al. [18]. Finite-differences ERT codes on structured grids have been mainly pursued by Klose et al [4–6,19–21]. In general, solving the ERT on structured grids is attractive because algorithm complexity is minimized as a result of the well-ordered nature of Cartesian grids.

Unstructured grids are generally seen to be superior to structured grids because they appear to be better at modeling arbitrarily shaped geometries. Structured grids, however, can also model geometries of arbitrary shape. Structured grids can resolve curved boundaries using the blocked-off region method and it has already been applied to light propagation problem using the simplified spherical harmonics approximation to the ERT by Klose et al. [21]. Modeling curved boundaries with structured grids often requires finely discretized grids (i.e. dense grids). Unfortunately, the computational cost of solving the ERT on dense grids can be prohibitively expensive. Performing calculations on a coarse grid can save computational time, but it can also lead to significant numerical error. The use of coarse grids can also result in increased modeling error. This occurs because coarse grids are often unable to accurately resolve the physical boundary, leading to inaccurate contributions from the boundary conditions. Additionally, grids that are too coarse at the boundary might not be able to accurately locate boundary sources and detectors. Fiber bundles are typically used to direct photons from a laser source to the medium. Similar sets of fibers are used to direct the escaping photons to a detector [22,23]. Thus, these sources and detectors are highly localized and a finely discretized grid is necessary to accurately describe their position on the surface. Furthermore, in the area of small animal fluorescence and bioluminescence imaging it is possible to have internal sources close to the tissue boundary and it is well known that these types of sources can lead to severe modeling errors 20]. Using a dense grid that can accurately capture boundary effects can reduce this source of error; however, this solution also leads to an increase in computation time.

The aim of this work is to reduce the computational effort required to solve the ERT on structured grids. Implementing an algorithm that solves the ERT on a grid that is refined only near boundaries can reduce computational cost. The resulting grid is a single grid with various levels of refinement - the grid is coarsely discretized in the interior and finely discretized near the boundary. These types of grids are called block-structured grids (BSGs) and are a subset of the more general adaptive mesh refinement techniques [16]. A single dense Cartesian grid is transformed into a relatively sparse grid. Solving the ERT on a BSG requires less computational effort than solving it on a single dense grid, as there are less grid nodes inside the computational domain. The area of Cartesian grids with local refinements has itself been a topic of research [16,24]. In particular, BSGs have been used in the field of computational fluid dynamics, especially in the study of shock-hydrodynamics [25] and general fluid flows [26]. Previous work has mainly focused on implementing finite volume methods on BSGs [27] that were adaptively refined on the interior of the medium. In our work we focus on refining the grid only near the tissue boundary,

In connection with the ERT only one group has used BSGs to date. In 1998, Jesse et al. [28] presented a method to solve the time-independent ERT on BSGs using a finite-volume method. Their code is limited to two-dimensional rectangular media with isotropic scattering and does not consider partially reflective boundary conditions. While useful for applications in nuclear physics and heat transfer, this algorithm is not suited for problems concerning light propagation in arbitrarily shaped media with highly anisotropic scattering and partially reflective boundary conditions. In the area of tissue optics, BSGs have been employed for solving the diffusion equation for fluorescent light propagation [29]. Their work is limited to treating rectangular geometries and has been developed as a finite-element method algorithm.

In this paper we go beyond the approaches previously presented and implement the first frequency-domain (FD) ERT with reflective boundary conditions on BSGs of arbitrary shape. The FD-ERT is particularly important in small animal imaging, where the diffusion approximation has limited validity, as light travels only a few mean-free path and non-diffusive boundary conditions dominate the solutions [7,18,30]. It is highly desirable to have a fast and accurate numerical solver for this and other similar applications (e.g. imaging of arthritic human finger joints [31]) that require solutions to the problem of light propagation in small domains.

The remainder of this paper is organized as follows. In Section 2 we present our BSG generating algorithm, review the multiple forms of the frequency domain ERT, discuss its finite difference discretization, and present our modification to these established methods for treating BSGs. We also present numerical phantoms used for simulations in this work. In section 3 we present our results, and we conclude with a discussion in Section 4.

## 2. Methods

The overall structure of our algorithms is shown in Fig. 1 . First, the boundary information is provided by giving the coordinates and normal vectors of the surface that encloses the medium under consideration. The code then determines the computational domain by filling the volume defined in the first step with structured grid points separated by user-defined spacings. This grid is transformed into a BSG using our grid generator, and the code then solves the light propagation problem using the discretized FD-ERT on this BSG. In the following sections we will describe each aspect of our code in more detail.

#### 2.1. Grid Generation

To generate the appropriate BSG we must first determine the active computational domain using information about the physical boundary provided as input (surface coordinates and normal vectors). This process, known as the blocked-off region (BOR) method, starts by defining a nominal domain in the form of a rectangular cuboid that completely encloses the arbitrarily shaped object defined by the boundary information [21]. The domain is then discretized on a Cartesian grid with a user specified spatial resolution along each axis (Δ*x*, Δ*y*, and Δ*z*). It is necessary that the resolution be small enough to capture all physical effects. In our work a grid spacing of 1/10 μ_{s} is generally sufficient. Next, the code finds all grid points within the discretized domain that lie within the volume enclosed by the boundary.

After discretizing Euclidean space with a fine grid, the BOR method defines objects of arbitrary shape by segmenting the grid into an “active” and an “inactive” region. The method proceeds as follows. First, grid points inside the nominal domain corresponding to the user provided surface coordinates are identified (called “boundary points”). These points are assigned to the active region. Next, the algorithm checks one grid point at a time and determines if it is a boundary point. When a boundary point is encountered, the algorithm uses the normal vector of the surface at that grid point to determine in which direction the “inside” and “outside” regions of the medium are, relative to the current grid point. All future non-boundary points in the “inside” region are labeled as “interior points” and are assigned to be part of the active domain (dark gray area, Fig. 2b
). Points in the “outside” region of the medium are labeled “exterior points,” are inactive, and are not used in future calculations (light gray area, Fig. 2b). We illustrate this process with the following example. Let the normal vector at boundary point (*i*,*j*,*k*) be ** n** =

*n*

_{x}**ê**+

_{1}*n*

_{y}**+**

*ê*_{2}*n*

_{z}**with**

*ê*_{3}*n*<0. Then, nodes (

_{x}*i*-1,

*j*,

*k*) and (

*i*+ 1,

*j*,

*k*) are classified as “exterior” and “interior” points, respectively. Furthermore, all grid points to the right of (

*i*+ 1,

*j*,

*k*) are also “interior” points (i.e. (

*i*+ 2,

*j*,

*k*), (

*i*+ 3,

*j*,

*k*), etc.). This labeling scheme continues until a new boundary point is encountered, at which point the process is restarted. The result of this subroutine is a Cartesian grid where all nodes are classified as exterior, boundary, or interior points. The resulting grid is used as input to the BSG generating subroutine. An overview of the major aspects of the BOR method is presented in Fig. 3 .

The BSG generator adjusts the grid inside the computationally active domain. As input for this step, the level of refinement needs to be specified by the user. The grid generating subroutine discretizes the computational domain with the finest expected spacing (Δ*x _{f}*) as provided by the user. A subroutine then superimposes the next coarsest grid (Δ

*x*= 2Δ

*x*) over the fine grid. Next, the algorithm removes coarse grid points on the boundary of the computational domain. In the final step, all fine grid points within the coarse grid are removed from the computational domain. This process is repeated to generate higher-order BSGs. The output of the BSG generator contains the following information for each grid point: 1) location: exterior, main boundary, or interior, 2) grid level, and 3) type of fine/coarse boundary. The term “fine/coarse boundary” refers to points in the “active” domain at the union of coarse and fine grid segments (Fig. 4a ). The C + + pseudo code for the algorithm generating a two-dimensional

_{f}*n-level*grid is as follows (finest grid = level 0, coarsest grid = level

*n*):

- 1. Generate grid hierarchy: (for
*m = 0*;*m < = n*;*m + +*)- a. Determine all level
*m*grid points in exterior, main boundary, or interior. - b. If
*m > 0*, for all interior grid points level*m*i. Delete all points level

*m*on the main numerical boundary.ii. Delete all points level

*m-1*interior to a cell of 4 level*m*points.iii. Delete all points level

*m-1*interior to a cell of 6 level*m*points.

- 2. Classify active grid points into four categories
- a. On main boundary
- b. On fine/coarse boundary (Fig. 4a)
i. Missing west neighbor

ii. Missing south neighbor

iii. Missing east neighbor

iv. Missing north neighbor

v. Coarse point (black dot)

c. Interior point, in active domain, and not on fine/coarse boundary.

d. Interior point not in active domain (black diamond, Fig. 4b).

The algorithm for generating three-dimensional BSGs is conceptually the same, however, the number of possible configurations of fine/coarse boundary points increases from 5 to 19 in step 2b of the pseudo code for the algorithm generating two-dimensional BSGs. To summarize, consider the example presented in Fig. 2. A cross-section through a mouse obtained from a magnetic resonance imaging (MRI) scan is shown in Fig. 2a. Discretization of this image by a dense grid, a 3-level BSG, and a coarse grid are shown in Fig. 2b-d, respectively. The fine grid and BSG approximate the physical boundary with the same accuracy; however, the fine grid is clearly denser. The coarse grid is less dense than the fine grid but it is a poor approximation of the true geometry. After comparing all three grids it is clear that BSGs can resolve physical boundaries as accurately as dense grids using but with significantly fewer points. These properties make BSGs attractive for solving the ERT at reduced computation cost with minimal loss in accuracy.

#### 2.2 Light Propagation Model

As a particular implementation of a light transport model on BSGs, this work focuses on solutions of the frequency-domain ERT. In this case, a boundary light source is intensity modulated, leading to propagation of so-called photon density waves in tissue. Photon density waves have been used in many applications to characterize optical properties of various scattering media including biomedical tissues [1,32]. We have previously developed a general mathematical framework for modeling these photon density waves [4].

To review, the general form of the FD-ERT and related boundary conditions are given by

A list of the associated variables and their meanings is given in Table 1 . While the radiance is an important quantity, it is the photon fluence rate given by

that is often the most important quantity in clinical medicine, as it describes the amount of energy per unit area at a given location (watt/area or photons s^{−1} cm^{−2}) [1,4]. A second important quantity in clinical applications is the partial current

This is the quantity that can be directly measured at a given boundary point [4,22,23,33].

The general form of the FD-ERT, Eqs. (1),2), can be adapted to account for fluorescence and bioluminescence effects which have become increasingly important in recent years [30,34–37]. In the case of fluorescence, one equation (*ERT 1*) is used to model the excitation field inside the medium due to a modulated boundary source. The fluorophore inside the medium is modeled as an internal source (*Q*). Fluorophore emission is a function, among other things, of the excitation field. Thus, a second ERT is used to model the fluorophore emission field (*ERT 2*). *ERT 1* is obtained from Eqs. (1),2) by defining the boundary source as

setting the internal source (*Q*) to zero, and defining the total attenuation coefficient as

*ERT 2* is obtained by setting the boundary source (*S*) to zero and defining the fluorescent source as a function of the excitation field. This relationship is given by

During emission, the total attenuation coefficient is only a function of the absorption and scattering parameters of the medium at the emission wavelength, and is given by

Bioluminescence can be modeled by a single ERT. Bioluminescence cannot be modulated; all information with regards to modulation frequency must be discarded. In this case, only internal sources exist and the governing equation is obtained from Eqs. (1),2) by setting boundary sources (*S*) to zero and defining the source of bioluminescence (*Q*). The total attenuation coefficient is given by Eq. (8).

#### 2.3 Discretization of Frequency Domain ERT

To discretize the FD-ERT, we implement a finite-differences upwind step-method for the spatial variable and a discrete-ordinates method (*S _{N}*) for the angular variable [4,20,38–41]. Here, we will first review how this is implemented on a single Cartesian grid. Subsequently, we will describe how this method is adapted to BSGs.

*2.3.1. Discretization on single grid:* The first step in discretizing the ERT is to use the discrete-ordinates method to replace the integral term with the *extended trapezoidal rule* [5]. This approximation is given by

where *k* is the ordinate number, *ψ _{k}* is the radiance in the

*k*ordinate (where

^{th}*k*denotes the

*k*discrete angle), and

^{th}*ω*is a predetermined ordinate weight with full level symmetry [40,41]. The integral term does not require special treatment for implementation on BSGs.

_{k}The partial derivative terms are discretized using the upwind-step method. The directional cosines of a given ordinate determine the upwind-direction for that particular ordinate. There are *eight* possible numerical schemes, one for each octant in the three dimensional Cartesian coordinate system. For example, when all directional cosines are positive, the upwind-step method requires an Euler step in the negative *x*-, *y*-, and *z*-axis. For this example, the discretization of Eqs. (1),2) is given by

The term *p _{kk’}* in Eq. (10) is the Henyey-Greenstein phase function and is given by

where *g* is the anisotropy factor. The radiance, *ψ*, can be solved from Eqs. (10),11) with any number of established algorithms. In this work we implement the source iteration technique (i.e. the matrix-free point-wise *Gauss-Seidel Method*). The fluence is given by

obtained by applying Eq. (9), the *extended trapezoidal rule*, to Eq. (3) [20,38,40–42].

*2.3.2. Discretization on block-structured grid*: To solve the FD-ERT on BSGs, the Euler step used in the step-method, Eq. (10), must be changed to a step of variable size and become dependent on the local grid. For example, Δ*x* becomes Δ*x _{ijl}* where

*ijl*denotes the current grid point and its value is determined by the size of the local grid. Implementing a finite differences numerical scheme with a variable Euler step in a matrix-free formulation on BSGs is a complicated endeavor and requires great care. The main difficulty arises when solving for the FD-ERT on mesh points on a fine/coarse grid boundary (i.e. points that straddle both a coarse grid and a fine grid region).

In general, the numerical stencil requires each interior grid point (*i*,*j*,*l*) to have four neighbors in two-dimensions and six neighbors in three-dimensions. However, mesh points on a fine/coarse boundary do not always have a full set of neighbors. This problem can be overcome by adjusting the numerical scheme for each individual fine/coarse boundary point, or by creating the missing point so that the normal scheme is applicable. In this work we consider the latter option. The missing point (virtual point) is constructed through interpolation using neighboring points. For illustration consider the two-dimensional example in Fig. 5a
. Here, five cases must be considered independently for a given octant. For example, it is clear from Eq. (10) that when all directional cosines are positive, points *ii* and *iii* will require creating a virtual point interior to the coarse grid. However, points *i* and *iv* have all neighbors necessary to complete the stencil. The fifth case (black dots in Fig. 5a) can be treated as points on the fine grid or on the coarse grid. By treating these boundary points as coarse grid points we assure that they will always have a complete set of neighbors and no further special treatment is necessary. The set of points used to create the virtual point varies according to the type of boundary point. The different types of boundary points can be reduced to five in two dimensions (Fig. 5a) and nineteen in three dimensions.

There are many ways to interpolate values for the virtual points needed to complete the differencing schemes. Possible interpolation methods include averaging of the nearest four neighbors, as well as bilinear, biquadratic, and bicubic interpolation [42]. Higher order interpolation schemes are typically preferred because they are able to better resolve the curvature of the solution field. Higher order methods, however, are computationally more expensive. As a result, the type of interpolation scheme to be used must be chosen based on accuracy and computational cost. In this work we implement the four-neighbor averaging scheme because it yields sufficiently accurate results. As an example, consider Fig. 5b. Consider solving the FD-ERT when all directional cosines are positive on the interior grid point (*i*,*j*,*l*), represented by the black triangle. It is clear from Eq. (10) that the solution at grid point (*i*-1,*j*,*l*), represented by a black dot, is necessary to solve the equation. However, that grid point does not exist and must be created by averaging the solution at the four neighboring grid points denoted by white dots.

#### 2.4. Numerical Phantoms

We compare the accuracy and computation speed of our BSG algorithm to a previously published single grid algorithm that also uses a combination of discrete ordinates and the finite differences upwind-step method to solve the FD-ERT [4]. We test the performance of the BSG algorithm on four distinct numerical phantoms (Fig. 6 ): (a) a 2-cm-diameter disk, (b) an anatomically accurate cross section through a mouse obtained from a MRI data set, (c) a cylinder that is 2cm in height and diameter, and (d) a homogenous three-dimensional mouse phantom obtained from an MRI data set. We solve the FD-ERT on these phantoms on single Cartesian grids, on both coarse and fine grids, and on BSGs. We evaluate the performance of the BSG algorithm by comparing these solutions.

The disk phantom is of interest because it represents an instance where a coarse grid cannot place a single grid point on the true phantom surface, whereas a finer grid can always better approximate the true boundary. This is a worst-case scenario for the coarse grid. The anatomically correct phantom is of interest because it is representative of the arbitrarily shaped geometries encountered in practice where boundaries can have convex and concave regions. Generating a mesh for the anatomically correct phantom and solving the light propagation problem is a good test for both the mesh-generating routine and the FD-ERT numerical solver.

As can be seen from Fig. 6a, the disk phantom has two embedded fluorescent probes. The properties of the fluorescent probes are listed in Table 2
. The frequency modulated boundary sources are defined on the top-left quadrant of the disk (represented by black dots). Thus, instead of defining a single boundary source we define a constant source area. This is important because the number of boundary points on a given model increases with decreasing Δ*x*. With this setup we ensure the number of photons injected into the phantom is independent of the number of boundary points. The source density is 8 × 10^{9} photons cm^{−2} sr^{−1}. The optical properties of the disk phantom are varied and they are summarized in Tables 2 and 3
.

Figure 6b shows the two-dimensional mouse-like phantom. As in the disk-shaped phantom, this phantom has two embedded fluorescent probes. The boundary sources for excitation are on the top-left quadrant and are defined as external boundary sources for reasons similar to those presented for the disk phantom. Source density is 8x10^{9} photons cm^{−2} sr^{−1}. The specific optical properties used are summarized in Tables 2 and 3.

The grid spacings used for simulations in two dimensions are Δ*x* = 2/256, 2/128, 2/64, 2/32 cm (where Δ*x* = Δ*y*). For each grid refinement, we determine the solution to the ERT on 1-, 2-, and 3-level BSGs and compare them to the benchmark solution. The benchmark solution for these simulations is the solution on the finest single grid (Δ*x* = 2/256 cm). In addition, we vary the optical parameters (μ_{a}, μ_{s}), modulation frequency (ω), and refractive index (*n*) of each phantom and analyze the performance of the algorithm for each case.

Simulations on three-dimensional phantoms (Fig. 6c,d) are carried out using parameters from case 4 in Table 3. The benchmark solutions for the cylinder and mouse phantoms are defined on grids with Δ*x* = 2/128 cm and Δ*x* = 0.1 cm, respectively (Δ*x* = Δ*y* = Δ*z*).

#### 2.5. Quantification of computation time and solution accuracy

Computation cost and accuracy of a solution is defined relative to its corresponding benchmark. Computation time is reported in seconds. Benchmark solutions for determining accuracy are obtained on grids with Δ*x* = 2/256 cm for all phantoms introduced in Section 2.4. The relative speed up (RSU) achieved with the BSG algorithm is used to compare the computational cost of a solution obtained on a BSG (*T _{BSG}*) and on a single fine grid (

*T*), where the solutions we compare are partial current measurements at the boundary as defined by Eq. (4). All simulations were performed on a computer with a 2.93 GHz Intel Core 2 Duo processor.

_{F}The accuracy of the partial current is of particular interest because optical tomographic reconstruction (absorption, fluorescence, as well as bioluminescence) is an optimization process requiring accurate boundary data as input to the reconstruction algorithm [4,5,9]. Accuracy is quantified by the mean percent error (MPE) relative to the benchmark solution and is given by

where *f _{i}^{a}* and

*f*are the approximate and benchmark solutions, respectively. Only partial currents are compared and the MPE is reported. This is done because cross-sectional fluence solutions on coarse grids cannot, in general, be directly compared to solutions computed on fine grids since the numerical boundary of each grid is unique. This point is perhaps best illustrated by referring to Fig. 2b-d.

_{i}^{b}## 3. Results

In this section we summarize the computation time required for each forward problem and the accuracy of the solution at the boundary using the MPE measure. As previously mentioned in Section 2.4, we focus on the five specific combinations of parameters presented in Tables 2 and 3. We note that throughout this section, the cited grid spacing (Δ*x*) always refers to the spacing of the finest grid (i.e. the grid near the boundaries). As an example, a 2-level BSG with grid spacing Δ*x* has an embedded section of coarse grid points whose spacing is 2Δ*x*. The terms 1L, 2L, and 3L in all results tables refer to the number of grid levels in the BSGs.

#### 3.1 Block Structured Grids

We begin our analysis by showing representative examples of BSGs generated by our algorithm. We present 1- and 2-level grids fitted to the two- and three-dimensional phantoms (Fig. 7 ). These results demonstrate that the BSG generating section of our algorithm works properly with arbitrarily shaped geometries. We have also verified that the algorithm works with less complicated geometries.

In addition, we present representative examples of fluence computed on both a single grid and a 2-level BSG (Fig. 8 ). The sample solutions were computed using the anatomically correct mouse cross-section phantom with two embedded fluorophores and optical properties corresponding to case 1 in Table 3. Note that the solutions are qualitatively similar. The grids on which the sample solutions were computed are shown in Fig. 7a,d.

During implementation of the BSG algorithm it became evident that the thickness of every layer of the BSG with distinct spacing should be at least 2 grid points in every direction. Figure 7d is an example of a 2-level BSG where the outer fine grid is at least 3 grid points thick in each direction. We found that without this restriction the algorithm for solving the FD-ERT becomes more complex because the grid becomes less structured and it is then necessary to introduce more elaborate interpolation schemes. The algorithm we have developed allows control of the thickness of the outermost fine grid. All results presented in the subsequent section were obtained on grids with an outermost fine grid thickness of at least 2 elements. The computation time required for generating a BSG was always less than 0.5% percent of the time required for solving the FD-ERT.

#### 3.2 Disk Phantom

To perform a quantitative error analysis we first studied the circular disk. In Fig. 9
we present examples of typical partial current measurements on the boundary of the medium, for excitation and fluorescence emission computed on singled grids (Fig. 9a-b) and on BSGs (Fig. 9c-d). Solutions computed on single grids were computed on grids with Δ*x* = 2/256, 2/128, and 2/64 cm. Solutions were computed on 1-, 2-, and 3- level BSGs with Δ*x* = 2/256 cm. It is clear that, compared to solutions on single coarse grids, the solution computed on BSGs better approximates the true solution to fluorescence excitation and emission. The solution to fluorescence emission is less accurate than the solution to the excitation problem. This occurs because the numerical error from the excitation solution propagates into the emission solution. Thus, the solution to emission computed on BSGs is more accurate than the solution to emission computed on single coarse grids because the error in the solution to excitation on BSGs is lower.

The results for computation time, mean percent error, and relative speed up achieved with the BSG algorithm are summarized in Table 4 . Relative to a single fine grid, the solution computed on a 2-level BSG is obtained 1.9 to 2.2 times faster, while a solution on a 3-level BSG is obtained 3.1 to 4.2 times faster.

As expected, the solutions computed on 2- and 3-levels BSGs are less accurate than the solutions computed on a single fine grid. The increase in error is relatively small. For example, in case 1, with Δ*x* = 2/64 cm, the MPE of the solution obtained from a single grid is 3.3%. The solution on a 2-level BSG has MPE of 7.4% and is obtained 2.1 times faster than the solution computed on a single grid. As was explained before, the interior of the 2-level BSG has Δ*x* = 2/32 cm. The MPE of the solution on a single grid with Δ*x* = 2/32 cm is 18.9%. Therefore, the solution computed on a 2-level BSG is twice as accurate. Results from the other cases are similar.

We note that the errors in the solutions computed on the coarse grid (Δ*x* = 2/32 cm) in cases 1 and 4 are similar, i.e. the refractive index mismatch at the boundary of the circular phantom does not have a significant impact on the error in the solution. This is unexpected.

#### 3.3 Small Animal Phantom

The results of our quantitative error and performance analysis using the mouse cross-section phantom are summarized in Table 5
. In general, the solutions to the excitation problem using a 2-level BSG are obtained about 2 times faster than the solution computed on a single fine grid with a relatively small increase in error. The notable exceptions to these observations are the results from 2-level BSGs with Δ*x* = 2/32 cm. In these cases the solution on 2-level BSGs is obtained only 1 time faster than the single grid solution and the increase in error is very large. This occurs because the interior of the BSG is very coarse (Δ*x* = 2/16 cm) causing the numerical error to be large.

Comparing case 1 and case 4 is particularly interesting because the results are significantly different while the only difference in the simulations is the use of refractive index mismatch in case 4. In case 1, when there is no refractive index mismatch at the phantom surface, the error in the solution computed on the coarsest grid is 25.3%. However, in case 4, this same error increases to 43.9% when the phantom is assigned a refractive index of 1.37 (thus creating a mismatch at the air/tissue boundary). This is evidence that a finer mesh is required when taking into account the refractive index mismatch at the tissue boundary.

As expected, the performance of the algorithm with the given set of grid discretization is heavily dependent on scattering and absorption parameters. For example, the error in the solution computed on coarse grids is particularly large when the absorption coefficient is large (case 2). However, the error is relatively small even on the coarse grid when the scattering coefficient is very small (case 3). We note that the performance of the algorithm is not dependent on modulation frequency (case 5).

#### 3.4 Three dimensional phantoms

Representative examples of solutions to the FD-ERT on the cylindrical and three-dimensional mouse phantoms are shown in Fig. 10
. Results from simulations on the cylindrical phantom are summarized in Table 6
. The error in the solution computed on a 2-level grid (Δ*x* = 2/128 cm) is only 0.28%, while the solution computed on a single coarse grid (Δ*x* = 2/64 cm) 2.94%. Thus, solving the FD-ERT on a 2-level grid instead of a coarse grid reduced the error by 2.66%. Similarly, the solution computed on a 3-level grid (Δ*x* = 2/128 cm) is 30.21%, while the error in the solution computed on the coarsest grid (Δ*x* = 2/32 cm) is 78.25%. The error in the solution is reduced by 48.29%. In addition, using BSGs reduces computation time. Solutions on 2- and 3- level grids are obtained 1.5 and 3.0 times faster than the solution on the fine grid, respectively.

Results from simulations on the three-dimensional mouse phantom are summarized in Table 7 . In this example the error in the solution computed on a 2-level BSG is 28.42%. The solution on the 2-level BSG was obtained 1.09 times faster than the single grid solution.

## 4. Discussions and Conclusion

The motivation for solving the FD-ERT on BSGs arises from the need to reduce computation time without sacrificing accuracy of numerical solutions to the FD-ERT on structured grids. The errors in the solutions to the light propagation problem obtained on structured grids come from (1) inherent truncation errors from approximating a continuous equation with a discrete numerical scheme and (2) a mismatch between the numerical and physical boundary. Finely discretized grids can accurately model complex boundaries but solving the ERT on these grids is computationally expensive. Solving the ERT on coarse grids requires less computational effort; however, the accuracy of the solution can be very poor. A compromise between fine grids and coarse grids is the use of BSGs as presented in this paper. The interior of a BSG is primarily a coarse grid, while the outer layers of the grid, which are closer to the physical boundary, are more finely discretized.

In this work we presented the first algorithm for solving the FD-ERT on BSGs. This is also the first algorithm that incorporates reflective boundary conditions and a subroutine for generating BSGs directly into the algorithm for solving the FD-ERT. The BSG generator uses boundary information to determine the computational domain, discretizes it with a fine grid, then it adaptively coarsens the grid up to a user defined level. The final computation grid is a union of fine and coarse grids, where the coarse grid is restricted to the interior and the grid spacing becomes smaller near the boundary. The ERT is solved on this grid with a combination of the upwind step method and the discrete ordinates approximation. We use the blocked-off region method to treat curved boundaries.

Solutions to the ERT computed on single structured grids are corrupted by inherent error of the numerical approximation to the continuous equation and error due to poorly resolved boundaries. The inherent numerical error arises from the first order upwind step method approximation to the derivative terms and is proportional to spatial discretization (Δ*x*). In addition, there is numerical error due to the *S _{N}* approximation to the integral terms (it decreases with increasing order of the

*S*method). Error due to poorly resolved boundaries arises when the single Cartesian grid does not accurately approximate the physical boundary. This error is particularly large when a coarse grid is used to approximate curved geometries.

_{N}The total error associated with solutions to the ERT on BSGs is similar to the errors associated with solutions on a single coarse grid; however, there are two fundamental differences: (1) boundary errors are reduced because refining the grid near the boundary more accurately captures boundary effects and (2) a new source of error is introduced by interpolating at boundaries between coarse and fine grid sections inside the computational domain. The inherent numerical error associated with using the upwind step method on BSGs is similar to the numerical error expected from using this scheme on a single coarse grid because the majority of the BSG is composed of coarse grid points.

The need for BSGs becomes apparent when geometries are arbitrarily shaped. Results from 2D and 3D mouse phantoms confirm that the error due to a poorly resolved boundary is significant, especially when the refractive index mismatch at the air/tissue interface is taken into account. In this case, we found that the error in the solution computed on a single coarse grid was 25.3% when the refractive index mismatch was not taken into account. However, the error increased to 43.9% when there was refractive index mismatch at the boundary.

Simulations on a disk and mouse phantom show that solutions on BSGs are always obtained faster than the corresponding fine grid solutions. Solutions on 2- and 3-level BSGs were obtained in about 1/3 and 1/4 the time it took to obtain the same solution on a single fine grid, respectively. In general, the speed up achieved by the algorithm was not affected by changes in optical properties, refractive index, or modulation frequency.

Through analysis of simulations on the phantoms we show that solutions on BSGs are significantly more accurate than solutions on single coarse grids. Increasing the refinement of the grid near the boundary decreased the overall error in the solution for all cases studied. The general trend can be summarized as follows: the MPE of partial current measurements from solutions computed on 2- and 3-level BSGs are reduced to about 1/3 and 1/6, respectively, from the MPE of a solution obtained on a single coarse grid.

Overall we find that solving the FD-ERT on BSGs yields an attractive algorithm for modeling the light propagation problem on geometries with arbitrary shape without using dense and finely discretized Cartesian grids. The algorithm provides a method to substantially reduce the error due to poorly resolved boundaries. Furthermore, solutions to the FD-ERT on BSG are obtained at a much lower computation cost compared to solutions computed on a single fine grid.

## Acknowledgments

The authors would like to thank Dr. Hyun K. Kim (Columbia University) for helpful discussion regarding unstructured mesh algorithms. This work was funded in part by grants from the National Cancer Institute (NCI) (5U54CA126513-03: Tumor Microenvironment Network - The role of inflammation and stroma in digestive cancer and 4R33CA118666: Small animal tomography system for green fluorescent protein imaging), the National Institute for Arthritis and Musculoskeletal and Skin Diseases (NIAMS-2R01 AR46255: Optical Tomographic Imaging of Joint Diseases), which are all part of the National Institutes of Health (NIH). A.D. Klose was supported in part by a grant from the National Center For Research Resources (Award Number UL1RR024156).

## References and links

**1. **A. J. Welch, and M. J. C. Van Gemert, *Optical-thermal response of laser-irradiated tissue*, (Plenum Press, New York, NY, 1995).

**2. **S.A. Prahl, M. Keijzer, S.L. Jacques, A.J. Welch, “A Monte Carlo model of light propagation in tissue,” in *Dosimetry of Laser Radiation in Medicine and Biology*,Vol. 5 of SPIE Institute Series (SPIE, 1989), pp. 102–111.

**3. **L. H. Wang, S. L. Jacques, and L. Q. Zheng, “MCML--Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Methods Programs Biomed. **47**(2), 131–146 (1995). [CrossRef] [PubMed]

**4. **A. D. Klose, “Radiative Transfer of Luminescence in Biological Tissue”, in *Light Scattering Reviews, Volume 4*, A.A. Kokhanovsky (Ed.), 293–345 (Springer, Berlin, 2009).

**5. **A. D. Klose and A. H. Hielscher, “Iterative reconstruction scheme for optical tomography based on the equation of radiative transfer,” Med. Phys. **26**(8), 1698–1707 (1999). [CrossRef] [PubMed]

**6. **A. D. Klose, U. Netz, J. Beuthan, and A. H. Hielscher, “Optical tomography using the time-independent equation of radiative transfer - Part 1: forward model,” J. Quant. Spectrosc. Radiat. Transf. **72**(5), 691–713 (2002). [CrossRef]

**7. **A. H. Hielscher, R. E. Alcouffe, and R. L. Barbour, “Comparison of finite-difference transport and diffusion calculations for photon migration in homogeneous and heterogeneous tissues,” Phys. Med. Biol. **43**(5), 1285–1302 (1998). [CrossRef] [PubMed]

**8. **K. Ren, G. S. Abdoulaev, G. Bal, and A. H. Hielscher, “Algorithm for solving the equation of radiative transfer in the frequency domain,” Opt. Lett. **29**(6), 578–580 (2004). [CrossRef] [PubMed]

**9. **H. K. Kim and A. H. Hielscher, “A PDE-constrained SQP algorithm for optical tomography based on the frequency-domain equation of radiative transfer,” Inverse Probl. **25**(1), 015010 (2009). [CrossRef]

**10. **O. Dorn, “A transport-backtransport method for optical tomography,” Inverse Probl. **14**(5), 1107–1130 (1998). [CrossRef]

**11. **M. B. Salah, F. Askri, and S. B. Nasrallah, “Unstructured control-volume ﬁnite element method for radiative heat transfer in a complex 2-D geometry,” Numer. Heat Transf. B **48**(5), 477–497 (2005). [CrossRef]

**12. **S. R. Arridge, M. Cope, and D. T. Delpy, “The theoretical basis for the determination of optical pathlengths in tissue: temporal and frequency analysis,” Phys. Med. Biol. **37**(7), 1531–1560 (1992). [CrossRef] [PubMed]

**13. **S. R. Arridge, M. Schweiger, M. Hiraoka, and D. T. Delpy, “A finite element approach for modeling photon transport in tissue,” Med. Phys. **20**(2), 299–309 (1993). [CrossRef] [PubMed]

**14. **M. S. Patterson, B. Chance, and B. C. Wilson, “Time resolved reflectance and transmittance for the non-invasive measurement of tissue optical properties,” Appl. Opt. **28**(12), 2331–2336 (1989). [CrossRef] [PubMed]

**15. **M. S. Patterson and B. W. Pogue, “Mathematical model for time-resolved and frequency domain fluorescence spectroscopy in biological tissues,” Appl. Opt. **33**(10), 1963–1974 (1994). [CrossRef]

**16. **V. D. Liseikin, *Grid generation methods*, Second Edition (Springer, Netherlands, 2010).

**17. **J. C. Rasmussen, A. Joshi, T. Pan, T. Wareing, J. McGhee, and E. M. Sevick-Muraca, “Radiative transport in fluorescence-enhanced frequency domain photon migration,” Med. Phys. **33**(12), 4685–4700 (2006). [CrossRef] [PubMed]

**18. **X. Gu, K. Ren, and A. H. Hielscher, “Frequency-domain sensitivity analysis for small imaging domains using the equation of radiative transfer,” Appl. Opt. **46**(10), 1624–1632 (2007). [CrossRef] [PubMed]

**19. **A. D. Klose and A. H. Hielscher, “Fluorescence tomography with simulated data based on the equation of radiative transfer,” Opt. Lett. **28**(12), 1019–1021 (2003). [CrossRef] [PubMed]

**20. **A. D. Klose, V. Ntziachristos, and A. H. Hielscher, “The inverse source problem based on the radiative transfer equation in optical molecular imaging,” J. Comput. Phys. **202**(1), 323–345 (2005). [CrossRef]

**21. **A. D. Klose, B. J. Beattie, H. Dehghani, L. Vider, C. Le, V. Ponomarev, and R. Blasberg, “In vivo bioluminescence tomography with a blocking-off finite-difference SP3 method and MRI/CT coregistration,” Med. Phys. **37**(1), 329–338 (2010). [CrossRef] [PubMed]

**22. **J. M. Lasker, J. M. Masciotti, M. Schoenecker, C. H. Schmitz, and A. H. Hielscher, “Digital-signal-processor-based dynamic imaging system for optical tomography,” Rev. Sci. Instrum. **78**(8), 083706 (2007). [CrossRef] [PubMed]

**23. **B. W. Pogue and G. Burke, “Fiber-optic bundle design for quantitative fluorescence measurement from tissue,” Appl. Opt. **37**(31), 7429–7436 (1998). [CrossRef] [PubMed]

**24. **R. B. Simpson, “Automatic local refinement for irregular rectangular meshes,” Int. J. Numer. Methods Eng. **14**(11), 1665–1678 (1979). [CrossRef]

**25. **M. J. Berger and J. Oliger, “Adaptive Mesh Refinement for Hyperbolic Partial Differential Equations,” J. Comput. Phys. **53**(3), 484–512 (1984). [CrossRef]

**26. **W. L. Chen, F. S. Lien, and M. A. Leschziner, “Local mesh refinement within a multi-block structured-grid scheme for genereal flows,” Comput. Methods Appl. Mech. Eng. **144**(3-4), 327–369 (1997). [CrossRef]

**27. **M. J. Berger and P. Colella, “Local adaptive mesh refinement for shock-hydrodynamics,” J. Comput. Phys. **82**(1), 64–84 (1989). [CrossRef]

**28. **J. P. Jessee, W. A. Fiveland, L. H. Howell, P. Colella, and R. B. Pember, “An Adaptive Mesh Reﬁnement Algorithm for the Radiative Transport Equation,” J. Comput. Phys. **139**(2), 380–398 (1998). [CrossRef]

**29. **A. Joshi, W. Bangerth, K. Hwang, J. C. Rasmussen, and E. M. Sevick-Muraca, “Fully adaptive FEM based fluorescence optical tomography from time-dependent measurements with area illumination and detection,” Med. Phys. **33**(5), 1299–1310 (2006). [CrossRef] [PubMed]

**30. **A. H. Hielscher, “Optical tomographic imaging of small animals,” Curr. Opin. Biotechnol. **16**(1), 79–88 (2005). [CrossRef] [PubMed]

**31. **A. K. Scheel, M. Backhaus, A. D. Klose, B. Moa-Anderson, U. J. Netz, K. G. Hermann, J. Beuthan, G. A. Müller, G. R. Burmester, and A. H. Hielscher, “First clinical evaluation of sagittal laser optical tomography for detection of synovitis in arthritic finger joints,” Ann. Rheum. Dis. **64**(2), 239–245 (2005). [CrossRef] [PubMed]

**32. **B. J. Tromberg, L. O. Svaasand, T. T. Tsay, and R. C. Haskell, “Properties of photon density waves in multiple-scattering media,” Appl. Opt. **32**(4), 607–616 (1993). [CrossRef] [PubMed]

**33. **U. J. Netz, J. Beuthan, and A. H. Hielscher, “Multipixel system for gigahertz frequency-domain optical imaging of finger joints,” Rev. Sci. Instrum. **79**(3), 034301 (2008). [CrossRef] [PubMed]

**34. **V. Ntziachristos, “Fluorescence molecular imaging,” Annu. Rev. Biomed. Eng. **8**(1), 1–33 (2006). [CrossRef] [PubMed]

**35. **M. J. Niedre, R. H. de Kleine, E. Aikawa, D. G. Kirsch, R. Weissleder, and V. Ntziachristos, “Early photon tomography allows fluorescence detection of lung carcinomas and disease progression in mice in vivo,” Proc. Natl. Acad. Sci. U.S.A. **105**(49), 19126–19131 (2008). [CrossRef] [PubMed]

**36. **A. Soubret, J. Ripoll, and V. Ntziachristos, “Accuracy of fluorescent tomography in the presence of heterogeneities: study of the normalized Born ratio,” IEEE Trans. Med. Imaging **24**(10), 1377–1386 (2005). [CrossRef] [PubMed]

**37. **O. Gheysens and F. M. Mottaghy, “Method of bioluminescence imaging for molecular imaging of physiological and pathological processes,” Methods **48**(2), 139–145 (2009). [CrossRef] [PubMed]

**38. **K. M. Case, and P. F. Zweifel, *Linear transport theory*, (Addison-Wesley, Reading, 1967).

**39. **J. J. Duderstadt, and W. R. Martin, *Transport theory*, (John Wiley, New York, 1979).

**40. **E. E. Lewis, and W. F. Miller, *Computational Methods of Neutron Transport*, (Wiley, New York, 1984).

**41. **B. G. Carlson, and K. D. Lathrop, “Transport theory - the method of discrete ordinates”, in *Computing Methods in Reactor Physics*, H. Greenspan et al, eds. (Gordon and Breach, New York, 1968, pp. 166–266).

**42. **K. E. Atkinson, *An Introduction to Numerical Analysis*, 2nd ed., (John Wiley & Sons, Canada, 1989).