Research Article

Vol. 37, No. 2 / February 2020 / Journal of the Optical Society of America A

231

Fast and robust reconstruction algorithm forfluorescence diffuse optical tomography
assuming a cuboid target
Chunlong Sun,1,2 Gen Nakamura,2 Goro Nishimura,3Manabu Machida5, *Yu Jiang,4 Jijun Liu,1 AND

1School of Mathematics, Southeast University, Nanjing 210096, China
Department of Mathematics, Hokkaido University, Sapporo 060-0810, Japan
3
Research Institute for Electronic Science, Hokkaido University, Sapporo 001-0020, Japan
4
School of Mathematics, Shanghai University of Finance and Economics, Shanghai 200433, China
5
Institute for Medical Photonics Research, Hamamatsu University School of Medicine, Hamamatsu 431-3192, Japan*Corresponding author: machida@hama-med.ac.jp
2

Received 16 September 2019; revised 4 December 2019; accepted 4 December 2019; posted 4 December 2019 (Doc. ID 377152);
published 10 January 2020

A fast algorithm for fluorescence diffuse optical tomography is proposed. The algorithm is robust against the choiceof initial guesses. We estimate the position of a fluorescent target by assuming a cuboid (rectangular parallelepiped)
for the fluorophore target. The proposed numerical algorithm is verified by a numerical experiment and an experiment with a meat phantom. The target position is reconstructed with a cuboid from measurements in the time
domain. Moreover, the long-time behavior of the emission light is investigated making use of the analytical solution
to the diffusion equation. © 2020 Optical Society of America
https://doi.org/10.1364/JOSAA.37.000231

1. INTRODUCTION
Fluorescence diffuse optical tomography (FDOT) is one typeof optical tomography that makes use of fluorescent light fromfluorophore. In FDOT, diffuse light from fluorophore, such as
indocyanine green (ICG), is detected on the boundary of biological tissue to obtain tomographic images [1,2]. FDOT can beformulated as an inverse source problem. As theoretical results,
in addition to the usual L 2 regularization, the L 1 regularization
was tested [3–6], and the total variation was considered [7].
An improvement in image quality was reported with the total
variation by making use of the Bregman distance [8]. In dataacquisition, the measured fluorescent light can be divided by themeasured excitation light to cancel unknown constants, which
was introduced as normalized measurements [9].
FDOT has been verified in vivo [10,11] and also in clinicalresearch for breast cancer [12–14]. Noncontact measurementswith a CCD camera were established, and fluorescence yield wasreconstructed for mice with lung tumors [15]. The resolution of
tomographic images is a function of the number of sources anddetectors. A large number of parallel channels can be employed
by a detection fiber array coupling onto a CCD camera, whichleads to fluorescence molecular tomography finding a target ofa few millimeters [16,17]. A tumor in a mouse was imaged by
the use of a conical mirror, by which light is collected by a CCD
1084-7529/20/020231-09 Journal © 2020 Optical Society of America

camera without attaching optical fibers to the mouse [18]. As
an alternative, the number of measured data can be increasedby time-dependent experiments [19]. The superiority of thetime-resolved approach over the continuous-wave approach in
FDOT was concluded [20], although the quality of images canbe improved even for a time-independent experiment if a large
data set is used [21]. In FDOT, usually it is not easy to know
precise values of absorption and scattering coefficients. The
reduction in the reconstruction errors due to this inaccuracy is
proposed with the Bayesian approximation error approach [22].
Among three types of FDOT (continuous wave, frequency
domain, and time domain), the measurement in time domain
has the most fruitful information. Although the governing
equation for near-infrared light in the frequency domain
reduces to the time-independent diffusion equation, the timeresolved curve for the time-domain FDOT is calculated fromthe time-dependent diffusion equation. Usually, the information carried by measurements in time domain is not fully
explored. By a transformation in time such as Fourier transform,
Laplace transform, Mellin transform, Mellin–Laplace transform, etc., the problem is reduced to an inverse problem for thetime-independent diffusion equation [23].
A straightforward way of solving the linear inverse problemof FDOT is the direct inversion with singular value decomposition (SVD). Although SVD is proven to be successful [21],232

Research Article

Vol. 37, No. 2 / February 2020 / Journal of the Optical Society of America A

the method is computationally expensive from the viewpointof memory and speed. An alternative approach is the use of
iterative methods. For example, the Levenberg–Marquardt(LM) method [24,25] has been successfully used for FDOTof breast cancer in the frequency domain [12]. By the regularization of total variation, edges in reconstructed images are notsmoothed. However, compared with the usual L 2 regularization, some care is necessary for the use of total variation. To
this end, the split Bregman iteration was developed [26,27]
and used in Ref. [7]. The three-dimensional FDOT in the time
domain with the time-dependent diffusion equation is still achallenging problem. In Ref. [19], reconstructed images for thetime-domain FDOT by the time-dependent diffusion equation were obtained in two dimensions by an iterative schemeof the Newton–Raphson inversion. Two-dimensional timedependent diffusion equations were solved for the FDOT by
the total light approach [28,29]. Although the original coupled diffusion equations are simplified by this approach, thecorresponding inverse problem becomes more difficult, since
the approach makes use of the weak nonlinearity of the FDOT
inverse problem, which is almost linear.
In this paper, we consider FDOT by measuring time-resolved
data from reflected light. Aiming at a fast and robust numericalalgorithm, we focus on the reconstruction of the position of
a fluorophore target, neglecting its precise shape. We assume
a cuboid or rectangular parallelepiped in the medium. Then
we try to find the cuboid that represents the true target. Dueto the assumption of the cuboid shape, the reconstruction of
our FDOT can be done by determining only several unknown
parameters. For the robustness and fast convergence of ourreconstruction scheme, we narrow the region of interest during
the reconstruction process. Before assuming a cuboid, we use
a cubic to have a good initial guess. The initial guess for the
cubic is chosen using the spatial knowledge of the intensity fromboundary measurements. The proposed reconstruction scheme
was tested both numerically and experimentally. We will consider the three-dimensional FDOT using the time-dependent
diffusion equation, which is the most computationally expensive FDOT. Using the proposed method, we show that thethree-dimensional FDOT in the time domain is feasible on alaptop computer (Apple MacBook Pro, 2.7 GHz Intel Core i5).
If we consider continuous-wave or frequency-domain measurements, two numerical integrals are removed in our formulation[see Um in Eq. (13)] and the computational time will be further
reduced. Furthermore, we study the asymptotic behavior of
the temporal profile of the emission light using the diffusion
equation.
The paper is organized as follows. In Section 2, we develop
the formulation of our FDOT and give an analytical formula forthe emission light. The proposed numerical scheme is describedin Section 3 with a numerical example of an ellipsoidal target. In
Section 4, we validate our numerical method with a beef sample
in which a fluorescent target is embedded. The detected emission light carries information on the target to be reconstructed.
In Section 5, we study the time dependence of the emission
light for a point fluorophore. Finally, Section 6 is devoted toconclusions.

2. FORMULATIONLet us suppose that a fluorescent target is embedded in biologicaltissue occupying the half space (−∞ < x < ∞, −∞ < y < ∞,0 < z < ∞). Let c be the speed of light in the medium. Letu e (r, t), u m (r, t) be the energy densities of the excitation
light and emission light, respectively. Here, r = (ρ, z), withρ = (x , y ) as the position and t as time. Boundary measurements are performed during 0 < t < T. We assume that the
reduced scattering coefficient µ0s and the absorption coefficientµa are constants everywhere in the medium. Although thewavelength of the emission light is longer than the wavelengthof the excitation light, and thus values of µa and µ0s are different
for the excitation light and emission light, we ignore this difference. Moreover, compared with µa , the absorption coefficientfor the fluorophore is negligible because it is nonzero only at
the position of the fluorophore. Then in the medium (z > 0),u e (r, t) and u m (r, t) obey the following diffusion equations(e.g., [19,20]):
1 ∂− D1 + µa u e = 0,
(1)
c ∂t


1 ∂
− D1 + µa u m = F ,(2)c ∂twhere D = 1/(3µ0s ), andF (r, t) =

n(r)
τ

Z t

e −(t−s )/τ u e (r, s )ds .

(3)0Here, n(r) is proportional to the fluorophore concentration,
whose proportionality constant is the product of the quantum
yield and absorption cross section, and τ is the fluorescence lifetime. We assume that u e = 0 and u m = 0 at t = 0. The sample
is illuminated at position rs = (ρ s , 0) by a pencil beam of thetemporal profile h(t) (t > 0) in the x −y plane. At z = 0, u e , u msatisfy the Robin boundary conditions as∂
u m + βu m = 0.
∂z
(4)
R1
1
The parameter β is given by β = 2D
(1 − 2 0 R(µ)µdµ)/(1 +R120 R(µ)µ dµ) with the Fresnel reflectance R(µ), which
depends on the refractive index of the medium [30]. Suppose
the out-going light is detected at rd on the boundary. The excitation light and emission light are detected through a response
function R as
Z t
Ue (rd , t; rs )=R(t − s )u e (rd , s )ds ,−∂u e + βu e = h(t)δ(ρ − ρ s ),
∂z

−

0

Um (rd , t; rs )=Z tR(t − s )u m (rd , s )ds ,(5)0where the function R is determined by the detector. These
Ue , Um correspond to experimentally measured lightexp
exp
Ue (rd , t; rs ), Um (rd , t; rs ).
0
Let G(r, r ; t) be the Green’s function that satisfies1 ∂− D1 + µa G(r, r0 ; t) = δ(r − r0 )δ(t), (6)c ∂tResearch Article

Vol. 37, No. 2 / February 2020 / Journal of the Optical Society of America A

with the initial condition G = 0 at t = 0 and the Robin bound∂ary condition − ∂zG + βG = 0. The Green’s function is
obtained as [31–34]

233

where
f 1 (ρ d , ρ s , t, s ; x 1 , x 2 , y 1 , y 2 ) =(x −x )2 +(y −y )2e −µa c t
− d s 4Dc t d se√
43 π 2 D2 t (t − s )sG(r, r0 ; t) = c (4π Dc t)−3/2 e −µa c t(x − x 0 )2 + (y − y 0 )2
× exp −
4Dc t!
g (z, z 0 ; t),× [ f 3 (x d , x s, t, s; x 2 ) − f 3 (x d , x s, t, s; x 1 )](7)

× [ f 3 (y d , y s, t, s; y 2 ) − f 3 (y d , y s, t, s; y 1 )] ,

(14)

and

where
g (z, z 0 ; t) =f 2 (t, s ; z1 , z2 ) =
0 2

exp −

(z + z )
4Dc t

!

0 2

+ exp −(z − z )4Dc t

× exp β(z + z ) + β Dc t erfc
2

0





!

We introduce the instrument response function as
Z tR(t − s )h(s )ds .
q (t) =

.

(8)

c

where we introduce
Z te −t /τ q (t − t 0 )dt 0 .
0

(12)

0

Ii =
0

Um (rd , t; rs ) =0.

The reconstruction of a fluorophore, which in general has different shapes, is considered. In our proposed inversion scheme,
we will reconstruct the fluorophore by identifying a cuboidn(r). The cuboid is expected to specify the location of the fluorophore. In this way, our FDOT is achieved by using only seven
unknown parameters: M, x 1 , x 2 , y 1 , y 2 , z1 , z2 . Although thenumber of unknown parameters is significantly fewer than othermethods that use voxels, the original linear inverse problem ofreconstructing n(r) from Um in (11) becomes nonlinear. The
reconstruction will be done with the LM method. Even though
there are only seven parameters, the choice of initial guesses isimportant. For the sake of robust reconstruction, we propose
to obtain a set of good initial guesses for the seven parameters
by narrowing the region of interest according to three steps
described below.
Let NSD be the number of source–detector pairs. The
positions of each source and detector are denoted by ris , rid(i = 1, ... , NSD ). At each detector, light is measured at Nt
temporal points.
We first look for a rough location of the target projected to
the x −y plane or a region of interest 0 on the boundary, underwhich the target is likely to be embedded, by observing
Z T

Equation (11) can be rewritten as

Q(s )



We can compute Um (rd , t; rs ) using (13).

(11)

0

(16)



3. IDENTIFICATION OF A CUBOID

0

Z t−s

t
4Dc (t − s )s(9)× G(r0 , rs ; s 0 )dr0 ds 0 ds ,

M

r
erf

where M > 0 is a constant. Then we obtain
Z t−sZZ tn(r0 )G(rd , r0 ; t − s − s 0 )Um (rd , t; rs ) = Q(s )Z t

s x d + (t − s )x sx−tf 3 (x d , x s , t, s ; x ) =The fluorophore target that we try to reconstruct may havedifferent shapes. In this paper, we focus on the position of thetarget and solve the inverse problem assuming a cuboid. Thatis, we will find the location of the target using a cuboid. In theformulation below, we assume that the target has the shape of a
cuboid.
Let c be a cuboid specified by x 1 , x 2 , y 1 , y 2 , z1 , z2 as x 1 <x < x 2 , y 1 < y < y 2 , z1 < z < z2 . We assume thatM, r ∈ c ,n(r) =(10)
0, r ∈/ c ,

D
Q(t) =
τ(15)

Here, we define00

g (0, z 0 ; t − s )g (z 0 , 0; s )dz 0 .

z1√− 2β π Dc t

z + z 0 + 2β Dc t
√
4Dc t

Z z2

f 1 (ρ d , ρ s , t − s , t ; x 1 , x 2 , y 1 , y 2 )0

× f 2 (t − s , t 0 ; z1 , z2 )dt 0 ds ,
(13)

exp

Um (rid , t; ris )dt,

i = 1, ... , NSD .

(17)

The observation time T = 3 ns for the numerical experiment inthis section and T = 10 ns for the meat phantom experiment in
Section 4. We call this step of setting 0 the topography process.
By using Ii with large values, we set the region 0.
Then we assume a cubic target whose location and size aredetermined by four parameters x 0 , y 0 , z0 , l , such that234

Vol. 37, No. 2 / February 2020 / Journal of the Optical Society of America A

lx1 = x0 − ,2lx2 = x0 + ,2lz1 = z0 − ,2lz2 = z0 + .
2ly1 = y0 − ,2ly2 = y0 + ,2(18)

Choosing the initial guess for x 0 , y 0 inside 0, we solve theinverse problem of determining x 0 , y 0 , z0 , l , M by the LM
method. We note that the choice of M is not strict comparedwith other parameters. The obtained cubic is denoted byacubic= (x 00 , y 00 , z00 , l 0 , M 0 ).
∗(19)

We refer to this step as the cubic tomography. Since we have only
five parameters instead of seven, there is a good stability for the
inverse problem in this step.
Finally, with the obtained values acubic
as the initial guess∗cuboida0, we solve the inverse problem of determining a cuboid
target. In this way, we can identify c for n(r) by obtaining
reconstructed values:
acuboid= (x 1∗ , x 2∗ , y 1∗ , y 2∗ , z1∗ , z2∗ , M ∗ ).
∗(20)

This last step is called the cuboid tomography. The algorithm issummarized as follows.
Step 1. (Topography process) Find 0 on the boundary.
Step 2. (Cubic tomography) By searching underneath 0,
obtain reconstructed values acubic.
∗
Step 3. (Cuboid tomography) Find acuboid, starting with∗acubic.
∗Below, we will demonstrate our algorithm with a numericalexperiment. We set τ = 0, h(t) = R(t) = δ(t). Let us set an
ellipsoidal target as{x /1.5 + y /3 + (z − 11) /1.5 ≤ 1},
2

2

2

2

2

2

(21)

Research Article
exp

where 1t = 6.67 ps, and t pi is the peak time of Um (rid , t; ris )(i = 1, ... , NSD ). Thus, we have NSD Nt (= 640) measured valexp
ues Um (rid , tki ; ris ) (i = 1, ... , NSD , k = 1, ... , Nt ). Let a be avector that contains unknown parameters to be reconstructed.
exp
We find a = a∗ by fitting Um to Um using the LM method
implemented on MATLAB [35,36].
(Step 1) We find I4 = I10 = I17 = I27 = 5.9 × 10−7 , whereas
I1 = I11 = I20 = I26 = 1.0 × 10−9 . Hence, we can set
0 = {−10 < x 0 < 10, −10 < y 0 < 10}.

(25)

At this moment, we have no knowledge about the depth atwhich the target is embedded. We suppose 0 < z0 < 30 together
with 0 < l < min(20, 2z0 ) and 0 < M < 10.
(Step 2) We set the initial values of x 0 , y 0 in 0. For example,
we can start from acubic= (x 0 , y 0 , z0 , l , M) = (2, 2, 5, 4, 0.1).
0
The reconstructed values are
acubic= (x 00 , y 00 , z00 , l 0 , M 0 ) = (0.0, 0.0, 11.24, 4.089, 0.0086).
∗(26)
We note that choosing x 0 , y 0 in 0 is important for fast convergence. If we pick a point outside 0 and set (x 0 , y 0 ) =(−15, −15), more than 100 iterations are needed, whereas
about 10 iterations are sufficient for x 0 , y 0 ∈ 0.
(Step 3) Now we give the initial guess acuboidfrom acubic
as0∗
acuboid= (x 1 , x 2 , y 1 , y 2 , z1 , z2 , M)
0

l0
l0
l0
l0
= x 00 − , x 00 + , y 00 − , y 00 + ,
2
2
2
2

l0
l0
z0 0 − , z0 0 + , M 0 .
2
2

(27)

Then assuming that the target is a cuboid, we obtainreconstructed values acuboidas∗where the unit of length is mm, with n(r) = 0.02 mm−1 . In the
numerical experiment, 5% noise is added, andexp

Um = Um (1 + 0.05ε),(22)

where ε is drawn from the standard Gaussian distribution. We
use the following NSD = 32 source–detector pairs p i = (ris ; rid )(i = 1, ... , NSD ) on the boundary:√p 4 j −3 = (ξ j , ζ j + 10 3; ξ j − 10, ζ j ),√p 4 j −2 = (ξ j , ζ j + 10 3; ξ j + 10, ζ j ),√p 4 j −1 = (ξ j , ζ j − 10 3; ξ j − 10, ζ j ),√p 4 j = (ξ j , ζ j − 10 3; ξ j + 10, ζ j )(23)

for j = 1, ... , 8. Here, (ξ1 , ζ1 ) = (−10, 10), (ξ2 , ζ2 ) =(−10, 0),(ξ3 , ζ3 ) = (−10, −10),(ξ4 , ζ4 ) = (0, −10),(ξ5 , ζ5 ) = (10, −10), (ξ6 , ζ6 ) = (10, 0), (ξ7 , ζ7 ) = (10, 10),and (ξ8 , ζ8 ) = (0, 10). We set Nt = 20 andtki = t pi + (k − 10)1t(k = 1, ... , Nt ),(24)

Fig. 1. Identification of the ellipsoidal target by a cuboid. The
obtained cuboid given by acuboid(blue) is shown with the ellipsoid∗(red).

Research Article

Vol. 37, No. 2 / February 2020 / Journal of the Optical Society of America A

235

Fig. 2. Same as Fig. 1 but cross sections are shown. The cuboid (blue) and ellipsoid (red) are shown, from the left, on the plane at z = 11, on thex −z plane, and on the y −z plane, respectively.

= (x 1∗ , x 2∗ , y 1∗ , y 2∗ , z1∗ , z2∗ , M ∗ )acuboid∗= (−1.074, 1.036, −2.146, 2.166, 9.908, 12.02, 0.029).

(28)
The reconstructed cuboid is shown in Fig. 1 with the ellipsoidal
target. The cuboid (blue) and ellipsoid (red) are shown in twodimensional planes in Fig. 2. We see that the center position of
the ellipsoidal target is correctly reconstructed, and the volumeof the target is well recovered by the cuboid. Furthermore, thereconstructed M ∗ is close to the true value 0.02 mm−1 .
Figure 2 shows that the reconstructed cuboid is the one thathas the largest area in the ellipse on each projected plane: x −y ,x −z, or y −z. In the x −y plane, the cross section of the targetisp an ellipse of area π(1.5)(3) ≈ 14.1. From the inequality(x 2 /1.52 )(y 2 /32 ) ≤ 12 (x 2 /1.52 + y 2 /32 ) = 12 , we see that

Figure 3 shows how the tube is placed in the meat phantom. The
tube has the shape of a cylinder of length 8 mm and diameter
2 mm. Time-dependent measurements were conducted using
a holder placed on the top of the beef sample as shown in Fig. 3.
Four optical fibers (two are for sources, and the other two are fordetectors) are attached to the holder. The motorized stage on
which the meat phantom is placed changes positions while theholder is fixed. A near-infrared pulsed laser at 780 nm is used toexcite the fluorescent target. See [37] for the experimental setup.
Since the source–detector separation is 2 cm, we can assume thehalf space.
Optical parameters for the meat phantom are µ0s =0.92 mm−1 , µa = 0.023 mm−1 . The refractive index is set to1.37. Moreover, τ = 0.6 ns [38]. By neglecting source–detector
pairs with tiny signals, we can use NSD = 16 source–detector
pairs p i = (ris ; rid ) (i = 1, ... , NSD ):√√√√p 1 = (−5 − 10 3, 0; −5, 10), p 2 = (−5 + 10 3, 0; −5, 10), p 3 = (−5 − 10 3, 5; −5, 15), p 4 = (−5 + 10 3, 5; −5, 15),√
p 5 = (−10 3, 0; 0, 10),√p 9 = (5 − 10 3, 0; 5, 10),√p 6 = (−10 3, 5; 0, 15),√
p 7 = (5 − 10 3, 5; 5, −5),√
p 10 = (5 − 10 3, −5; 5, 5),√
p 13 = (−15 + 10 3, 0; −15, 10),√p 8 = (5 − 10 3, 5; 5, 15),√
p 11 = (−10 3, −5; 0, 5),√
p 14 = (−15 + 10 3, 5; −15, −5),√
p 12 = (−10 + 10 3, 0; −10, 10),√p 15 = (−15 + 10 3, 5; −15, 15),√
p 16 = (−10 + 10 3, 5; −10, 15),(29)
the maximum area of the rectangle inside the ellipse is nine,since the area is (2x )(2y ) ≤ 4(1.5 × 3)/2 = 9. The area ofthe cross section of the reconstructed cuboid (Fig. 2 (Left))is (x 2∗ − x 1∗ )(y 1∗ − y 2∗ ) = 9.098. Similarly, in the x −z plane[Fig. 2 (middle)] and y −z plane [Fig. 2 (right)], the cross sections of the cuboid have the largest areas in the corresponding
ellipses.
We emphasize that the narrowing process 0 → acubic→0
cuboida0
is essential. If the initial guess acubic
for
step
2
is
used
as
the
0
initial guess for the cuboid tomography, the iteration of the LM
method does not converge except for some special cases.
4. BEEF EXPERIMENT
Let us reconstruct a fluorescent tube embedded in beef with ourreconstruction scheme. The tube contains 1 µM ICG solution.

where the unit of coordinates is mm. The background
fluorescence is subtracted from the signal at each source–
detector pair. Since√ the signal from the source–detector pair
p 0 = (−15 − 10 3, 5; −15, 15) contains little contribution
from the target fluorescence, we used the emission light fromp 0 as the background fluorescence. Sources and detectors areshown in Fig. 4. Red small disks show sources and blue circles
show detectors. The source–detector pairs are denoted by pinklines. Moreover, the green rectangle shows the position of thereconstructed tube in (32) projected on the x −y plane. We
set Nt = 20 and 1t = 6.1 ps. At the ith source–detector pair(i = 1, ... , NSD ), measured times used for reconstruction are
tki = t pi + (k − 11)1t (k = 1, ... , Nt ), where t pi is the peak
exptime of Um (rid , t; ris ).
Since measured values are large in this region, suchasI5 = 3.3 × 104 ,I6 = 4.0 × 104 ,I12 = 2.9 × 104 ,236

Research Article

Vol. 37, No. 2 / February 2020 / Journal of the Optical Society of America A

Fig. 3. (Left) A tube of ICG was placed inside the beef. (Middle) The tube was embedded at a depth of about 16 mm. (Right) Boundary measurements were performed using optical fibers attached to a holder on the top of the beef.

y

p16

p8

p15
p0
p3

p6

p1

p5
p11

p4

p9

p14

p12
p10

p2

x

p13
p7
5mm

Fig. 4. Schematic figure of measurements on the boundary. Red small disks show sources, and blue circles show detectors. The source–detector
pairs are denoted by pink lines. The green rectangle shows the position of the reconstructed tube projected on the x −y plane.

I16 = 3.4 × 104 , first we set (step 1)
0 = {−10 < x < 0, 5 < y < 20}.

(30)

Although the value of n(r), i.e., M for the cube and cuboid, isa parameter to be reconstructed, in this beef experiment, M isdetermined only up to a constant that comes from the propertyof the instrument. We set (x 0 , y 0 , z0 ) = (−5, 10, 7), l = 2, andobtain (step 2)(x 00 , y 00 , z00 ) = (−4.12, 7.72, 17.25),

l 0 = 3.92.

(31)

With the above values as the initial guess, we use the iterative
method once again (step 3). The reconstructed values are
obtained as
x 1∗ = −5.16,

x 2∗ = −3.11,z1∗ = 16.05,

z2∗ = 16.34.

y 1∗ = 3.83,

The exact position of the embedded tube in beef is notknown. If the tube is parallel to the y axis as the green rectangle in Fig. 4, the reconstructed diameter of the tube is
x 2∗ − x 1∗ = 2.05 mm, and the reconstructed length of the
tube is y 2∗ − y 1∗ = 8.2 mm, where the true diameter and length
are 2 mm and 8 mm, respectively. The reconstructed depth
position is estimated as (z1∗ + z2∗ )/2 = 16.2 mm, which isconsistent with the depth shown in Fig. 3 (middle). Comparedwith these reconstructed values, the reconstructed thicknessz2∗ − z1∗ = 0.29 mm is different from the diameter of the tube.
In light of the fact that the ellipsoid is correctly reconstructed in
Section 3, the reason for the short reconstructed thickness may
be attributed to experimental conditions such as the number
and positions of source–detector pairs and measurement noise.

y 2∗ = 12.03,
5. LONG-TIME BEHAVIOR(32)

The position of the cuboid is what we expected. In particular,
we see that the orientation of the cylinder must be almost parallel to the y axis as shown in Fig. 4. Thus, the position of thefluorescent tube is reconstructed with the proposed numerical
scheme.

The measured emission light Um has information on the fluorophore target. Let us investigate the long-time behavior of
Um . We will see below that the fluorescence lifetime τ can be
estimated by the long-time behavior of Um .
For simplicity, we consider a point targetn(r) = δ(r − rc ).

(33)

Research Article

Vol. 37, No. 2 / February 2020 / Journal of the Optical Society of America A

That is, c is a point rc . According to (11), the emission lightUm (rd , t; rs ) can be calculated as
Z t
1
Um (rd , t; rs ) =Q(t − s )e −µa c s
16(π D)3 c 0
Z s0×ξ(s , s 0 )e ζ (s ,s ) ds 0 ds .
(34)
0Here we introduced
1[(s − s 0 )s 0 ]3/2"
p× 1 − β π Dc (s − s 0 )Wξ(s , s 0 ) =zc + 2β Dc (s − s 0 )p
4Dc (s − s 0 )

!#




√
zc + 2β Dc s 0
× 1 − β π Dc s 0 W,√
4Dc s 0(35)

237

Since we are interested in the long-time behavior ofUm (rd , t; rs ), we take the Laplace transform for t and takeonly small values of the Laplace variable into account. Then forlarge t, we have1|rs − rc | + |rd − rc |Um (rd , t; rs ) ∝ exp −t−.
√τ
2c µa D
(43)
Therefore, we obtain−1∂τ = − ln Um (rd , t; rs )
.
(44)
∂t
In Sections 3 and 4, we have developed the reconstruction
algorithm using a pre-determined fluorescence lifetime, which
is not always known a priori. Kumar et al. pointed out that thefluorescence lifetime can be estimated from the asymptotic
behavior of the temporal profile of the emission light [39,40].
Here, we showed that their formula also holds true in thepresence of the boundary.

andζ (s , s 0 ) = −|rd − rc |2|rs − rc |2−,
4Dc (s − s 0 )4Dc s 0(36)

where we write2

erfc(x ) = e −x W(x ).
√We note that W(x ) ∼ 1/(x π ) for large x , andzc + 2β Dc t p= 2βzc .
√t∈(0,T)
4Dc t
min

(37)

(38)

Let us suppose
zc 

1
.
2β

(39)

Using the method of steepest descent, we obtainUm (rd , t; rs ) ≈zc2(|rs − rc | + |rd − rc |)3√|rs − rc ||rd − rc |8(π D)5/2 cZ te −µa c s(|rs − rc | + |rd − rc |)2×
Q(t − s ) 3/2 exp −
s
4Dc s
0
×

1
1
ds ,
C 0 zc + 2β Dc |rs − rc |s C 0 zc + 2β Dc |rd − rc |s(40)

where C 0 = |rs − rc | + |rd − rc |. Below, we put h(t) = R(t) =δ(t) for simplicity.
We note that
Q(t) =

D −t/τe.
τ(41)

By considering large t, we assume the following relations:t  τ,

τ

zc
.
2β Dc

(42)

6. CONCLUSION
By giving up the reconstruction of the detailed shape of thetarget, we can estimate the position of the target by reconstructing only several unknown parameters. Even for these several
parameters, the LM method is not stable or converges slowly
unless good initial guesses are used [41,42]. Thus, the proposed
procedure of narrowing target domains as 0 → acubic→ acuboid00
is important for the iterative method to work.
In general, the fluorescence lifetime is also an unknown
parameter. By the analysis in Section 5, we found that the lifetime can be estimated by observing the long-time behaviorof Um (rd , t; rs ) in the case of a point fluorophore. Kumar
et al. obtained the formula (44) for an infinite homogeneous
medium, and its validity was checked by Monte Carlo simulation [39]. This has been developed as the asymptotic time
domain approach, in which resolution below the point spreadfunction can be obtained [40,43]. In this paper, the formula(44) was derived in the presence of the boundary.
In Ref. [44], the domain is assumed to contain simply connected subdomains with constant optical parameters. In thispaper, the cuboid target with constant optical parameters freely
moves in the domain during the search process. The computational advantage of assuming a cuboid is that the solution Um is
obtained in a simple expression, as shown in (13).
The proposed algorithm was developed in the half-space. The
ideas of assuming a simple shape for the target and of decreasing
the number of unknown parameters step by step can be applied
to more general cases where the diffusion equation must be
numerically solved by the finite difference method or finite
element method [45]. In addition, the use of the LM method
is not essential; the proposed algorithm also works with other
iterative schemes such as the conjugate gradient method and thedamped Gauss–Newton method.
In this paper, the algorithm was formulated assuming a single
target. The generalization of the method for multiple targets isstraightforward, at least if the number of targets is known.

238

Vol. 37, No. 2 / February 2020 / Journal of the Optical Society of America A

Funding. Japan Society for the Promotion of Science(15H05740, 15K21766, 17H02081, 17K05572, 18K03438,
19K04421); National Natural Science Foundation of China(11971104, 11531005, 11871149, 11971121, 91730304);
Postgraduate Research & Practice Innovation Program of
Jiangsu Province (KYCX18_0051).
Acknowledgment. The authors appreciate the A3 foresight program: Modeling and Computation of Applied Inverse
Problems by the Japan Society for the Promotion of Science(JSPS). C. Sun is supported by the Postgraduate Research &Practice Innovation Program of Jiangsu Province, and NSFC.
G. Nakamura is supported by a Grant-in-Aid for Scientific
Research of JSPS. G. Nishimura is supported by a Grand-in-Aid
for Scientific Research. Y. Jiang is supported by the National
Natural Science Foundation of China. J. Liu is supported
by NSFC. M. Machida is supported by a Grant-in-Aid for
Scientific Research of JSPS. He also acknowledges the support
by HUSM Grant-in-Aid.
Disclosures. The authors declare no conflicts of interest.
REFERENCES
1. V. Ntziachristos, C. Bremer, and R. Weissleder, “Fluorescence imaging with near-infrared light: new technological advances that enablein vivo molecular imaging,” Eur. Radiol. 13, 195–208 (2003).
2. V. Ntziachristos, J. Ripoll, L. V. Wang, and R. Weissleder, “Looking
and listening to light: the evolution of whole-body photonic imaging,”
Nat. Biotechnol. 23, 313–320 (2005).
3. J.-C. Baritaux, K. Hassler, and M. Unser, “An efficient numericalmethod for general Lp regularization in fluorescence molecular
tomography,” IEEE Trans. Med. Imag. 29, 1075–1087 (2010).
4. D. Han, J. Tian, S. Zhu, J. Feng, C. Qin, B. Zhang, and X. Yang, “Afast reconstruction algorithm for fluorescence molecular tomography
with sparsity regularization,” Opt. Express 18, 8630–8646 (2010).
5. D. Han, X. Yang, K. Liu, C. Qin, B. Zhang, X. Ma, and J. Tian, “Efficientreconstruction method for L1 regularization in fluorescence molecular tomography,” Appl. Opt. 49, 6930–6937 (2010).
6. P. Mohajerani, A. A. Eftekhar, J. Huang, and A. Adibi, “Optimal sparse
solution for fluorescent diffuse optical tomography: theory and phantom experimental results,” Appl. Opt. 46, 1679–1685 (2007).
7. A. Behrooz, H.-M. Zhou, A. A. Eftekhar, and A. Adibi, “Total variation
regularization for 3D reconstruction in fluorescence tomography:
experimental phantom studies,” Appl. Opt. 51, 8216–8227 (2012).
8. J. F. P.-J. Abascal, J. Chamorro-Servent, J. Aguirre, S. Arridge, T.
Correia, J. Ripoll, J. J. Vaquero, and M. Desco, “Fluorescence diffuse
optical tomography using the split Bregman method,” Med. Phys. 38,
6275–6284 (2011).
9. V. Ntziachristos and R. Weissleder, “Experimental three-dimensional
fluorescence reconstruction of diffuse media by use of a normalized
Born approximation,” Opt. Lett. 26, 893–895 (2001).
10. A. Koenig, L. Hervé, V. Josserand, M. Berger, J. Boutet, A. Da Silva,
J.-M. Dinten, P. Peltié, J.-L. Coll, and P. Rizo, “In vivo mice lung tumor
follow-up with fluorescence diffuse optical tomography,” J. Biomed.
Opt. 13, 011008 (2008).
11. V. Ntziachristos, C. H. Tung, C. Bremer, and R. Weissleder,“Fluorescence molecular tomography resolves protease activity
in vivo,” Nat. Med. 8, 757–761 (2002).
12. A. Corlu, R. Choe, T. Durduran, M. A. Rosen, M. Schweiger, S. R.
Arridge, M. D. Schnall, and A. G. Yodh, “Three-dimensional invivo fluorescence diffuse optical tomography of breast cancer inhumans,” Opt. Express 15, 6696–6716 (2007).
13. S. J. Erickson, S. L. Martinez, J. DeCerce, A. Romero, L. Caldera,and A. Godavarty, “Three-dimensional fluorescence tomography of

Research Article

human breast tissues in vivo using a hand-held optical imager,” Phys.
Med. Biol. 58, 1563–1579 (2013).
14. S. van de Ven, A. Wiethoff, T. Nielsen, B. Brendel, M. van der Voort, R.
Nachabe, M. Van der Mark, M. Van Beek, L. Bakker, L. Fels, S. Elias,
P. Luijten, and W. Mali, “A novel fluorescent imaging agent for diffuse
optical tomography of the breast: first clinical experience in patients,”
Mol. Imag. Biol. 12, 343–348 (2010).
15. L. Hervé, A. Koenig, A. Da Silva, M. Berger, J. Boutet, J. M. Dinten, P.
Peltié, and P. Rizo, “Noncontact fluorescence diffuse optical tomography of heterogeneous media,” Appl. Opt. 46, 4896–4906 (2007).
16. V. Ntziachristos and R. Weissleder, “Charge-coupled-device based
scanner for tomography of fluorescent near-infrared probes in turbid
media,” Med. Phys. 29, 803–809 (2002).
17. F. Stuker, J. Ripoll, and M. Rudin, “Fluorescence molecular tomography: principles and potential for pharmaceutical research,”
Pharmaceutics 3, 229–274 (2011).
18. C. Li, G. S. Mitchell, J. Dutta, S. Ahn, R. M. Leahy, and S. R. Cherry,“A three-dimensional multispectral fluorescence optical tomography
imaging system for small animals based on a conical mirror design,”
Opt. Express 17, 7571–7585 (2009).
19. F. Gao, H. Zhao, L. Zhang, Y. Tanikawa, A. Marjono, and Y. Yamada,“A self-normalized, full time-resolved method for fluorescence diffuse
optical tomography,” Opt. Express 16, 13104–13121 (2008).
20. N. Ducros, C. D’Andrea, A. Bassi, and F. Peyrin, “Fluorescence diffuse optical tomography: time-resolved versus continuous-wave inthe reflectance configuration,” IRBM 32, 243–250 (2011).
21. G. Y. Panasyuk, Z.-M. Wang, J. C. Schotland, and V. A. Markel,“Fluorescent optical tomography with large data sets,” Opt. Lett. 33,
1744–1746 (2008).
22. M. Mozumder, T. Tarvainen, S. Arridge, J. P. Kaipio, C. D’Andrea, andV. Kolehmainen, “Approximate marginalization of absorption andscattering in fluorescence diffuse optical tomography,” Inverse Probl.
Imaging 10, 227–246 (2016).
23. S. R. Arridge, “Optical tomography in medical imaging,” Inverse
Probl. 15, R41–R93 (1999).
24. K. Levenberg, “A method for the solution of certain non-linear problems in least squares,” Quart. Appl. Math. 2, 164–168 (1944).
25. D. W. Marquardt, “An algorithm for least-squares estimation of nonlinear parameters,” J. Soc. Ind. Appl. Math. 11, 431–441 (1963).
26. T. Goldstein and S. Osher, “The split Bregman method forL1-regularized problems,” SIAM J. Imaging Sci. 2, 323–343 (2009).
27. J. F. Cai, S. Osher, and Z. Shen, “Split Bregman methods and frame
based image restoration,” SIAM J. Multiscale Model. Simul. 8, 337–
369 (2009).
28. A. Marjono, A. Yano, S. Okawa, F. Gao, and Y. Yamada, “Total lightapproach of time-domain fluorescence diffuse optical tomography,”Opt. Express 16, 15268–15285 (2008).
29. S. Okawa, A. Yano, K. Uchida, Y. Mitsui, M. Yoshida, M. Takekoshi,
A. Marjono, F. Gao, Y. Hoshi, I. Kida, K. Masamoto, and Y. Yamada,“Phantom and mouse experiments of time-domain fluorescence
tomography using total light approach,” Biomed. Opt. Express 4,
635–651 (2013).
30. A. Ishimaru, Wave Propagation and Scattering in Random Media(Academic, 1978), Vol. 1.
31. H. S. Carslaw and J. C. Jaeger, Conduction of Heat in Solids (Oxford
University, 1959).
32. K. Yosida and S. Ito, Functional Analysis and Differential Equations(Iwanami, 1976).
33. A. H. Hielscher, S. L. Jacques, L. Wang, and F. K. Tittel, “The influence of boundary conditions on the accuracy of diffusion theory intime-resolved reflectance spectroscopy of biological tissues,” Phys.
Med. Biol. 40, 1957–1975 (1995).
34. M. Machida and G. Nakamura, “Born series for the photon diffusion
equation perturbing the Robin boundary condition, to appear in J.
Math. Phys.
35. R. Fletcher, “A modified Marquardt subroutine for nonlinear leastsquares,” Technical Report AERE-R 6799 (1971).
36. J. J. Moré, “The Levenberg-Marquardt algorithm: implementationand theory,” in Numerical Analysis, G. A. Watson, ed., Vol. 630 ofLecture Notes in Mathematics (Springer, 1977).

Research Article

Vol. 37, No. 2 / February 2020 / Journal of the Optical Society of America A

37. G. Nishimura, “Contrast improvement in indocyanine green fluorescence sensing in thick tissue using a time-gating method,” Biomed.
Opt. Express 10, 1234–1249 (2019).
38. G. Nishimura, K. Awasthi, and D. Furukawa, “Fluorescence lifetime
measurements in heterogeneous scattering medium,” J. Biomed.
Opt. 21, 075013 (2016).
39. A. T. N. Kumar, J. Skoch, B. J. Bacskai, D. A. Boas, and A. K. Dunn,“Fluorescence-lifetime-based tomography for turbid media,” Opt.
Lett. 30, 3347–3349 (2005).
40. A. T. N. Kumar, S. B. Raymond, G. Boverman, D. A. Boas, and B. J.
Bacskai, “Time resolved fluorescence tomography of turbid media
based on lifetime contrast,” Opt. Express 14, 12255–12270 (2006).
41. M. Hanke, “A regularizing Levenberg-Marquardt scheme, with applications to inverse groundwater filtration problems,” Inverse Probl. 13,
79–96 (1997).

239

42. Y. Jiang and G. Nakamura, “Convergence of Levenberg-Marquardt
method for the inverse problem with an interior measurement,”
J. Inverse Ill-Posed Probl. 27, 195–215 (2018).
43. W. L. Rice, S. Hou, and A. T. N. Kumar, “Resolution below the point
spread function for diffuse optical imaging using fluorescence lifetime
multiplexing,” Opt. Lett. 38, 2038–2040 (2013).
44. A. Zacharopoulos, M. Schweiger, V. Kolehmainen, and S. Arridge,“3D shape based reconstruction of experimental data in diffuse
optical tomography,” Opt. Express 17, 18940–18956 (2009).
45. Q. Zhu, H. Dehghani, K. M. Tichauer, R. W. Holt, K. Vishwanath, F.
Leblond, and B. W. Pogue, “A three-dimensional finite element modeland image reconstruction algorithm for time-domain fluorescenceimaging in highly scattering media,” Phys. Med. Biol. 56, 7419–7434(2011).