Photoacoustic imaging is based on the detection of acoustic waves that are generated by thermoelastic expansion induced by optical absorption of laser light. The majority of photoacoustic imaging systems can be classified as either microscopy systems or computed tomography (CT) systems.1 In general, the microscopy-based systems do not require a reconstruction step as the lateral resolution is set by either a focused transducer or focused laser illumination. CT systems can be further subdivided into scanning systems (also called synthetic aperture) and array systems. From a mathematical point of view, synthetic apertures and arrays require the same reconstruction algorithm, although arrays tend to have smaller individual transducers so the directivity can be more uniform over the imaging area. Transducer arrays have a speed advantage as they do not need to be mechanically scanned, which can be important in many applications. In addition, the array-based imaging systems have the potential to be integrated with ultrasound imaging systems, and indeed many of the arrays in use in photoacoustic systems are taken from ultrasound systems.2–4 In this paper we consider a linear transducer array that is similar to the ones used in current photoacoustic and ultrasound imaging systems.
There are several reconstruction algorithms used in linear array imaging systems. Filtered backprojection is based on an analytical inversion of the solution to the photoacoustic wave equation.5 It is theoretically exact for a continuous, infinite linear transducer array, and when the object to be imaged is acoustically homogeneous and has no acoustic absorption. Unfiltered (or simple) backprojection, although not theoretically exact, is numerically equivalent to the delay-and-sum beamformer method employed in ultrasound imaging systems, and hence is very readily implemented in combined ultrasound-photoacoustic imaging systems. The Fourier method of image reconstruction is also based on inverting the solution to the photoacoustic wave equation, but this inversion is performed in the frequency domain.6 For sufficiently small datasets or given sufficiently powerful computer hardware, all of these simple reconstruction algorithms can be used in real-time imaging systems. There are more sophisticated reconstruction algorithms that can correct for acoustic inhomogeneity and absorption such as iterative methods.7 In general, these iterative algorithms produce high-quality images; unfortunately, using current computers they are not capable of real-time imaging.
The minimum-norm least-squares (MNLS) solution to reconstructing a photoacoustic image is given by the pseudoinverse of the system matrix, which relates the initial pressure distribution to the measured pressure wave. In this paper, we calculate the pseudoinverse using a singular-value decomposition (SVD). The pseudoinverse calculation involves dividing by the nonzero singular values of the system matrix. In practice, we have to regularize the problem by throwing out small singular values that would otherwise amplify the noise. Where to set the regularization threshold depends on the specific signal and noise levels. The noise in an imaging system can generally be divided into object-dependent and object-independent noise. In many imaging systems, such as x-ray CT, the object-dependent noise varies strongly across detectors. In photoacoustics, however, the object-dependent noise is relatively uniform across transducer elements (as long as the object is not too close to the transducer). The object-independent noise (such as electronic noise) is also uniform across the transducer elements. The regularization level still depends on the relative strength of the signal and noise, but it is possible to calculate a few pseudoinverse matrices based on different singular-value truncations, store these in memory, and then decide in real time which one to use based on image quality. Due to the more uniform noise properties, we expect this method to work relatively better in photoacoustics than, say, x-ray CT, which would benefit from a more advanced regularization. Using this method the reconstruction step consists only of a matrix-vector multiplication, which is very fast. In a recent publication, the system matrix of a photoacoustic imaging system was measured experimentally, and the corresponding pseudoinverse was calculated ahead of time and used for real-time imaging.8 This method has the advantage of including all the relevant acoustic properties in the system matrix, but the disadvantages of requiring the experimental measurement of the system matrix and the inevitable inclusion of system noise in the measurement. In this paper we develop a model for a two-dimensional (2-D) system matrix and use that to calculate the pseudoinverse. We also develop a three-dimensional (3-D) forward model to simulate a photoacoustic imaging system.
A natural question then arises as how best to select an algorithm to optimize image quality. In a previous work,9 we performed a preliminary comparison between the filtered backprojection and pseudoinverse reconstruction algorithms. In that paper, we looked at various metrics such as resolution, streak length, and ideal observer signal-to-noise ration (SNR) for a point source object; however, these metrics can be misleading. The resolution and streak length tell nothing of the noise properties of the algorithms, and the ideal observer SNR only measures how much information is lost for each algorithm in going from the raw data to the reconstructed image, which does not necessarily correlate with image quality as judged by human observer task performance.
In this paper we demonstrate the use in photoacoustics of so-called task-based metrics of image quality that seek to define image quality based not on physical measures but on the utility of an image for performing a predefined task. Such task-based image-quality assessment has gained widespread attention in nuclear medicine and x-ray imaging modalities.10 While the gold standard of task-based assessment is human observer studies using the receiver operating characteristic (ROC) methodology, these are tedious and labor-intensive. As an alternative, various model observers have been proposed that seek to emulate human observer performance in well-defined tasks, such as the detection of a known signal in a known or random background. The channelized Hotelling observer (CHO) has been shown to predict human performance particularly well.10 It entails integrating images over a set of frequency-domain channels that mimic the frequency response of the human visual system and then applying a linear discriminant to the output. To assess the image quality of the reconstruction algorithms in this study, we use the CHO with two sets of channels that have been shown to predict the performance of human observers: constant 11 and difference of Gaussian (DOG) channels.12
The photoacoustic equation can be written as113 Fig. 1 and consists of a linear transducer array located on the axis with a cylindrical lens on its face and consisting of individual transducers, an object to be imaged in the half-space , and a light source above the object pointing in the positive direction. The cylindrical lens rejects acoustic waves that do not originate near the imaging plane. We assume uniform light illumination on the imaging plane to avoid modeling photon propagation. With this geometry, we can separate the heating function spectrum into , where is the absorption coefficient and is the irradiance spectrum. The average pressure on each individual transducer is the pressure integrated over and divided by the area, where the width in the direction is . The cylindrical acoustic lens imparts a phase of to the wave, where is the acoustic path length of the lens (refractive index times total thickness of the lens) and is the focal length. The lens has a height of and the transducer has a transfer function denoted by . This gives the measured voltage spectrum as
In this paper we consider a 32-element transducer with cylindrical lens focal length of 6 mm, height of 3 mm, individual element width of 0.15 mm, and center-to-center element spacing of 0.15 mm. The transducer transfer function is Gaussian with a center frequency of 3.5 MHz and 50% bandwidth, and we assume the laser pulse is sufficiently short so that , where is a constant. Note that the transducer spacing satisfies the Nyquist sampling criterion. For the forward model, we use the full 3-D expression given by Eq. (3). For the system matrix used to calculate the pseudoinverse, we use a 2-D approximation of Eq. (3) with . For the discrete version, we use the notation , where is a vector representing the voltage measured by the transducers, is the system matrix, and is a vector representing the initial pressure distribution of the object being imaged. We use a basis with a pixel size of 37.5 μm. The temporal sampling rate is 100 ns with 71 samples per transducer, which produces a system matrix of size .
We performed an SVD on the system matrix, which was then used to calculate the pseudoinverse. The SVD factors the system matrix into , where and are unitary matrices, the dagger represents conjugate transpose, and is a diagonal matrix consisting of the singular values in decreasing order. The pseudoinverse is then given by , where is a diagonal matrix consisting of the reciprocal of singular values above a threshold. The thresholding of the singular values serves as a regularization. Note that this method is also referred to as truncated SVD inversion. Figure 2 shows a plot of the magnitude of the singular values. The image, which is an estimation of the initial pressure , is then given by .
Filtered Backprojection and Simple Backprojection
The filtered backprojection formula5 can be written as14 A Hamming window was applied along the transducer elements for apodization. The simple backprojection algorithm is equivalent to the filtered backprojection algorithm except the pressure is not time-integrated or filtered with the Ram-Lak filter. Note that this filtered backprojection formula is built on a foundation developed by Norton in the context of reflection tomography.15 An alternative algorithm has been developed by Burgholzer et al. but has not been implemented in the present study.16
The Fourier transform-based formula6 can be written as
For the filtered backprojection, simple backprojection, and Fourier based methods we deconvolve the transducer response using
Figures 3 and 4 show the image reconstruction of a point source and a 0.5-mm cube, respectively, for various singular-value truncations as well as the backprojection, filtered backprojection, and Fourier-based algorithms with different regularization parameters. Uncorrelated Gaussian noise is added to before reconstruction to simulate electronic noise. The standard deviation of the noise is set to be equal to the peak signal level divided by 10. There is no background for these images. Figure 5 shows the image reconstruction of a lumpy background17 for the different algorithms. For each lumpy background image, there is a corresponding image with lumpy background and point source signal. All of the images in Fig. 5 contain the same background; in the next section the number of lumps and their locations are treated as random variables. The lumpy background is only present in a volume represented by a 0.5 mm cube at the center of the object. Uncorrelated Gaussian noise is also added to before reconstruction, where, again, the standard deviation is set to the peak signal level divided by 10.
As mentioned in the Introduction, the CHO has been shown to predict human observer performance for various tasks. The choice of channels depends on the specific task to be performed. In this paper we use two sets of channels: constant 11 and DOG channels.12 Constant channels have been shown to work well at predicting human performance when the noise is additive and Gaussian. DOG channels work well for additive Gaussian noise as well as for a lumpy background. The channel profiles are shown in Fig. 6.
In order to quantify image quality, we compute the channelized Hotelling SNR given by117). For the uniform background tasks, uncorrelated Gaussian noise is added to before reconstruction, as discussed in the previous section. These tasks are referred to as signal known exactly, background known exactly (SKE/BKE). For the lumpy background17 task, the number of lumps and their locations are treated as random variables. The number of lumps, , is modeled as a Gaussian random variable with a mean and variance of 1250. That number was selected using , where is the number of discrete object points in a 0.5-mm cube in our simulations. The location of the lumps is modeled as a uniform random variable over the volume represented by a 0.5-mm cube at the center of the object. Each lump is modeled as a point source with a magnitude that is 100 times lower than the point source signal magnitude. Uncorrelated Gaussian noise is also added to before reconstruction for this task. This task is referred to as SKE in a lumpy background and is a good model of tumor contrast agent imaging where the signal represents uptake of the contrast agent by the tumor and the lumpy background represents contrast agent uptake by the surrounding tissue.
We define the CHO efficiency as
In Figs. 7Fig. 8 through 9, we show the CHO efficiency for the three tasks. For the SKE/BKE tasks, the efficiency is highest when the regularization is chosen so that the noise is lowest even though the resolution is worse as seen in Figs. 3 and 4. The performance of the algorithms is similar for the SKE/BKE tasks except that the Fourier-based algorithm performs best for the identification of the cube as determined by the difference of Gaussian channels. However, the SKE/BKE task is generally not a realistic model of actual imaging systems. For the SKE in a lumpy background task, which is more realistic, the method based on the pseudoinverse performs best according to both sets of channels. Also notice that for this task the optimal regularization is not necessarily to choose the lowest possible noise amplification. It is also interesting that for the three tasks studied in this paper, the simple backprojection performs as well as the filtered backprojection.
In this paper, we have used task-based image quality metrics to compare image reconstruction algorithms in photoacoustic tomography. The study was not intended to be the final word on algorithm selection but rather to introduce the use of task-based assessment and model observers to the field. As the field grows and matures, it could be valuable to deploy such tools for optimizing systems and algorithms for various tasks.
We have shown that for the case of a linear transducer array and an SKE in a lumpy background task, the image reconstruction algorithm based on the pseudoinverse performs best according to the CHOs we look at in this paper. Future work is planned to study the validity of the channels used in this paper by performing human observer experiments. Future work is also planned to study tasks that include signal variability. It should be noted that in this paper we assume uniform light illumination on the imaging plane. Depending on the specific system and task to be performed this may not be a realistic assumption, and in this case the laser irradiation geometry and photon propagation should be included in the forward model.
Another interesting conclusion of this study is that the simple backprojection algorithm performs as well as the filtered backprojection algorithm for the tasks we looked at. Since the simple backprojection algorithm is numerically equivalent to the delay-and-sum beamformer reconstruction employed in existing ultrasound imaging systems, combined ultrasound-photoacoustic imaging systems may be able to use this algorithm without loss of performance. This would be advantageous in terms of system design. The type of analysis presented in this paper could be used to determine if this conclusion is valid for the specific task the combined ultrasound-photoacoustic imaging system is being built for.
This work was supported by ACS Research Scholar Grant RSG-08-117-01-CCE to Patrick La Rivière.