## 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,^{7} 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.^{8}

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).

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

## Eq. 1

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} q(y_d,t) = \int _\mathbb {R}p_\alpha (0,y_d,z,t)dz. \end{equation} \end{document} $$q({y}_{d},t)={\int}_{\mathbb{R}}{p}_{\alpha}(0,{y}_{d},z,t)dz.$$## Eq. 2

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} q_{0 \alpha }(x,y) = \int _{\mathbb {R}} p_{0\alpha } (x,y,z) dz \end{equation} \end{document} $${q}_{0\alpha}(x,y)={\int}_{\mathbb{R}}{p}_{0\alpha}(x,y,z)dz$$*x*,

*y*) that are parallel to the

*z*−axis. Since it would be equivalent to fix the object and rotate the family of lines, instead the projections (

*q*

_{0α}) provide 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\}$ , see Fig. 2.

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

## Eq. 3

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} Q_0(k_x,k_y)=\frac{2ck_y}{\text{sign}(k_y)\sqrt{k_x^2+k_y^2}}Q\left(k_x,c\cdot \text{sign}(k_y)\sqrt{k_x^2+k_y^2}\right), \end{equation} \end{document} $${Q}_{0}({k}_{x},{k}_{y})=\frac{2c{k}_{y}}{\text{sign}\left({k}_{y}\right)\sqrt{{k}_{x}^{2}+{k}_{y}^{2}}}Q\left({k}_{x},c\xb7\text{sign}\left({k}_{y}\right)\sqrt{{k}_{x}^{2}+{k}_{y}^{2}}\right),$$*Q*(

*k*

_{x}, ω) is the Fourier transform of the measured data

*q*(

*x*

_{d}, ·,

*t*) and

*Q*

_{0}(

*k*

_{x},

*k*

_{y}) is the Fourier transform of the projections

*q*

_{0}(

*x*,

*y*). In order to derive this, the dispersion relation [TeX:] \documentclass[12pt]{minimal}\begin{document}$k_y=\sqrt{(\omega /c)^2-k_x^2})$\end{document} ${k}_{y}=\sqrt{{(\omega /c)}^{2}-{k}_{x}^{2}})$ was used.

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)

## Eq. 4

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} T[\mathbf g](\omega _k):=\sum _{n=0}^{N-1}e^{-i\omega _kn2\pi /N}g_n, \end{equation} \end{document} $$T\left[\mathbf{g}\right]\left({\omega}_{k}\right):=\sum _{n=0}^{N-1}{e}^{-i{\omega}_{k}n2\pi /N}{g}_{n},$$_{k}=

*k*ω as is usually done. Therefore, this sum allows the nodes [TeX:] \documentclass[12pt]{minimal}\begin{document}$(\omega _k)_{k=N/2}^{N/2-1}$\end{document} ${\left({\omega}_{k}\right)}_{k=N/2}^{N/2-1}$ to be sampled on an arbitrary grid. Evaluating the sums in Eq. 4 requires [TeX:] \documentclass[12pt]{minimal}\begin{document}$\mathcal O(N^2)$\end{document} $\mathcal{O}\left({N}^{2}\right)$ operations, which of course does not lead to a fast algorithm. When using the NUFFT algorithm as presented in Ref. 4, we have two choices:

apply standard FFT to the data, then do the inverse transform on a nonuniform grid (via the NUFFT),

apply NUFFT to the data at nonequidistant grid points in order to get nonequidistant results in frequency domain, then do the equidistant inverse Fourier transform.

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

## Eq. 5

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} \sum _{n=0}^{N-1} e^{-i \omega n 2\pi /N} g_n = \sum _{j \in \mathbb {Z}} e^{- i \pi (\omega - j/c)} \hat{\Psi }(\omega - j/c) \hat{G}_j \,, \end{equation} \end{document} $$\sum _{n=0}^{N-1}{e}^{-i\omega n2\pi /N}{g}_{n}=\sum _{j\in \mathbb{Z}}{e}^{-i\pi (\omega -j/c)}\widehat{\Psi}(\omega -j/c){\widehat{G}}_{j}\phantom{\rule{0.16em}{0ex}},$$## Eq. 6

[TeX:] \documentclass[12pt]{minimal}\begin{document} \begin{equation} \hat{G}_j := \frac{c}{2\pi } \left[ \sum _{n =0}^{N-1} \frac{g_n e^{-i j n 2\pi /(Nc)}}{ \Psi (n 2\pi /N - \pi)} \right] \,, \qquad j \in \mathbb {Z}\,, \end{equation} \end{document} $${\widehat{G}}_{j}:=\frac{c}{2\pi}\left[\sum _{n=0}^{N-1}\frac{{g}_{n}{e}^{-ijn2\pi /\left(Nc\right)}}{\Psi (n2\pi /N-\pi )}\right]\phantom{\rule{0.16em}{0ex}},\phantom{\rule{2.em}{0ex}}j\in \mathbb{Z}\phantom{\rule{0.16em}{0ex}},$$*c*> 1 is an oversampling factor, α < π(2

*c*− 1) is the window size, and [TeX:] \documentclass[12pt]{minimal}\begin{document}$(g_n)_{n=0}^{N-1} \in \mathbb {C}^N$\end{document} ${\left({g}_{n}\right)}_{n=0}^{N-1}\in {\mathbb{C}}^{N}$ is the vector containing the data. Equation 5 expresses the Fourier transform of the signal

*g*

_{n}evaluated at an arbitrary frequency ω.

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

## 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:

sutures in order to test whether the use of the algorithm influences the image resolution limit

mouse mammary-carcinoma

zebrafish embryos 2 days-post-fertilization.

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.^{6} Even better is a box-scan, consisting of three linear scans (in the directions upward, forward, and downward the sample).^{10} 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.

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}\mathrm{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\mathrm{exp}\left[-\frac{4\mathrm{log}\left(2\right){(x-{x}_{0})}^{2}}{{x}_{\mathrm{FHWM}}^{2}}\right].$$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.

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{\u27e8{S}_{s}\u27e9-\u27e8{S}_{b}\u27e9}{\sigma},$$*S*

_{s}⟩ is the mean value of the (suture) signal, ⟨

*S*

_{b}⟩ is the mean value of the background noise, and σ its standard deviation.

^{11}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
^{12} and Parajuli and Doppler.^{13} Mice were euthanized by CO_{2} 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.

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

## 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.

Standard FFT algorithm | |
---|---|

Pro | Con |

No pre-computations | |

Distortions and shadow-artifacts | |

biased towards detector array | |

NUFFT algorithm | |

Pro | Con |

Pre-computations necessary | |

Reduction of artifacts | |

homogeneous reconstruction area |

## Acknowledgments

This work has been supported by the Austrian Science Fund, Project Nos. S 10502-N20, S 10505-N20, S 10508-N20, and S 10509-N20. We thank Sonja Töchterle from the University of Innsbruck for her help with the zebrafish samples.

## References

**,” Med. Phys., 35 (12), 5758 –5767 (2008). https://doi.org/10.1118/1.3013698 Google Scholar**

*Prospects of photoacoustic tomography***,” IEEE Trans. Med. Imaging, 28 (11), 1727 –1735 (2009). https://doi.org/10.1109/TMI.2009.2022623 Google Scholar**

*A reconstruction algorithm for photoacoustic imaging based on the nonuniform FFT***,” SIAM J. Sci. Comput. (USA), 14 (6), 1368 –1393 (1993). https://doi.org/10.1137/0914081 Google Scholar**

*Fast Fourier transforms for nonequispaced data***,” J. Fourier Analy. Appl., 9 (5), 431450 (2003). https://doi.org/10.1007/s00041-003-0021-1 Google Scholar**

*Non-equispaced fast Fourier transforms with applications to tomography***,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control, 52 (9), 1577 –1583 (2005). https://doi.org/10.1109/TUFFC.2005.1516030 Google Scholar**

*Thermoacoustic tomography with integrating area and line detectors***,” Inverse Prob. Eng., 23 S81 –S64 (2007). https://doi.org/10.1088/0266-5611/23/6/S07 Google Scholar**

*Experimental evaluation of reconstruction algorithms for limited view photoacoustic tomography with line detectors***,” Proc. SPIE, 7177 71770T (2009). https://doi.org/10.1117/12.808873 Google Scholar**

*Comparison of optical and piezoelectric integrating line detectors***,” J. App. Phy., 105 (10), 102026 (2009). https://doi.org/10.1063/1.3116133 Google Scholar**

*Characterization of integrating ultrasound detectors for photoacoustic tomography***,” IEEE Trans. Med. Imag., 21 (7), 823828 (2002). https://doi.org/10.1109/TMI.2002.801172 Google Scholar**

*Exact frequency-domain reconstruction for thermoacoustic tomography I: Planar geometry***,” J. Biomed. Opt., 16 036007 (2011). https://doi.org/10.1117/1.3556720 Google Scholar**

*Photoacoustic tomography of**ex vivo*mouse hearts with myocardial infarction**,” Opt. Express, 17 (24), 21414 –21426 (2009). https://doi.org/10.1364/OE.17.021414 Google Scholar**

*Multispectral optoacoustic tomography (MSOT) scanner for whole-body small animal imaging***,” Proc. Nati. Acad. Sci. U.S.A., 89 10578 –10582 (1992). https://doi.org/10.1073/pnas.89.22.10578 Google Scholar**

*Expression of the neu protooncogene in the mammary epithelium of transgenic mice induces metastatic disease***,” In Vitro Cell Dev. Biol: Anim., 45 442 –450 (2009). https://doi.org/10.1007/s11626-009-9212-7 Google Scholar**

*Precision-cut slice cultures of tumors from MMTV-neu mice for the study of the ex vivo response to cytokines and cytotoxic drugs***,” Opt. Lett., 31 (6), 781 –783 (2006). https://doi.org/10.1364/OL.31.000781 Google Scholar**

*Image reconstruction in optoacoustic tomography for dispersive acoustic media*