## 1.

## Introduction

Photoacoustic tomography is an emerging imaging technique that combines the contrast of optical absorption with the resolution of ultrasound images (see, for instance, Ref. 1). In experiments, an object is irradiated by a short-pulsed laser beam. Depending on the absorption properties of the material, some of the pulse energy is absorbed and converted into heat. This leads to a thermoelastic expansion, which causes a pressure rise, resulting in an ultrasonic wave called photoacoustic signal. The signal is detected by an array of ultrasound transducers outside the object. Using this signal, the initial pressure is reconstructed, offering a three-dimensional (3-D) image proportional to the amount of absorbed energy at each position. This is the imaging parameter of photoacoustics.

Common measurement setups rely on small ultrasound sensors, which are arranged “uniformly” along simple geometries, such as planes, spheres, or cylinders (see, for instance, Refs. 12.3.4.–5). A nonequispaced arrangement of transducers aligned on a spherical array has already been used by Xiang et al.^{6} Here, we investigate photoacoustic reconstructions from ultrasound signals recorded at “not necessarily equispaced” positions on a planar surface. In other words, we use an irregular sensor point arrangement, where sensor points are denser toward the center. This is done in order to maximize image quality, when the number of sensor points is a limiting factor. Our current setup allows to acquire data from 50 sensor points each second. As a consequence, data acquisition of a typical sample requires several minutes. Reducing this acquisition time is a crucial step in advancing photoacoustic tomography toward clinical and preclinical application. This approach can also be seen as new way for dealing with the limited view problem, where deficiencies are caused by a small detection region.

For the planar arrangement of point-like detectors, there exist several approaches for reconstruction, including numerical algorithms based on filtered backprojection formulas and time-reversal algorithms (see, for instance, Refs. 3 and 78.–9).

The suggested algorithm in the present work realizes a Fourier inversion formula [see Eq. (1) below] using a fast Fourier transform (FFT) method suitable for nonequispaced grid points and frequencies.^{10}11.^{–}^{12} Such methods have been designed for evaluation of Fourier transforms at nonequispaced points in frequency domain, or nonequispaced data points in spatial, respectively, temporal domain. In the present realization introduced in Ref. 13, the prior is called nonequispaced range nonuniform FFT (NER-NUFFT), whereas the latter is called nonequispaced data nonuniform FFT (NED-NUFFT). Both NUFFT methods have proven to achieve high accuracy and simultaneously reach the computational efficiency of conventional FFT computations on regular grids.^{13}

For the reconstruction, we propose a novel combination of NED- and NER-NUFFT, which we call NEDNER-NUFFT, based on the following considerations:

1. The discretization of the analytic inversion formula [Eq. (1)] contains evaluations at nonequidistant sample points in frequency domain.

2. In addition, and this comes from the motivation of this paper, we consider evaluation at nonuniform sampling points.

The first issue can be solved by a NER-NUFFT implementation: For two-dimensional (2-D) photoacoustic inversion with “uniformly” placed sensors on a “measurement line,” such an implementation has been considered in Ref. 14. Furthermore, this method was used for biological photoacoustic imaging in Ref. 15. In both papers, the imaging was realized in 2-D due to the use of integrating line detectors.^{16}^{,}^{17} In this paper, we will analyze the NER-NUFFT in a 3-D imaging setup with point sensors for the first time. The second issue is solved by employing the NED-NUFFT;^{13} thus, the name NEDNER-NUFFT for the combined reconstruction algorithm.

The outline of this work is as follows: In Sec. 2, we outline the basics of the Fourier reconstruction approach by presenting the underlying photoacoustic model. We state the Fourier domain reconstruction formula [Eq. (1)] in a continuous setting. Moreover, we figure out two options for its discretization. We point out the necessity of a fast and accurate algorithm for computing the occurring discrete Fourier transforms (DFTs) with nonuniform sampling points. In Sec. 3, we briefly explain the idea behind the NUFFT. We state the NER-NUFFT (Sec. 3.1) and NED-NUFFT (Sec. 3.2) formulas in the form we need it to realize the reconstruction on a nonequispaced grid. In Sec. 4, we introduce the 3-D experimental setup.

The sections thereafter describe the realized experiments. In Sec. 5, we compare the NER-NUFFT with conventional FFT reconstruction for synthetic data in 3-D. For the real data comparisons, we add a time-reversal reconstruction. Section 6 explains how we choose and implement the nonequispaced sensor placement. In Sec. 7, we turn to the NEDNER-NUFFT in 2-D with simulated data, in order to test different sensor arrangements in an easily controllable environment. In Sec. 8, we interpolate an irregular equisteradian sensor arrangement data from experimentally acquired datasets. We apply our NEDNER-NUFFT approach to the nonuniform data and quantitatively compare the reconstructions to regular grid reconstructions. We conclude with a summary of the results in Sec. 9, where we also discuss the benefits and limitations of the presented methods.

## 2.

## Numerical Realization of a Photoacoustic Inversion Formula

Let $U\subset {\mathbb{R}}^{d}$ be an open domain in ${\mathbb{R}}^{d}$, and $\mathrm{\Gamma}$ a $d-1$ dimensional hyperplane not intersecting $U$. Mathematically, photoacoustic imaging consists in solving the operator equation:

where $f$ is a function with compact support in $U$ and $\mathbf{Q}[f]$ is the trace on $\mathrm{\Gamma}\times (0,\infty )$ of the solution of the equation:An explicit inversion formula for $Q$ in terms of the Fourier transforms of $f$ and $g:=\mathbf{Q}[f]$ has been first formulated by Norton and Linzer^{18} and introduced to photoacoustics by Köstli et al.^{19} Let $(\mathit{x},y)\in {\mathbb{R}}^{d-1}\times {\mathbb{R}}^{+}$. Assume without loss of generality (by choice of proper basis) that $\mathrm{\Gamma}$ is the hyperplane described by $y=0$. Then, the reconstruction reads as follows:

## (1)

$$\mathbf{F}[f](\mathit{K})=\frac{2{K}_{y}}{\kappa (\mathit{K})}\mathbf{F}[\mathbf{Q}f][{\mathit{K}}_{\mathit{x}},\kappa (\mathit{K})],$$For the numerical realization, these three steps have to be realized in discrete form. We hereby largely follow the description in Ref. 14. In addition, we give the formula accordingly for the case of measurements taken on nonequispaced grid points.

Let us denote evaluations of a function $\varphi $ at sampling points $({\mathit{x}}_{m},{y}_{n})\in {(-X/2,X/2)}^{d-1}\times (0,Y)$ by

For convenience, we will modify this notation in case of evaluations on an equispaced Cartesian grid. We define the $d$-dimensional grid

In frequency domain, we have to sample symmetrically with respect to ${K}_{y}$. Therefore, we also introduce the interval:

Since we will have to deal with evaluations that are partially in-grid, partially not necessarily in-grid, we will also use combinations of Eqs. (2) and (3). In this paper, we will make use of discretizations of the source function $f$, the data function $g$ and their Fourier transforms $\widehat{f}$ and $\widehat{g}$, respectively.Let, in the following,

## (4)

$${\widehat{f}}_{\mathit{j},l}\approx \frac{2l}{{\kappa}_{\mathit{j},l}}\sum _{n\in {\mathrm{G}}_{y}}{\mathrm{e}}^{-2\pi \mathrm{i}\text{\hspace{0.17em}}{\kappa}_{\mathit{j},l}n/{N}_{y}}\phantom{\rule{0ex}{0ex}}\xb7\sum _{\mathit{m}\in {\mathbf{G}}_{x}}{\mathrm{e}}^{-2\pi \mathrm{i}(\mathit{j}\xb7\mathit{m}+\mathrm{ln})/{N}_{x}^{d-1}}{g}_{\mathit{m},n},$$This is the formula from Ref. 14.

## Remark 1

Note that we use the interval notation for the integer multi-indices for notational convenience. Moreover, we also choose the length of the Fourier transforms to be equal to ${N}_{x}$ in the first $d-1$ dimensions, respectively. This could be generalized without changes in practice.

Now, we assume to sample $g$ at $M$, not necessarily uniform, points ${\mathit{x}}_{m}\in {(-X/2,X/2)}^{d-1}$: then,

## (5)

$${\widehat{f}}_{\mathit{j},l}\approx \frac{2l}{{\kappa}_{\mathit{j},l}}\sum _{n\in {\mathrm{G}}_{y}}{e}^{-2\pi \mathrm{i}{\kappa}_{\mathit{j},l}n/{N}_{y}}\phantom{\rule{0ex}{0ex}}\xb7\sum _{m=1}^{M}\frac{{h}_{m}}{{\mathrm{\Delta}}_{x}^{d-1}}{e}^{-2\pi \mathrm{i}(\mathit{j}\xb7{\mathit{x}}_{m})/M}{g}_{m,n}.$$Equation (5) can be interpreted as follows: once we have computed the Fourier transform of the data and evaluated the Fourier transform at nonequidistant points with respect to the third coordinate, we obtain the (standard, equispaced) Fourier coefficients of $f$. The image can then be obtained by applying standard FFT techniques.

The straightforward evaluation of the sums on the righthand side of Eq. (5) would lead to a computational complexity of order ${N}_{y}^{2}\times {M}^{2}$. Usually, this is improved by the use of FFT methods, which have the drawback that they need both the data and evaluation grid to be equispaced in each coordinate. This means that if we want to compute Eq. (5) efficiently, we have to interpolate both in domain and frequency space. A simple way of doing that is by using polynomial interpolation. It is used for photoacoustic reconstruction purposes, for instance, in the “k-wave” toolbox for MATLAB.^{20} Unfortunately, this kind of interpolation seems to be suboptimal for Fourier interpolation with respect to both accuracy and computational costs.^{13}^{,}^{21}

A regularized inverse k-space interpolation has already been shown to yield better reconstruction results.^{22} The superiority of applying the NUFFT, compared to linear interpolation, has been shown theoretically and computationally by Haltmeier et al.^{14}

## 3.

## Nonuniform Fast Fourier Transform

This section is devoted to the brief explanation of the theory and the applicability of the nonuniform Fourier transform, where we explain both the NER-NUFFT (Sec. 3.1) and the NED-NUFFT (Sec. 3.2) in the form (and spatial dimensions) we utilize them afterward.

The NEDNER-NUFFT algorithm used for implementing [Eq. (5)] essentially (up to scaling factors) consists of the following steps:

1. Compute a $d-1$ dimensional NED-NUFFT in the $\mathit{x}$-coordinates due to our detector placement.

2. Compute a one-dimensional NER-NUFFT in the ${K}_{y}$-coordinate as indicated by the reconstruction formula [Eq. (5)].

3. Compute an equispaced $d$-dim inverse FFT to obtain a $d$-dimensional picture of the initial pressure distribution.

Such an algorithm has been introduced in Ref. 11. In contrast to the present approach that is based on Ref. 13, it uses Gaussian bells as window functions. The presented error estimate could be improved in Ref. 12. For a qualitative comparison of the most common available methods, see Ref. 13.

## 3.1.

### Nonequispaced Range Nonuniform Fast Fourier Transform Case

With the NER-NUFFT, it is possible to efficiently evaluate the DFT at nonequispaced positions in frequency domain.

To this end, we introduce the one-dimensional DFT, evaluated at nonequispaced grid points ${\kappa}_{l}\in \mathbb{R}$:

## (6)

$${\widehat{\varphi}}_{l}=\sum _{n\in {\mathrm{G}}_{y}}{\varphi}_{n}{\mathrm{e}}^{-2\pi \mathrm{i}{\kappa}_{l}n/N},\phantom{\rule[-0.0ex]{2em}{0.0ex}}l=1,\dots ,M.$$In order to find an efficient algorithm for evaluation of Eq. (6), we use a window function $\mathrm{\Psi}$, an oversampling factor $c>1$ and a parameter $c<\alpha <\pi (2c-1)$ that satisfy:

1. $\mathrm{\Psi}$ is continuous inside some finite interval $[-\alpha ,\alpha ]$ and has its support in this interval and

2. $\mathrm{\Psi}$ is positive in the interval $[-\pi ,\pi ]$.

Then (see Refs. 13 and 14), we have the following representation for the Fourier modes occurring in Eq. (6):

## (7)

$${e}^{-\mathrm{i}x\theta}=\frac{c}{\sqrt{2\pi}\mathrm{\Psi}(\theta )}\sum _{k\in \mathbb{Z}}\widehat{\mathrm{\Psi}}(x-k/c){\mathrm{e}}^{-\mathrm{i}k\theta /c},\phantom{\rule[-0.0ex]{2em}{0.0ex}}|\theta |\le \pi .$$## (8)

$${\widehat{\varphi}}_{l}\approx \sum _{k=-K+1}^{K}{\widehat{\mathrm{\Psi}}}_{l,k}\sum _{n\in {\mathrm{G}}_{y}}\frac{{\varphi}_{n}}{{\mathrm{\Psi}}_{n}}{e}^{-2\pi i\mathit{ln}/cN},\phantom{\rule{0ex}{0ex}}l=1,\dots ,M.$$## (9)

$${\mathrm{\Psi}}_{n}:=\mathrm{\Psi}(2\pi n/{N}_{y}-\pi ),{\widehat{\mathrm{\Psi}}}_{l,k}:=\frac{c}{\sqrt{2\pi}}{\mathrm{e}}^{-\mathrm{i}\pi [{\kappa}_{l}-({\mu}_{l,k})]}\widehat{\mathrm{\Psi}}[{\kappa}_{l}-({\mu}_{l,k})],$$The choice of $\mathrm{\Psi}$ is made in accordance with the assumptions above, so we need $\mathrm{\Psi}$ to have compact support. Furthermore, to make the approximation in Eq. (8) reasonable, its Fourier transform $\widehat{\mathrm{\Psi}}$ needs to be concentrated as much as possible in $[-K,K]$. In practice, a common choice for $\mathrm{\Psi}$ is the Kaiser–Bessel function, which fulfills the needed conditions, and its Fourier transform is analytically computable.

## 3.2.

### Nonequispaced Data Nonuniform Fast Fourier Transform Case

A second major aim of the present work is to handle data measured at nonequispaced acquisition points ${\mathit{x}}_{m}$ in an efficient and accurate way. Therefore, we introduce the nonequispaced data, $d-1$ dimensional DFT

## (10)

$${\widehat{\varphi}}_{\mathit{j}}=\sum _{m=1}^{M}{\varphi}_{m}{\mathrm{e}}^{-2\pi \mathrm{i}(\mathit{j}\xb7{\mathit{x}}_{m})/N},\phantom{\rule{0ex}{0ex}}\mathit{j}\in {\mathbf{G}}_{x}.$$^{13}as described in Sec. 3.1. The representation [Eq. (7)] is here used for each entry of $\mathit{j}$ and inserted (with now setting $\theta =2\pi n/N$) into formula [Eq. (10)], which leads to

## (11)

$${\widehat{\varphi}}_{\mathit{j}}\approx \frac{1}{{\mathrm{\Psi}}_{\mathit{j}}}\sum _{m=1}^{M}\sum _{\mathit{k}\in {\{-K,\dots ,K-1\}}^{d-1}}{\varphi}_{m}{\widehat{\mathrm{\Psi}}}_{\mathit{j},\mathit{k}}\phantom{\rule{0ex}{0ex}}\xb7{\mathrm{e}}^{-2\pi \mathrm{i}(\mathit{j}\xb7{\mathit{\mu}}_{m,\mathit{k}})/cM},$$Further remarks on the implementation of the NED- and NER-NUFFT, as well as a summary about the properties of the Kaiser–Bessel function and its Fourier transform can be found in Refs. 13 and 14.

## 4.

## Experimental Setup

Before we turn to the evaluation of the algorithm, we describe the photoacoustic setup. A detailed explanation and characterization of the working principles of our setup can be found in Ref. 23. It consists of a Fabry–Pérot (FP) polymer film sensor for interrogation,^{24}^{,}^{25} a 50-Hz pulsed laser source and a subsequent optical parametric oscillator, which emits optical pulses. These pulses have a very narrow bandwidth and can be tuned within the visible and near-infrared range. The optical pulses propagate through an optical fiber. When the light is emitted, it diverges and impinges upon a sample. Some of this light is absorbed and partially converted into heat. This leads to a pressure rise generating a photoacoustic wave, which is then recorded via the FP-sensor head. The sensor head consists of an $\sim 38\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ thick polymer (Parylene C), which is sandwiched between two dichroic dielectric coatings. These dichroic mirrors have a noteworthy transmission characteristic. Light from 600 to 1200 nm can pass the mirrors largely unattenuated, whereas the reflectivity from 1500 to 1650 nm (sensor interrogation band) is about 95%.^{23} The acoustic pressure of the incident photoacoustic wave produces a change in the optical thickness of the polymer film. A focused continuous wave laser, operating within the interrogation band, can now determine the change of thickness at the interrogation point via FP interferometry. The frequency response of this specific setup of up to 100 MHz has been analytically predicted, based on a model used in Ref. 25 and experimentally confirmed.^{23} There is a linear rolloff reaching zero at 57.9 MHz, with a subsequent rise.

## 5.

## Comparison of the Nonequispaced Range Nonuniform Fast Fourier Transform Reconstruction with Fast Fourier Transform and Time Reversal

In this section, several reconstruction methods will be compared for regular grids. This will be done with synthetic data as well as with experimental data. All CPU-based reconstructions are carried out with a workstation PC (Quad Core at 3.6 GHz). All parameters that are not unique to the reconstruction method are left equal.

## 5.1.

### Synthetic Data

For the comparison of different implementations of the FFT-based reconstruction, we conduct a forward simulation of a solid sphere on a $200\times 200\times 100$ computational grid, using the “k-wave 1.1” MATLAB library. The maximum intensity projections of the $xy$ and the $xz$ plane of all reconstructions are shown in Fig. 1.

To obtain the closest possible numerical reconstruction of the Fourier-based inversion formula, we directly evaluate the righthand side of Eq. (4), and subsequently invert a 3-D equispaced Fourier transform by applying the conventional 3-D (inverse) FFT. We call this reconstruction “direct FT.” It serves as ground truth for computing the correlation coefficient (Appendix B).

For the NUFFT reconstruction, the temporal frequency oversampling factor in Eq. (8) is set to $c=2$ and the interpolation length is set to $K=2$. In the linear interpolation FFT case, we use both $c=1$ and $c=2$. The FFT reconstruction with $c=1$ was conducted via the k-wave toolbox, which does not provide oversampling options out of the box. However, the oversampling still can be achieved in a computationally not optimal way by adding zeros at the end of the data term in the time dimension. After reconstruction, the temporal dimension translates into the $z$ axis. In the $xy$ dimension, no oversampling is performed.

The correlation coefficient and the computational time of the methods can be found in Table 1. The errors indicate the superiority of the NUFFT reconstruction in comparison to linear interpolation, with a comparable computational effort (Table 1). The results also show that in the FFT case, an artificial oversampling in the temporal frequency dimension is highly recommended.

## Table 1

The first column compares computational times. In the last four columns, the difference to a full correlation with the direct FT reconstruction method is given in %, for the three-dimensional data and the three maximum intensity projections.

Time | (Correlation-100) in % | ||||
---|---|---|---|---|---|

(s) | 3-D | xy | xz | yz | |

NER ($c=2$) | 59 | 0.005 | 0.0003 | 0.001 | 0.001 |

FFT ($c=2$) | 56 | 3.457 | 0.54 | 0.45 | 0.45 |

FFT ($c=1$) | 53 | 14.00 | 0.65 | 0.81 | 0.81 |

## 5.2.

### Experimentally Acquired Data

For an overall qualitative assessment, two datasets, of a 3.5- and a 5-day-old chick embryo, are used.^{26} This corresponds to the development stages HH21 and HH27 of the Hamburger–Hamilton (HH) criterion.^{27} The data are sampled with a spatial step size of $60\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$, covering an area of $1.008\times 1.008\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{mm}}^{2}$ (3.5 days) and $1.02\times 1.02\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{mm}}^{2}$ (5 days) and a time step size of 16 ns, corresponding to a maximum frequency of 31.25 MHz. To avoid aliasing, the signal is low pass filtered to the maximal spatial frequency of 25.3 MHz. A full reconstruction with the NER-NUFFT for the 5-day-old embryo is shown in Fig. 6.

The time reversal reconstruction is performed via the k-wave toolbox. The spatial upsampling factor in the $x$, $y$ direction is set to 2. For time reversal, this is realized by linearly interpolating the sensor data to a finer grid, whereas for the FFT-based reconstructions, zero padding in the Fourier domain is performed.

The oversampling factor for the FFT reconstructions in the time domain is $c=2$. The number of time steps used for the reconstruction covers more than twice the depth range of the visible objects and is 280 for the HH21 and 320 for the HH27 embryo.

In Table 2, a comparison for the computational time is shown. For the NUFFT case, $\mathrm{\Psi}$ as defined in Eq. (9) can be precomputed, which roughly halves reconstruction time in subsequent reconstructions using the same discretization, as has been already reported in Ref. 15. Moreover, the computation time improves by a factor of 200 when using FFT-based reconstructions instead of time reversal.

## Table 2

Comparison of the computational effort of two chick embryo data sets, with different reconstruction methods.

Time (s) | ||
---|---|---|

Embryo stage27 | HH21 | HH27 |

NER-NUFFT ($k=2$, $c=2$) | 21 | 24 |

NER-NUFFT with precomputed $\mathrm{\Psi}$ | 13 | 14 |

FFT with linear interpolation ($c=2$) | 20 | 23 |

Time reversal | 7236 | 7659 |

The relative Tenenbaum sharpness (Appendix B) for the 3-D data of the 3.5-day-old chick embryo was slightly better for the NER-NUFFT reconstruction (44.2) than for time reversal (43.0) and FFT with linear interpolation (41.1). A comparison of clippings of the maximum intensity projection (MIP) in the $xz$ plane is shown in Fig. 2. The time reversal reconstruction seems smoothed compared to the FFT reconstructions, which is probably a result of the different spatial upsampling modalities.

In the bottom graphs of Fig. 2, the cumulative reconstructed signal for each layer is plotted, as fraction of the NER-NUFFT cumulative signal. The additional falloff for the FFT with linear interpolation has been determined by a line fit. For the 3.5-day-old embryo, it was 10.6% per mm and 13.7% per mm for the 5-day-old embryo. While it intuitively makes sense that the $z$-axis is primarily affected by errors introduced by a suboptimal implementation of Eq. (1), this problem needs further research to be fully understood.

## 6.

## Nonequispaced Sensor Placement

The current setups allow data acquisition at just one single sensor point for each laser pulse excitation. Since our laser is operating at 50 Hz, data recording of a typical sample requires several minutes. In an effort to reduce data acqusition time, we try to maximize the image quality for a given number of acquisition points and a given region of interest.

Our newly implemented NEDNER-NUFFT is ideal for dealing with nonequispaced positioned sensors, as error analyses for the NED- and the NER-NUFFT indicate.^{13} This newly gained flexibility of sensor positioning offers many possibilities to enhance the image quality compared to a rectangular grid.

Also any nonequispaced grids that may arise from a specific experimental setup can be efficiently computed via the NEDNER-NUFFT approach.

## 6.1.

### Equiangular and Equisteradian Projections

In this article, we use the NEDNER-NUFFT to tackle the limited view or limited aperture problem, for the case of a limited number of available detectors, which can be placed discretionary on a planar surface. To understand the limited view problem, it is helpful to define a detection region. According to Ref. 28, this is the region which is enclosed by the normal lines from the edges of the sensor. Mathematically speaking, the wavefront propagates on straight lines in the direction of the singularity [Ref. 29, Chapter VIII]. As a consequence, the reconstruction is locally stable if the straight line through the normal to the object boundary passes through the detector surface.^{30} Therefore, certain edges are invisible to the detector, as depicted in Fig. 3. One approach to overcome this problem experimentally has been made by enclosing the target in a reverberant cavity.^{31} In addition, a lot of effort has been made to enhance reconstruction techniques in order to deal with the limited view problem.^{28}^{,}^{32}33.34.35.^{–}^{36}

Our approach to deal with this problem is different. It takes into account that in many cases, the limiting factor is the number of sensor points and the limited view is a consequence of this constraint. We use an irregular grid arrangement that is dense close to a center of interest and becomes sparser the further away the sampling points are located. We realize this by means of an equiangular, or equisteradian sensor arrangement, where for a given point of interest, each unit angle or steradian gets assigned one sensor point. This arrangement can also be seen as a mock hemispherical detector.

For the equiangular sensor arrangement, a point of interest is chosen. Each line, connecting a sensor point with the point of interest, encloses a fixed angle to its adjacent line. In this sense, we mimic a circular sensor array on a straight line. The position of the sensor points is pictured on top of the third image in Fig. 4.

The obvious expansion of an equiangular projection to 3-D is the equisteradian projection. Here, we face a problem analogous to the problem of placing equispaced points on a 3-D sphere and then projecting the points, from the center of the sphere, onto a 2-D plane outside the sphere (the detector plane).

The algorithm used for this projection is explained in detail in Appendix A. Our input variables are the diameter of the detection region, which we define as the diameter of the disc where the sensor points are located, the distance of the center of interest from the sensor plane $r$ and the desired number of acquisition points. In the top left section of Figs. 7 and 8, the sensor arrangements are depicted.

## 6.2.

### Weighting Term

To determine the weighting term ${h}_{m}$ in [Eq. (5)] for 3-D, we introduce a function that describes the density of equidistant points per unit area ${\rho}_{p}$. In our specific case, ${\rho}_{p}$ describes the density on a sphere around a center of interest. Further we assume that ${\rho}_{p}$ is spherically symmetric and decreases quadratically with the distance from the center of interest $r:{\rho}_{p,s}\propto 1/{r}^{2}$. We now define ${\rho}_{p,m}$ for a plane positioned at distance ${r}_{0}$ from the center of interest. In this case, ${\rho}_{p,s}(r)$ attenuates by a factor of $\mathrm{sin}\text{\hspace{0.17em}}\alpha $, where $\alpha =\mathrm{arcsin}({r}_{0}/r)$ is the angle of incidence. Hence, ${\rho}_{p,m}\propto {r}_{0}/{r}^{3}$. This yields a weighting term of

Analogously, we can derive ${h}_{m}$ for 2-D:For the application of this method to the FP setup, it is noteworthy that there is a frequency dependency on sensitivity which itself depends on the angle of incidence. These characteristics have been extensively discussed in Ref. 37.

## 7.

## Application of the Nonequispaced Data and Nonequispaced Range Nonuniform Fast Fourier Transform Case with Synthetic Data in Two Dimensions

A tree phantom, designed by Brian Hurshman and licensed under CC BY 3.0,^{38} is chosen for the 2-D computational experiments on a grid with $x=1024$ $z=256$ points. A forward simulation is conducted via k-wave 1.1.^{20} The forward simulation of the k-wave toolbox is based on a first order k-space model. A perfectly matched layer of 64 grid points is added. Also, white noise is added to obtain an signal to noise ratio of 30 dB.

In Fig. 4, our computational phantom is shown at the top. For each reconstruction, a subset of 32 out of the 1024 possible sensor positions was chosen. In Fig. 4, their positions are marked at the top of each reconstructed image. For the equispaced sensor arrangements, we let the distance between two adjacent sensor points sweep from 1 to 32, corresponding to a detection region sweep from 32 to 1024. The sensor points are always centered in the $x$-axis.

To compare the different reconstruction methods, we use the correlation coefficient and the Tenenbaum sharpness. These quality measures are explained in Appendix B.

We apply the correlation coefficient only within the region of interest marked by the white circle in Fig. 4. The Tenenbaum sharpness was calculated on the smallest rectangle, containing all pixels within the circle. The results are shown in Fig. 5.

The Tenenbaum sharpness for the equiangular sensor placement is 23001, which is above all values for the equispaced arrangements. The correlation coefficient is 0.913 compared to 0.849, for the best equispaced arrangement. In other words, the equiangular arrangement is 42.3% closer to a full correlation than any equispaced grid.

In Fig. 4, the competing reconstructions are compared. While the crown of the tree is depicted quite well for the equispaced reconstruction, the trunk of the tree is barely visible. This is owed to the limited view of the detection region. As the equispaced interval and the detection region increase, the trunk becomes visible, but at the cost of the crown’s quality. In the equiangular arrangement, a tradeoff between these two effects is achieved. Additionally, the weighting term for the outmost sensors is 17 times the weighting term for the sensor point closest to the middle. This amplifies the occurrence of artifacts, particularly outside the region of interest.

The bottom image in Fig. 4 shows the equiangular sensor arrangement, where the missing sensor points are polynomially interpolated to an equispaced grid and a NER-NUFFT reconstruction is applied afterward. The interpolation is conducted for every time step from our subset to all 1024 sensor points. The correlation coefficient for this outcome was 0.772 while the sharpness measure is 15654. This outcome exemplifies the clear superiority of the NUFFT to conventional FFT reconstruction when dealing with irregular grids.

## 8.

## Application of the Nonequispaced Data and Nonequispaced Range Nonuniform Fast Fourier Transform Case with Experimental Data in Three Dimensions

We will now examine if the positive effects of the NEDNER-NUFFT reconstruction with nonequispaced detectors are transferred to 3-D data. For these comparisons, we use the datasets of the two chick embryos already presented in Sec. 5.2. By the use of polynomial interpolation for each time step, we map this data to discretionarily placed points on the acquisition plane. Thus, sensor data is obtained for regular and irregular grids with arbitrary step sizes. The sensor positions are indicated by the red dots in Figs. 7 and 8.

This procedure allows us to use a full reconstruction as a ground truth and thus ensures a quantitative quality control via the correlation coefficient (Appendix B). We are also safe from any experimental errors that could be introduced between measurements. A drawback is that we can only interpolate to step sizes $\geqq 60\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ without loss of information. Therefore, the presented images are always made with rather few sensor points and naturally of a lower quality. However, we want to emphasize that this is a result of our experimental procedure.

For all comparison reconstructions, the NEDNER-NUFFT has been used for practical reasons. While the NER-NUFFT cannot deal with nonrectangular grids, it is equivalent to the NEDNER-NUFFT for rectangular regular grids. Using the NEDNER-NUFFT, the spacing of the computational grid can be chosen freely. It corresponds to the width of the Kaiser–Bessel function for interpolation [see Eq. (11)]. If the computational grid is much finer than the local sensor point density, a strong signal close to the sensor surface will produce high intensity spots with an intensity distribution according to the Kaiser–Bessel function, instead of a homogeneous area. Making the computational grid coarser than the sensor point density produces a more blurry reconstruction with a reduced lateral resolution. The computational grid therefore is chosen as fine as possible without reducing the lateral resolution.

We use the two chick embryo datasets to extract the irregular sensor data via layerwise polynomial interpolation. A clipping of the $\mathrm{MIP}xz$ of reconstruction of both chick embryos is shown in Fig. 2. A full NER-NUFFT reconstruction of the 5-day-old chick embryo is shown in Fig. 6. For the comparisons, we define a region of interest (ROI) in the form of a cylinder with a height to diameter ratio of 8:9. The area where the sensor points are located for a given reconstruction will be called the “detection region.” The proportions between the ROI and the detection region are the same for all measurements as depicted in Fig. 8.

In order to avoid spatial aliasing, the time data have been lowpass filtered with a cutoff frequency according to ${F}_{\text{cutoff}}={c}_{\text{sound}}/2dx$, where ${c}_{\text{sound}}$ is the sound speed and $dx$ is the step size. In the equisteradian grid, the (locally varying) stepsize $dx$ has been defined as the distance to the nearest neighboring point.

We now conduct a fair comparison between the equisteradian sensor arrangement described in Appendix A and regular grid arrangements for the given ROI. This is done by maximizing image fidelity, while always using (approximately) the same number of sensor points. The comparisons are undertaken for three different ROIs. This is done to show that the advantages of the equisteradian arrangement are not confined to a single case, but rather consistent for different features and volume sizes. All selected ROIs need to have a detection region, that is, fully covered by the underlying data set.

The results are shown in Figs. 2, 7, and 8. The figures are organized in a similar manner and depict different reconstructions via MIPs. On the top left segment, the equisteradian grid arrangement is shown. The top right segment shows the ground truth: a NEDNER-NUFFT reconstruction using all original sensor points that are placed within the detection region. In the segments below, the regular grid reconstructions are shown. On the left, the detection region coincides with the region of interest, on the right, it coincides with the detection region. Figures 7 and 8 also show the sensor point placement. While in the small detection region, the reconstructions have a good resolution, edges are blurred and certain features are invisible. As the detection region increases, these features appear, at the cost of reduced overall resolution. The equisteradian grid arrangement has a rather high resolution toward the center, while still displaying the mentioned features.

All the above figures contain four graphs, which depict the correlation coefficient for the three MIPs and the volume data. The detection region for the regular grids is increased and shown on the $x$-axis, while the number of acquisition points stays constant. The computational grid was chosen to be as small as possible but greater than the step size and has been increased in steps of $30\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ in order to have a consistent $\mu \mathrm{m}/\text{pixel}$ spacing in the reconstruction.

The results show that the equisteradian arrangement consistently produces reconstructions that outperform every regular grid arrangement. It provides a good combination of a large detection region and high resolution at the region of interest. This is demonstrated by the juxtaposition of different reconstructions and also confirmed via the correlation coefficient.

## 9.

## Conclusions

## 9.1.

### Summary and Results

We computationally implemented a 3-D nonuniform FFT photoacoustic image reconstruction, called NER-NUFFT to efficiently deal with the nonequispaced Fourier transform evaluations arising in the reconstruction formula.

In the computational results, it could be shown that the NER-NUFFT is much closer (more than 100 times in the test piece) to perfect correlation than the FFT reconstruction with linear interpolation.

We then used real datasets for comparison, recorded with a FP-planar sensor setup^{23} and included the k-wave time reversal algorithm. Regarding reconstruction time, the results of Haltmeier et al.^{14} and Schulze et al.^{15} could be confirmed for 3-D, where the FFT with linear interpolation performs similar to NER-NUFFT. Additionally, the NER-NUFFT reconstruction time could be significantly reduced (almost halved), if the values of the interpolation functions $\mathrm{\Psi}$ and $\widehat{\mathrm{\Psi}}$ had already been precomputed for the chosen discretization. The time reversal computation took more than 300 times longer on a CPU, than any FFT-based reconstruction. Concerning image quality, the NER-NUFFT and time reversal reconstruction perform on a very similar level, while the conventional FFT method fails to correctly image the depth-dependent intensity falloff. While this falloff is almost synchronous for time reversal and NER-NUFFT, there was an additional intensity drop of about 10% per mm in the linear interpolation FFT-based reconstruction.

The second application of the NUFFT approach concerned the applicability of irregular grid arrangements, which were new in photoacoustic tomography. In fact, this was done by implementing the NEDNER-NUFFT. Our goal was to maximize image quality in a given region of interest, using a limited number of sensor points. To do this, we developed an equiangular sensor placement for 2-D and an equisteradian placement in 3-D, which assigns one sensor point to each angle/steradian for a given center of interest.

For the 2-D simulations, we showed that this arrangement enhances the image quality for a given region of interest and a confined number of sensor points in comparison to regular grids.

In 3-D, we used the aforementioned chick embryo data and reconstructed with an interpolated subset of the original sensor data. We thus conducted a fair comparison between regular grid arrangements and the equisteradian arrangement with a limited number of sensor points for three ROIs. While the volume of these regions ranged from 1.7 to $13.7\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{mm}}^{3}$, the shape always remained a cylinder with a height-to-diameter ratio of $8:9$.

For our ROIs, the correlation of the equisteradian arrangement to the full reconstruction was consistently higher than any regular grid arrangement, using an almost equal number of sensor points.

## 9.2.

### Discussion

For the case of regular sampled grids, the results of Haltmeier et al.^{14} were confirmed for 3-D in the synthetic data experiments. The synthetic data results further show the importance of using a zero-pad factor of at least 2 in the time domain, when using FFT-based reconstruction methods. In the case of real data, the main identifiable difference was the additional intensity drop for greater depth of the FFT reconstruction in comparison to the other two methods. The great computational advantage of using FFT-based reconstructions makes it the most suitable method for most cases in a planar sensor geometry setup. There was no detectable difference in the reconstruction quality between the NER-NUFFT and time reversal. From our point of view, the use of the NER-NUFFT therefore seems to be especially useful in the case of high resolution imaging in relatively deep-lying regions.

The NEDNER-NUFFT implementation allowed to efficiently reconstruct data from nonequispaced sensor points. This is used to extend the primary application of a planar sensor surface, recording images over a large area, by the possibility to image a well defined region of interest with a shorter acquisition time. Thus, we designed a sensor mask to better image small regions in larger depths at a fixed number of sensor points, using a design that projects an equispaced hemispherical detector geometry onto a planar sensor surface. For this case, our sensor arrangement produced consistently better reconstructions than any regular grid, because it allowed to maintain a high resolution within our ROI, while still capturing features that could only be detected outside the region of interest. In our test examples, the NEDNER-NUFFT further enhances the image quality in deep regions while maintaining a reasonable computational effort.

In comparison to a real hemispherical detector, there is an increase of acoustic attenuation. This is countered by greater accessibility, scalability, and flexibility on the planar detector. The ROI not necessarily needs to fit into a spherical shape, the size of the region of interest just has an upper bound and the number of acquisition points is limited by the measurement time.

As an outlook, we mention that the case where the field of view is much larger than the imaging depth has not been investigated in this paper. For this case, a similar approach of expanding the field of view by nonequispaced sensor point placement is possible. This could mitigate the image degradation toward the boundaries of the detection region. However, the achievable benefit of such a method would decrease with an increase of the ratio of the detection region area to the detection region boundary and the maximum imaging depth.

## Appendices

## Appendix A:

### Algorithm for Equisteradian Sensor Arrangement

In our algorithm, the diameter of the detection region and the distance of the center of interest from the sensor plane is defined. The number of sensor points $N$ will be rounded to the next convenient value.

Our point of interest is placed at $z={r}_{0}$, centered at a square $xy$ grid. The point of interest is the center of a spherical coordinate system, with the polar angle $\theta =0$ at the $z$-axis towards the $xy$-grid.

First, we determine the steradian $\mathrm{\Omega}$ of the spherical cap from the point of interest that projects onto the acquisition point plane via

This leads to a unit steradian $\omega =\mathrm{\Omega}/N$ with $N$ being the number of sensors one would like to record the signal with. The sphere cap is then subdivided into slices $k$ which satisfy the condition: where ${\theta}_{1}$ encloses exactly one unit steradian $\omega $ and ${j}_{k}$ has to be a power of two, in order to guarantee some symmetry. The value of ${j}_{k}$ doubles, when ${r}_{s}>1.8\xb7{r}_{k}$, where ${r}_{s}$ is the chord length between two points on $k$ and ${r}_{k}$ is the distance to the closest point on $k-1$. These values are chosen in order to approximate local equidistance between acquisition points on the sensor surface.The azimuthal angles for a slice $k$ are calculated according to

with $i=0,\dots ,j-1$, where stems from the former slice $k-1$. The sensor points are now placed on the $xy$-plane at the position indicated by the spherical angular coordinates:## Appendix B:

### Quality Measures

In the case where a ground truth image is available, we choose the correlation coefficient $\rho $, which is a measure of the linear dependence between two images ${U}_{1}$ and ${U}_{2}$. Its range is $[-\mathrm{1,1}]$. A correlation coefficient close to 1 indicates linear dependence.^{39} It is defined via the variance $\mathrm{Var}({U}_{i})$ of each image and the covariance $\mathrm{Cov}({U}_{1},{U}_{2})$ of the two images:

## (12)

$$\rho ({U}_{1},{U}_{2})=\frac{\mathrm{Cov}({U}_{1},{U}_{2})}{\sqrt{\mathrm{Var}({U}_{1})\mathrm{Var}({U}_{2})}}.$$We decided not to use the widely applied ${L}^{p}$ distance measure because it is a morphological distance, meaning it defines the distance between two images by the distance between their level sets. Therefore, two identical linearly dependent images can have a correlation coefficient of 1 and still a huge ${L}^{p}$ distance. This can be dealt with by normalizing the data as in Refs. 28 and 35. We choose the correlation coefficient instead, because in experimentally acquired data, single high intensity artifacts can occur, which would have a disproportionately large effect on the normalized ${L}^{p}$ distance.

In case of experimentally collected data, there are only a few methods available for the comparison of different reconstruction methods. A possible way for measuring sharpness is obtained from a measure for the high frequency content of the image.^{40} Out of the plethora of published focus functions, we select the Tenenbaum function, because of its robustness to noise:

Like the ${L}^{2}$ norm and unlike the correlation coefficient, the Tenenbaum function is an extensive measure, meaning it increases with image dimensions. Therefore, we normalized it to ${\overline{F}}_{\text{Tenenbaum}}={F}_{\text{Tenenbaum}}/N$, where $N$ is the number of elements in $U$.

## Acknowledgments

This work is supported by the Medical University of Vienna, the European projects FAMOS (FP7 ICT 317744) and FUN OCT (FP7 HEALTH 201880), Macular Vision Research Foundation (MVRF, USA), Austrian Science Fund (FWF), Project P26687-N25 (Interdisciplinary Coupled Physics Imaging), and the Christian Doppler Society (Christian Doppler Laboratory “Laser development and their application in medicine”). We further want to thank Barbara Maurer and Wolfgang J. Weninger from the Center for Anatomy and Cell Biology at the Medical University of Vienna for providing us with the chick embryos.