Photoacoustic computed tomography (PACT) is an emerging imaging modality that combines the high optical contrast of blood-rich structures with the high spatial resolution of ultrasound detection.1–3 In PACT, the biological tissue of interest is irradiated with a short laser pulse. Under the condition of thermal confinement, the absorption of the optical energy results in the generation of pressure waves via the photoacoustic (PA) effect.4,5 These pressure waves are subsequently detected using broadband ultrasound transducers. The image reconstruction problem in PACT is estimating the absorbed optical energy density from the measured PA signals. Such an image may be of great importance for preclinical and clinical use as there exists a strong correlation between optical absorption and the pathological condition of tissue.3,6,7
The majority of PACT reconstruction methods8,9 are based on idealized imaging models that assume an acoustically homogeneous medium. However, these assumptions are violated in many applications of PACT. For example, in small animal imaging applications, bone and/or gas pockets represent strong sources of acoustic heterogeneity. When the spatially variant acoustic properties of the medium are not accounted for in the imaging model, the reconstructed images may contain significant distortions and artifacts.10,11 However, in practice, it is challenging to estimate these acoustic properties accurately.12
To circumvent the need for detailed information regarding an object’s acoustic properties as well as to suppress artifacts due to the presence of acoustic heterogeneity, we previously proposed a half-time-based reconstruction method.13 The half-time-based reconstruction method exploits data redundancies14 to uniquely and stably reconstruct images from measurement data that are uniformly truncated with respect to the temporal coordinate. Although the half-time-based reconstruction method mitigates acoustic heterogeneity-induced artifacts, the reconstructed images can still contain significant distortions. Moreover, half-time-based methods do not employ any a priori information about the location of strong acoustic heterogeneities in the reconstruction.
In addition to truncation-based strategies, work has also been conducted on incorporating statistical information about the object to mitigate artifacts due to acoustic heterogeneities.12 In that approach, a priori information about the acoustic properties of the object is utilized to probabilistically weigh the tomographic contribution of each detector to a pixel in the reconstructed image.12 Such a statistical approach was shown to have merit for mitigating artifacts due to weak acoustic heterogeneities. However, in the presence of strong heterogeneities, the approach has not been demonstrated to be effective.12
In this work, a PACT reconstruction method that is based on a variable data truncation (VDT) approach is proposed. This method represents a generalization of the half-time reconstruction method. The VDT-based reconstruction method employs a priori information about the location of the isolated acoustic heterogeneities but does not require information regarding its acoustic properties. In the VDT-based method, the degree of temporal truncation is dependent on the location of the isolated acoustically heterogeneous region relative to the ultrasound transducer positions. Due to the adaptive nature of the temporal truncation, artifacts arising from the acoustic heterogeneities can, in some cases, be more effectively suppressed in images reconstructed by this method as compared with images reconstructed by the half-time-based method.
Background: Imaging Models and Reconstruction Methods for Photoacoustic Computed Tomography
Photoacoustic Wavefield Propagation: Continuous and Discrete Formulation
Let denote the photoacoustically induced pressure wavefield at location and time . Additionally, let denote the initial pressure distribution generated within the object due to the PA effect and denote the vector-valued acoustic particle velocity. In our formulation of the PACT imaging model, the object and the surrounding medium are assumed to possess homogeneous and lossless acoustic properties. Let denote the medium’s uniform speed of sound (SOS) value and denote the distribution of the uniform ambient density value.
In practice, the detected pressure wavefield is discretized temporally and spatially at specific transducer locations. Let denote the discretized pressure signal. Unless stated otherwise, throughout the study, we shall neglect the effects of the acousto-electric impulse response as well as the spatial impulse response of the transducer.17–19 A continuous-to-discrete PACT imaging model20,21 can be generally expressed as
To obtain a discrete-to-discrete (D-D) imaging model for use in numerically simulating PA wavefield propagation, a finite-dimensional approximate representation of the object function needs to be introduced. Also, we will assume that the reconstruction method estimates a 2-D slice of the object in plane with the circular aperture . Thus, the finite-dimensional representation of the 2-D slice of the object function is given as22 etc.
Thus, given , , and , a D-D imaging model is given as23-based interpolation is performed. The goal of image reconstruction in a discrete setting is to determine an estimate of using the measured data .
Full-Time-Based Reconstruction Methods
In this section, the salient features of the full-time-based backprojection (BP) method and the full-time-based iterative reconstruction algorithm are discussed. In the full-time-based methods, the images are reconstructed from full-time data. The full-time data refer to the data recorded at all transducers for time , where . Here, is the time it takes for pressure waves to propagate from the center of the circular measurement aperture to the ’th transducer in an acoustically homogeneous medium, where .
A variety of BP-type PACT image reconstruction methods24–27 have been developed based on the continuous imaging model described by Eq. (4). In the presence of discretization, these methods provide an approximate estimate of the object function. The finite-dimensional estimate of the object function obtained using the BP method is given as2828 the BP term in Eq. (9) does not contain the temporal derivative of pressure. This allows us to mitigate the impact of noise in the reconstructed image as the derivative term amplifies the contribution of noise in the measured data to the reconstructed image. Moreover, the BP formula in Eq. (9) assumes a full spherical detector surface.28 Since the detector geometry used in our applications is circular, the BP reconstruction formula presented in Eq. (9) is only approximately applicable. Nevertheless, the reconstruction formula presented in Eq. (9) has been used empirically to yield reconstructed images of the initial pressure distribution.29,30 However, a quantitatively accurate 2-D filtered BP reconstruction formula for circular detector geometry has been established.26,27
Optimization-based image reconstruction algorithms
Since BP-type reconstruction methods are based on the discretization of an analytical reconstruction formula, they possess several limitations. For instance, BP reconstruction methods require the measured data to be densely sampled on an aperture that encloses the object. Moreover, analytical reconstruction methods ignore measurement noise. Hence, they can yield reconstructed images that have suboptimal trade-offs between image variance and spatial resolution. The use of iterative PACT image reconstruction algorithms can circumvent these shortcomings.18
Commonly employed iterative PACT reconstruction algorithms seek to compute penalized least squares estimates by solving an optimization problem of the form
Half-Time-Based Reconstruction Methods
It has been shown that images reconstructed using full-time data can contain significant distortions when acoustic variations in the object are not accurately modeled.11,14,31 To address this problem, an image reconstruction method based on a data truncation strategy known as the half-time method13,14 was previously proposed. In the half-time method, all transducer measurements are temporally truncated at a specific delay time , where . For a circular measurement aperture, it can be shown that .
In half-time-based reconstruction methods, the data vector is truncated at a constant delay time . This truncation process can be described by a truncation matrix that acts on to produce a truncated data vector as
Details regarding the construction of truncation matrix for the half-time method are described next.
Define the matrix as
By replacing the full-time data vector with a half-time data vector in a full-time reconstruction method, a half-time reconstruction method is readily obtained. For example, a half-time BP method is expressed as
Image Reconstruction Based on Variable Truncation Methods
Even though half-time-based reconstruction methods mitigate acoustic heterogeneity-induced artifacts, the reconstructed image can, in some cases, still contain significant artifacts. In scenarios where an acoustically heterogeneous region is localized and its support is relatively small compared with the area of the reconstruction region, the half-time truncation method can be readily extended to a VDT method. Unlike the half-time method, the temporal truncation strategy employed in the VDT method utilizes a priori information about the location of the acoustically heterogeneous region. The VDT method is applicable when the acoustic heterogeneity is isolated in location from the features that are likely to be imaged and when the support of the acoustically heterogeneous region is small relative to the area of the reconstruction region. In small animal imaging applications, bone, spinal column, and gas pockets represent isolated acoustic heterogeneities.
In the VDT method, the measurement data recorded at each transducer are temporally truncated based on the distance between the corresponding transducer and the nearest isolated acoustic heterogeneity. As a result, all data that have been strongly influenced by an acoustic heterogeneity are discarded and not employed for image reconstruction.
As depicted in Fig. 1(a), the pressure signals that have been distorted by traveling through the acoustic heterogeneity are not truncated in the half-time measurement data. When a reconstruction method employing such data does not account for this distortion, it can be a source of significant artifacts in the reconstructed image. On the other hand, the VDT measurement data do not contain pressure signals that are reflected by or transmitted through the acoustic heterogeneity [see Fig. 1(b)]. Hence, the artifacts introduced by these distorted pressure signals are eliminated in the images reconstructed using VDT-based methods. Additionally, for the transducer shown in Fig. 1(b), the VDT-based method uses more temporal data than the half-time-based reconstruction method. In the subsequent section, we will describe the construction of the truncation matrix for the VDT method.
Construction of Variable Data Truncation Matrix
Similar to the half-time-based method, the truncation process in the VDT-based method can be described by a truncation matrix that acts on the data vector . Thus, replacing the truncation matrix in Eq. (14) with , we haveFig. 2). The nearest grid point to the ’th transducer at that corresponds to an acoustic heterogeneity is given as
Variable Data Truncation Reconstruction Methods
VDT-based reconstruction methods can be formed readily by replacing the full-time data vector in a full-time reconstruction method by its truncated version. Thus, the VDT-based BP method to reconstruct is given as
Computer-simulation studies were conducted to compare the performance of the VDT- and half-time-based reconstruction methods.
The -space pseudospectral method32,33 for numerically solving the PA wave equation was implemented in the MATLAB® -wave toolbox.34 This toolbox was employed to compute the action of forward operator . A 2-D circular scanning geometry consisting of 512 transducers that were evenly distributed in a circle of radius 50 mm was employed. The numerical phantoms employed in this study to generate forward data are shown in Fig. 3. The acoustic heterogeneity employed in the simulation, shown in Fig. 3(b), imitated an air void ( and ) while the background medium consisted of water ( and ). The initial pressure phantom in Fig. 3(a) consisted of two line absorbers placed perpendicular to one another, operator .
Assuming ideal point-like transducers, the simulated pressure data corresponding to the numerical phantoms were calculated using the -space pseudospectral method for the scanning geometry described. A grid with a pitch of was employed to simulate the pressure data. A total of 5200 temporal samples were computed at each transducer location with a time step of . In addition, the generated forward data on the sensors were contaminated with additive white Gaussian noise to achieve a signal-to-noise ratio of 10.
For both the iterative and BP reconstruction methods, a constant SOS of was assumed. For the VDT-based reconstruction methods, based on the known location of the acoustic heterogeneity and the assumed , the truncation matrix was established as prescribed in Eq. (20). The truncation matrix was utilized to compute the truncated data vector . Similarly, the truncation matrix and the truncated pressure data vector were also computed for use with the half-time-based reconstruction methods. Furthermore, the fast iterative shrinkage/thresholding method (FISTA)35,36 was implemented to solve the optimization problems in Eqs. (18) and (24).
The images reconstructed using the half-time- and VDT-based methods are shown in Fig. 4. In the images reconstructed using the half-time-based methods, the artifacts due to the acoustic heterogeneity (air void) are marked by white arrows in Figs. 4(a) and 4(c). The same artifacts, however, are mitigated in the images reconstructed using VDT-based methods, which are shown in Figs. 4(b) and 4(d). In addition, differences exist between the images reconstructed using the BP and iterative methods. This difference is due to the action of the TV-penalty term in the iterative algorithm, which smoothes the image while preserving its edges. The mitigation of the artifacts due to the air void can be quantified by calculating the root-mean-squared error (RMSE) between the reconstructed pressure distribution and the true pressure distribution. The RMSEs between the reconstructed pressure distribution and the original pressure distribution are shown in Table 1. The RMSE values for the images reconstructed using the VDT-based method are lower than the RMSE values for the images reconstructed using the half-time-based methods for both the BP and iterative methods.
The RMSE between the reconstructed pressure distribution and the initial pressure distribution.
|Half-time based||VDT based|
The performances of the half-time- and VDT-based reconstruction methods were also investigated using experimental data. The experimental data were acquired in the Optical Imaging Laboratory at Washington University in St. Louis.
The 2-D PACT system employing a ring array composed of 512 transducer elements distributed in a radius 50 mm was employed in the experimental studies. A short laser pulse (DLS9050, Continuum, 5- to 9-ns pulse width) with a repetition rate of 50 Hz at a wavelength of 532 nm was used to irradiate a sample located in the center of the measurement system. The generated acoustic signals were detected by transducers with a center frequency of 5 MHz and a bandwidth > 90%. The data recorded by the transducers were digitized at a sampling rate of 40 MHz.
Experimental Phantom Studies
The experimental agar phantom, shown in Fig. 5, consisted of two linear optical absorbers (human hair) placed approximately perpendicular to one another with an air void insert in the lower right quadrant.
The phantom was placed at the center of the circular transducer array and PA signals were acquired for 2000 time points with a sampling rate of 40 MHz. Since the VDT-based reconstruction methods require a priori information about the acoustically heterogenous region, the location of the air void relative to the transducers needed to be determined. This was accomplished through visual inspection of the image shown in Fig. 5. For both the reconstruction methods, the constant SOS of the homogeneous background was assumed to be . Additionally, for the iterative reconstruction algorithm, the experimentally obtained data were preprocessed prior to reconstruction. The preprocessing of the data involved temporally upsampling by a factor of 4 and filtering with a Hann-window low-pass filter with a cutoff frequency of 10 MHz. The preprocessing was done to avoid any issues with the numerical stability of the wave equation solver.34 Furthermore, the FISTA35,36 was implemented to solve the optimization problem in Eqs. (18) and (24).
Images reconstructed using the half-time- and VDT-based BP and iterative reconstruction methods are shown in Fig. 6. In the images reconstructed using half-time-based approaches, the artifacts due to the presence of an air void are marked with red arrows. Also, the location of the acoustically heterogeneous region (air void) that was extracted from Fig. 5 is delineated by a white boundary, as shown in Figs. 6(b) and 6(d).
From Fig. 6, the artifacts due to acoustic heterogeneity are found to be mitigated in the images reconstructed using VDT-based reconstruction methods. The artifacts are present in the lower right quadrant, as shown in Figs. 6(a) and 6(c). In addition, we also observe that the VDT-based methods are robust with regards to errors in the estimation of the boundary of the acoustically heterogeneous region. Thus, even with such inaccuracies in estimating the location and shape of air void, the VDT-based methods achieve effective mitigation of artifacts due to acoustic heterogeneities.
The second set of experimental data was acquired in vivo from a mouse trunk.
All experimental animal procedures were carried out in conformity with laboratory animal protocols approved by the Animal Studies committee of Washington University in St. Louis. In small animal imaging applications, the major sources of acoustic heterogeneity are bones, spinal column, and gas voids. These heterogeneities perturb the acoustic wavefield causing severe distortions in the reconstructed images.
For the purposes of this study, the VDT-based reconstruction methods were employed to mitigate artifacts only due to the spinal column. To obtain a priori information about the location of the acoustically heterogeneous region, we developed a semiautomatic segmentation method. The semiautomatic segmentation method estimated the boundary of the spinal column from the images reconstructed using the half-time-based BP method. As the VDT-based reconstruction methods are robust to errors in the estimation of boundary of the acoustically homogeneous region, the segmentation method did not need to be exact. The boundary of the spinal cord that was extracted by the segmentation method is depicted by the white heart-shaped contours in Figs. 7(c) and 7(d). We can observe that the segmentation method overestimated the area of the spinal column. Thus, all of the pressure signals reflected off or transmitted through the spinal column are truncated.
For both methods, the constant SOS of the homogeneous background was assumed to be . To select the optimal SOS value of the background, we scanned over a range of values and compared the image quality of the images reconstructed using the half-time-based BP method. The value that gave the best image quality was picked as the optimal SOS of the homogenous background. Similar to the experimental agar phantom data set, the data obtained from the mouse trunk were preprocessed prior to applying the iterative reconstruction algorithm. The preprocessing involved temporally upsampling by a factor of 4 and filtering using a Hann-window low-pass filter with a cutoff frequency of 10 MHz. Furthermore, the FISTA35,36 was implemented to solve the optimization problems in Eqs. (18) and (24).
Images reconstructed using the half-time- and VDT-based BP and iterative reconstruction methods are shown in Fig. 7. The acoustically heterogeneous region is delineated by heart-shaped contours in the images reconstructed using VDT-based methods. To compare the performance of the VDT- and the half-time-based reconstruction methods, we analyzed the differences in the regions enclosed by white rectangular boxes in Fig. 7. The zoomed-in delineated regions of Fig. 7 are shown in Fig. 8.
In Fig. 8, vessels of interest are marked by white arrows. In the images reconstructed using half-time-based methods, even though we observe the bifurcation of the lower part of the main vessel, the top part of the main vessel is obscured [see Figs. 8(a) and 8(c)]. However, in images reconstructed using VDT-based reconstruction methods, as shown in Figs. 8(b) and 8(d), in addition to the bifurcation of the main vessel, the top part of the vessel along with a vessel branching off from it is visible. The vessel branch and the main vessel are marked by white arrows in Figs. 8(b) and 8(d).
When using half-time-based reconstruction methods, the pressure signals reflected off and transmitted through the acoustic heterogeneity (spinal column) can interfere with PA signals coming from the blood vessels, ultimately obscuring the blood vessels in the reconstructed images. However, with the VDT-based reconstruction methods, the pressure signals reflected or transmitted through the acoustic heterogeneity are truncated. Thus, certain features obscured in the images reconstructed using half-time-based methods can be clearly seen in the images reconstructed using VDT-based methods.
A VDT approach to PACT image reconstruction was proposed and investigated. The performance of VDT- and half-time-based reconstruction methods were compared using simulated and experimental data sets. The PACT reconstruction methods employed, the iterative and BP methods, were modified to include data truncation strategies in their formulation. From the results, the artifacts, due to acoustic heterogeneities, were mitigated much more effectively in images reconstructed using VDT-based methods as compared with the images reconstructed using half-time-based methods. This improvement in image quality was significant in the reconstructed mouse trunk images, where we found that structures obscured in the images reconstructed using half-time-based methods were visible in the images reconstructed using VDT-based methods.
However, if the size of the heterogeneity is large or if there are multiple heterogeneities, the VDT-based approach can lead to excessive temporal truncation of the data. For such cases, the VDT-based methods may not perform better than the half-time-based methods. The problem of multiple heterogeneities represents a topic for further study.
L. V. W. has a financial interest in Microphotoacoustics, Inc., which, however, did not support this work.
This work was supported in part by National Institutes of Health award Nos. CA1744601, EB01696301, and 5T32EB01485505.
Joemini Poudel received his BASc degree in biomedical engineering from Simon Fraser University in Vancouver, Canada. He is currently a PhD student in biomedical engineering at Washington University in St. Louis (WUSTL). His research interest involves application of reconstruction algorithms for photoacoustic computed tomography (PACT).
Thomas P. Matthews is a PhD student in the lab of Dr. Mark Anastasio. He has received a MS degree in biomedical engineering from WUSTL, and BE degree in engineering and an AB degree in physics from Dartmouth College. His research interests include image reconstruction methods in ultrasound computed tomography and PACT, with a focus on optimization-based methods, accurate physics-based models of imaging systems, and high-performance scientific computing.
Lei Li received his BS and MS degrees from Harbin Institute of Technology, China, in 2010 and 2012, respectively. He is working as a graduate research assistant under the tutelage of Dr. Lihong Wang at Washington University. His current research focuses on photoacoustic microscopy and PACT, especially on improving photoacoustic imaging speed and applying it to brain functional imaging and small-animal whole-body imaging.
Mark A. Anastasio received his PhD from the University of Chicago in 2001. He is currently a professor of biomedical engineering at WUSTL. His research interests include tomographic image reconstruction, imaging physics, and the development of computed biomedical imaging systems. He is an internationally recognized expert in the fields of imaging science, phase-contrast x-ray imaging, and photoacoustic computed tomography.
Lihong V. Wang received his PhD at Rice University, Houston, Texas. Currently, he is holding the Gene K. Beare distinguished professorship of biomedical engineering at Washington University in St. Louis. He has published 400 peer-reviewed journal articles and delivered 400 keynote, plenary, or invited talks. His google scholar h-index and citations have reached 110 and over 50,000, respectively.