Conventional photoacoustic computed tomography (PACT) image reconstruction methods assume that the object and surrounding medium are described by a constant speed-of-sound (SOS) value. In order to accurately recover fine structures, SOS heterogeneities should be quantified and compensated for during PACT reconstruction. To address this problem, several groups have proposed hybrid systems that combine PACT with ultrasound computed tomography (USCT). In such systems, a SOS map is reconstructed first via USCT. Consequently, this SOS map is employed to inform the PACT reconstruction method. Additionally, the SOS map can provide structural information regarding tissue, which is complementary to the functional information from the PACT image. We propose a paradigm shift in the way that images are reconstructed in hybrid PACT-USCT imaging. Inspired by our observation that information about the SOS distribution is encoded in PACT measurements, we propose to jointly reconstruct the absorbed optical energy density and SOS distributions from a combined set of USCT and PACT measurements, thereby reducing the two reconstruction problems into one. This innovative approach has several advantages over conventional approaches in which PACT and USCT images are reconstructed independently: (1) Variations in the SOS will automatically be accounted for, optimizing PACT image quality; (2) The reconstructed PACT and USCT images will possess minimal systematic artifacts because errors in the imaging models will be optimally balanced during the joint reconstruction; (3) Due to the exploitation of information regarding the SOS distribution in the full-view PACT data, our approach will permit high-resolution reconstruction of the SOS distribution from sparse array data

Anastasio, Mark A.

Sheng, Qiwei

Wang, Kun

Wang, Lihong V.

Xia, Jun

Zhu, Liren

Caltech Authors - Main

PROCEEDINGS OF SPIE
SPIEDigitalLibrary.org/conference-proceedings-of-spie
Photoacoustic computed tomography
without accurate ultrasonic
transducer responses
Qiwei  Sheng, Kun  Wang, Jun  Xia, Liren  Zhu, Lihong V.
Wang, et al.
Qiwei  Sheng, Kun  Wang, Jun  Xia, Liren  Zhu, Lihong V. Wang, Mark A.
Anastasio, "Photoacoustic computed tomography without accurate ultrasonic
transducer responses," Proc. SPIE 9323, Photons Plus Ultrasound: Imaging
and Sensing 2015, 932313 (11 March 2015); doi: 10.1117/12.2081056
Event: SPIE BiOS, 2015, San Francisco, California, United States
Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 9/4/2018  Terms of Use: https://www.spiedigitallibrary.org/terms-of-use
Photoacoustic computed tomography without accurate
ultrasonic transducer responses
Qiwei Shenga, Kun Wanga, Jun Xiab, Liren Zhua, Lihong V. Wanga, and Mark A. Anastasioa
a Department of Biomedical Engineering, Washington University in St. Louis,
St. Louis, MO 63130
b Department of Biomedical Engineering, The University at Buffalo, Buffalo, NY 14260
ABSTRACT
Photoacoustic computed tomography (PACT) holds great promise for biomedical imaging. When the imag-
ing system employs conventional piezoelectric ultrasonic transducers, the ideal photoacoustic (PA) signals are
degraded by the transducers’ acousto-electric impulse responses (EIRs) during the measurement process. If un-
accounted for, this can degrade the accuracy of the reconstructed image. In principle, the effect of the EIRs on
the measured PA signals can be ameliorated via deconvolution; images can be reconstructed subsequently by
application of a reconstruction method that assumes an idealized EIR. Alternatively, the effect of the EIR can be
incorporated into an imaging model and implicitly compensated for during reconstruction. In either case, the ef-
ficacy of the correction can be limited by errors in the assumed EIRs. In this work, a joint optimization approach
to PACT image reconstruction is proposed for mitigating errors in reconstructed images that are caused by use
of an inaccurate EIR. A variable projection method is employed to refine the measured EIR during the process of
reconstructing the sought-after absorbed optical energy density. Computer-simulation and experimental studies
are conducted to investigate the numerical properties of the method and demonstrate its value for mitigating
image distortions and enhancing the visibility of fine structures.
1. INTRODUCTION
In PACT,1–5 tissues are illuminated by short laser pulses to generate internal broadband photoacoustic
wavefields via the photoacoustic effect.1,2 From the measurements of the propagated acoustic wavefields outside
the object, an image reconstruction method can be employed to estimate the absorbed optical energy density
within the tissue.
A variety of PACT imaging systems have been developed. When the system employs piezoelectric transducers,
the PA signals at the transducer locations are convolved with the transducers’ acousto-electric impulse responses
(EIRs) to generate the output digital signals. If unaccounted for, this degradation of the measurement data will
result in a modulation of the spatial frequency components of the estimated absorbed optical energy density
distribution.6 In principle, if the EIR is known accurately, the effect of the EIRs on the PA signals can be
removed via deconvolution or compensated for implicitly during image reconstruction by use of an imaging
model incorporating the EIR.7–9 However, acquiring an accurate measurement of the EIR can be difficult.10
Furthermore, there are typically hundreds of individual elements with unique EIRs of their own in an transducer
array, and it can be an arduous task to characterize each EIR.
The purpose of this study is to develop a joint optimization approach to simultaneously estimate the absorbed
optical energy density and the EIR of the transducer array. To achieve this, a variable projection method11–13
is employed to refine the measured EIR during the process of reconstruction.
(Send correspondence to:)
Mark Anastasio: E-mail: anastasio@wustl.edu
Photons Plus Ultrasound: Imaging and Sensing 2015, edited by Alexander A. Oraevsky, Lihong V. Wang
Proc. of SPIE Vol. 9323, 932313 · © 2015 SPIE · CCC code: 1605-7422/15/$18
doi: 10.1117/12.2081056
Proc. of SPIE Vol. 9323  932313-1
Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 9/4/2018
Terms of Use: https://www.spiedigitallibrary.org/terms-of-use
2. BACKGROUND
A discrete-to-discrete imaging model1,2, 5, 14,15 for PACT can be established as
u = H(h)θ, (1)
where θ is a finite-dimensional representation of the sought-after object, i.e., absorbed optical energy density
distribution, u denotes the measured data vector, h denotes the vector of EIR samples, and H denotes a system
matrix that maps the absorbed optical energy density distribution to the measured data vector. Generally,
the number of unknowns in h is significantly smaller than the number of unknowns in θ. The notation H(h)
indicates the dependence of H on h. A conventional penalized least squares image reconstruction method seeks
the solution
θˆ = argmin
θ≥0
‖u−H(h)θ‖2 + λR1(θ). (2)
Here R1(θ) denotes the penalty term, whose impact is controlled by the regularization parameter λ.
A standard gradient descent algorithm solves Eqn. (2) iteratively by updating the estimate of θ by taking a
small step along the direction opposite the gradient of the cost function. The step size is commonly determined
by using a line search method. A conventional iterative reconstruction method based upon a gradient descent
algorithm is described in Alg. 1. In Alg. 1, Pc denotes the operator that projects all negative values to 0, γk is
Algorithm 1 Conventional iterative reconstruction algorithm for estimation of θ
1: h← experimentally-measured EIR.
2: k ← 0 {k is the number of algorithm iteration}
3: while stopping criterion is not satisfied do
4: k ← k + 1
5: θk ← Pc
{
θk−1 − γk∇θφ(θk−1,h)
}
6: end while
7: θˆ ← θk.
the step size, ∇θ is the gradient operator with respect to θ, and φ(θ,h) = ‖u−H(h)θ‖2+λR1(θ). Note that the
EIR h is fixed in the conventional iterative reconstruct method, and therefore errors in h could be propagated
into the sought-after object θ.
3. JOINT RECONSTRUCTION
In order to mitigate errors in the reconstructed images, we formulate the image reconstruction problem as a
joint numerical optimization problem
(θˆ, hˆ) = argmin
θ≥0,h
ϕ(θ,h), (3)
where the cost function ϕ(θ,h) is defined as
ϕ(θ,h) ≡ ‖u−H(h)θ‖2 + λR1(θ) + αR2(h). (4)
Here, R1(θ) and R2(h) represent penalty terms, whose impacts are controlled by the regularization parameters
λ and α, respectively.
4. VARIABLE PROJECTION METHOD
A variable projection method is employed to reformulate the minimization problem (3). The variable projec-
tion method exploits the separable structure of the image model (1) by taking advantage of the fact that Eqn. (1)
is bilinear with respect to (θ,h) and h contains relatively few unknowns compared to θ such that a closed form
solution is possible when θ is given. Instead of explicitly separating the variables θ and h as in coordinate
Proc. of SPIE Vol. 9323  932313-2
Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 9/4/2018
Terms of Use: https://www.spiedigitallibrary.org/terms-of-use
descent, the variable projection method implicitly eliminates the variable h to obtain a reduced problem that
depends only on θ
min
θ≥0
ϕ(θ,h(θ)), (5)
where h(θ) is a solution of
min
h
ϕ(θ,h). (6)
To compute the gradient of the reduced cost function for the gradient descent algorithm, it is useful to note that
∇θϕ(θ,h(θ)) = ∇θϕ(θ,h)
∣∣∣
h=h(θ)
, (7)
which simplifies the gradient calculation; the gradient computation prescribed by Eqn. (7) is identical to that
employed by the standard gradient descent methods for penalized least squares image reconstruction method. A
complete description of the variable projection method is given in Alg. 2.
Algorithm 2 Variable projection (VP) algorithm for joint estimation of θ and h
1: h0 ← experimentally-measured EIR.
2: θ0 ← arg minθ≥0 ϕ(θ,h0)
3: k ← 0 {k is the number of algorithm iteration}
4: while stopping criterion is not satisfied do
5: k ← k + 1
6: hk ← arg minh ϕ(θk−1,h)
7: θk ← Pc
{
θk−1 − γk∇θϕ(θk−1,hk)
}
8: end while
9: θˆ ← θk and hˆ← hk.
5. DESCRIPTION OF COMPUTER-SIMULATION STUDIES
Computer-simulation studies were conducted to demonstrate the use of the proposed joint reconstruction
approach and the variable projection (VP) method.
5.1 Generation of pressure data
A two-dimensional (2D) numerical phantom was employed to represent the absorbed optical energy density
distribution as shown in Fig. 1a. The phantom was contained in a region of 22.0 × 22.0 mm2 and consisted of
 
 
0
0.2
0.4
0.6
0.8
1
(a) Phantom
0 0.4 0.8 1.2 1.6−0.5
0
0.5
Time (µ s)
 
 
Measured EIR 1
Measured EIR 2
(b) EIRs
Figure 1: (a) A absorbed energy density phantom. (b) True EIR and inaccurate EIR.
six disks with various radii and absorbed optical energy densities. We employed a 2D ring-shaped transducer
array of radius 25 mm. The transducer array consisted of 128 uniformly-distributed detecting elements each of
Proc. of SPIE Vol. 9323  932313-3
Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 9/4/2018
Terms of Use: https://www.spiedigitallibrary.org/terms-of-use
aw
which sampled at 20 MHz. The speed-of-sound was assumed to be constant and set at c0 = 1.5 mm/µs. At each
transducer, a data vector was formed by use of the C-D imaging model,9 and therefore no inverse crime was
committed. We calculated 600 equally spaced temporal samples over the interval [10, 25) µs. Subsequently, the
noiseless data vector uq was obtained by convolving the pressure data with EIR-1 in Fig. 1b. To the noise-free
data set, we added 3% Gaussian white noise to simulate noise-contaminated measurements.
5.2 Image reconstruction methods
From the simulated noiseless and noisy data, images were reconstructed by solving the minimization problem
in Eqn. (3) by use of Algorithm 2. Conventional quadratic smoothness penalties were employed:
R1(θ) =
N−1∑
n=0
∑
k∈Nn
([θ]n − [θ]k)2 , (8)
R2(h) =
I−2∑
i=0
([h]i+1 − [h]i)2, (9)
where Nn is the index set of four neighboring pixels of the n-th pixel.
The reconstructed image was represented by 440× 440 pixels. The initial guess of the EIR employed in the
VP algorithm was different than the EIR that was assumed when generating the simulated data. This served
to simulate a situation in which an experimentally measured EIR contained errors. Algorithm 2 was terminated
after 500 iterations since it was observed that the changes in the reconstructed images with more iterations were
negligible. EIR-2 in Fig. 1b was employed to initialize the VP algorithm. For comparison, we also reconstructed
images by use of a conventional iterative image reconstruction method that considered the EIR to be fixed. The
accuracy of the reconstructed images was quantified by their root-mean-square error (RMSE) from the numerical
phantom.
6. NUMERICAL RESULTS
The images reconstructed by use of the conventional iterative method and the VP method from noise-free
data are shown in Fig. 2a and Fig. 2b, respectively, and their profiles are plotted in Fig. 2c. As confirmed by
 
 
0
0.1
0.2
0.3
0.4
0.5
X
(a)
 
 
0
0.2
0.4
0.6
0.8
1
X
(b)
−11 0 11
0
0.2
0.4
0.6
0.8
1
X (mm)
A
 
 
VPM
CIM
Phantom
(c)
Figure 2: (a) and (b) Images reconstructed from noiseless data using the conventional iterative method (CIM)
with λ = 1.0 × 10−4 and the VP method (VPM) with λ = 1.0 × 10−4 and α = 1000, respectively. (c) Image
profiles through noiseless images seen in Fig. 2 corresponding to the reconstructions with the VP method (solid
line), with the conventional iterative method (dashed line), and phantom (dotted line). The locations of the
profiles are indicated by the “X” arrows in the Fig. 2a, and 2b, respectively.
comparing with the numerical phantom in Fig. 1a, the use of an inaccurate EIR can result in strong artifacts
and distortions in images reconstructed by use of the conventional methods (RMSE = 0.1552). As revealed by
Fig. 2b and the profiles in 2c, the VP algorithm yielded an image with fewer artifacts and distortions, and the
Proc. of SPIE Vol. 9323  932313-4
Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 9/4/2018
Terms of Use: https://www.spiedigitallibrary.org/terms-of-use
0 0.4 0.8 1.2 1.6−0.5
0
0.5
Time (µ s)
 
 
Recovered
True EIR
Initial EIR
Figure 3: Recovered EIR
image fidelity was improved as reflected by the reduced RMSE of 0.0084. The accuracy of the recovered EIR
was also improved as indicated by Fig. 3.
Figure 4 shows the images reconstructed from the noisy data and their profiles. The RMSEs are 0.1601 and
0.0283 for the conventional iterative method and VP method, respectively. These preliminary results suggest
that the proposed joint reconstruction method is numerically stable.
 
 
0
0.1
0.2
0.3
0.4
0.5
X
(a)
 
 
0.2
0.4
0.6
0.8
1
X
(b)
−11 0 11
0
0.2
0.4
0.6
0.8
1
X (mm)
A
 
 
VPM
CIM
Phantom
(c)
Figure 4: Images reconstructed from noise data using (a) the conventional iterative method with λ = 5.0× 10−3
and (b) the VP algorithm with λ = 5.0× 10−3 and α = 10000. (c) Image profiles through noiseless images seen
in Fig. 4 corresponding to the reconstructions with the VP method (solid line), with the conventional iterative
method (dashed line), and phantom (dotted line). The locations of the profiles are indicated by the “X” arrows
in the Fig. 4a, and 4b, respectively.
7. SUMMARY
It is known that neglecting or employing a contaminated EIR in PACT imaging reconstruction can result
in distortions and artifacts in the reconstructed images. However, accurate measuring or modeling EIR can be
difficult. In this study, we proposed a joint reconstruction approach for PACT that mitigates artifacts in the
reconstructed images caused by use of an inaccurate EIR. A nonlinear least squares minimization problem was
formulated, which exploited the separable nature of the imaging model, and a VP method was employed to solve
the minimization problem. The results demonstrate that the joint reconstruction approach for estimating both
the EIR and the absorbed optical energy density can increase the fidelity of the reconstructed images.
ACKNOWLEDGMENTS
This work was supported in part by NIH awards CA167446 and EB01696301.
Proc. of SPIE Vol. 9323  932313-5
Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 9/4/2018
Terms of Use: https://www.spiedigitallibrary.org/terms-of-use
REFERENCES
1. Oraevsky, A. A. and Karabutov, A. A., [Optoacoustic Tomography ], Biomedical Photonics Handbook, CRC
Press (2003).
2. Wang, L. V., “Tutorial on photoacoustic microscopy and computed tomography,” Selected Topics in Quan-
tum Electronics, IEEE Journal of 14, 171–179 (Jan 2008).
3. Kruger, R. A., Reinecke, D. R., and Kruger, G. A., “Thermoacoustic computed tomographytechnical con-
siderations,” Medical Physics 26(9) (1999).
4. Cox, B. T., Arridge, S. R., Ko¨stli, K. P., and Beard, P. C., “Two-dimensional quantitative photoacoustic
image reconstruction of absorption distributions in scattering media by use of a simple iterative method,”
Applied Optics 45, 1866–1875 (Mar 2006).
5. Wang, K. and Anastasio, M. A., “Photoacoustic and thermoacoustic tomography: Image formation princi-
ples,” in [Handbook of Mathematical Methods in Imaging ], Scherzer, O., ed., Springer (2011).
6. Anastasio, M. A., Zhang, J., Modgil, D., and La Riviere, P., “Application of inverse source concepts to
photoacoustic tomography,” Inverse Problems 23, S21–S35 (2007).
7. 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).
8. Rosenthal, A., Ntziachristos, V., and Razansky, D., “Model-based optoacoustic inversion with arbitrary-
shape detectors,” Medical Physics 38(7) (2011).
9. 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).
10. Rosenthal, A., Ntziachristos, V., and Razansky, D., “Optoacoustic methods for frequency calibration of
ultrasonic sensors,” Ultrasonics, Ferroelectrics and Frequency Control, IEEE Transactions on 58, 316–326
(February 2011).
11. Chung, J. and Nagy, J., “An efficient iterative approach for large-scale separable nonlinear inverse problems,”
SIAM Journal on Scientific Computing 31(6), 4654–4674 (2010).
12. Aravkin, A. Y. and van Leeuwen, T., “Estimating nuisance parameters in inverse problems,” Inverse Prob-
lems 28(11), 115016 (2012).
13. Golub, G. and Pereyra, V., “The differentiation of pseudo-inverses and nonlinear least squares problems
whose variables separate,” SIAM Journal on Numerical Analysis 10(2), 413–432 (1973).
14. Xu, M. and Wang, L. V., “Photoacoustic imaging in biomedicine,” Review of Scientific Instruments 77(4),
– (2006).
15. Wang, L. V., ed., [Photoacoustic Imaging and Spectroscopy ], Optical Science and Engineering, CRC Press
(March 2009).
Proc. of SPIE Vol. 9323  932313-6
Downloaded From: https://www.spiedigitallibrary.org/conference-proceedings-of-spie on 9/4/2018
Terms of Use: https://www.spiedigitallibrary.org/terms-of-use


Photoacoustic computed tomography without accurate ultrasonic transducer responses

https://authors.library.caltech.edu/89448/1/932313.pdf

Photoacoustic computed tomography without accurate ultrasonic transducer responses

Abstract

Similar works

Full text

Available Versions

Caltech Authors - Main