">
IOP P
UBLISHING
P
HYSICS IN
M
EDICINE AND
B
IOLOGY
Phys. Med. Biol. 53 (2008) 39213942
doi:10.1088/0031-9155/53/14/013
Fast iterative image reconstruction methods for fully
3D multispectral bioluminescence tomography
Sangtae Ahn
1
, Abhijit J Chaudhari
1,
3
, Felix Darvas
1,
4
,
Charles A Bouman
2
and Richard M Leahy
1
1
Signal and Image Processing Institute, University of Southern California, Los Angeles,
CA 90089, USA
2
School of Electrical and Computer Engineering, Purdue University, West Lafayette, IN 47907,
USA
E-mail:
leahy@sipi.usc.edu
Received 11 August 2007, in nal form 3 June 2008
Published 1 July 2008
Online at
stacks.iop.org/PMB/53/3921
Abstract
We investigate fast iterative image reconstruction methods for fully 3D
multispectral bioluminescence tomography for applications in small animal
imaging. Our forward model uses a diffusion approximation for optically
inhomogeneous tissue, which we solve using a nite element method (FEM).
We examine two approaches to incorporating the forward model into the
solution of the inverse problem. In a conventional direct calculation approach
one computes the full forward model by repeated solution of the FEM problem,
once for each potential source location. We describe an alternative on-the-y
approach where one does not explicitly solve for the full forward model.
Instead, the solution to the forward problem is included implicitly in the
formulation of the inverse problem, and the FEM problem is solved at each
iteration for the current image estimate. We evaluate the convergence speeds of
several representative iterative algorithms. We compare the computation cost
of those two approaches, concluding that the on-the-y approach can lead to
substantial reductions in total cost when combined with a rapidly converging
iterative algorithm.
(Some gures in this article are in colour only in the electronic version)
1. Introduction
Bioluminescence tomography is an in vivo imaging technique that localizes and quanties
bioluminescent sources in a small animal. It has recently gained a great deal of attention as
3
Present address: Department of Biomedical Engineering, University of California, Davis, CA 95616, USA.
4
Present address: Department of Neurological Surgery, University of Washington, Seattle, WA 98105, USA.
0031-9155/08/143921+22$30.00
© 2008 Institute of Physics and Engineering in Medicine
Printed in the UK
3921 3922
S Ahn et al
a promising means for macroscopic in vivo imaging of gene expression and other molecular
and cellular-level processes (Contag and Bachmann
2002
, Ntziachristos et al
2005
). The
objective of bioluminescence tomography is to reconstruct the three-dimensional (3D) spatial
distribution of bioluminescent sources inside the animal from images of the light emitted
through the animal surface. Collecting measurement data in multiple spectral bands and
from more than one view, so that the entire animal surface is imaged, facilitates tomographic
reconstruction of the 3D optical source distribution (Chaudhari et al
2005
).
We have set up a bioluminescence tomographic imaging system using mirrors to collect
four independent views (Chaudhari et al
2005
). Multispectral data were acquired to improve
localization (Chaudhari et al
2005
, Dehghani et al
2006
, Cong and Wang
2006
, Wang et al
2006b
, Kuo et al
2007
, Allard et al
2007
, Lv et al
2007
). We then used an FEM solution of
the diffusion equation to construct realistic forward models. We have successfully localized
bioluminescent sources with this system, using the TOAST software to compute the forward
model (Arridge et al
1993
). However, TOAST is primarily intended to solve the diffuse optical
tomography problem and is not well adapted for bioluminescence tomography applications;
the forward model computation cost is high.
Fully 3D multispectral bioluminescence tomography is a computationally challenging
inverse problem because (1) FEM-based forward models, which allow inhomogeneous tissue
and realistic geometry, require the inversion of a very large matrix to solve the forward problem,
(2) the use of multispectral data, which helps localize deep sources, increases the problem size
by the number of spectral bins and (3) the inverse problem is intrinsically ill-posed due to the
nature of the photon diffusion process and the limited information contained in data collected
only on the animal surface.
In this paper we investigate various numerical techniques in order to minimize the cost
of computing inverse solutions. The inverse problem is equivalent to a search for the source
distribution which best predicts the measured data while also satisfying an a priori assumption
regarding the smoothness of the source distribution. In order to iteratively solve the inverse
problem, a forward model which maps a source to the data domain must be computed
repeatedly.
We explore two different approaches to incorporating the forward solutions
into iterative algorithms for the inverse problem: (1) a direct calculation approach where the
forward problem is solved in advance for every source location and then a forward solution
is computed as a linear combination of these precomputed solutions, and (2) an on-the-y
approach where the FEM matrix inversion problem is solved as needed at each iteration using
the current estimate of the source conguration. In this way, the full forward model need
never be computed. We evaluate those approaches, combined with several different iterative
algorithms, to determine the combination of forward and inverse algorithm that will minimize
the computational cost.
Our goal in this paper is to introduce and evaluate methods that can rapidly reconstruct
3D bioluminescent images while still using a realistic FEM forward model. To the best of
our knowledge, the on-the-y approach is novel to this application and can substantially
reduce the total reconstruction time, as shown in section
4
.
The iterative algorithms
we consider in this paper include the class of incremental gradient or ordered subset
methods, which have been widely explored for applications in nuclear medicine imaging
(Hudson and Larkin
1994
, Browne and De Pierro
1996
, Ahn and Fessler
2003
) but
have not been investigated in the optical imaging literature. The algebraic reconstruction
technique (ART) (Gordon et al
1970
), which is a special case of the incremental gradient
algorithms, was used for diffuse optical tomography (Arridge and Schweiger
1998
). Our
comparison of reconstruction times for the different methods described above should help
guide the selection of reconstruction and forward modeling algorithms for developers Fast image reconstruction methods for bioluminescence tomography
3923
of 3D bioluminescent and uorescent imaging systems when using FEM-based forward
models.
In section
2
we formulate the forward problem based on the diffusion equation using FEM
models and describe approaches to incorporating the FEM solvers into iterative algorithms.
In section
3
we dene the inverse problem using regularized least squares (RLS) and describe
various iterative algorithms for computing the inverse solution. Finally, in section
4
we evaluate
different reconstruction methods in terms of computation times using in vivo bioluminescent
imaging data from a mouse with an implanted brain tumor.
2. Forward problem
In the forward problem one needs to predict the photon ux from the animal surface for a
given bioluminescent source distribution using known optical scatter and absorption properties
within the animal. Since iterative solutions to the inverse problem require multiple solutions
to the forward problem, a computationally efcient forward solution is important.
2.1. Forward model
Photon propagation in turbid media can be described by the radiative transfer or Boltzmann
transport equation (Arridge
1999
).
While this model can be used in bioluminescence
tomography (Klose
2007
), the high computation cost and detailed knowledge required of
the medias optical properties limit its use in practice. In contrast, the diffusion approximation
under the assumption of isotropic scattering has been used extensively (Arridge et al
1993
,
Schweiger et al
1995
, Jiang
1998
, Arridge
1999
, Hielscher et al
1999
, Gu et al
2004
, Wang
et al
2004
, Cong et al
2005
, Slavine et al
2006
, Dehghani et al
2006
, Comsa et al
2006
,
Soloviev
2007
, Kuo et al
2007
). The diffusion approximation is reasonably accurate in soft
tissue in the near-infrared region where scattering dominates absorption (Arridge
1999
, Shen
et al
2007
). We, therefore, focus on the forward model based on the following steady-state
diffusion equation (Arridge
1999
),
{ · (r, ) + µ
a
(r, )
}(r, ) = q(r, ),
for
r
,
(1)
subject to a Robin boundary condition
(r, )
+ 2(r, )G
n(r)
· (r, ) = 0
for
r

(2)
where q represents the bioluminescent source distribution, denotes the photon density,
r
R
3
denotes the location vector,
R
3
is the animal volume,
is the animal surface, n is a unit vector pointed outwardly normal to the surface , is the wavelength and
(r, )
= 1/[3{µ
a
(r, )
+ µ
s
(r, )
}], with µ
a
and µ
s
being the absorption and reduced
scattering coefcients, respectively. In the boundary condition (
2
), G is a parameter modeling
internal reection at the boundary and can be computed as G
= (1 + )/(1 ), where 1.4399n
2
int
+ 0.7099n
1
int
+ 0.6681 + 0.0636n
int
with n
int
being the refractive index
of tissue (Schweiger et al
1995
).
The measurable photon ux at r

is given by
m(r, )
= (r, )/(2G). The mapping from the source distribution q in the volume to
the photon ux m through the surface can be obtained by solving (
1
) for given q.
Approaches to solving the diffusion equation (
1
) inc