**128×128 pixels**. In order to compare the algorithms, we use channelized Hotelling observers that predict the detection ability of human observers. We use two sets of channels: constant

*Q [/**i] and difference of Gaussian spatial frequency channels. We look at three tasks, identification of a point source in a uniform background, identification of a 0.5-mm cube in a uniform background, and identification of a point source in a lumpy background. For the lumpy background task, which is the most realistic of the tasks, the method based on the pseudoinverse performs best according to both sets of channels.*

*
*
Petschke and La Rivière: Comparison of photoacoustic image reconstruction algorithms using the channelized Hotelling observer where ${v}_{s}$ is the speed of sound, $p(\mathbf{r},t)$ is the pressure, $\beta $ is the thermal coefficient of volume expansion, ${C}_{P}$ is the specific heat at constant pressure, and $H(\mathbf{r},t)$ is the heating function. The Green’s function solution for the pressure spectrum can be written aswhere $k=\omega /{v}_{s}$. The geometry we are considering is shown in Fig. 1 and consists of a linear transducer array located on the $x$ axis with a cylindrical lens on its face and consisting of $N$ individual transducers, an object to be imaged in the half-space ${z}_{0}>0$, and a light source above the object pointing in the positive ${y}_{0}$ direction. The cylindrical lens rejects acoustic waves that do not originate near the imaging $({x}_{0}-{z}_{0})$ 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 $\tilde{H}({\mathbf{r}}_{\mathbf{0}},\omega )=\mu ({\mathbf{r}}_{\mathbf{0}})\tilde{I}(\omega )$, where $\mu ({\mathbf{r}}_{\mathbf{0}})$ is the absorption coefficient and $\tilde{I}(\omega )$ is the irradiance spectrum. The average pressure on each individual transducer is the pressure integrated over $x$ and $y$ divided by the area, where the width in the $x$ direction is $w$. The cylindrical acoustic lens imparts a phase of ${e}^{-ikd}{e}^{ik{y}^{2}/(2f)}$ to the wave, where $d$ is the acoustic path length of the lens (refractive index times total thickness of the lens) and $f$ is the focal length. The lens has a height of $h$ and the transducer has a transfer function denoted by $\tilde{T}(\omega )$. This gives the measured voltage spectrum as where $n$ denotes each individual transducer and runs from $-(N-1)/2$ to $(N-1)/2$. The middle transducer has $n=0$, and the center of the middle transducer is located at $x=0$. The extra ${z}_{0}/|\mathbf{r}-{\mathbf{r}}_{\mathbf{0}}|$ is an obliquity factor.where the brackets indicate expectation values, ${\widehat{\mathbf{f}}}_{\mathbf{c}\mathbf{h}}={\mathbf{U}}_{\mathbf{c}\mathbf{h}}\widehat{\mathbf{f}}$, ${\mathbf{U}}_{\mathbf{c}\mathbf{h}}$ is a matrix that represents the spatial frequency filtering and summing of the channels, and ${\mathbf{K}}_{\mathbf{c}\mathbf{h}}$ is the covariance matrix of ${\widehat{\mathbf{f}}}_{\mathbf{c}\mathbf{h}}$. Note that the covariance matrix is only a $5\times 5$ matrix for the constant $Q$ channels and a $3\times 3$ matrix for the DOG channels, which means the number of simulations required to make the matrix full rank is greatly reduced as compared to the full image covariance matrix. In this study we use 500 simulations to construct Eq. (7). For the uniform background tasks, uncorrelated Gaussian noise is added to $\mathbf{g}$ 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 background

## 1.

## Introduction

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}3.^{–}^{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 $Q$^{11} and difference of Gaussian (DOG) channels.^{12}

## 2.

## System Matrix

The photoacoustic equation can be written as^{1}

## (1)

$$({\nabla}^{2}-\frac{1}{{v}_{s}^{2}}\frac{{\partial}^{2}}{\partial {t}^{2}})p(\mathbf{r},t)=-\frac{\beta}{{C}_{P}}\frac{\partial H(\mathbf{r},t)}{\partial t},$$^{13}

## (2)

$$\tilde{p}(\mathbf{r},\omega )=\frac{i\omega \beta}{4\pi {C}_{P}}{\int}_{{V}_{0}}{\mathrm{d}}^{3}{\mathbf{r}}_{\mathbf{0}}\tilde{H}({\mathbf{r}}_{\mathbf{0}},\omega )\frac{{e}^{-ik|\mathbf{r}-{\mathbf{r}}_{\mathbf{0}}|}}{|\mathbf{r}-{\mathbf{r}}_{\mathbf{0}}|},$$## (3)

$${\tilde{V}}_{n}(\omega )=\frac{i\omega \beta \tilde{I}(\omega )\tilde{T}(\omega ){e}^{-ikd}}{4\pi hw{C}_{P}}{\int}_{-w/2+nw}^{w/2+nw}{\int}_{-h/2}^{h/2}\mathrm{d}x\mathrm{d}y{e}^{ik{y}^{2}/(2f)}\phantom{\rule{0ex}{0ex}}{\int}_{{V}_{0}}{\mathrm{d}}^{3}{\mathbf{r}}_{\mathbf{0}}\mu ({\mathbf{r}}_{\mathbf{0}}){z}_{0}\frac{{e}^{-ik|\mathbf{r}-{\mathbf{r}}_{\mathbf{0}}|}}{{|\mathbf{r}-{\mathbf{r}}_{\mathbf{0}}|}^{2}},$$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 $\tilde{I}(\omega )\tilde{T}(\omega )={I}_{0}\tilde{T}(\omega )$, where ${I}_{0}$ 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 $y={y}_{0}=0$. For the discrete version, we use the notation $\mathbf{g}=\mathbf{H}\mathbf{f}$, where $\mathbf{g}$ is a vector representing the voltage measured by the transducers, $\mathbf{H}$ is the system matrix, and $\mathbf{f}$ is a vector representing the initial pressure distribution of the object being imaged. We use a $128\times 128\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{pixel}$ 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 $2272\times 16384$.

## 3.

## Image Reconstruction

## 3.1.

### Pseudoinverse

We performed an SVD on the system matrix, which was then used to calculate the pseudoinverse. The SVD factors the system matrix into $\mathbf{H}=\mathbf{U}\mathbf{S}{\mathbf{V}}^{\u2020}$, where $\mathbf{U}$ and ${\mathbf{V}}^{\u2020}$ are unitary matrices, the dagger represents conjugate transpose, and $\mathbf{S}$ is a diagonal matrix consisting of the singular values in decreasing order. The pseudoinverse is then given by ${\mathbf{H}}^{+}=\mathbf{V}{\mathbf{S}}^{+}{\mathbf{U}}^{\u2020}$, where ${\mathbf{S}}^{+}$ 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 $\mathbf{f}$, is then given by $\widehat{\mathbf{f}}={\mathbf{H}}^{+}\mathbf{g}$.

## 3.2.

### Filtered Backprojection and Simple Backprojection

The filtered backprojection formula^{5} can be written as

^{14}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}

## 3.3.

### Fourier Method

The Fourier transform-based formula^{6} can be written as

## 3.4.

### Deconvolution

For the filtered backprojection, simple backprojection, and Fourier based methods we deconvolve the transducer response using

where the asterisk indicates complex conjugation and $\zeta $ is a regularization parameter that varies from ${10}^{-6}$ to 10. For small values of $\zeta $ the deconvolution is the inverse of the transducer response and noise amplification occurs. When $\zeta $ is 10, $\tilde{D}(\omega )\approx {\tilde{T}}^{*}(\omega )$, which is a matched filter. The matched filter lowers the noise at the expense of resolution. For all of the reconstructed images we perform envelope detection by taking the absolute value of the analytic function of the image.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 $\mathbf{g}$ 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 background^{17} 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 $\mathbf{g}$ before reconstruction, where, again, the standard deviation is set to the peak signal level divided by 10.

## 4.

## Image Quality

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 $Q$^{11} and DOG channels.^{12} Constant $Q$ 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 by^{11}

## (7)

$${\mathrm{SNR}}_{\mathbf{c}\mathbf{h}}^{2}=\langle {\widehat{\mathbf{f}}}_{\mathbf{c}\mathbf{h}}^{\u2020}\rangle {\mathbf{K}}_{\mathbf{c}\mathbf{h}}^{-\mathbf{1}}\langle {\widehat{\mathbf{f}}}_{\mathbf{c}\mathbf{h}}\rangle ,$$^{17}task, the number of lumps and their locations are treated as random variables. The number of lumps, $L$, is modeled as a Gaussian random variable with a mean and variance of 1250. That number was selected using $\langle L\rangle /N=1/100$, where $N=1.25\times {10}^{5}$ 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 $\mathbf{g}$ 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

where ${\mathrm{SNR}}_{\text{data}}^{2}$ is given by## (9)

$${\mathrm{SNR}}_{\text{data}}^{2}=\langle {\mathbf{g}}^{\u2020}\rangle {\mathbf{K}}_{\mathbf{g}}^{-1}\langle \mathbf{g}\rangle .$$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.

## 5.

## Conclusion

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.

## Acknowledgments

This work was supported by ACS Research Scholar Grant RSG-08-117-01-CCE to Patrick La Rivière.

## References

C. LiL. V. Wang, “Photoacoustic tomography and sensing in biomedicine,” Phys. Med. Biol. 54(19), R59–R96 (2009).PHMBA70031-9155http://dx.doi.org/10.1088/0031-9155/54/19/R01Google Scholar

J. J. Niederhauseret al., “Combined ultrasound and optoacoustic system for real-time high-contrast vascular imaging *in vivo*,” IEEE Trans. Med. Imaging 24(4), 436–440 (2005).ITMID40278-0062http://dx.doi.org/10.1109/TMI.2004.843199Google Scholar

D. W. Yanget al., “Fast full-view photoacoustic imaging by combined scanning with a linear transducer array,” Opt. Express 15(23), 15566–15575 (2007).OPEXFF1094-4087http://dx.doi.org/10.1364/OE.15.015566Google Scholar

R. Olafssonet al., “Real-time, contrast enhanced photoacoustic imaging of cancer in a mouse window chamber,” Opt. Express 18(18), 18625–18632 (2010).OPEXFF1094-4087http://dx.doi.org/10.1364/OE.18.018625Google Scholar

D. ModgilP. J. La Rivière, “Implementation and comparison of reconstruction algorithms for two-dimensional optoacoustic tomography using a linear array,” J. Biomed. Opt. 14(4), 044023 (2009).JBOPFO1083-3668http://dx.doi.org/10.1117/1.3194293Google Scholar

K. P. KöstliP. C. Beard, “Two-dimensional photoacoustic imaging by use of Fourier-transform image reconstruction and a detector with an anisotropic response,” Appl. Opt. 42(10), 1899–1908 (2003).APOPAI0003-6935http://dx.doi.org/10.1364/AO.42.001899Google Scholar

J. Zhanget al., “Simultaneous reconstruction of speed-of-sound and optical absorption properties in photoacoustic tomography via a time-domain iterative algorithm,” Proc. SPIE 6856, 68561F (2008).PSISDG0277-786Xhttp://dx.doi.org/10.1117/12.764171Google Scholar

M. B. Roumeliotiset al., “Singular value decomposition analysis of a photoacoustic imaging system and 3D imaging at 0.7 FPS,” Opt. Express 19(14), 13405–13417 (2011).OPEXFF1094-4087http://dx.doi.org/10.1364/OE.19.013405Google Scholar

A. PetschkeP. J. La Rivière, “Photoacoustic image reconstruction using the pseudoinverse of the system matrix with the potential for real time imaging,” Proc. SPIE 8223, 82233J (2012).PSISDG0277-786Xhttp://dx.doi.org/10.1117/12.909145Google Scholar

H. H. BarrettK. J. Myers, Foundations of Image Science, pp. 913–1000, Wiley, Hoboken, NJ (2004).Google Scholar

K. J. MyersH. H. Barrett, “Addition of a channel mechanism to the ideal-observer model,” J. Opt. Soc. Am. A 4(12), 2447–2457 (1987).JOAOD60740-3232http://dx.doi.org/10.1364/JOSAA.4.002447Google Scholar

C. K. AbbeyH. H. BarrettD. W. Wilson, “Observer signal-to-noise ratios for the ML-EM algorithm,” Proc. SPIE 2712, 47–58 (1996).PSISDG0277-786Xhttp://dx.doi.org/10.1117/12.236860Google Scholar

A. PetschkeP. J. La Rivière, “Comparison of intensity-modulated continuous-wave lasers with a chirped modulation frequency to pulsed lasers for photoacoustic imaging applications,” Biomed. Opt. Express 1(4), 1188–1195 (2010).BOEICL2156-7085http://dx.doi.org/10.1364/BOE.1.001188Google Scholar

R. S. C. Cobbold, Foundations of Biomedical Ultrasound, pp. 492–607, Oxford, New York, (2006).Google Scholar

S. J. Norton, “Reconstruction of a reflectivity field from line integrals over circular paths,” J. Acoust. Soc. Am. 67(3), 853–863 (1980).JASMAN0001-4966http://dx.doi.org/10.1121/1.383964Google Scholar

P. Burgholzeret al., “Temporal back-projection algorithms for photoacoustic tomography with integrating line detectors,” Inverse Probl. 23(6), S65–S80 (2007).INPEEY0266-5611http://dx.doi.org/10.1088/0266-5611/23/6/S06Google Scholar

J. P. RollandH. H. Barrett, “Effect of random background inhomogeneity on observer detection performance,” J. Opt. Soc. Am. A 9(5), 649–658 (1992).JOAOD60740-3232http://dx.doi.org/10.1364/JOSAA.9.000649Google Scholar

*
*Citation
Download Citation
Adam Petschke,
Patrick J. La Rivière,
"Comparison of photoacoustic image reconstruction algorithms using the channelized Hotelling observer," Journal of Biomedical Optics 18(2), 026009 (5 February 2013). http://dx.doi.org/10.1117/1.JBO.18.2.026009
Submission: Received ; Accepted