Non-Equispaced Grid Sampling in Photoacoustics with a Non-Uniform FFT

To obtain the initial pressure from the collected data on a planar sensor arrangement in photoacoustic tomography, there exists an exact analytic frequency domain reconstruction formula. An efficient realization of this formula needs to cope with the evaluation of the data's Fourier transform on a non-equispaced mesh. In this paper, we use the non-uniform fast Fourier transform to handle this issue and show its feasibility in 3D experiments with real and synthetic data. This is done in comparison to the standard approach that uses linear, polynomial or nearest neighbor interpolation. Moreover, we investigate the effect and the utility of flexible sensor location to make optimal use of a limited number of sensor points. The computational realization is accomplished by the use of a multi-dimensional non-uniform fast Fourier algorithm, where non-uniform data sampling is performed both in frequency and spatial domain. Examples with synthetic and real data show that both approaches improve image quality.


Introduction
Photoacoustic tomography is an emerging imaging technique that combines the contrast of optical absorption with the resolution of ultrasound images (see for instance 27 ).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 3D 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 2,27,[29][30][31] ).A nonequispaced arrangement of transducers aligned on a spherical array has already been used by. 28ere 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 towards the center.This is done in order to maximize image quality, when the number of sensor points is the limiting factor.This approach can be used for dealing with the limited view problem where deficiencies are caused by a small detection region and is motivated by the capabilities and requirements of our experimental setup.
For the planar arrangement of point-like detectors there exist several approaches for reconstruction, including numerical algorithms based on filtered back-projection formulas and time-reversal algorithms (see for instance 17,29,32,33 ).
The suggested algorithm in the present work realizes a Fourier inversion formula (see (1) below) using the non-uniform fast Fourier transform (NUFFT).This method has been designed for evaluation of Fourier transforms at non-equispaced points in frequency domain, or non-equispaced data points in spatial, respectively temporal domain.The prior is called NER-NUFFT (nonequispaced range non-uniform FFT), whereas the latter is called NED-NUFFT (non-equispaced data non-uniform FFT).Both algorithms have been introduced in. 9Both NUFFT methods have proven to achieve high accuracy and simultaneously reach the computational efficiency of conventional FFT computations on regular grids. 9or 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 (1) contains evaluations at non-equidistant sample points in frequency domain.
2. In addition, and this comes from the motivation of this paper, we consider evaluation at non-uniform sampling points.
The first issue can be solved by a NER-NUFFT implementation: For 2D photoacoustic inversion with uniformly placed sensors on a measurement line, such an implementation has been considered in. 12Furthermore, this method was used for biological photoacoustic imaging in. 23In both papers the imaging was realized in 2D due to the use of integrating line detectors. 5,21 n this paper we will analyze the NER-NUFFT in a 3D imaging setup with point sensors for the first time.The second issue is solved by employing the NED-NUFFT. 9Thus the name NEDNER-NUFFT for the combined reconstruction algorithm.
The outline of this work is as follows: In section 2 we outline the basics of the Fourier reconstruction approach by presenting the underlying photoacoustic model.We state the Fourier domain reconstruction formula (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 with non-uniform sampling points.In section 3 we briefly explain the idea behind the NUFFT.We state the NER-NUFFT (subsection 3.1) and NED-NUFFT (subsection 3.2) formulas in the form we need it to realize the reconstruction on a non-equispaced grid.In section 4 we introduce the 3D experimental setup.
The sections thereafter describe the realized experiments.In section 5 we compare the NER-NUFFT with conventional FFT reconstruction for synthetic data in 3D.For the real data comparisons we add a time reversal reconstruction.Section 6 explains how we choose and implement the non-equispaced sensor placement.In section 7 we turn to the NEDNER-NUFFT in 2D with simulated data, in order to test different sensor arrangements in an easily controllable environment.In section 8 we interpolate an irregular equi-steradian sensor arrangement data from experimentally acquired data-sets.We apply our NEDNER-NUFFT approach to the non-uniform data and quantitatively compare the reconstructions to regular grid reconstructions.We conclude with a summary of the results in section 9, where we also discuss the benefits and limitations of the presented methods.

Numerical Realization of a Photoacoustic Inversion Formula
Let U ⊂ R d be an open domain in R d , and Γ 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 Q[f ] is the trace on Γ × (0, ∞) of the solution of the equation In other words, the photoacoustic imaging problem consists in identifying the initial source f from measurement data g = p| Γ×(0,∞) .
An explicit inversion formula for Q in terms of the Fourier transforms of f and g := Q[f ] has been first formulated by 20 and introduced to photoacoustics by. 16Let (x, y) ∈ R d−1 ×R + .Assume without loss of generality (by choice of proper basis) that Γ is the hyperplane described by y = 0. Then the reconstruction reads as follows: where F denotes the d-dimensional Fourier transform: Here, the variables x, K x are in R d−1 , whereas y, K y ∈ R.
For the numerical realization these three steps have to be realized in discrete form: We denote evaluations of a function ϕ at sampling points For convenience, we will modify this notation in case of evaluations on an equispaced Cartesian grid.We define the d-dimensional grid and assume our sampling points to be located on m∆ x , n∆ y , where and write ϕ m,n = ϕ(m∆ x , n∆ y ) , where ∆ x := X/N x resp.∆ y := Y /N y are the occurring step sizes.
In frequency domain, we have to sample symmetrically with respect to K y .Therefore, we also introduce the interval G Ky := {−N y /2, . . ., N y /2 − 1}.
Since we will have to deal with evaluations that are partially in-grid, partially not necessarily ingrid, we will also use combinations of ( 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 f resp.ĝ.Let in the following denote the d-dimensional discrete Fourier transform with respect to space and time.By discretizing formula (1) via Riemann sums it follows where This is the formula from. 12emark 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 x m ∈ (−X/2, X/2) d−1 : Then, fj,l ≈ 2l κ j,l n∈Gy e −2πiκ j,l n/Ny The term h m represents the area of the detector surface around x m and has to fulfill Note that the original formula (4) can be received from ( 5) by choosing {x m } to contain all points on the grid ∆ x G x .Formula (5) can be interpreted as follows: Once we have computed the Fourier transform of the data and evaluated the Fourier transform at non-equidistant 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 right hand side of (5) would lead to a computational complexity of order N 2 y × 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 (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. 25nfortunately, this kind of interpolation seems to be sub-optimal for Fourier-interpolation with respect to both accuracy and computational costs 9,34 A regularized inverse k-space interpolation has already been shown to yield better reconstruction results. 15The superiority of applying the NUFFT, compared to linear interpolation, has been shown theoretically and computationally by. 12 3 The non-uniform fast Fourier transform (NUFFT) 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 (subsection 3.1) and the NED-NUFFT (subsection 3.2) in the form (and spatial dimensions) we utilize them afterwards.
The NEDNER-NUFFT algorithm used for implementing (5) essentially (up to scaling factors) consists of the following steps: 1. Compute a d − 1 dimensional NED-NUFFT in the 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 (5).
3. Compute an equispaced d-dim inverse FFT to obtain a d dimensional picture of the initial pressure distribution.

The non-equispaced range (NER-NUFFT) case
With the NER-NUFFT (non-equispaced range -non-uniform FFT) it is possible to efficiently evaluate the discrete Fourier transform at non-equispaced positions in frequency domain.
To this end, we introduce the one dimensional discrete Fourier transform, evaluated at nonequispaced grid points κ l ∈ R: In order to find an efficient algorithm for evaluation of (6), we use a window function Ψ, an oversampling factor c > 1 and a parameter c < α < π(2c − 1) that satisfy: 1. Ψ is continuous inside some finite interval [−α, α] and has its support in this interval and Then (see 9,12 ) we have the following representation for the Fourier modes occurring in ( 6): By assumption, both Ψ and Ψ are concentrated around 0. So we approximate the sum over all k ∈ Z by the sum over the 2K integers k that are closest to κ l + k.By choosing θ = 2πn/N − π and inserting (7) in (6), we obtain Here K denotes the interpolation length and where µ l,k is the nearest integer (i.e. the nearest equispaced grid point) to κ l + k.
The choice of Ψ is made in accordance with the assumptions above, so we need Ψ to have compact support.Furthermore, to make the approximation in (8) reasonable, its Fourier transform Ψ needs to be concentrated as much as possible in [−K, K].In practice, a common choice for Ψ is the Kaiser-Bessel function, which fulfills the needed conditions, and its Fourier transform is analytically computable.

The non-equispaced data (NED-NUFFT) case
A second major aim of the present work is to handle data measured at non-equispaced acquisition points x m in an efficient and accurate way.Therefore we introduce the non-equispaced data, The theory for the NED-NUFFT is largely analogous to the NER-NUFFT 9 as described in Subsection 3.1.The representation ( 7) is here used for each entry of j and inserted (with now setting θ = 2πn/N ) into formula (10), which leads to where the entries in µ m,k are the nearest integers to x m + k.Here we have used the abbreviations for the needed evaluations of Ψ and Ψ.
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. 9,12 he 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. 36It consists of a Fabry Pérot (FP) polymer film sensor for interrogation, 3,4 a 50 Hz pulsed laser source and a subsequent optical parametric oscillator (OPO) 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 approximately 38 µ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%. 36The 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 4 and experimentally confirmed. 36There is a linear roll-off reaching zero at 57.9 MHz, with a subsequent rise.

Comparison of the NER-NUFFT reconstruction with FFT 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 @ 3.6 GHz).All parameters that are not unique to the reconstruction method are left equal.

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 • 200 • 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 figure 1.
To obtain the closest possible numerical reconstruction of the Fourier based inversion formula, we directly evaluate the right hand side of formula (4), and subsequently invert a 3D equispaced Fourier transform by applying the conventional 3D (inverse) FFT.We call this reconstruction direct FT.It serves as ground truth for computing the correlation coefficient (appendix .2).
For the NUFFT reconstruction, the temporal frequency oversampling factor in ( 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 doesn't 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   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 3D data and the 3 maximum intensity projections.
time (s) embryo stage 13 HH21 HH27 NER-NUFFT (k=2,c=2)  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.

Experimentally acquired data
For an overall qualitative assessment two data sets, of a 3.5, and a 5 day old chick embryo, are used. 18This corresponds to the development stages HH21 and HH27 of the Hamburger & Hamilton (HH) criterion. 13The data are sampled with a spatial step size of 60 µm, covering an area of 1.008 • 1.008mm 2 (3.5 days) and 1.02 • 1.02mm 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 figure 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, Ψ as defined in (9) can be precomputed, which roughly halves reconstruction time in subsequent reconstructions using the same discretization, as has been already reported in. 23Moreover, the computation time improves by a factor of 200 when using FFT-based reconstructions instead of time reversal.x-axis (mm) x-axis (mm) x-axis (mm) x-axis (mm) x-axis (mm) x-axis (mm) All reconstructions are normalized, so their maximum value is 1.On the bottom graphs, the cumulative signal for each z-axis layer is plotted as a fraction of the corresponding NER-NUFFT layer.While the depth dependent signal between the time reversal reconstruction and the NER-NUFFT roughly remains the same, for the FFT with linear interpolation a fall-off can be observed.The linear fit suggests a reduction of 10.6 % per mm for the 3.5 days old chick embryo and 13.7 % per mm for the 5 days old data.
The relative Tenenbaum sharpness (appendix .2) for the 3D data of the 3.5 days 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 figure 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 figure 2 the cumulative reconstructed signal for each layer is plotted, as fraction of the NER-NUFFT cumulative signal.The additional fall-off 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 zaxis is primarily affected by errors introduced by a sub-optimal implementation of equation 1, this problem needs further research to be fully understood.

Non-Equispaced 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.Reducing this acquisition time is a crucial step in advancing photoacoustic tomography towards clinical and preclinical application.Therefore, 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 non-equispaced positioned sensors, as error analyses for the NED-and the NER-NUFFT indicate. 9This newly gained flexibility of sensor positioning offers many possibilities to enhance the image quality compared to a rectangular grid.
Also any non-equispaced grids that may arise from a specific experimental setup can be efficiently computed via the NEDNER-NUFFT approach.

Equi-angular and equi-steradian 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, 35 this is the region which is enclosed by the normal lines from the edges of the sensor.Mathematically speaking, the wave front propagates on straight lines in the direction of the singularity [14, 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. 19Therefore certain edges are invisible to the detector, as depicted in figure 3.One approach to overcome this problem experimentally has been made by enclosing the target in a reverberant cavity. 6In addition, a lot of effort has been made to enhance reconstruction techniques in order to deal with the limited view problem. 1,8,10,24,26,35 Our proach 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 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 equi-angular, or equi-steradian 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 equi-angular 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 figure 4.
The obvious expansion of an equi-angular projection to 3D is the equi-steradian projection.Here we face a problem analogous to the problem of placing equispaced points on a 3D sphere and then projecting the points, from the center of the sphere, onto a 2D plane outside the sphere (the detector plane).
The algorithm used for this projection is explained in detail in appendix .1.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 figures 7 and 8 the sensor arrangements are depicted.

Weighting term
To determine the weighting term h m in (5) for 3D we introduce a function that describes the density of equidistant points per unit area ρ p .In our specific case, ρ p describes the density on a sphere around a center of interest.Further we assume that ρ p is spherically symmetric and decreases quadratically with the distance from the center of interest r: ρ p,s ∝ 1/r 2 .We now define ρ p,m for a plane positioned at distance r 0 from the center of interest.In this case ρ p,s (r) attenuates by a factor of sin α, where α = arcsin(r 0 /r) is the angle of incidence.Hence ρ p,m ∝ r 0 /r 3 .This yields a weighting term of: h m (r) ∝ r 3 Analogously we can derive h m for 2D: 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. 7

Application of the NEDNER-NUFFT with Synthetic Data in 2D
A tree phantom, designed by Brian Hurshman and licensed under CC BY 3.01 , is chosen for the 2 dimensional computational experiments on a grid with x = 1024 z = 256 points.A forward simulation is conducted via k-wave 1.1. 25The forward simulation of the k-wave toolbox is based on a first order k-space model.A PML (perfectly matched layer) of 64 grid points is added.Also white noise is added to obtain an SNR (signal to noise ratio) of 30 dB.
In figure 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 figure 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 .2.
We apply the correlation coefficient only within the region of interest marked by the white circle in figure 4. The Tenenbaum sharpness was calculated on the smallest rectangle, containing all pixels within the circle.The results are shown in figure 5.
The Tenenbaum sharpness for the equi-angular 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 equi-angular arrangement is 42.3 % closer to a full correlation than any equispaced grid.
In figure 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 equi-angular arrangement a trade off 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 figure 4 shows the equi-angular sensor arrangement, where the missing sensor points are polynomially interpolated to an equispaced grid and a NER-NUFFT reconstruction is applied afterwards.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.

Application of the NEDNER-NUFFT with Experimental Data in 3D
We will now examine if the positive effects of the NEDNER-NUFFT reconstruction with nonequispaced detectors are transferred to 3D data.For these comparisons we use the data sets of the two chick embryos already presented in section 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 figures 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 .2).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 60µ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 non rectangular 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 (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.8, the blue ROI in figure 7.For the red ROI the focus point for the equi-steradian arrangement is shown and the detection region, which marks the area where sensor points are located.The distance r of the 'focus point' from the plane governs the size of the detection region and the ROI, according to the proportions shown.Within the detection region all sensor points occupy the same steradian from the focus point's perspective.
We use the two chick embryo data sets to extract the irregular sensor data via layer-wise polynomial interpolation.A clipping of the MIPxz of reconstruction of both chick embryos is shown in figure 2. A full NER-NUFFT reconstruction of the 5 day old chick embryo is shown in figure 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 is the same for all measurements as depicted in 8.
In order to avoid spatial aliasing the time data have been low pass filtered with a cut off frequency according to F cutof f = c sound /2dx where c sound is the sound speed and dx the step size.In the equi-steradian 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 equi-steradian sensor arrangement described in appendix .1 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 figures 7, 8 and 8.The figures are organized in a similar manner and depict different reconstructions via maximum intensity projections (MIPs).On the top left segment, the equi-steradian 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 is 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 equi-steradian grid arrangement has a rather high resolution towards 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µm in order to have a consistent µm/pixel spacing in the reconstruction.
The results show that the equi-steradian 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.

Summary and results
We computationally implemented a 3D non-uniform FFT photoacoustic image reconstruction, called NER-NUFFT (non equispaced range-non uniform FFT) 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 data sets for comparison, recorded with a FP-planar sensor setup 36 and included the k-wave time reversal algorithm.Regarding reconstruction time the results of 23 and 12 could be confirmed for 3D, 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 Ψ and Ψ had already been pre-computed 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 fall-off.While this fall-off 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 (non equispaced data-NER-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 equi-angular sensor placement for 2D and an equi-steradian placement in 3D, which assigns one sensor point to each angle/steradian for a given center of interest.
For the 2D 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 3D 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 equi-steradian arrangement with a limited number of sensor points for three regions of interest.While the volume of these regions ranged from 1.7 to 13.7 mm 3 the shape always remained a cylinder with a height to diameter ratio of 8:9.
For our regions of interest, the correlation of the equi-steradian arrangement to the full reconstruction, was consistently higher than any regular grid arrangement, using an almost equal number of sensor points.

Discussion
For the case of regular sampled grids the results of 12 where confirmed for 3D 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 non-equispaced 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 region of interest, 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 region of interest 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 non-equispaced sensor point placement is possible.This could mitigate the image degradation towards 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.
.1 Appendix: Algorithm for equi-steradian 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 θ = 0 at the z-axis towards the xy-grid.
First we determine the steradian Ω of the spherical cap from the point of interest, that projects onto the acquisition point plane via Ω = 2π (1 − cos (θ max )) .
This leads to a unit steradian ω = Ω/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 ω j k = 2π (cos (θ k−1 ) − cos (θ k )) , where θ 1 encloses exactly one unit steradian ω 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 • 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: (pol, az) = ((θ k + θ k+1 ) /2, ϕ i,k ) .2 Appendix: Quality measures In the case where a ground truth image is available, we choose the correlation coefficient ρ, which is a measure of the linear dependence between two images U 1 and U 2 .Its range is [−1, 1].A correlation coefficient close to 1 indicates linear dependence. 22It is defined via the variance Var (U i ) of each image and the covariance Cov (U 1 , U 2 ) of the two images: 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. 24,35 e 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. 11Out of the plethora of published focus functions we select the Tenenbaum function, because of its robustness to noise: with g as the Sobel operator: 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 F Tenenbaum = F Tenenbaum /N , where N is the number of elements in U .

Fig 1
Fig 1 Maximum intensity projections (MIPs) in the xy and xz plane of different reconstructions, for a solid sphere.The direct Fourier Transform (top) serves as ground truth.c denominates the upsampling factor in the time domain and k the interpolation width of the Kaiser-Bessel function in (4).

Fig 2
Fig 2The top shows clippings of the MIPs along the y-axis for the two chick embryos.The three reconstruction methods from top to bottom are: FFT with linear interpolation, NER-NUFFT and time reversal.All reconstructions are normalized, so their maximum value is 1.On the bottom graphs, the cumulative signal for each z-axis layer is plotted as a fraction of the corresponding NER-NUFFT layer.While the depth dependent signal between the time reversal reconstruction and the NER-NUFFT roughly remains the same, for the FFT with linear interpolation a fall-off can be observed.The linear fit suggests a reduction of 10.6 % per mm for the 3.5 days old chick embryo and 13.7 % per mm for the 5 days old data.

Fig 3
Fig 3 Depiction of the limited view problem.Edges whose normal vectors cannot intersect with the sensor surface are invisible to the sensor.The invisible edges are the coarsely dotted lines.The detection region is marked by a grey background.The finely dotted lines are used to construct the invisible edges.Edges perpendicular to the sensor surface are invisible for a plane sensor.

Fig 4
Fig 4 Various reconstructions of a tree phantom (top) with different sensor arrangements.All sensor arrangements are confined to 32 sensor points.The sensor positions are indicated as white rectangles on the top of the images.The second image shows the best (see figure 5) equispaced sensor arrangement, with a distance of 13 points between each sensor.The third image shows the NEDNER-NUFFT reconstruction with equi-angular arranged sensor positions.The bottom image shows the same sensor arrangement, but all omitted sensor points are polynomially interpolated and afterwards a NER-NUFFT reconstruction was conducted.

Fig 5
Fig 5 Correlation coefficient and Tenenbaum sharpness for equispaced sensor arrangements with intervals between the sensor points reaching from 1 to 32.The maximum of the correlation coefficient is at 13.The corresponding reconstruction is shown in figure 4. The straight lines indicate the results for the equi-angular projection.

Fig 6
Fig 6 MIPs (maximum intensity projections) of a full NER-NUFFT reconstruction of a 5 day old chick embryo (HH27), cropped along the y-axis.Two cylindrical ROIs (regions of interest) are indicated, each by a circle and two rectangles.The red ROI is discussed in figure 8, the blue ROI in figure 7.For the red ROI the focus point for the equi-steradian arrangement is shown and the detection region, which marks the area where sensor points are located.The distance r of the 'focus point' from the plane governs the size of the detection region and the ROI, according to the proportions shown.Within the detection region all sensor points occupy the same steradian from the focus point's perspective.

Fig 9
Fig 9 Reconstruction comparisons of the center region of the 3.5 day old chick embryo (HH21).The depiction is analog to figure 7, without the MIPyz and the sensor placement.

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