Computational Realization of a Non-Equidistant Grid Sampling in Photoacoustics with a Non-Uniform FFT

Computational Realization of a Non-Equidistant Grid Sampling in Photoacoustics with a Non-Uniform FFT

Abstract

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. This is done in comparison to the standard approach that uses polynomial interpolation. Moreover, we investigate the effect and the utility of flexible sensor location on the quality of photoacoustic image reconstruction. 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. We show that with appropriate sampling the imaging quality can be significantly improved. Reconstructions with synthetic and real data show the superiority of this method. Keywords: Image reconstruction, Photoacoustics, non-uniform FFT

1Introduction

Photoacoustic tomography is an emerging imaging technique that combines the good contrast of optical absorption with the resolution of ultrasound images (see for instance [19]). In experiments an object is irradiated by a short-pulsed laser beam. Depending on the absorption properties of the material, some light 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 pressure distribution at the time of the laser excitation 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 covering the specimen of interest (see for instance [21] below).

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 [20]).

The suggested algorithm in the present work realizes a Fourier inversion formula (see 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 (non-equispaced range-non uniform FFT), whereas the latter is called NED-NUFFT (non-equispaced data-non uniform FFT). Both algorithms have been introduced in [7]. Both NUFFT methods haven proven to achieve high accuracy and simultaneously reach the computational efficiency of conventional FFT computations on regular grids [7]. This work investigates photoacoustic reconstructions from ultrasound signals recorded at non-equispaced positions on a planar surface. To the best of our knowledge, this is a novel research question in Photoacoustics, where the use of regular grids is the common choice. 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 formulas (see ) contains evaluation 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 [9]. This method was used for biological photoacoustic imaging in [15]. In both papers the imaging could be realized in 2D because integrating line detectors [4] have been used for data recording. In this paper, however, the focus is on 3D imaging, because measurements are taken with point sensors. Experimentally, we show the applicability and superiority of the (NED)NER-NUFFT reconstruction formula in three spatial dimensions, compared with standard interpolation FFT reconstruction. To be precise, in this paper we conduct three dimensional imaging implemented using NEDNER-NUFFT, with ultrasound detectors aligned non-uniformly on a measurement plane. To easily assess the effect of a given arrangement, also 2D numerical simulations have been conducted, to support the argumentation. We quantitatively compare the results with other computational imaging methods: As a reference we use the -wave toolbox [17] with a standard FFT implementation of the inversion algorithm. The NEDNER-NUFFT yields an improvement of the lateral and axial resolution (the latter even by a factor two).

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 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 discuss the 3D experimental setup. The NER-NUFFT is compared with conventional FFT reconstruction. A test chart is used to quantify resolution improvements in comparison to the k-wave FFT reconstruction with linear interpolation. In axial direction this improvement was about 170% while reducing the reconstruction time by roughly 35%.

In Section 5 we then turn to the NEDNER-NUFFT in 2D with simulated data, in order to test different sensor arrangements in an easily controllable environment. An equiangular arrangement turns out to yield an improvement of over 40% compared to the best choice of equispaced sensor arrangement. Furthermore, we use the insights gained from the 2D simulations to develop an equi-steradian sensor arrangement for our 3D measurements. We apply our NEDNER-NUFFT approach on these data and quantitatively compare the outcomes with reconstructions from equispaced data obtained by the NER-NUFFT approach. Our results show a significant improvement of the already superior NER-NUFFT.

2Numerical Realization of a Photoacoustic Inversion Formula

Let be an open domain in , and a dimensional hyperplane not intersecting . Mathematically, photoacoustic imaging consists in solving the operator equation

where is a function with compact support in and is the trace on of the solution of the equation

In other words, the photoacoustic imaging problem consists in identifying the initial source from measurement data .

An explicit inversion formula for in terms of the Fourier transforms of and has been found in [25]. Let . Assume without loss of generality (by choice of proper basis) that is the hyperplane described by . Then the reconstruction reads as follows:

where denotes the -dimensional Fourier transform:

and

Here, the variables are in , whereas .

For the numerical realization these three steps have to be realized in discrete form: We denote evaluations of a function at sampling points by

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

and assume our sampling points to be located on , where

and write

where resp. are the occurring step sizes.

In frequency domain, we have to sample symmetrically with respect to . 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 and . In this paper, we will make use of discretizations of the source function , the data function and their Fourier transforms resp. .

Let in the following

denote the -dimensional discrete Fourier transform with respect to space and time. By discretizing formula via Riemann sums it follows

where

This is the formula from [9].

Now, we assume to sample at , not necessarily uniform, points : Then,

The term represents the area of the detector surface around and has to fulfil . Note that the original formula can be received from by choosing to contain all points on the grid .

Formula 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 . The image can then be obtained by applying standard FFT techniques.

The straightforward evaluation of the sums on the right hand side of would lead to a computational complexity of order . 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 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 [17]. Unfortunately, this kind of interpolation seems to be sub-optimal for Fourier-interpolation with respect to both accuracy and computational costs [7]

A regularized inverse k-space interpolation has already been shown to yield better reconstruction results [10]. The superiority of applying the NUFFT, compared to linear interpolation, has been shown theoretically and computationally by [9].

3The non-uniform fast Fourier transform

This section is devoted to the brief explanation of the theory and the applicability of the non-uniform 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 utilise them afterwards.

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

  1. Compute a dimensional NED-NUFFT in the -coordinates due to our detector placement.

  2. Compute a one-dimensional NER-NUFFT in the -coordinate as indicated by the reconstruction formula .

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

3.1The 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 non-equispaced grid points :

In order to find an efficient algorithm for evaluation of , we use a window function , an oversampling factor and a parameter that satisfy:

  1. is continuous inside some finite interval and has its support in this interval and

  2. is positive in the interval .

Then (see [7]) we have the following representation for the Fourier modes occurring in (Equation 7):

By assumption, both and are concentrated around . So we approximate the sum over all by the sum over the integers that are closest to . By choosing and inserting in , we obtain

Here denotes the interpolation length and

where is the nearest integer (i.e. the nearest equispaced grid point) to .

The choice of is made in accordance with the assumptions above, so we need to have compact support. Furthermore, to make the approximation in (Equation 9) reasonable, its Fourier transform needs to be concentrated as much as possible in . In practice, a common choice for is the Kaiser-Bessel function, which fulfils the needed conditions, and its Fourier transform is analytically computable.

3.2The non-equispaced data (NED-NUFFT) case

A second major aim of the present work is to handle data measured at non-equispaced acquisition points in an efficient and accurate way. Therefore we introduce the non-equispaced data, dimensional DFT

The theory for the NED-NUFFT is largely analogous to the NER-NUFFT [7] as described in subSection 3.1. The representation is here used for each entry of and inserted (with now setting ) into formula , which leads to

where the entries in are the nearest integers to . 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 summery about the properties of the Kaiser-Bessel function and its Fourier transform can be found in [7].

4Comparison of NER-NUFFT and k-wave FFT

Before we turn to the evaluation of the algorithm we describe the photoacoustic setup. Our device consists of a FP (Fabry Pérot) polymer film sensor for interrogation [2]. A pulsed laser source and a subsequent optical parametric oscillator (OPO) provide optical pulses. These pulses have a very narrow bandwidth and can be tuned within the visible and near infrared spectrum. The optical pulses are then transmitted via an optical fibre. When the light is emitted it diverges and impinges upon a sample with homogeneous fluence, thus generating a photoacoustic signal. This signal is then recorded via the FP-sensor head. The sensor head consists of an approximately thick polymer (Parylene C) which is sandwiched between two dichroic dielectric coatings. These dichroic mirrors have a noteworthy transmission characteristic. Light from to can pass the mirrors largely unabated, whereas the reflectivity from to (sensor interrogation band) is about 95% [27]. The incident photoacoustic wave produces a linear 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.

We choose two 2D targets for comparison, a star and a USAF (US Air Force) resolution test chart. Both targets are made of glass with a vacuum-deposited durable chromium coating. The star target has 72 sectors on a pattern diameter of 5 mm, with an unresolved core diameter of .

The targets are positioned in parallel to the sensor surface at a distance of about , and water is used as coupling medium between the target and the sensor.

The interpolation length for the NER-NUFFT reconstruction is . The computational times are shown in Tab. ? showing that the linearly interpolated FFT is about 30 % slower than the NER-NUFFT.

Segment of the MIP (maximum intensity projection) in the z-axis of a star sample, reconstructed using FFT and linear interpolation (top left) and using the NUFFT (top right) with a square plotted around the center. The intensity of the reconstructed image along the sidelines of the square is plotted, for both reconstructions. The purple line indicates the frequency of the star sample.
Segment of the MIP (maximum intensity projection) in the -axis of a star sample, reconstructed using FFT and linear interpolation (top left) and using the NUFFT (top right) with a square plotted around the center. The intensity of the reconstructed image along the sidelines of the square is plotted, for both reconstructions. The purple line indicates the frequency of the star sample.

A segment of the reconstructed star target is shown in Fig. ?. The intensity is plotted against the sides of an imaginary square () placed around the center of the star phantom. It is clearly visible, that the FFT reconstruction is not able to represent the line pairs, when the density exceeds , corresponding to a resolution of , whereas they are still largely visible for the NER-NUFFT reconstruction.

Comparison between the NED-NUFFT and FFT reconstruction for a USAF-chart and a comparison of computational times for both phantoms. The improvement in percent was caluclated by:
NUFFT FFT Improvement
FWHM axial LSF
FWHM lateral LSF
Time: Star target 140 s 189 s
Time: USAF chart 298 s 384 s

For a quantitative comparison of the resolution we use the USAF chart. It is recorded on an area of sensor points corresponding to with a grid spacing of and a time resolution of . As yet there is no standardized procedure to measure the resolution of a photoacoustic imaging system. We proceed similar to [27], by fitting a line spread function (LSF) and an edge spread function (ESF) to the intensities of our reconstructed data. For the LSF to be meaningful, its source has to approximate a spatial delta function. This is the case in the -axis, since the chrome coating of the USAF target is just about thick.

Segment of the xy-MIP of an USAF chart reconstruction conducted with FFT (left) and NER-NUFFT (right). On the bottom images the data-sets used for the quantitative resolution are marked. The black square depicts the 49 xy-coordinates used for the LSF fit along the z-axis for the axial resolution. The data for a single point is shown in the bottom right inlay. The white lines show the intensity fit for the lateral resolution. The bottom left inlay shows the data for a single white line.
Segment of the -MIP of an USAF chart reconstruction conducted with FFT (left) and NER-NUFFT (right). On the bottom images the data-sets used for the quantitative resolution are marked. The black square depicts the 49 -coordinates used for the LSF fit along the -axis for the axial resolution. The data for a single point is shown in the bottom right inlay. The white lines show the intensity fit for the lateral resolution. The bottom left inlay shows the data for a single white line.

We fit a Cauchy-Lorentz distribution to the -axis of our reconstructed data for 49 adjacent -coordinates. Their positions are marked as a black square within the white square, depicted in the bottom images of Fig. ?. The reconstruction in the -axis for a single point is shown in the bottom right inlay of Fig. ?, for 8 points, covering a distance of in the -direction, around the maximum intensity. A fit of the Lorentz distribution is shown for both reconstruction methods. The FWHM (full width, have maximum) of the Lorentz distribution,

is the parameter . The output is the intensity in dependence from the -axis, and and are fitting parameters.

The average and standard deviation of for both of the 49 datasets are shown in Table ?. The line spread function FWHM of the FFT-reconstruction turns out to be more than twice as big as the one of the NER-NUFFT reconstruction.

For the lateral resolution, there is no target approximating a delta function, so we have to use the ESF instead:

The here is the FWHM of the associated LSF, and , and are fitting parameters. The ESF requires a step function as source, of which our target provides plenty. We choose the long side of 12 bars, marked by white lines in the bottom images of Fig. ?, for this fit. The data for a particular line is shown in the bottom left inlay. We omitted all datasets, where only one point marked the transition from low to high intensity, rendering the edge fit unreliable and resulting in unrealistic improvements of our new method well over 100%. Finally we averaged over 15 edges. The results are shown in Table ?. While the deviation between different edges is quite high, the improvement of our NER-NUFFT reconstruction for the 15 evaluated edges ranged from to .

5Non-equidistant grid sampling

The current setups allow data acquisition at just one single sensor point for each laser pulse excitation. Since our laser is operating at data recording of a typical sample requires several minutes. Reducing this acquisition time is a crucial step in advancing Photoacoustics towards clinical and preclinical application. Therefore, in this work we try to maximize the image quality for a given number of acquisition points. We are able to do this, because our newly implemented NEDNER-NUFFT is ideal for dealing with non-equispaced positioned sensors. This newly gained flexibility of sensor positioning offers many possibilities to enhance the image quality, compared against a rectangular grid. For instance a hexagonal grid was found to yield an efficiency of compared with for an exact reconstruction of a wave-number limited function [13].

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

Here, we will use it to tackle the limited view problem. Many papers deal with the limited view problem, when reconstructing images [6]. 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 largely a consequence of this limitation.

In our approach we therefore use a 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 equi-steradian sensor arrangement, where for a given point of interest each unit angle or steradian gets assigned one sensor point.

Depiction of the limited view problem. Edges whose normal vector 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 always invisible for a plane sensor.
Depiction of the limited view problem. Edges whose normal vector 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 always invisible for a plane sensor.

To understand the limited view problem, it is helpful to define a detection region. According to [26], this is the region which is enclosed by the normal lines from the edges of the sensor. Pressure waves always travel in the direction of the normal vector of the boundary of the expanding object. Therefore certain boundaries are invisible to the detector as depicted in Fig. ?.

5.1Equiangular and equi-steradian projection sensor mask

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 that sense we mimic a circular sensor array on a straight line. The position of the sensor points can be seen on top of the third image in Fig. ?.

The obvious expansion of an equiangular projection to 3D is an equi-steradian projection. This problem is analogous to the problem of placing equispaced points on a 3D sphere and then projecting the points, from the center of the sphere, through the points, onto a 2D plane outside the sphere.

We developed an algorithm for this problem, which is explained in detail in Appendix Section 10. Our input variables are the grid size, the distance of the center of interest from the sensor plane and the desired number of acquisition points, which will be rounded to the next convenient value.

A sensor arrangement with 1625 points on a grid is shown in the top left image in Fig. ?.

5.2Weighting term

To determine the weighting term in Equation 6 for 3D we introduce a function that describes the density of equidistant points per unit area . In our specific case, describes the density on a sphere around a center of interest. Further we assume that is spherically symmetric and decreases quadratically with the distance from the center of interest : . We now define for a plane positioned at distance from the center of interest. In this case attenuates by a factor of , where is the angle of incidence. Hence . By applying the spacing of the regular grid . This yields a weighting term of:

Analogously we can derive for 2D:

We applied a normalization after the reconstruction to all measurements.

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, which has been extensively discussed in [5]. The angle of incidence for our specific setup is . At this angle, the frequency components around get attenuated by more than . Below the frequency response remains quite stable (attenuation below for the measurement angles occurring in our setup.

6Computational assessment of different non-equispaced grid arrangements in two dimensions to tackle the limited view problem

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 Fig.) equispaced sensor arrangement, with a distance of 13 points between each sensor. The third image shows the NEDNER-NUFFT reconstruction with equiangular arranged sensor positions. The bottom image shows the same sensor arrangement, but all omitted sensor points are linearly interpolated and afterwards a NER-NUFFT reconstruction was conducted.
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 Fig.) equispaced sensor arrangement, with a distance of 13 points between each sensor. The third image shows the NEDNER-NUFFT reconstruction with equiangular arranged sensor positions. The bottom image shows the same sensor arrangement, but all omitted sensor points are linearly interpolated and afterwards a NER-NUFFT reconstruction was conducted.

A tree phantom, designed by Brian Hurshman and licensed under CC BY 3.0 1, is chosen for the 2 dimensional computational experiments on a grid with points. A forward simulation is conducted via k-wave [17]. The forward simulation of the k-wave toolbox is based on a first order k-space model. A PML (perfectly matched layer) of 64 gridpoints is added, as well as of noise.

In Fig. ? 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. ? 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. The sensor points where always centered in the -axis.

To compare the different reconstruction methods we used the correlation coefficient and the Tenenbaum sharpness. These quality measures are explained in Appendix Section 11.

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

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 Fig. . The straight lines indicate the results for the equiangular projection.
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 Fig. . The straight lines indicate the results for the equiangular projection.

The Tenenbaum sharpness for the equiangular sensor placement was 23001, which is above all values for the equispaced arrangements. The correlation coefficient was 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. ? 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 due to the limited view problem. When the equispaced interval increases, the tree becomes visible, but at the cost of the crown’s quality. In the equiangular 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 artefacts, particularly outside of our region of interest.

The bottom image in Fig. ? shows the equiangular sensor arrangement, reconstructed in a conventional manner. The missing sensor points are interpolated to an equispaced grid, and a NER-NUFFT reconstruction is applied afterwards. We conducted a linear interpolation from our subset to all 1024 sensor points. The correlation coefficient for this outcome was 0.7348 while the sharpness measure was 21474. This outcome exemplifies the clear superiority of the NUFFT to conventional FFT reconstruction when dealing with non-equispaced grids.

73D application of the NED-NER-NUFFT with real data

MIPs of all 3 planes for the NER-NUFFT reconstruction of a yarn phantom.
MIPs of all 3 planes for the NER-NUFFT reconstruction of a yarn phantom.

For a qualitative assessment of our new sensor arrangement we need a 3D phantom. We choose a yarn which we record on a rectangular grid with sensor points, with a grid spacing of and time sampling of Hence an area of is covered. As coupling medium water is used, in which the yarn is fully immersed.

To determine the utility of non-equispaced grid sampling, we follow a certain routine. First we acquire a densely sampled dataset. Then we use a very small subset of the initially collected sensor data, to test different sensor arrangements. Therefore we can always use the complete reconstruction as our model standard together with the quality measurements explained in Appendix Section 11.

An upsampling factor of 2 was used for all reconstructions, hence the reconstructed image of the MIP for the plane consists of pixels. The complete reconstruction with the NER-NUFFT took 154 seconds. Maximum intensity projections (MIP) of this full reconstruction for all axis are shown in Fig. ?. The MIP in the plane is our model standard, for comparison with the other reconstructions.

The sensor placement for the cirular arrangement is shown on the top left image, comprising 1625 sensor points. On the top right an equispaced sensor arrangement with 41\times41=1689 sensor points is displayed. The intervall between two adjacent sensor points is 5 for this configuration. The bottom images show MIPs of the NED-NER-NUFFT reconstruction only using the sensor points shown above.
The sensor placement for the cirular arrangement is shown on the top left image, comprising 1625 sensor points. On the top right an equispaced sensor arrangement with sensor points is displayed. The intervall between two adjacent sensor points is for this configuration. The bottom images show MIPs of the NED-NER-NUFFT reconstruction only using the sensor points shown above.

For the equi-steradian sensor mask we choose our center of interest right in the center of the -MIP where the little knot can be seen, off the sensor surface. The resulting sensor mask, including the reconstruction is shown in Fig. ?, it consists of sensor points (or 3.18 % of the initial number of sensor points). The weighting term, accounting for sensor sparsity is 9.7 times higher for the outermost sensor point, than for centermost sensor points in the -plane. The reconstruction for this arrangement took 134 seconds.

We compared this arrangement to rectangular grids, which all had sensor points (or 3.29 % of the initial number of sensor points) but varying distances between two adjacent points. The grid with an interval between sensor points of 5 is shown on the top left in Fig. ?.

Correlation coefficient calculated on different number of pixels, for 3 different reconstructed MIPs. Pixels of interest are added according to the intensity of the corresponding pixels of the model standard image (Fig. ). In the inlay the correlation coefficient mask is shown for 20087 points which corresponds to all pixels with at least 4% of the maximum value of the model standard and is the value used in Fig. .
Correlation coefficient calculated on different number of pixels, for 3 different reconstructed MIPs. Pixels of interest are added according to the intensity of the corresponding pixels of the model standard image (Fig. ). In the inlay the correlation coefficient mask is shown for 20087 points which corresponds to all pixels with at least 4% of the maximum value of the model standard and is the value used in Fig. .

In Fig. ? the equi-steradian grid is compared with 2 equispaced grids. To get a more precise measure of the correlation between the reconstructed images and the model standard, we calculated the correlation coefficient only within a region of interest. The region of interest is firstly confined by a centered disc, whose boundary is shown as a dotted circle in the inlay in Fig. ?. The -axis shows the number of pixels of interest used to calculate the correlation coefficient within this disc. These pixels are increased according to the intensity of the corresponding pixels of the model’s standard MIP. Fig. ? demonstrates that the correlation coefficient for the equi-steradian arrangement, always remains better, within the depicted disc, when competing against the two strongest equispaced grids.

Correlation Coefficient calculated on pixels of interest with at least 4% of the maximum value of the model standard, for an increasing diameter of the confining disc. The straight line at 8.4\,\mathrm{mm} corresponds to the straight line in Fig. . The other straight lines indicate the side length of the square of a particular equispaced sensor point arrangement.
Correlation Coefficient calculated on pixels of interest with at least 4% of the maximum value of the model standard, for an increasing diameter of the confining disc. The straight line at corresponds to the straight line in Fig. . The other straight lines indicate the side length of the square of a particular equispaced sensor point arrangement.

In Fig. ? the 4 equispaced arrangements, with intervals reaching from 2 to 5 are compared to the equi-steradian grid. Here the pixels of interest are set to a threshold of at least 4% of the maximum value, while the diameter is increased. The dotted straight lines indicate the side length of the square of a particular equispaced sensor point arrangement. As expected the correlation coefficient for this equispaced sensor arrangement starts to decline around that threshold.

The correlation coefficient for the equi-steradian grid starts to fall behind the equispaced grid of interval 5 towards the end. This is expected, since the equi-steradian grid is only meant to give better results for a region of interest around the center. This is very clearly the case. Between a diameter of 2 to 8 mm the correlation coefficient is on average 25.14 % closer to a value of than its strongest equispaced contender of interval 4, and 30.87 % better than the interval 5 grid. At a diameter of the correlation coefficient for the equi-steradian grid was 0.960 compared to 0.947 for the interval 4 grid and 0.944 for the interval 5 grid. This is 24.8 % and 28.6% closer to full correlation.

Tenenbaum Sharpness calculated on pixels of interest with at least 4% of the maximum value of the model standard. The Tenenbaum sharpness is calculated on the smallest rectangle, that contains all pixels of interest.
Tenenbaum Sharpness calculated on pixels of interest with at least 4% of the maximum value of the model standard. The Tenenbaum sharpness is calculated on the smallest rectangle, that contains all pixels of interest.

Fig. ? shows the normalized Tenenbaum sharpness. The Tenenbaum sharpness, unlike the correlation coefficient, cannot be calculated on non-adjacent grid points, therefore it has been calculated on the smallest rectangle, that contains all pixels of interest. The equi-steradian has the highest sharpness for most values, with the interval 4 grid being very slightly better around a diameter of 3-4 mm. There is a drop of the Tenenbaum sharpness towards the end.

8Discussion and Conclusion

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 non-equispaced Fourier transform evaluations arising in the reconstruction formula. This method was compared with the k-wave implemented FFT reconstruction, which uses a polynomial interpolation. The two reconstruction methods where compared using 2D targets. The lateral resolution showed an improvement of which is in good agreement with the illustrative results for the star target. The axial resolution showed an improvement of . The computation time was about less, for the NER-NUFFT than the linearly interpolated FFT reconstruction. In conclusion the NER-NUFFT reconstruction proved to be unequivocally superior to conventional linear interpolation FFT reconstruction methods.

We further implemented the NED-NER-NUFFT (non equispaced data-NER-NUFFT), which allowed us to efficiently reconstruct from data recorded at non-equispaced placed sensor points. This newly gained flexibility was used to tackle the limited view problem, by placing sensors more sparsely further away from the center of interest. We developed an equiangular 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. In the 2D computational simulations we showed that this arrangement significantly enhances image quality in comparison to regular grids.

In 3D we conducted experiments, where a yarn phantom was recorded. The maximum intensity projection (MIP) of the full reconstruction was compared to MIPs of reconstructions that only used about 3% of the original data. Within our region of interest, the correlation of our image was 0.96, which is 24.8% closer to full correlation than the best equispaced arrangement, reconstructed from slightly more sensor points.

The sensor placement to tackle the limited view problem, combined with the NED-NER-NUFFT gives significantly better results for an object located at the center of a bigger sensor surface. This result was confirmed in the 2D simulation as well as for real data in 3D.

Acknowledgment

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

10Appendix: Algorithm for equi-steradian sensor arrangement

In our algorithm, the grid size and the distance of the center of interest from the sensor plane is defined. The number of sensor points will be rounded to the next convenient value.

Our point of interest is placed at , centered at a square grid. The point of interest is the center of a spherical coordinate system, with the polar angle at the -axis towards the -grid.

First we determine the steradian of the spherical cap from the point of interest, that projects onto the acquistion point plane via

This leads to a unit steradian with being the number of sensors one would like to record the signal with. The sphere cap is then subdivided into slices which satisfy the condition

where encloses exactly one unit steradian and has to be a power of two, in order to guarantee some symmetry. The value of doubles, when , where is the chord length between two points on and is the distance to the closest point on . These values are chosen in order to guarantee We apply some restrictions to approximate equidistance between acquisition points on the sensor surface.

The azimuthal angles for a slice are calulated according to:

with , where

stems from the former slice . The sensor points are now placed on the -plane at the position indicated by the spherical angular coordinates:

11Appendix: Quality measures

The correlation coefficient is a measure of the linear dependence between two images and .

Its range is and a correlation coefficient close to 1 indicates linear dependence [14]. It is defined via the variance, of each image and the covariance, of the two images:

We did not choose the widely used distance measure because it is a morphological distance measure, 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 and still a huge distance. Normalizing the images only mitigates this problem, because single, high intensity artifacts or a varying intensity over the whole image can greatly alter the distance.

While the correlation coefficient is a good measure for the overall similarity between two images it does not include any sharpness measure. Hence blurred edges are punished very little, in comparison to slight variations of homogeneous areas. To address this shortcoming we chose a sharpness measure or focus function as a second quality criterion. Sharpness measures are obtained from some measure of the high frequency content of an image [8]. They have also been used to select the best sound speed in photoacoustic image reconstruction [18]. Out of the plethora of published focus functions we select the Tenenbaum function, because of its robustness to noise:

with as the Sobel operator:

Like the norm and unlike the correlation coefficient the Tenenbaum function is an extensive measure, meaning it increases with image dimensions. Therefore we normalized it to , where is the number of elements in .

Footnotes

  1. http://thenounproject.com/term/tree/16622/

References

  1. Biomedical photoacoustic imaging.
    P. Beard. Interface Focus, 1:602–631, 2011.
  2. Two-dimensional ultrasound receive array using an angle-tuned fabry-perot polymer film sensor for transducer field characterization and transmission ultrasound imaging.
    P.C. Beard. IEEE Trans. Ultrason., Ferroeletr., Freq. Control, 52(6):1002–1012, June 2005.
  3. Transduction mechanisms of the fabry-perot polymer film sensing concept for wideband ultrasound detection.
    P.C. Beard, F. Perennes, and T.N. Mills. IEEE Trans. Ultrason., Ferroeletr., Freq. Control, 46(6):1575–1582, November 1999.
  4. Thermoacoustic tomography with integrating area and line detectors.
    P. Burgholzer, C. Hofer, G. Paltauf, M. Haltmeier, and O. Scherzer. IEEE Trans. Ultrason., Ferroeletr., Freq. Control, 52(9):1577–1583, September 2005.
  5. The frequency-dependent directivity of a planar fabry-perot polymer film ultrasound sensor.
    B.T. Cox and P.C. Beard. IEEE Trans. Ultrason., Ferroeletr., Freq. Control, 54(2):394–404, 2007.
  6. Influence of limited-view scanning on depth imaging of photoacoustic tomography.
    W. Dan, C. Tao, X.-J. Liu, and X.-D. Wang. Chinese Physics B, 21(1):014301, 2012.
  7. Non-equispaced fast Fourier transforms with applications to tomography.
    K. Fourmont. J. Fourier Anal. Appl., 9(5):431–450, 2003.
  8. A comparison of different focus functions for use in autofocus algorithms.
    F. C. A. Groen, I. T. Young, and G. Ligthart. Cytometry, 6(2):81–91, 1985.
  9. A reconstruction algorithm for photoacoustic imaging based on the nonuniform FFT.
    M. Haltmeier, O. Scherzer, and G. Zangerl. IEEE Trans. Med. Imag., 28(11):1727–1735, November 2009.
  10. Fourier reconstruction in optoacoustic imaging using truncated regularized inverse k-space interpolation.
    M. Jaeger, S. Schüpbach, A. Gertsch, M. Kitz, and M. Frenz. Inverse Probl., 23:S51–S63, 2007.
  11. Mathematics of thermoacoustic tomography.
    P. Kuchment and L. Kunyansky. European J. Appl. Math., 19:191–224, 2008.
  12. Photoacoustic tomography with integrating area and line detectors.
    G. Paltauf, R. Nuster, M. Haltmeier, and P. Burgholzer. In L. V. Wang, editor, Photoacoustic Imaging and Spectroscopy, Optical Science and Engineering, pages 251–263. CRC Press, Boca Raton, FL, 2009.
  13. Sampling and reconstruction of wave-number-limited functions in -dimensional Euclidean spaces.
    D. P. Petersen and D. Middleton. Inform.and Control, 5:279–323, 1962.
  14. Handbook of Mathematical Methods in Imaging.
    O. Scherzer, editor. Springer, New York, 2011.
  15. On the use of frequency-domain reconstruction algorithms for photoacoustic imaging.
    R. Schulze, G. Zangerl, M. Holotta, D. Meyer, F. Handle, R. Nuster, G. Paltauf, and O. Scherzer. J. Biomed. Opt., 16(8):086002, 2011.
  16. Reconstruction of high quality photoacoustic tomography with a limited-view scanning.
    C. Tao and X. Liu. Opt. Express, 18(3):2760–2766, Feb 2010.
  17. K-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wace fields.
    B. E. Treeby and B. T. Cox. J. Biomed. Opt., 15:021314, 2010.
  18. Automatic sound speed selection in photoacoustic image reconstruction using an autofocus approach.
    B. E. Treeby, T. K. Varslot, E. Z. Zhang, J. G. Laufer, and P. C. Beard. Biomed. Opt. Express, 16(9):090501–090501–3, 2011.
  19. Photoacoustic Imaging and Spectroscopy.
    L. V. Wang, editor. Optical Science and Engineering. CRC Press, Boca Raton, 2009.
  20. Exact frequency-domain reconstruction for thermoacoustic tomography–I: Planar geometry.
    M. Xu and L. V. Wang. IEEE Trans. Med. Imag., 21:823–828, 2002.
  21. Time-domain reconstruction for thermoacoustic tomography in a spherical geometry.
    M. Xu and L. V. Wang. IEEE Trans. Med. Imag., 21(7):814–822, 2002.
  22. Analytic explanation of spatial resolution related to bandwidth and detector aperture size in thermoacoustic or photoacoustic reconstruction.
    M. Xu and L. V. Wang. Phys. Rev. E, 67(5):0566051–05660515, 2003.
  23. Universal back-projection algorithm for photoacoustic computed tomography.
    M. Xu and L. V. Wang. Phys. Rev. E, 71(1), 2005.
  24. Photoacoustic imaging in biomedicine.
    M. Xu and L. V. Wang. Rev. Sci. Instruments, 77(4), 2006.
  25. Exact frequency-domain reconstrcution for thermoacoustic tomography — I: Planar geometry.
    Y. Xu, D. Feng, and L. V. Wang. IEEE Trans. Med. Imag., 21(7):823–828, 2002.
  26. Reconstructions in limited-view thermoacoustic tomography.
    Y. Xu, L. V. Wang, G. Ambartsoumian, and P. Kuchment. Med. Phys., 31(4):724–733, 2004.
  27. Backward-mode multiwavelength photoacoustic scanner using a planar Fabry-Perot polymer film ultrasound sensor for high-resolution three-dimensional imaging of biological tissues.
    E. Zhang, J. Laufer, and P. Beard. App. Opt., 47:561–577, 2008.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
""
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
   
Add comment
Cancel
Loading ...
22063
This is a comment super asjknd jkasnjk adsnkj
Upvote
Downvote
""
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Submit
Cancel

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test
Test description