A Reconstruction Algorithm for Photoacoustic Imaging based on the Nonuniform FFT

Photoacoustic imaging (PAI) is a novel promising tool for visualizing light absorbing structures in an optically scattering medium, which carry valuable information for medical diagnostics. It is based on the generation of acoustic waves by illuminating an object with pulses of non-ionizing electromagnetic radiation, and combines the high contrast of pure optical and the high resolution of ultrasonic imaging. The method has demonstrated great promise for a variety of biomedical applications, such as imaging of animals [1, 2], early cancer diagnostics [3, 4], and imaging of vasculature [5, 6].

When an object is illuminated with short pulses of non-ionizing electromagnetic radiation, it absorbs a fraction of energy and heats up. This in turn induces acoustic (pressure) waves, that are recorded with acoustic detectors outside of the object. Other than in conventional ultrasound imaging, where the source of acoustic waves is an external transducer, in PAI the source is the imaged object itself. The frequency bandwidth of the recorded signals is therefore generally broad and depends on the size and the shape of illuminated structures.

Let $C_{0}^{\infty}(\mathbb{H})$ denote the space of smooth functions with bounded support in the half space $\mathbb{H}:=\mathbb{R}^{d-1}\times(0,\infty)$, $d\geq 2$. Consider the initial value problem

with $f\in C_{0}^{\infty}(\mathbb{H})$. Here $\Delta$ denotes the Laplacian in $\mathbb{R}^{d}$ and $\partial_{t}$ is the derivative with respect to $t$. We write $\mathbf{x}=(x,y)$, $x\in\mathbb{R}^{d-1}$, $y\in\mathbb{R}$, and define the operator $\operatorname{\mathbf{Q}}:C_{0}^{\infty}(\mathbb{H})\to C^{\infty}(\mathbb{R}^{d})$ by

Photoacoustic imaging for planar recording geometry is concerned with reconstructing $f\in C_{0}^{\infty}(\mathbb{H})$ from incomplete and possibly erroneous knowledge of $\operatorname{\mathbf{Q}}f$. Of practical interest are the cases $d=2$ and $d=3$, see [30, 31, 32, 33].

The discrete Fourier transform of a vector $\mathbf{g}=(g_{n})_{n=0}^{N-1}\in\mathbb{C}^{N}$ with respect to the nodes $\boldsymbol{\omega}=(\omega_{k})_{k=-N/2}^{N/2-1}$ (with $N$ even) is defined by

Direct evaluation of the $N$ sums in (4) requires $\mathcal{O}(N^{2})$ operations. Using the classical fast Fourier transform (FFT) this effort can be reduced to $\mathcal{O}(N\log N)$ operations. However, application of the classical FFT is restricted to the case of evenly spaced nodes $\omega_{k}=k$, $k=-N/2,\dots,N/2-1$.

The one-dimensional nonuniform FFT (see [18, 19, 20, 29, 21, 22, 23]) is an approximate but highly accurate method for evaluating (4) at arbitrary nodes $\omega_{k}$, $k=-N/2,\dots,N/2-1$, in $\mathcal{O}(N\log N)$ operations.

In this section we apply the one-dimensional nonuniform FFT to photoacoustic imaging. Throughout the following we restrict our attention to two dimensions, noting that the general case $d\geq 2$ can be treated in an analogous manner.

Assume that $f$ is a smooth function that vanishes outside $(0,X)^{2}$, and set $g:=w_{\rm cut}\operatorname{\mathbf{Q}}f$, where $w_{\rm cut}$ is as in (2). Fourier reconstruction names an implementation of (3), that uses discrete data

with $(m,n)\in\{0,\dots,N-1\}^{2}$ and reconstructs an approximation

with $(m,n)\in\{0,\dots,N-1\}^{2}$. Here $f^{\dagger}$ is defined by (3), $N$ is an even number and $\Delta_{\rm samp}:=X/N$. In the Appendix we show that the sampling in (10), (9) is sufficiently fine, provided that $\Delta_{\rm samp}\leq\pi/\Omega$, where $\Omega$ is the essential bandwidth of $f$.

Discretizing (3) with the trapezoidal rule gives

where

One notices that the inner sums in (11),

can be exactly evaluated with $N$ one-dimensional FFTs, and the outer sums

can be approximately evaluated with $N$ one-dimensional nonuniform FFTs. Denoting the resulting approximation by $\hat{g}_{k,l}\simeq\hat{g}_{k}(\omega_{k,l})$ and setting

we finally find

with the inverse two-dimensional FFT algorithm.

The nonuniform FFT based reconstruction algorithm is summarized in Algorithm 2. Its numerical complexity can easily be estimated. Evaluating (12) requires $N\mathcal{O}(N\log N)$ operations ($N$ one-dimensional FFTs), evaluating (13) requires $N\mathcal{O}(N\log N)$ operations ($N$ non-uniform FFTs), and (15) is evaluated with the inverse two-dimensional FFT in $\mathcal{O}(N^{2}\log N)$ operations. Therefore the overall complexity of Algorithm 1 is $\mathcal{O}(N^{2}\log N)$.

In the next section we numerically compare Algorithm 2 with standard Fourier algorithms presented in the literature [41, 13], which all differ in the way how the sums in (13) are evaluated:

1. 1.

   Direct Fourier algorithm. Equation (13) cannot be evaluated with the classical FFT algorithm because the nodes $\omega_{k,l}$ are non-equispaced. The most simple way to evaluate (13) is with direct summation. Because there are $N^{2}$ such sums, direct Fourier reconstruction requires $\mathcal{O}(N^{3})$ operations. Consequently it does not lead to a fast algorithm. However, since (13) is evaluated exactly, it is optimally suited to evaluate the image quality of reconstructions with fast Fourier algorithms.

2. 2.

   Interpolation based algorithm. A fast and simple alternative to direct Fourier reconstruction is as follows: Choose an oversampling factor $c\geq 1$ and exactly evaluate

   $$\hat{g}_{k}(\omega):=\Delta_{\rm samp}\sum_{n=0}^{N-1}e^{-i\omega n2\pi/N}\tilde{g}_{k,n}\,,$$

   at the uniformly spaced nodes $\omega=\Delta_{\rm samp}j/c$, $j\in\{0,\dots,Nc-1\}$, with the one-dimensional FFT algorithm. In a next step, linear interpolation is used to find approximate values $\hat{g}_{k,l}\simeq\hat{g}_{k}(\omega_{k,l})$, see [13]. Evaluating $\hat{g}_{k,l}$ with linear interpolation requires $\mathcal{O}(N^{2})$ operation and therefore the overall numerical effort of linear interpolation based Fourier reconstruction is $\mathcal{O}(N^{2}\log N)$.

   Algorithms using nearest neighbor interpolation instead of the linear one have the same numerical complexity and have also been applied to PAI (see, e.g., [42]). Higher order polynomial interpolation has been applied in [43] for a cubic recording geometry.

3. 3.

   Truncated $\operatorname{sinc}$ reconstruction. If the function $\Psi$ in Algorithm 2 is chosen as the characteristic function of the interval $[-c\pi,c\pi]$, $c\geq 1$, then the nonuniform fast Fourier transform reduces to the truncated $\operatorname{sinc}$ interpolation considered in [41]. However, due to the slow decay of $\operatorname{sinc}(\omega)$, truncation will introduce a non-negligible error in the reconstructed image (see Remark 3).

The Fourier algorithms are also be compared with a numerical implementation of the back-projection formula

where $r=\sqrt{(x-x^{\prime})^{2}+y^{2}}$ denotes the distance between the detector location $(x^{\prime},0)$ and the reconstruction point $(x,y)$. Equation (16) has been obtained in [8] by applying the method of descent to the three-dimensional universal back-projection formula discovered by Xu and Wang [10]. Again (16) gives an exact reconstruction only if it is applied to complete data $(\operatorname{\mathbf{Q}}f)(x,t)$, $(x,t)\in\mathbb{R}^{2}$. In the numerical experiments the back-projection formula is applied to the partial data $w_{\rm cut}\operatorname{\mathbf{Q}}f$, see (2), and implemented with $\mathcal{O}(N^{3})$ operation counts as described in [8, Section 3.3].

In the following we numerically compare the proposed nonuniform FFT based algorithm with standard Fourier algorithm and the back projection algorithm based on (16).

The cutoff function $w_{\rm cut}$ is constructed by convolution of

with the characteristic function of $[0,X]^{2}$, where $\epsilon$ is a small parameter and $C_{\epsilon}$ is chosen in such a way that $\int_{\mathbb{R}^{2}}\varphi_{\epsilon}(x,t)dxdt=1$. Typically, $\epsilon$ is chosen as a “small” multiple of the sampling step size $\Delta_{\rm samp}=X/N$.

In all numerical experiments we take $X=1$, and $N=512$. The window width $\alpha$ is chosen to be slightly smaller than $\pi(2c-1)$, where $c$ is the oversampling factor that determines the accuracy of the Fourier reconstruction algorithms.

We presented a novel fast Fourier reconstruction algorithm for photoacoustic imaging using a limited planar detector array. The proposed algorithm is based on the nonuniform FFT. Theoretical investigation as well as numerical simulations show that our algorithm produces better images than existing Fourier algorithms with a similar numerical complexity. Moreover the proposed algorithm has been shown to be stable against data perturbations.

[Sampling and Resolution]

Let $f$ be smooth function that vanishes outside $(0,X)^{d}$, and define $g$, $f^{\dagger}$ by (2), (3). We further assume that $\operatorname{\mathbf{F}}w_{\rm cut}$ is concentrated around zero and, that $f$ is essentially bandlimited with essential bandwidth $\Omega$, in the sense that $(\operatorname{\mathbf{F}}f)(\mathbf{K})$ is negligible for $|\mathbf{K}|\geq\Omega$. Note that since $f$ has bounded support, $\operatorname{\mathbf{F}}f$ cannot vanish exactly on $\{|\mathbf{K}|\geq\Omega\}$.

• •

  Sampling of $g$. Equation (1) implies that

  $$(\operatorname{\mathbf{F}}g)(K_{x},\omega)=\left(\operatorname{\mathbf{F}}w_{\rm cut}\right)\ast\left(\operatorname{\mathbf{F}}\operatorname{\mathbf{Q}}f\right)(K_{x},\omega)\,,$$

  (17)

  with

  $$(\operatorname{\mathbf{F}}\operatorname{\mathbf{Q}}f)(K_{x},\omega)=\frac{2\omega(\operatorname{\mathbf{F}}f)\left(K_{x},\operatorname{sign}(\omega)\sqrt{\omega^{2}-|K_{x}|^{2}}\right)}{\operatorname{sign}(\omega)\sqrt{\omega^{2}-|K_{x}|^{2}}}$$

  if $|\omega|>|K_{x}|^{2}$ and $(\operatorname{\mathbf{F}}\operatorname{\mathbf{Q}}f)(K_{x},\omega)=0$ otherwise. The assumption that $f$ has essential bandwidth $\Omega$ and equation (17) imply that $(\operatorname{\mathbf{F}}g)(K_{x},\omega)$ is negligible outside the set

  $$\{(K_{x},\omega):|K_{x}|\leq|\omega|\leq\Omega\}\subset(-\Omega,\Omega)^{d}\,.$$

  Now Shannon’s sampling theorem [39, 40] states that $g$ is sufficiently fine sampled if the step size in $x$ and in $t$ satisfies the Nyquist condition $\Delta_{\rm samp}\leq\pi/\Omega$.

• •

  Sampling of $f^{\dagger}$. Similar considerations as above again show that $f^{\dagger}$ is essentially bandlimited with essential bandwidth $\Omega$. Shannon’s sampling theorem implies that $f^{\dagger}$ can be reliable reconstructed from discrete samples taken with step size $\Delta_{\rm samp}\leq\pi/\Omega$.

If $f$ has essential bandwidth larger than $\Omega$, the function $g$ has to be filtered with a low pass-filter before sampling. Otherwise, sampling introduces aliasing artifacts in the reconstructed image.

Theoretically, the resolution (at least of the visible parts) can be increased ad infinity by simply decreasing the sampling size $\Delta_{\rm samp}$. In practical applications, several other factors such as the bandwidth of the ultrasound detection system limit the bandwidth of the data, and therefore the resolution of reconstructed images [33]. This, however, also guarantees that in practice a moderate sampling step size $\Delta_{\rm samp}$ gives correct sampling without aliasing.