PROCEEDINGS OF SPIE
SPIEDigitalLibrary.org/conference-proceedings-of-spie
Image reconstruction in
photoacoustic tomography with
heterogeneous media using an
iterative method
Chao Huang, Kun Wang, Lihong V. Wang, Mark A.
Anastasio
Chao Huang, Kun Wang, Lihong V. Wang, Mark A. Anastasio, "Image
reconstruction in photoacoustic tomography with heterogeneous media using
an iterative method," Proc. SPIE 8581, Photons Plus Ultrasound: Imaging and
Sensing 2013, 858106 (4 March 2013); doi: 10.1117/12.2007072
Event: SPIE BiOS, 2013, San Francisco, California, United States
Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 9/5/2018 Terms of Use: https://www.spiedigitallibrary.org/terms-of-use
Image reconstruction in photoacoustic tomography with
h
eterogeneous media using an iterative method
Chao Huang, Kun Wang, Lihong V. Wang, and Mark A. Anastasio
Department of Biomedical Engineering, Washington University in St. Louis,
St. Louis, MO 63130
ABSTRACT
There remains an urgent need to develop effective photoacoustic computed tomography (PACT) image recon-
struction methods for use with acoustically inhomogeneous media. Transcranial PACT brain imaging is an im-
portant example of an emerging imaging application that would benefit greatly from this. Existing approaches
to PACT image reconstruction in acoustically heterogeneous media are limited to weakly varying media, are
computationally burdensome, and/or make impractical assumptions regarding the measurement geometry. In
this work, we develop and investigate a full-wave approach to iterative image reconstruction in PACT for media
possessing inhomogeneous speed-of-sound and mass density distributions. A key contribution of the work is the
formulation of a procedure to implement a matched discrete forward and backprojection operator pair, which
facilitates the application of a wide range of modern iterative image reconstruction algorithms. This presents
the opportunity to employ application-specific regularization methods to mitigate image artifacts due to mea-
surement data incompleteness and noise. Our results establish that the proposed image reconstruction method
can effectively compensate for acoustic aberration and reduces artifacts in the reconstructed image.
Keywords:
Photoacoustic tomography, optoacoustic tomography, thermoacoustic tomography, iterative image
reconstruction, acoustic heterogeneity, k-space pseudospectral method
1. INTRODUCTION
Conventional image reconstruction algorithms in photoacoustic computed tomography (PACT) are based on the
assumption that the acoustic mdeium is homogeneous.
1–5
However, this assumption is not warranted in many
biomedical applications of PACT. For example, in transcranial PACT of primates, the presence of the skull will
introduce strong acoustic heterogeneities, which could severely distort the PA wavefields.
6,7
In these applications
of PACT, if the reconstruction algorithm does not account for the acoustic heterogeneities, the reconstructed
images can contain significant artifacts.
Several image reconstruction methods have been developed to compensate for the effects of smooth variations
in a medium’s speed-of-sound (SOS) distribution.
8–10
These methods are based on a geometrical acoustics (GA)
approximation, which utilizes the Eikonal equation to model acoustic wave front propagation. However, the GA
approximation is based on the assumption that the length scale of the SOS variation is much greater than the
acoustic wavelength, which can be violated when the media process strong acoustic heterogeneities.
Full-wave approaches to PACT image reconstruction have also been proposed.
11–13
While these methods
are based on solutions to the exact PA wave equation, which permits a broader domain of applicability, they
also possess certain practical limitations. For example, finite element methods (FEMs) have intensive com-
putational burden, which is especially problematic for three-dimensional (3D) applications of PACT. Although
time-reversal (TR) methods are mathematically exact in their continuous forms in 3D homogeneous media,
12
they
are predicated on the assumption that the measurement surface encloses the object, which is often impractical
in biomedical applications of PACT.
In this work, we develop and investigate a full-wave approach to iterative image reconstruction in PACT with
media possessing heterogeneous SOS and mass density distributions. We construct a discrete imaging model for
(Send correspondence to Mark A. Anastasio)
M
ark A. Anastasio: E-mail: anastasio@seas.wustl.edu, Telephone: 1 314 935 7208
Photons Plus Ultrasound: Imaging and Sensing 2013, edited by Alexander A. Oraevsky, Lihong V. Wang,
Proc. of SPIE Vol. 8581, 858106 · © 2013 SPIE · CCC code: 1605-7422/13/$18 · doi: 10.1117/12.2007072
Proc. of SPIE Vol. 8581 858106-1
Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 9/5/2018
Terms of Use: https://www.spiedigitallibrary.org/terms-of-use
PACT that is based on the exact photoacoustic PA wave equation
, and establish a procedure to implement a
matched forward and backprojection operator pair associated with the discrete imaging model. This permits
application of a wide-range of iterative image reconstruction algorithms that can mitigate the effects of data
incompleteness and noise. The forward and backprojection operators are based on the k-space pseudospectral
method,
14
which possesses significant computational advantages over finite-difference and finite-element methods.
The developed reconstruction methodology is investigated by use of computer simulation studies.
2. BACKGROUND
In a heterogeneous lossless fluid medium, the propagation of the thermoacoustically-induced pressure wavefield
p
(
r
, t
) can be modeled by the following three coupled equations
13,15
∂
∂t
u
(
r
,
t
) =
−
1
ρ
0
(
r
)
∇
p
(
r
,
t
)
,
(1)
∂
∂t
ρ
(
r
,
t
) =
−
ρ
0
(
r
)
∇
u
(
r
, t
)
,
(2)
p
(
r
, t
) =
c
0
(
r
)
2
ρ
(
r
, t
)
,
(3)
subject to the initial conditions:
p
0
(
r
)
≡
p
(
r
, t
)
|
t
=0
= Γ(
r
)
A
(
r
)
,
u
(
r
, t
)
|
t
=0
= 0
,
(4)
where
u
(
r
, t
) is the acoustic particle velocity,
c
0
(
r
) is the medium’s SOS distribution,
A
(
r
) is the absorbed optical
energy density distribution, Γ(
r
) is the Grueneisen parameter, and
ρ
(
r
, t
) and
ρ
0
(
r
) denote the distributions of
the medium’s acoustic and ambient densities, respectively.
To compute the numerical solutions of the coupled equations (1), (2), (3), the k-space pseudospectral method
can be utilized to discretize the coupled equations, and the discrete form is
13
u
i
m
+1
=
u
i
m
+
Φ
i
p
m
,
(5)
ρ
i
m
+1
=
ρ
i
m
+
Ψ
i
u
i
m
+1
,
(6)
p
m
+1
=
C
3
X
i
=1
ρ
i
m
+1
,
(7)
where
p
m
≡
(
p
(
r
1
, m
∆
t
)
,
, p
(
r
N
, m
∆
t
))
T
lexicographically ordered vector representation of sampled values of
p
(
r
, t
) at time step
m
∆
t
(
m
= 1
,
, M
) and at vertices of a 3D Cartesian grid
r
n
(
n
= 1
,
, N
).
M
is the
total number of time steps, and
N
≡
N
1
N
2
N
3
is the toal number of vertices, where
N
i
the number of vertices
along the
i
-th dimension (
i
= 1
,
2
,
3).
u
i
m
,
ρ
i
m
are defined in similar manner as
p
m
for
u
(
r
, t
) and
ρ
(
r
, t
) in the
i
-th dimension, respectively. The diagonal matrix
C
≡
diag(
c
2
0
(
r
1
)
,
, c
2
0
(
r
N
)) represents the sampled vaules of
squared SOS distribution
c
0
(
r
)
2
. The matrices
Φ
i
and
Ψ
i
are defined as
Φ
i
p
m
≡−
∆
t
Q
−
1
F
−
1
{
j
K
i
◦
κ
◦
F
{
p
m
}}
,
(8)
Ψ
i
u
i
m
+1
≡−
∆
t
QF
−
1
{
j
K
i
◦
κ
◦
F
{
u
i
m
+1
}}
,
(9)
where
j
is the imaginary unit,
Q
≡
diag(
ρ
0
(
r
1
)
,
, ρ
0
(
r
N
)) is the sampled vaules of the ambient density
ρ
0
(
r
),
◦
is the Hadamard product,
F
is the operator that first converts a
N
×
1 lexicographically ordered vector into
a
N
1
×
N
2
×
N
3
3D matrix then performs the 3D forward discrete Fourier transform (DFT), and
F
−
1
is the
operator that first carries out the 3D inverse DFT then convert the resulting
N
1
×
N
2
×
N
3
3D matrix into a
N
×
1 lexicographically ordered vector.
K
i
is a 3D matrix with elements given by
K
1
n
1
n
2
n
3
= 2
π
n
1
−
1
L
1
,
K
2
n
1
n
2
n
3
=
2
π
n
2
−
1
L
2
,
K
3
n
1
n
2
n
3
=
2
π
n
3
−
1
L
3
,
(
10)
where
n
i
= 1
,
, N
i
(
i
= 1
,
2
,
3), and
L
i
denotes the length of the spatial grid in the
i
-th dimension.
The 3D matrix
κ
= sinc(
1
2
∆
t
c
min
K
) is the k-space operator, where sinc(
x
) =
sin(
x
)
x
,
c
m
in
is the minimum
of
c
0
(
r
),
K
is a 3D matrix defined as
K
≡
q
P
3
i
=
1
K
i
◦
K
i
, and the sinc function and square root function are
both element-wise operations.
Proc. of SPIE Vol. 8581 858106-2
Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 9/5/2018
Terms of Use: https://www.spiedigitallibrary.org/terms-of-use
3. DISCRETE IMAGING MODEL AND ITERATIVE IMAGE RECONSTRUCTIO
N
The k-space pseudospectral method for numerically solving the photoacoustic wave equation described in Section
2 can be employed to establish a discrete PACT imaging model. The resulting discrete imaging operator
H
, also
known as the system matrix, and its adjoint or backprojection operator
H
†
can be employed with an iterative
algorithm to accomplish the PACT image reconstruction task.
3.1 Discrete imaging model
Equations (5), (6) and (7) can be described by a single matrix equation to determine the updated wavefield
variables after a time step ∆
t
as
v
m
+1
=
Wv
m
,
(11)
where
v
m
= (
u
1
m
,
u
2
m
,
u
3
m
,
ρ
1
m
,
ρ
2
m
,
ρ
3
m
,
p
m
)
T
is a 7
N
×
1 vector containing all the wavefield variables at the time
step
m
∆
t
. The 7
N
×
7
N
propagator matrix
W
is defined as
W
≡
I
N
×
N
0
N
×
N
0
N
×
N
0
N
×
N
0
N
×
N
0
N
×
N
Φ
1
0
N
×
N
I
N
×
N
0
N
×
N
0
N
×
N
0
N
×
N
0
N
×
N
Φ
2
0
N
×
N
0
N
×
N
I
N
×
N
0
N
×
N
0
N
×
N
0
N
×
N
Φ
3
Ψ
1
0
N
×
N
0
N
×
N
I
N
×
N
0
N
×
N
0
N
×
N
Ψ
1
Φ
1
0
N
×
N
Ψ
2
0
N
×
N
0
N
×
N
I
N
×
N
0
N
×
N
Ψ
2
Φ
2
0
N
×
N
0
N
×
N
Ψ
3
0
N
×
N
0
N
×
N
I
N
×
N
Ψ
3
Φ
3
CΨ
1
CΨ
2
CΨ
3
C
C
C
D
,
(12)
where
D
≡
C
P
3
i
=1
Ψ
i
Φ
i
,
I
N
×
N
is the
N
×
N
identity matrix, and
0
N
×
N
is the
N
×
N
zero matrix.
The wavefield quantities can be propagated forward in time from
t
= 0 to
t
= (
M
−
1)∆
t
as
v
0
v
1
.
.
.
v
M
−
1
=
T
M
−
1
T
1
v
0
0
7
N
×
1
.
.
.
0
7
N
×
1
,
(13)
where the 7
N M
×
7
N M
matrices
T
m
(
m
= 1
,
, M
−
1) are defined in terms of
W
as
T
m
=
I
7
N
×
7
N
0
7
N
×
7
N
.
.
.
.
.
.
.
.
.
0
7
N
×
7
N
I
7
N
×
7
N
0
7
N
×
7
N
W
0
(
m
+1)
·
7
N
×
(
M
−
m
)
·
7
N
0
(
M
−
m
−
1)
·
7
N
×
m
·
7
N
0
(
M
−
m
−
1)
·
7
N
×
(
M
−
m
)
·
7
N
,
(14)
with
W
residing between the (7
N
(
m
−
1)+1)-th to 7
N m
-th rows and the (7
N m
+1)-th to 7
N
(
m
+1)-th columns
of
T
m
.
From the equation of state in Eqn. (3) and initial conditions Eqn. (4), the vector (
v
0
,
0
,
,
0
)
T
can be
computed from the initial pressure distribution
p
0
as
v
0
0
7
N
×
1
.
.
.
0
7
N
×
1
=
T
0
p
0
,
(15)
where
T
0
≡
(
τ
,
0
7
N
×
N
,
,
0
7
N
×
N
)
T
,
(16)
τ
≡
(
0
N
×
N
,
0
N
×
N
,
0
N
×
N
,
1
3
C
−
1
,
1
3
C
−
1
,
1
3
C
−
1
,
I
N
×
N
)
T
,
(
17)
Proc. of SPIE Vol. 8581 858106-3
Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 9/5/2018
Terms of Use: https://www.spiedigitallibrary.org/terms-of-use
In general, the transducer locations
r
d
l
a
t which the PA data
ˆ
p
are measured do not coincide with the vertices
of the Cartesian grid at which the values of the propagated field quantities are computed. The measured PA
data
ˆ
p
can be related to the computed field quantities via an interpolation operation as
ˆp
=
S
v
0
v
1
.
.
.
v
M
−
1
,
(18)
where
ˆ
p
≡
ˆ
p
0
ˆ
p
1
.
.
.
ˆ
p
M
−
1
,
ˆ
p
m
≡
p
(
r
d
1
, m
∆
t
)
.
.
.
p
(
r
d
L
, m
∆
t
)
,
(19)
and
S
≡
Θ 0
L
×
7
N
0
L
×
7
N
0
L
×
7
N
Θ
0
L
×
7
N
.
.
.
.
.
.
.
.
.
.
.
.
0
L
×
7
N
0
L
×
7
N
Θ
.
(20)
Here,
Θ
≡
[
s
1
,
,
s
L
]
T
, where
s
l
(
l
= 1
,
, L
) is a 1
×
7
N
row vector in which all elements are zeros except
the 4 corresponding to acoustic pressure at 4 grid nodes
r
l,
1
,
r
l,
2
,
r
l,
3
,
r
l,
4
that are nearest to the transducer
location
r
d
l
. In other words, these 4 entries are interpolation coefficients to compute the acoustic pressure at the
l
-th transducer, and their values are given by the barycentric coordinates of
r
d
l
with respect to
r
l,
1
,
r
l,
2
,
r
l,
3
,
r
l,
4
,
which are determined by Delaunay triangulation.
16
By use of Eqns. (13), (15), and (18), we obtains the discrete PACT imaging model
ˆ
p
=
ST
M
−
1
T
1
T
0
p
0
=
Hp
0
,
(21)
where the explicit form of the system matrix is identified as
H
≡
ST
M
−
1
T
1
T
0
.
3.2 Recontruction method
Upon the establishment of the discrete PACT imaging model, the image reconstruction task, which is to estimate
the initial pressure distribution
p
0
from knowledge of the measured data
ˆ
p
, can be accomplished by iteratively
solving Eqn. (21) with an appropriate regularization. This can be achieved by seeking solutions of an optimization
problem, which, in this work, was designed as a penalized least squares (PLS) cost function that contained a
total variation (TV) penalty term
ˆ
p
0
= arg min
p
0
≥
0
k
ˆp
−
Hp
0
k
2
+
λ
|
p
0
|
TV
,
(22)
where
λ
is the regularization parameter, and a non-negativity constraint was employed. For the 3D case, the
TV-norm is defined as
|
p
0
|
TV
=
N
X
n
=1
q
([
p
0
]
n
−
[
p
0
]
n
−
1
)
2
+
([
p
0
]
n
−
[
p
0
]
n
−
2
)
2
+ ([
p
0
]
n
−
[
p
0
]
n
−
3
)
2
,
(23)
where [
p
0
]
n
denotes the
n
-th grid node, and [
p
0
]
n
−
1
, [
p
0
]
n
−
2
,
[
p
0
]
n
−
3
are neighboring nodes before the
n
-th node
along the first, second and third dimension, respectively. The fast iterative shrinkage/thresholding algorithm
(FISTA)
17,18
was employed to solve Eqn. (22).
Proc. of SPIE Vol. 8581 858106-4
Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 9/5/2018
Terms of Use: https://www.spiedigitallibrary.org/terms-of-use
When iterative algorithms, including FISTA, are employed to
minimize a PLS cost function, the action of
the operators
H
and its adjoint
H
†
must be computed.
19
In this work, the action of
H
was computed by use of
the MATLAB k-Wave toolbox.
20
It can be verified the action of
H
†
can be computed as
v
M
−
1
=
Θ
T
ˆ
p
M
−
1
,
(24)
v
m
−
1
=
Θ
T
ˆ
p
m
−
1
+
W
T
v
m
, m
=
M
−
1
,
,
1
(25)
p
bp
=
τ
T
v
0
.
(26)
4. NUMERICAL STUDIES AND RESULTS
Numerical studies were conducted to assess the effectiveness and robustness of the proposed discrete imaging
model in studies of iterative image reconstruction from incomplete data sets in two-dimensional (2D) PACT.
Note that, although the above discrete imaging model is based on the 3D PA wave equation, the 2D formulation
is contained as a special case. The performance of the iterative reconstruction method was compared to an
existing TR reconstruction method.
13
4.1 Numerical studies
The low contrast disc phantom shown in Fig. 1-(a) was employed to investigate the robustness of the iterative
reconstruction methods with respect to noise in the measurement data and different types of data incompleteness
resulting from three differenet scanning geometries. A ‘full-view’ scanning geometry utilized 180 transducers that
were evenly distributed on a circle of radius 40 mm. A ‘few-view’ scanning geometry utilized 60 transducers that
were equally distributed on the circle. Finally, a ‘limited-view’ scanning geometry utilized 90 transducers that
were evenly located on a semi-circle of radius 40 mm. At each transducer location, a total of 20,000 temporal
samples of simulated pressure data were computed at time step ∆
t
= 30 ns, all of which were employed by the
TR image reconstruction method. However, only the first 1,500 temporal samples were employed by the iterative
reconstruction method. The simulated data were then contaminated by 3% (with respect to maximum value of
noiseless data) additive white Gaussian noise (AWGN). Figure 1-(b) and (c) display the SOS and density maps
used in the simulation studies, which were deduced from the X-ray CT data of a monkey skull.
7
All simulation
studies were computed on a Cartesian grid of 512
×
512 pixels with a pitch of 0.2 mm, and the regularization
parameter
λ
was empirically selected to have a value of 0.001.
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
(a)
Speed of sound (m/s)
1600
1800
2000
2200
2400
2600
2800
(b)
Mass density (kg/m
3
)
1200
1400
1600
1800
2000
(c)
F
igure 1. Panel (a) is the disc numerical phantoms employed to represent
p
0
in the simulation studies. Panel (b) and (c)
are SOS and density maps deduced from the X-ray CT data of a monkey skull.
4.2 Numerical results
The reconstructed images corresponding to the three scanning geometries are displayed in Figs. 2 - 4. In each
figure, the results in the top row correspond to use of the TR reconstruction method, while the bottom row
shows the corresponding results obtained by use of the iterative method. The profiles shown in each figure
are along the ‘Y’-axis indicated by the arrow in Fig. 2-(a). The red lines and blue lines corresdpond to profiles
through the phantom and reconstructed images, respectively. With all of three scanning geometries, the iterative
results contain lower noise levels than the TR results, which suggests that the iterative method is more robust
Proc. of SPIE Vol. 8581 858106-5
Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 9/5/2018
Terms of Use: https://www.spiedigitallibrary.org/terms-of-use
to the noise in the measurement. With the few-view and the limi
ted-view scanning geometries, the images
reconstructed from the iterative method contain fewer distortions and artifacts than the TR results. Also, the
values of the images reconstructed from the iterative method are much closer to the values of the phantom than
those produced by the TR method. These results indicate that the iterative method is more robust to the data
incompleteness.
Y
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
(a)
0
10
20
30
0
0.5
1
Position Y (mm)
p
0
Phantom
TR
(b)
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
(c)
0
10
20
30
0
0.5
1
Position Y (mm)
p
0
Phantom
Iterative
(d)
F
igure 2. Panel (a) and (c) are reconstructed images corresponding to the full-view scanning geometry by use of the TR
method and iterative method, respectively. Panel (b) and (d) are the corresponding profiles.
5. SUMMARY
To compensate for the effects of acoustic heterogeneities in PACT, we proposed and investigated an iterative
image reconstruction method that was based on the exact PA wave equation. We established a discrete imaging
model based on the k-space pseudospectral method, and provided a procedure of implementing the forward
and backprojection operators associated with the discrete imaging model. The matched operator pair was
employed in an iterative image reconstruction algorithm that sought to minimize a TV-regularized PLS cost
function. The developed reconstruction methodology was investigated by computer simulation studies, and the
results demonstrated that the reconstruction methodology can effectively reduce noise level and mitigate image
artifacts due to data incompleteness.
ACKNOWLEDGMENTS
This work was supported in part by NIH awards EB010049 and CA167446.
REFERENCES
1. Finch, D., Patch, S., and Rakesh, “Determining a function from its mean values over a family of spheres,”
SIAM Journal of Mathematical Analysis
35
, 1213–1240 (2004).
Proc. of SPIE Vol. 8581 858106-6
Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 9/5/2018
Terms of Use: https://www.spiedigitallibrary.org/terms-of-use
0
0.5
1
1.5
2
(a)
0
10
20
30
0
0.5
1
Position Y (mm)
p
0
Phantom
TR
(b)
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
(c)
0
10
20
30
0
0.5
1
Position Y (mm)
p
0
Phantom
Iterative
(d)
F
igure 3. Panel (a) and (c) are reconstructed images corresponding to the few-view scanning geometry by use of the TR
method and iterative method, respectively. Panel (b) and (d) are the corresponding profiles.
2. Xu, Y. and Wang, L. V., “Universal back-projection algorithm for photoacoustic computed tomography,”
Physical Review E
71
(016706) (2005).
3. Kunyansky, L. A., “Explicit inversion formulae for the spherical mean radon transform,”
Inverse Prob-
lems
23
, 373–383 (2007).
4. Wang, K., Ermilov, S., Su, R., Brecht, H.-P., Oraevsky, A., and Anastasio, M., “An imaging model incorpo-
rating ultrasonic transducer properties for three-dimensional optoacoustic tomography,”
Medical Imaging,
IEEE Transactions on
30
, 203 –214 (feb. 2011).
5. Wang, K. and Anastasio, M. A., “A simple fourier transform-based reconstruction formula for photoacous-
tic computed tomography with a circular or spherical measurement geometry,”
Physics in Medicine and
Biology
57
(23), N493 (2012).
6. Jin, X., Li, C., and Wang, L. V., “Effects of acoustic heterogeneities on transcranial brain imaging with
microwave-induced thermoacoustic tomography,”
Medical Physics
35
(7), 3205–3214 (2008).
7. Huang, C., Nie, L., Schoonover, R. W., Guo, Z., Schirra, C. O., Anastasio, M. A., and Wang, L. V.,
“Aberration correction for transcranial photoacoustic tomography of primates employing adjunct image
data,”
Journal of Biomedical Optics
17
(6), 066016 (2012).
8. Xu, Y. and Wang, L., “Effects of acoustic heterogeneity in breast thermoacoustic tomography,”
Ultrasonics,
Ferroelectrics and Frequency Control, IEEE Transactions on
50
, 1134 –1146 (sept. 2003).
9. Modgil, D., Anastasio, M. A., and La Rivire, P. J., “Image reconstruction in photoacoustic tomography with
variable speed of sound using a higher-order geometrical acoustics approximation,”
Journal of Biomedical
Optics
15
(2), 021308–021308–9 (2010).
10. Jose, J., Willemink, R. G. H., Resink, S., Piras, D., van Hespen, J. C. G., Slump, C. H., Steenbergen, W.,
van Leeuwen, T. G., and Manohar, S., “Passive element enriched photoacoustic computed tomography (per
Proc. of SPIE Vol. 8581 858106-7
Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 9/5/2018
Terms of Use: https://www.spiedigitallibrary.org/terms-of-use
−0.2
0
0.2
0.4
0.6
0.8
1
1.2
(a)
0
10
20
30
0
0.5
1
Position Y (mm)
p
0
Phantom
TR
(b)
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
(c)
0
10
20
30
0
0.5
1
Position Y (mm)
p
0
Phantom
Iterative
(d)
F
igure 4. Panel (a) and (c) are reconstructed images corresponding to the limited-view scanning geometry by use of the
TR method and iterative method, respectively. Panel (b) and (d) are the corresponding profiles.
pact) for simultaneous imaging of acoustic propagation properties and light absorption,”
Opt. Express
19
,
2093–2104 (Jan 2011).
11. Yao, L. and Jiang, H., “Enhancing finite element-based photoacoustic tomography using total variation
minimization,”
Appl. Opt.
50
, 5031–5041 (Sep 2011).
12. Hristova, Y., Kuchment, P., and Nguyen, L., “Reconstruction and time reversal in thermoacoustic tomog-
raphy in acoustically homogeneous and inhomogeneous media,”
Inverse Problems
24
(5), 055006 (2008).
13. Treeby, B. E., Zhang, E. Z., and Cox, B. T., “Photoacoustic tomography in absorbing acoustic media using
time reversal,”
Inverse Problems
26
(11), 115003 (2010).
14. Cox, B. T., Kara, S., Arridge, S. R., and Beard, P. C., “k-space propagation models for acoustically
heterogeneous media: Application to biomedical photoacoustics,”
The Journal of the Acoustical Society of
America
121
(6), 3453–3464 (2007).
15. Morse, P. M., “Theoretical acoustics,” in [
Theoretical Acoustics
], Princeton University Press (1987).
16. Lee, D. T. and Schachter, B. J., “Two algorithms for constructing a delaunay triangulation,”
International
Journal of Parallel Programming
9
, 219–242 (1980). 10.1007/BF00977785.
17. Beck, A. and Teboulle, M., “Fast gradient-based algorithms for constrained total variation image denoising
and deblurring problems,”
Image Processing, IEEE Transactions on
18
, 2419 –2434 (nov. 2009).
18. Wang, K., Su, R., Oraevsky, A. A., and Anastasio, M. A., “Investigation of iterative image reconstruction
in three-dimensional optoacoustic tomography,”
Physics in Medicine and Biology
57
(17), 5399 (2012).
19. Fessler, J. A., “Penalized weighted least-squares reconstruction for positron emission tomography,”
IEEE
Transactions on Medical Imaging
13
, 290–300 (1994).
20. Treeby, B. and Cox, B., “k-wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic
wave fields,”
Journal of Biomedical Optics
15
, 021314 (2010).
Proc. of SPIE Vol. 8581 858106-8
Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 9/5/2018
Terms of Use: https://www.spiedigitallibrary.org/terms-of-use