1.

Introduction

Photoacoustic imaging (PAI) is cross-sectional or three dimensional (3D) imaging based on the photoacoustic effect. PAI combines the advantages of optical imaging contrast with ultrasonic spatial resolution. In Ref. 1 the reader can find an overview of current state-of-the-art imaging modalities.

The combination of accurate detectors and efficient reconstruction algorithms is crucial for a successfully operating PAI system. Frequency domain algorithms are computationally efficient because they usually rely on an implementation with the fast Fourier transform (FFT). The standard FFT algorithm has the major drawback that it requires a uniform sampling of data points, while the backprojection formula for the photoacoustic inverse problem, described below, requires a nonuniform sampling. Several proposals exist to overcome this problem. In this work we focus on a novel inversion scheme based on the nonuniform FFT (NUFFT) as described in Ref. 2. The NUFFT has been introduced in Ref. 3, and the combination with the inverse Radon transform for solving tomographic problems was described in Ref. 4.

This paper is organized as follows: In Sec. 2 we give an overview of our experimental setup from which the data was obtained. Then we review the mathematics of the NUFFT reconstruction algorithm in Sec. 3. The most relevant section is Sec. 4, where we investigate how imaging properties are affected by the use of the novel algorithm. There, we show its accuracy when applied to clinically and biologically relevant examples. Finally, in Sec. 5 we conclude our work by a final discussion of the investigated algorithms.

2.

Experimental Setup

In this section we give a short description of our experimental setup used to acquire the data for our studies. Our setup is equipped with an integrating line detector for acquiring the measurement data (see Ref. 5 for a description of the setup). Integrating line detectors increase the image quality with respect to the resolution and occurring reconstruction artifacts, especially if a free-beam interferometric setup is used.^{6, 7} The integrating detector used in the presented studies is realized with a free-beam Mach–Zehnder interferometer (MZI). Optical beams as part of an interferometer have omnidirectional response,⁷ but need focusing in order to achieve high temporal and spatial resolution. Our focus was well below 60 μm along the focal length. Finally, the actual resolution limit is determined by the sampling density, the involved electronics, and penetration depth. Typical values are 80 to 120 μm.⁸

The setup is shown in Fig. 1. The emerging ultrasound waves propagate through a tank filled with water which acts as a coupling medium between the sample and the detector. The pressure wave leads to a change in the refractive index of the water and therefore causes a phase shift in the measurement beam of the MZI. The interferometer is actively stabilized by a piezocontrolled mirror, which counteracts slow oscillations as well as long-term drift. Since the optical components have to be stabilized, it is necessary to move the sample relative to the beam. As a pump laser, we used an optical parametric oscillator pumped by an Nd:YAG laser system (Surelite II-10, Continuum, Santa Clara, California).

Fig. 1

Schematic of experimental setup. Ultrasound waves are detected with a free-beam MZI. BS: beam splitter, M: mirror, M+PT: piezo-stabilized mirror, L: lens, BPF: band-pass filter, PD: photodiodes, HPF: high-pass filter.

Data acquisition for a complete 3D dataset consists of two steps: First, the sample is moved to a fixed detector position and rotated about 360 deg. The photoacoustic transient is recorded for each laser pulse at equidistant angular positions. Having recorded the signals for all angles and detector positions, a 3D dataset can be generated as described below.

3.

Mathematical Review

In order to describe the data acquisition process outlined above, we assume that a Cartesian coordinate system is placed beside an object that is rotated by an angle α. Furthermore, we assume that the linear integrating detector moves along the y −axis as described in Fig. 2. If p _α denotes the pressure that is emitted by the rotated object, described by p _0α, the measurement data are given by

Fig. 2

Left: Pressure integrals over lines parallel to the z −axis running through (0, y _d, 0) are collected from an object that has been rotated by an angle α. Right: Rotating the parallel family of lines instead of the object yields the Radon transform in each plane [TeX:] \documentclass[12pt]{minimal}\begin{document}$\mathbb {R}\times \lbrace y\rbrace$\end{document} $\mathbb{R}\times \left\{y\right\}$ .

Inverting the two-dimensional Radon transform in each plane gives the desired 3D reconstruction. In the following we will focus on the reconstruction of a single projection and thus will omit α in the following.

The first FFT-based inversion scheme for a similar photoacoustic inverse problem using point detectors was proposed in Ref. 9. For our case the backprojection identity reads

It seems a straightforward step to apply the FFT to the data, scale it back in frequency space using Eq. 3, and apply the FFT again. However, for equally spaced samples of the frequency, the use of this dispersion relation implies that k _y are sampled nonequidistant. In other words, the backprojection in frequency domain via the dispersion relation requires nonuniform sampling of the frequency data points used in the inverse transform. In order to apply FFT algorithms, which assume sampling on an equidistant grid, the data points have to be resampled via interpolation before transforming back to position space. As discussed in all the works cited above, interpolation in Fourier space causes considerable artifacts.

One can avoid this problem by using the discrete Fourier transform (DFT)

The first case is called nonequispaced data and the second case is called nonequispaced results. As discussed in the reference, the second option is easier to implement and provides a faster algorithm.

In the following we will summarize the results obtained in Ref. 2. Let [TeX:] \documentclass[12pt]{minimal}\begin{document}$\Psi : \mathbb {R}\rightarrow \mathbb {R}$\end{document} $\Psi :\mathbb{R}\to \mathbb{R}$ be an appropriate function that satisfies:

and let [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\Psi }:=\mathcal F\lbrace \Psi \rbrace$\end{document} $\widehat{\Psi }:=\mathcal{F}\left\{\Psi \right\}$ be its Fourier transform. Ψ is called window function. With this definition, we are ready to explain the NUFFT algorithm, which is based on the identity

The window function [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hat{\Psi }$\end{document} $\widehat{\Psi }$ has to be chosen such that it is concentrated around zero (and decaying rapidly away from zero). Then in the summation on the right-hand side of Eq. 5 only a few terms have to be taken into account. This is the main approximation in NUFFT-based algorithms. According to the above references, the Kaiser–Bessel window function is the best choice such that the approximation error is small. It is also important to note that we impose [TeX:] \documentclass[12pt]{minimal}\begin{document}$cN\in \mathbb {N}$\end{document} $cN\in \mathbb{N}$ , which allows us to expand the summation in Eq. 6 to [0, cN − 1]. Then, the equation represents an oversampled discrete Fourier transform. This allows us to apply the standard FFT in this step of the reconstruction algorithm simply by appending (c − 1)N zeros to the data vector g _n. We found that c can be chosen smaller than the typical value of 1.5 that has been used in Ref. 4. We chose c = 38/32 for all our reconstruction examples. Increasing the oversampling does not improve the image quality significantly.

For further analysis we refer the interested reader to the publication by Haltmeier ²

4.

Applications

The improvements achieved by NUFFT reconstruction are demonstrated below using several biological and clinical examples. All data were generated with the PAI system described above. We examined:

Only detector array geometries which fulfill the conditions of the theoretical derivations are stable. Unstable reconstructions suffer from artifacts in the projection images. Mainly, pressure values below zero are generated in the final image and appear either as shadow-like artifacts or an inhomogeneous background. For the application of FFT-based back projection algorithms as described above, this has the implication that the detector array positions must lie on an infinite straight line, as Eq. 3 has been derived under this condition. Of course, infinite detector arrays do not exist in reality, even if approximating by an array length much larger than the sample size produces some artifacts. This is known as the limited view problem. Therefore, it was suggested to make two linear scans perpendicular to each other such that the detector array encloses the object with a triangle.⁶ Even better is a box-scan, consisting of three linear scans (in the directions upward, forward, and downward the sample).¹⁰ The two or three reconstructed images are then superimposed. We used the box scan detector configuration in all examples discussed below. See also Fig. 3.

Fig. 3

Sketch of our scanning geometries. Depicted are an L-scan geometry (left) and a box-scan geometry (right). The shaded gray areas represent the set of coordinates where a stable reconstruction is possible.

All source codes were written in MATLAB (MathWorks, Natick, Massachusets), and the reconstructions were processed on a standard 64-bit personal computer. As discussed in Ref. 2, the overall complexity of the algorithm is [TeX:] \documentclass[12pt]{minimal}\begin{document}${\mathcal O}(cN^2\log cN)$\end{document} $\mathcal{O}\left(c{N}^2\log cN\right)$ , provided that the window functions Ψ are precomputed. The computational effort for the window functions depends on the choice of the oversampling c and the number of terms retained in the sum 6. Let M, a small integer number, denote the number of terms retained in the sum. Then the complexity of the precomputations is [TeX:] \documentclass[12pt]{minimal}\begin{document}${\mathcal O}(cNM)$\end{document} $\mathcal{O}\left(cNM\right)$ . For example, computation of a 300×300 px image using the FFT routine took 0.115 s. The NUFFT routine took 0.120 s without, and 0.252 s with precomputations. Notice that without precomputations, runtimes are comparable.

As described in the following paragraphs, we could show that the image quality obtained via the NUFFT algorithm is superior to the FFT-reconstructed images.

4.1.

Resolution and Image Quality

As an important remark we point out that the experimental resolution limit should not be affected by the use of different algorithms. We test the resolution limit in the two images to outrule the possibility that either Fourier domain algorithm causes smearing artifacts.

We measured clinical sutures with a diameter of 50 μm. This is below the resolution limit of our system, and, therefore, we can obtain a resolution limit by fitting a Gaussian distribution to the gray value profile of the sutures. In order to resolve the Gaussian, the pixel (voxel) size has to be significantly smaller than the anticipated spatial resolution limit. In this case, the pixel size was 27 μm. Commonly, one uses full-width-half-maximum (FWHM) of the Gaussian to quantize the resolution performance of an imaging system. Therefore, the fitting function ρ(x) reads

Eq. 7

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} \rho (x)=A\exp \left[-\frac{4\log (2)(x-x_0)^2}{x_{\rm FHWM}^2}\right]. \end{equation} \end{document} $$\rho \left(x\right)=A\exp \left[-\frac{4\log \left(2\right){(x-x_0)}^2}{x_{\mathrm{FHWM}}^2}\right].$$

The signals of two sutures lying side by side are depicted in Fig. 4. A curve fit has been applied to the normalized FFT and NUFFT gray value profiles. The fit to the FFT data gives 125 and 138 μm for the left and right suture, respectively, giving an average of 131.5 μm at FWHM. The fit to the NUFFT data gives 121 and 131 μm for the left and right suture, respectively, giving an average of 126 μm at FWHM. As the suture does not imitate a perfect point-spread function, i.e., a Dirac δ-function, the actual resolution limit is below the above values. Note that subtraction of the suture diameter as was done for example in Ref. 11 does not give a correct resolution performance of an imaging system in general since the finite suture diameter and detector size lead to a deconvolution problem.

Fig. 4

Normalized gray value profiles and Gaussian fits of two sutures (diameter 50 μm). On average, the resolution limits are 81.5 μm (FFT algorithm) and 76 μm (NUFFT algorithm). For a confidence level of 95%, the goodness of the fits is given by R ² = 0.871 (FFT algorithm) and R ² = 0.887 (NUFFT algorithm). We also show the FWHM values as horizontal lines (NUFFT: dashed line; FFT: dashed-dotted line).

Differences in resolutions using different reconstruction algorithms typically range about 10 μm as can be seen in Ref. 6, page S90]. Using 95% confidence bounds gives a confidence interval of 114 to 150 μm for the FFT data and 108 to 144 μm for the NUFFT data. We conclude that the resolution limit is not affected significantly by the use of either algorithm as expected.

In order to investigate the image quality, we used a phantom consisting of sutures lying close to each other. The sutures had diameters 20, 50 and 70 μm. The images are displayed in Fig. 5. Here, we observe the important fact that images reconstructed with the FFT algorithm tend to be biased toward the detector (located at the upper edge of the images). The bias is due to two reasons in combination: First, it is never possible to record all pressure data. On the other hand, the reconstruction algorithm assumes time data from zero to infinity. The timescale of the oscilloscope is always chosen as the best compromise between resolution and best approximation of an infinite time axis. Second, the FFT algorithm neglects the nonuniform sampling of data points in the reconstruction process. Together, this causes the aforementioned bias. The error varies with the frequencies occurring in the reconstruction process and cannot be given analytically. As a consequence, the 50 μm suture, which was tied to a ribbon, is not resolved completely by the FFT algorithm. The lower part of the ribbon is missing.

Fig. 5

Sutures of diameters 20, 50, and 70 μm reconstructed (a) with the FFT algorithm and (b) with the NUFFT algorithm. Observe that the suture bound to a ribbon (diameter 50 μm) is clearly visible in the NUFFT image, whereas its lower part is missing in the FFT image.

Another important aspect is that the shadow-like artifacts caused by the FFT algorithm are eliminated by the use of the NUFFT algorithm. Clearly, some of these artifacts in the FFT reconstructions are also caused by the limited view problem. Thus, we observe that the NUFFT algorithm is more impervious to artificial reflections that correspond to the limited view problem.

We also determined the signal-to-noise ratio (SNR) in both images. As reference, we chose the signal maximum of the 70 μm suture. We obtained an SNR of ∼43 for the FFT image and ∼42 for the NUFFT image. This clearly shows that the SNR itself is not improved by the use of the NUFFT algorithm. In image analysis the contrast-to-noise ratio (CNR) is important also. It is given by

Eq. 8

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} \text{CNR} = \frac{\langle S_s\rangle -\langle S_b\rangle }{\sigma }, \end{equation} \end{document} $$\text{CNR}=\frac{⟨S_s⟩-⟨S_b⟩}{\sigma },$$

where ⟨S _s⟩ is the mean value of the (suture) signal, ⟨S _b⟩ is the mean value of the background noise, and σ its standard deviation.¹¹ The resulting CNR’s were ∼34 for the FFT algorithm and ∼38 for the NUFFT algorithm. The quantities are again comparable, but it seems that using the NUFFT algorithm improves the CNR slightly.

We also computed the local noise power spectrum for suitable regions of interest in our images. We found that the use of the different algorithms again causes no significant differences.

4.2.

Test on Biological Samples

Let us begin by explaining the methods used to obtain our samples. Tumor bearing mouse mammary tumor virus (MMTV)-neu mice were obtained from the Central Laboratory Animal Facility of the Innsbruck Medical University (Innsbruck, Austria). Female individuals serve as a model system for breast cancer and develop one or more tumors along the mammary ridges at an age of about 5 to 7 months. A detailed description of the mouse strain can be found in Guy ¹² and Parajuli and Doppler.¹³ Mice were euthanized by CO₂ and sacrificed by cervical dislocation. The tumors were excised and transferred into formaldehyde immediately for conservation. For measurements, the tumors were embedded into agar and provided with a stick to enable mounting on the translational stage. All animals were treated in accordance with the Austrian animal welfare law and animal experiment act.

For the results, see Fig. 6. Notice that the blurring is significantly reduced by the use of the NUFFT algorithm. The detector positions lie along the top, left, and right margins of the image. Since the blurring only occurs along the vertical axis, we conclude that the limited view effect caused the blurring in the FFT image as the detection curve does not enclose the bottom of the image. We observe again that the NUFFT algorithm does not seem to be affected by this problem so severely.

Fig. 6

Box-scan projection images of an ex vivo mouse carcinoma. (a) FFT reconstruction. (b) NUFFT reconstruction. The background in the FFT image appears blurred, whereas in the NUFFT image it is more homogeneous.

A main motivation for improving the reconstruction algorithm was to increase resolution of small structures, such as zebrafish embryos. In the first place, the weak signals obtained from 2 day old zebrafish embryos could not be resolved in accurate detail (Fig. 7). The images in the first row represent the first sample, embedded in agar and measured with a box side length of 19 mm. The images in the second row represent the second sample, which was measured with a length of 13 mm because the agar cylinder could be reduced in diameter for this sample. We see the consequences of two important effects: The high frequencies caused by the small structures suffer from ultrasound attenuation in water. This is a well-known effect.¹⁴ Therefore, reducing the box side length and thus reducing the probe-detector distance considerably increased the image quality. Again, the homogeneity of the NUFFT images is better than that of the regular FFT images. The most important observation is that the embryos further from the detector (lowest regions in the images) appear blurred and distorted in the FFT images. In the final configuration, i.e., 13 mm box length and use of the NUFFT algorithm, the reconstructed image corresponds very well with the microscopy image.

Fig. 7

Box-scan of 2 day old zebrafish embryos. Two samples: sample 1 measured with box length 19 mm, sample 2 with 13 mm. (a) Microscope image of sample 1, (b) reconstructed with FFT algorithm, (c) reconstructed with NUFFT algorithm. (d) Microscope image of sample 2, (e) reconstructed with FFT algorithm, (f) reconstructed with NUFFT algorithm. Zebrafish embryos were grown as described in Ref. 15. Scanning was performed with paraformaldehyde-fixed embryos [(a)–(c)] or life embryos [(d)–(f)] embedded in 1.5% low melt agarose.

5.

Conclusion and Outlook

In summary, we have shown that the nonuniform FFT reconstruction algorithm indeed provides superior image quality when compared to the standard FFT algorithm. This concerns artifacts, as well as background gray values, which appear more homogeneous for the NUFFT images. We found that resolution and SNR (CNR) are not significantly improved by the use of the NUFFT over the standard FFT. We observed that the most capital improvement of the NUFFT reconstruction is its capability to reliably resolve image details that are farther away from the detector than about half the total image distance. In Table 1 we give the pros and cons for the two Fourier algorithms in summary.

Table 1

Pros and cons of the investigated frequency domain algorithms.