# A Framework for Directional and Higher-Order

Reconstruction in Photoacoustic Tomography

###### Abstract

Photoacoustic tomography is a hybrid imaging technique that combines high optical tissue contrast with high ultrasound resolution. Direct reconstruction methods such as filtered backprojection, time reversal and least squares suffer from curved line artefacts and blurring, especially in case of limited angles or strong noise. In recent years, there has been great interest in regularised iterative methods. These methods employ prior knowledge on the image to provide higher quality reconstructions. However, easy comparisons between regularisers and their properties are limited, since many tomography implementations heavily rely on the specific regulariser chosen. To overcome this bottleneck, we present a modular reconstruction framework for photoacoustic tomography. It enables easy comparisons between regularisers with different properties, e.g. nonlinear, higher-order or directional. We solve the underlying minimisation problem with an efficient first-order primal-dual algorithm. Convergence rates are optimised by choosing an operator dependent preconditioning strategy. Our reconstruction methods are tested on challenging 2D synthetic and experimental data sets. They outperform direct reconstruction approaches for strong noise levels and limited angle measurements, offering immediate benefits in terms of acquisition time and quality. This work provides a basic platform for the investigation of future advanced regularisation methods in photoacoustic tomography.

keywords: photoacoustic tomography, variational image reconstruction, total generalised variation, directional regularisation, compressive sampling, convex optimisation.

## 1 Introduction

Photoacoustic tomography (PAT), also known as optoacoustic tomography, is an intrinsically hybrid imaging technique that combines the high spectroscopic contrast of tissue constituents to light, with the high resolution of ultrasound imaging techniques. Tissue is illuminated with nanosecond laser pulses, causing some of the optical energy to be absorbed and converted into heat. This leads to thermoelastic expansion and creates ultrasound waves that are detected at the boundary. These detected signals are employed to reconstruct the spatial distribution of absorbed optical energy inside the tissue. Common methods applied are filtered backprojection (FBP) [22, 14, 40, 39] or time reversal [18, 34]. In this work, we develop a general variational framework that incorporates prior knowledge on the reconstruction, offering higher quality reconstructions than direct methods and providing robustness against noise and compressive sampling.

The optical absorption coefficient in tissue varies due to the relative presence of various components such as haemoglobin, oxy-haemoglobin, melanin, lipids and water. These chromophores have specific spectral signatures, and the use of multi-wavelength PAT can potentially reveal molecular specific information. This is directly linked to tissue physiopathology and function and has diagnostic potential. Photoacoustics is being researched for applications in various fields of biomedicine [43] such as brain imaging in small animals [9], breast cancer imaging in humans [17, 33], and imaging inflamed human synovial joints in rheumatoid arthritis [7, 36]. In most of these applications, detection of signals from haemoglobin and oxy-haemoglobin enables the visualisation of blood vessels. Several disorders are characterised by a dysregulation of blood vessel function, but also with uncontrolled creation of blood vessels or angiogenesis.

The standard methods of filtered backprojection (FBP) and time reversal suffer from curved line artefacts, especially when the noise level is high or the placement of detectors for measurements (sampling) is coarse [28]. Recently, there has been intense interest in solving the PAT reconstruction problem iteratively with a specific focus on regularised reconstruction. The regularisers used for these reconstructions have many different forms, depending on the prior assumptions on the image. Total variation (TV) regularisation is a powerful tool for noise removal [30], and generates the desired images with sharp edges. In [2], TV-regularised reconstructions on simulated data are presented, where an analytical adjoint operator for the k-space forward operator [34] is derived to model wave propagation. In [1] a similar reconstruction method is applied to compressively sampled simulated and experimental data. Although from the theory it is not obvious [27], the TV basis appears to work well in the area of compressive sampling for the reconstruction of an absorbed energy density map with abrupt edges [12]. Other works that combine total variation with compressive sampling in PAT can be found in [42] and [24].

In case of a heterogeneous fluence rate, such as an exponentially decaying one, a piecewise constant prior on the reconstruction is not realistic. Higher-order total variation methods, such as total generalised variation (TGV) [5] are better suited to deal with such data. The use of TGV shows great promise for image reconstruction in tomographic settings, such as MRI [20], CT [41, 25], ultrasound waveform tomography [23] and diffusion tensor imaging [35]. Image denoising with TGV as a post-processing step has made its way into optical coherence tomography [10] and optical diffraction tomography [21]. However, solving these optical reconstruction problems with TGV regularisation within the reconstruction algorithm is still an open research topic. A recent report [15] shows TGV-regularised reconstructions in PAT on vascular data. Here, measurements are taken by making snapshots of the acoustic waves using a CCD camera.

Vascular structures show a strong anisotropy, that can be promoted by using directional wavelets [19] or curvelets [6] as ‘building blocks’ for the reconstruction. Provost and Lesage [27] have shown that wavelets and curvelets are very sparse representations of the measurement data, and can help to recover anisotropic features in the reconstruction. This has been confirmed in [12] through new phantom and in vivo experiments.

To our knowledge, a general reconstruction framework for PAT, in which regularisers can be easily exchanged, does not exist. In this work, we give such a variational framework. Its applicability is demonstrated by using it for three reconstruction models, which contain TV, TGV and directional wavelet regularisers. While others based their methods on specific measurement geometries, our framework makes no assumption on the geometry, requires no specific physics modelling and can be applied to the PAT reconstruction problem in a 2- or 3-dimensional setting. We develop not only a TV, but also a TGV regularised reconstruction model, which can deal with the more realistic assumption of a decaying heterogeneous fluence rate. A reconstruction model using directional wavelet regularisation is developed for images containing anisotropic vascular structures. The numerical implementation with a primal-dual algorithm is a modular one: the regularisation or the data-fidelity can be chosen differently without having to change the structure. The algorithm has a convergence rate that is at least linear and that can be accelerated [26] for specific choices of the regulariser. The performance of the framework is demonstrated on synthetic data, as well as on experimental data. Robustness against noise and compressive sampling is shown.

The remainder of this paper is organised as follows: in section 2 the variational method for the regularised reconstruction of PAT is derived. After writing the reconstruction problem as a saddle-point problem, the numerical implementation with the first-order primal-dual algorithm PDHGM is explained in section 3. In section 4, the experimental setup that is employed to generate our results is explained. Moreover, the specific forward model that is used for obtaining the synthetic data is derived. The digital and experimental phantoms that are created as test cases for our reconstruction framework are shown in section 5. After an extensive analysis of the method on challenging synthetic data in section 6 and experimental data in section 7, we conclude and discuss the results in section 8.

## 2 Variational methods using regularisation

We consider the image reconstruction problem of finding an estimate for the absorbed energy density map from given (preprocessed) measurements . Here , where is the dimension of the space in which must lie and is the space-time domain of the detection surface between times and . Mathematically, the inverse problem can be formulated as the solution to the equation

(1) |

where is a bounded linear operator. This operator can be described in many ways, since there are multiple methods that model the forward process in photoacoustic imaging. For instance, it can be described as a k-space method that models linear acoustic propagation in time [34]. In our simulations and experiments, we use the spherical mean operator [22] to model the forward process.

Instead of using a direct reconstruction method such as FBP, we solve the inverse problem in a variational setting. We consider the following minimisation problem:

(2) |

The first term is the data fidelity term, which ensures that the reconstructed image, after applying the forward operator , is close to the noisy data . The -norm has been chosen because the photoacoustic measurements are predominantly affected by noise from the system electronics and thermal noise from the transducers, which are known to be of additive Gaussian type [42]. The second term is an weighted regularisation term which becomes large when our prior assumptions on are violated. Note that this is a very general framework: as long as is convex, we can put it in our numerical framework, which is explained in section 3.

Note that if we do not have any prior assumption on the data and we choose , we obtain for (2) the least squares (LS) reconstruction. The LS reconstruction can be evaluated by using a gradient descent algorithm or, if is a small enough matrix operator, by solving directly. However, since LS is very sensitive to noise in the data , Tikhonov regularised least squares (LS-T) is often used: the term in (2) then takes the form . This can again be evaluated directly by solving . In section 6, our models will be compared with LS-T.

### 2.1 Regularisation with total (generalised) variation

In this work we use the concept of total generalised variation (TGV), as introduced by Bredies et al. [5]. TGV is a generalisation of TV, in which higher order derivatives are considered instead of only first order derivatives in TV. Bredies et al. proposed the following definition:

(3) |

With the choice of , the desired order of regularity can be obtained, while the parameter gives a weight on every regularity level. By choosing and , the definition of TGV coincides with the definition of TV

(4) |

In this work, we consider and as regularisers. In our applications, we work with a discrete image , on which a gradient or derivative is always numerically defined. Therefore, if we assume or , the expressions for and respectively can be simplified. After substituting these expressions in (2), we obtain the following minimisation problems:

(5) | |||

(6) |

where and is the symmetrised gradient [5]. In the TV-regularised functional (5), the parameter determines the smoothness of the solution. In the TGV-regularised functional (6) one can choose the influence of first order and second order smoothness by choosing combinations of and . The effect of different combinations of and on the reconstruction of one-dimensional TV and TGV eigenfunctions was analysed in [3]. From now on, when TGV is mentioned, we mean the second order method .

### 2.2 Regularisation with directional wavelets

Blood vessel visualisation is an important aspect for the application of photoacoustic imaging in biomedicine. An image containing vascular structures can be sparsely represented by only a small number of elements in some basis or in a so called dictionary. There is a big amount of possibilities to represent this anisotropic data. A rough distinction can be made [29] between analytic transforms, such as the Gabor and wavelet transform; analytic dictionaries, such as curvelets and contourlets; and trained dictionaries, such as sparse PCA and K-SVD.

In this work, we make use of an enhancement of the discrete wavelet transform, but as our framework is very general, it is easy to replace it with a transform of dictionary that suits the data best. The discrete wavelet transform is in essence a one-dimensional transform and the obvious extension to a multi-dimensional setting is not satisfactory, since it lacks directionality. The dual-tree complex wavelet transform [19] does have directionally selective filters and is thus better able to represent multi-dimensional directional data. We propose the following minimisation problem

(7) |

where is the complex dual-tree wavelet transform and is the corresponding wavelet domain. Note the close similarity in structure between (5) and (7). Because of this similarity, we will not elaborate on the wavelet case in the numerical treatment of our models, but only mention the small changes that have to be made to the TV implementation to obtain the wavelet implementation.

## 3 Numerical Framework

In section 2 we considered the convex, possibly non-smooth minimisation problem (2). TV and TGV are defined via a maximisation problem over their dual variable(s), which can be seen in (3) and (4). Substituting the definition of TV/TGV in the minimisation problem changes it into a saddle-point problem, which can be solved through the use of primal-dual algorithms. These primal-dual algorithms can also be used for our wavelet regulariser or other regularisers that are not defined via their dual variables: with the use of Fenchel duality, the specific problem may be rewritten to its dual maximisation problem or its saddle-point problem.

### 3.1 Description of the saddle-point problem

With the use of Fenchel duality, we can rewrite the general minimisation problem (cf. (2))

(8) |

into the general saddle-point problem

(9) |

where is the convex conjugate of . We obtain our TV, TGV and wavelet models by describing the functionals and and the operator appropriately, which is done in (3.1) and (11).

(10) | ||||

(11) | ||||

Here we made a splitting and , in the same manner as in [31]. Note that we use tildes on the operators that appear in their convex conjugate form in the primal-dual problem. The convex conjugate, also known as Fenchel–Legendre transform, is defined as

(12) |

### 3.2 Proximal operators for TV, TGV and wavelet model

A key tool for the use of primal-dual algorithms is the so called proximal splitting method. The proximal operator for a scaled functional is defined as

(13) |

which gives a small value of , while keeping close to . Because the right hand side consists of a strongly convex quadratic part and a convex functional , it is strongly convex and thus (13) has a unique solution.

As can be seen in section 3.3, the saddle-point problems (3.1) and (11) now have the correct structure to be numerically implemented. Therefore, the proximal operators for all and are needed, which are derived below. The convex conjugate for is

(14) |

After checking the first and second order optimality conditions, we find that maximises (14). We fill this in and calculate the proximal operator

(15) |

The convex conjugate for is

(16) |

Here we choose the pointwise 2-norm of , which means for . This choice has the effect that the TV-eigenfunction will be a -dimensional sphere instead of a -dimensional cube (1-norm) or a -dimensional diamond (-norm) [11, Theorem 4.1]. For the wavelet model, is a scalar function instead of a vector function, which means that all pointwise -norms are equal. The proximal operator for is calculated

(17) |

Similarly, the proximal operator for reads

(18) |

It is easily seen that and .

### 3.3 First-order primal dual algorithm

We make use of the modified primal-dual hybrid gradient algorithm (PDHGM) as proposed by Chambolle and Pock [8], which can be seen as a generalisation to the PDHG algorithm [44]. PDHGM has the advantage that it has a direct solution at every step. Furthermore it can be shown that PDHGM has a convergence rate of at least and can even be , where is the number of steps in the algorithm [8, 26].

In section 3.1 we have split the functionals and in multiple parts. Therefore, we also need to split the proximal operators and evaluate them separately. A similar algorithm to [31, Algorithm 4] is used.

Algorithm 1 shows the implementation of the TV model. We perform two separate dual steps, which are both used as input for the primal step after that. The implementation of the wavelet model is obtained by changing and div to and respectively.

Algorithm 2 shows the implementation of the TGV model, which is built in the same spirit as the TV algorithm.

#### 3.3.1 Discretisation

The measurement data is sampled in both space and time: the detectors are modelled as point detectors at specific locations and measure the pressure at sampled time instances with interval . The image that we wish to reconstruct is also a discretised one. In case these discretisation do not perfectly match (which is often the case), we might get significant errors. To avoid these errors, we use the interpolation procedure as explained in [39, Chapter 5].

Since we are using an iterative procedure, it is computationally expensive to evaluate the discrete integral operator at every step. Therefore, we calculate the discretised matrix once and use it iteratively. Instead of analytically determining and then discretising it, we use the adjoint matrix . Accordingly, we also use matrix versions for the operators and . For the complex dual-tree wavelet transform and its adjoint (which in this case is equivalent to its inverse), we make use of the implementation by Cai and Li^{2}^{2}2available at http://eeweb.poly.edu/iselesni/WaveletSoftware/dt2D.html.

#### 3.3.2 Algorithm parameter selection

In the work of Chambolle and Pock [8], it is shown that the PDHGM algorithm always converges if . Here is the combined operator as defined in (3.1) and (11). A uniform bound for and might be undesired if the separate operators and , or have very different norms: all the proximal steps have a small step size, while this might not be needed for the steps related to only one of these operators.

In [26] this problem has been solved for matrix operators. Instead of using one value for and in the PDGHM algorithm, one can use diagonal matrices and that satisfy . With this choice, the sequence generated by the algorithm weakly converges to an optimal solution of the saddle-point problem.

For our specific operators, it appears that , but they show a very similar structure within each operator. Therefore, we only choose two sets of two operators, since they have shown to give a smoother convergence plot than with the parameter choice in [26], while the convergence rate is similar. We choose one set of two parameters for the operator and one for the combination of other operators. More precisely, we choose such that and . Here , or in the TV or wavelet model and is the norm of the combined other operators in the TGV model.

It is easily shown that the bound is satisfied by this choice if we write as a concatenation of the matrix version of all operators. For the TGV model, we have

and the diagonal matrices

This gives us the estimate

(19) |

A similar estimate can be given for the TV and wavelet model.

#### 3.3.3 Algorithm validation

The PDHGM algorithm with its specific discretisation and parameter choices is validated by looking at two criteria. Firstly, the duality gap is taken into consideration. The duality gap gives the (non-negative) difference between the primal and the dual functional. As the duality gap approaches zero, the solution gets asymptotically close to the desired solution for the primal-dual problem [4]. Hence the algorithm is validated by checking if the duality gap approaches zero. Moreover, a stopping criterion for the algorithm can be constructed based on the duality gap. In the formulation of (9), the duality gap reads

(20) |

The second criterion is the residual of the variables, defined as

(21) |

where any variable can be filled in instead of . The residual for whenever , which means that the sequence converges.

## 4 Experimental setup

The proposed methods are defined for solving the general problem . Here is the desired reconstruction, is (preprocessed) data and is some operator that maps to . Because of this generic structure, the methods can be applied to reconstruct images for many photoacoustic tomographs with many different forward models.

In our experimental setup, we make use of the finger joint imager as specified in [37]. A side-illumination laser delivers pulses of optical energy with a wavelength of 532 nm. For the recording of pressure waves a 1D piezoelectric detector array with 64 elements in a half-circle is used. The detector array has a central frequency of 7.5 MHz. Both the detectors and fibre bundles are simultaneously rotatable over 360 degrees. They can also be translated in order to image multiple slices. The detectors have a narrow focus (0.6 mm plane thickness) in one dimension, making it suitable for 2D slice based imaging. A schematic overview of the experimental setup is shown in Figure 1. For a more extensive explanation of the setup, we refer to [37].

### 4.1 Forward model

Based on the work of [22] and [38], we write our forward model as the projection of the absorbed energy distribution over spherical surfaces. Under the assumption of a homogeneous speed of sound , we obtain the following expression for the pressure measured by the detector at location :

(22) |

where is the thermal expansion coefficient, is the specific heat, is the illumination profile of the laser, is the impulse response of the detectors and is the absorbed energy. Instead of finding an explicit expression for and , we make use of the calibration measurement of an approximate delta pulse [38] located at , i.e. . After finding an analytic expression for in terms of and [39], we obtain

(23) |

where is the retarded time .

### 4.2 Preprocessing for reconstruction

We can write (23) concisely as

(24) |

We left the dependencies on , and out for an uncluttered notation; the convolution with is performed in the time domain. Recall that we are solving the problem , where in this specific setup, is defined as in (24). In order to get our preprocessed data , we first solve the deconvolution problem

(25) | |||

(26) |

After applying the convolution theorem to (25), we obtain , which for is equivalent to

(27) |

Unfortunately, the right hand side of (27) is undefined when . Moreover any noise on can have a big influence on the expression. Therefore, we construct a regularised deconvolution filter

such that

(28) |

Here is a small parameter such that it suppresses noise on , but leaves higher parts of it intact. By consecutively applying (28) and (26), we obtain our preprocessed measurement data for (1).

## 5 Test objects

In order to test our reconstruction framework with its regularisers, both digital and experimental phantoms have been created. For any of the three regularisers a specific digital phantom has been created, that fulfils the prior assumption on the image to be reconstructed. An experimental phantom with a vascular structure has been created to test the effectiveness of the different regularisers on real data.

### 5.1 Digital phantoms

The first digital phantom is shown in Figure 1(a). For this phantom, it is assumed that the fluence rate in the illuminated object is homogeneous. This means that sharp changes in the absorbed energy are expected and it thus consists of piecewise constant values. This image consists of both large and small objects with a variety of shapes, among which squares, discs and elongated rectangles. The larger discs could resemble dermis in the finger, while the small discs and rectangles could resemble blood vessels respectively perpendicular and in the 2D-plane of focus. Because of the sharp discontinuities, we will reconstruct this image with the TV method and compare this with FBP and LS-T reconstructions.

The second digital phantom is shown in Figure 1(b). For this phantom, it is assumed that the fluence rate is decaying from the edge of the disc to the middle of the disc. It resembles tissue that absorbs enough light such that the absorbed energy decays as we are deeper inside the tissue. Since the fluence rate is assumed heterogeneous, the TGV method will be used for the reconstruction.

The third digital phantom, shown in Figure 1(c), is a vascular structure, which has been obtained by preprocessing part of a retinal image from the DRIVE database [32]. Such a structure can be expected when blood vessels lie in the plane of focus. Moreover, such a phantom is similar to 3D vascular structures that can be expected in photoacoustic breast imaging.

### 5.2 Synthetic data acquisition

The acquisition of the preprocessed data from the measured pressure has a strong dependency on the calibration measurement , as can be seen in section 4.2. If we want to simulate the acquisition of , we will also have to make use of a calibration measurement , as can be seen in (23). In order to avoid an inverse crime, we will use a different calibration measurement for the forward model as for the reconstruction. Moreover, the discretisation of the ground truth image will be different from the discretisation of the reconstructed image .

### 5.3 Experimental phantoms

In order to experimentally test the capability of the regularisers, we developed a phantom with absorbers that resemble a vascular structure. The vessel-shaped absorbers (filaments) were constructed of sodium alginate (SA) gel carrying iron oxide nanoparticles to impart absorption. For the latter we used commercially available superparamagnetic iron oxide (SPIO) nanoparticles (Endorem - Guerbet, Villepinte, France). A dilute Endorem dispersion was mixed with SA solution in distilled water to arrive at a final 2% (w/v) of SA solution. The filaments were fabricated by extruding the SA solution through a syringe with a 30g needle and allowing to fall into 0.7 M calcium chloride (CaCl) solution. SA undergoes gel formation in the presence of Ca ions in water. The resulting gel sinks to the bottom and hardens for 15 minutes. Finally, the filaments are isolated and washed three times with distilled water to ensure removal of residual Ca.

The cylindrical phantom (diameter 24 mm) was made of Agar gel with Intralipid to provide optical scattering. To create the gel, first 3% (v/v) Agar was dissolved in water by boiling it in the microwave until a clear solution had been formed. Next, the temperature of the mixture was decreased to 40C under continuous stirring. A 3% (v/v) aqueous solution of 20% stock Intralipid was added drop-wise with stirring. This provides a reduced scattering coefficient (s’) of 9.7 cm at 532 nm.

To prepare the phantom, the Agar solution was poured halfway into a tube as mould and allowed to harden. A certain number of filaments were then laid on the surface of the stiff Agar gel (Figure 2(a)). The whole was then topped up with more Agar solution by carefully pouring. When the absorbing agar solution had set, the phantom was removed from the mould.

As the reconstructions will later show, some of the filaments moved during the pouring of the Agar solution. For later comparisons, we have registered (aligned) the photograph of the experimental phantom with the reconstruction. This has been achieved with an affine and B-spline grid based registration, for which a toolbox is available on Matlab Central^{3}^{3}3Available at https://nl.mathworks.com/matlabcentral/fileexchange/20057-b-spline-grid--image-and-point-based-registration. The registered image is shown in Figure 2(b). The registered image has been segmented to enable future analysis. The segmented image is shown in Figure 2(c). The thin lines in Figure 2(b) are registration artefacts and have been removed. The segmentation serves as a (quasi) ground truth reference for validating the morphology of our experimental reconstructions.

### 5.4 Quality measures

For the synthetic phantoms, a digital ground truth image is available to compare the reconstructions with. Besides a visual comparison, we make use of the peak signal-to-noise ratio (PSNR), defined as

where is the reconstructed image, the ground truth image and the number of pixels in the image.

A quality measure based on image intensities is not possible for the experimental reconstructions, since no digital ground truth image is available. However, when interested in the geometry and location of the vascular structure, only the segmentation of the reconstruction is of importance. This reconstruction segmentation can be compared with the ground truth segmentation (Figure 2(c)). We follow this idea and make use of the receiver operating characteristic (ROC) curve. The ROC curve is obtained by plotting the false positive rate (one minus specificity) against the true positive rate (sensitivity) of the thresholded reconstruction for various threshold values. The true positive rate gives the fraction of pixels inside the ground truth segmentation that are correctly classified as such by the reconstruction segmentation. The false positive rate gives the fraction of pixels outside the ground truth segmentation that are incorrectly classified as such by the reconstruction segmentation. As the threshold for the reconstruction varies, more pixels will be correctly segmented, at the cost of incorrectly segmented pixels. An ideal situation would be one where the true positive rate increases, without changing the false positive rate, i.e. an almost vertical line, followed by an almost horizontal line.

## 6 Synthetic results

Our reconstruction framework has been tested on both synthetic and experimental data. It is compared with two standard direct reconstruction methods, namely filtered backprojection (FBP) and Tikhonov-regularised least squares (LS-T), cf. section 2. The specific FBP algorithm that was used is explained in [39]. In the synthetic case, for the LS-T, TV, TGV and wavelet methods, the regularisation parameters are chosen such that the PSNR between the ground truth and the reconstruction are best. In case of data with additive Gaussian noise, the optimal regularisation parameters can be found with an L-curve method [16]. Robustness against noise and compressive sampling will be shown.

### 6.1 Robustness under compressive sampling

Synthetic preprocessed measurement data have been obtained according to (23) with a uniform sampling of 64 dectors over a 172 degree array. Six rotations of 60 degrees have been performed, giving us an almost uniform sampling in a total of 384 detection locations. For the first comparison, different reconstruction methods under uniform compressive sampling are considered. It is unclear if uniform compressive sampling gives the best reconstruction quality with respect to the different reconstruction methods (see e.g. Haber et al. on experimental design [13]). However, it intuitively makes sense to place the detectors uniformly. Moreover, since the detector need some space, they are often placed in such a fashion.

In Figure 4 the results of three reconstruction methods (FBP, LS-T, TV) are shown for coarse sampling (16 detectors), moderate sampling (64 detectors) and fine sampling (384 detectors). All methods give a good reconstruction in case of high sampling, although every method has its own reconstruction bias: TV gives contrast loss in smaller structures, LS-T gives blurred structure edges, whereas FBP gives curved line artefact in the whole image. The TV method is able to keep important features for a longer time as we sample with less detectors: with moderate sampling, it is still possible to detect the minor contrast changes in the upper right square, whereas this is not possible in the other reconstructions. With coarse sampling, the TV reconstruction still gives sharp edges and structures with the right intensity, while this is not the case in the other reconstructions.

For the second digital phantom, we again compare with FBP and LS-T. The reconstructions using a uniform sampling of 16 detectors are shown in Figure 5. The TGV method gives the desired smooth intensity within the disc, while keeping the sharp discontinuities. Moreover, it does not show the curved line artefacts that are visible in the FBP reconstruction.

In Figure 6, the reconstructions of the third phantom using a uniform sampling of 32 detectors have been shown. FBP performs poorly for this vascular data set, since it is not completely clear which parts of the reconstruction are vascular structures and which are curved line artefacts. The LS-T and wavelet reconstructions show less pronounced curved line artefacts than the FBP reconstruction. The wavelet reconstruction additionally has a higher intensity within the vessels and shows a smoother background. It is striking that regularisation with wavelets is not as effective in removing the curved line artefacts as previous regularisation with TV or TGV, as can be seen around and in Figure 6. The probable reason for this is that not only the vascular structure, but also the curved line artefacts can be sparsely represented in the directional wavelet basis. Note that this is only the case in a 2D setting, since in 3D, a backprojection artefact looks like part of the surface of a sphere instead of a curved line. Directional wavelets might be much more effective in removing these 3D artefacts, as long as a suitable transform is used that can sparsely represent anisotropic vascular structures.

In Figure 7, a plot of the PSNR values for different reconstruction methods under compressive sampling is shown. Both the TV and the TGV method perform better than FBP and LS-T. The wavelet method gives a minor improvement under compressive sampling, probably due to the problem that artefacts are too similar to the vascular structure. No noise has been added to the data for this comparison.

### 6.2 Robustness against noise

As explained in section 2, we expect mainly additive Gaussian noise, because of the system electronics and thermal noise from the transducers. After the simulated measured pressure has been obtained, Gaussian noise with zero mean and standard deviation was added to the measured pressure. In Figure 8, a comparison between reconstruction methods for the noisy data has been shown. The TV, TGV and wavelet method outperform FBP and LS-T for both low and high noise levels. The results were obtained with a uniform sampling of 192 detectors.

## 7 Experimental results

For the acquisition of experimental data, the cylindrical phantom was scanned in six rotations, giving a total uniform sampling of 384 detectors. We compare the TV, TGV and wavelet reconstruction methods with FBP by using both the full data (384 detectors) and using 16.7% compressively sampled data (64 detectors) for the reconstructions. The LS-T reconstruction is omitted, since under low noise and middle to high sampling, FBP has shown to outperform LS-T for the synthetic case. The regularisation parameters for the regularised reconstruction have been chosen similar to the parameters found in the synthetic tests, with a small deviation such that the reconstruction is visually best.

In Figure 9 it can be seen that FBP gives a sharp reconstruction with the expected curved line artefacts in the background. Due to these artefacts, the intensity of the anisotropic structure is not uniform along the structure and undesired low intensity gaps are created. All three regularisation methods give a much smoother reconstruction, both in the background and in the anisotropic structure. Moreover, the small gaps of the FBP reconstruction are filled. The methods give the expected result: a piecewise constant reconstruction for TV, a piecewise linear reconstruction for TGV and an anisotropy enhanced reconstruction for wavelet. Since the dual-tree wavelet approach makes use of smooth elongated structure as it’s ‘building blocks’, this reconstruction looks very similar to the TGV reconstruction, although the anisotropic structure shows a slightly sharper boundary in the wavelet reconstruction.

From the compressive sampling reconstruction (Figure 10) we can draw similar conclusions: the overall reconstructions are smoother and many gaps that are apparent in the FBP reconstruction are filled in the other reconstructions. If we focus on the region below and between and , it is clear that FBP gives strong curved line artefacts, while with the other methods these artefacts are clearly reduced. In these reconstructions, there is much less similarity between the TGV and wavelet results. Due to the use of direction elements, the wavelet reconstruction presents a much better connected structure than the TGV reconstruction.

As can be seen in Figure 10(a), the ROC curves of the regularised reconstructions show an almost ideal behaviour. The ROC curve of the FBP reconstruction lies much lower, which shows that this method is a poor choice when it comes to thresholding capability. In Figure 10(b), we see a similar relation between the reconstruction methods, although all curves lie lower than in Figure 10(a). It is interesting to see that all regularised reconstruction methods give a similar ROC curve. TGV and wavelet almost completely overlap, while TV is a bit higher in the steep part, but continues slightly lower in the second part of the graph. This tells us that either method will help us to correctly segment the photoacoustic reconstruction for anisotropic vascular structures.

## 8 Summary and outlook

In this work we have derived a general variational framework for regularised PAT reconstruction where no specific measurement geometry or physics modelling is assumed. The framework can be applied to a PAT reconstruction problem in a 2- or 3-dimensional setting. The primal-dual implementation is a modular one and the optimisation parameters are chosen such that the algorithm is efficient. Specific choices for the regulariser made in this paper are TV, second order TGV and the complex dual-tree wavelet transform. These regularisers can respectively deal with the prior assumptions of either sharp discontinuities in combination with a homogeneous and heterogeneous fluence rate, or anisotropic structures. The methods have been compared on both synthetic and experimental data with two widely used direct reconstruction methods.

We have shown that the modularity of the variational framework enables easy exchange of regularisers, which arise from different prior assumptions on the reconstruction. Using this framework, one can choose the terms within as desired, without having to change the structure of the algorithm.

Furthermore it was shown that our TV and TGV methods perform better than direct reconstruction methods: they are better able to handle both low and high noise levels and give better reconstructions under uniform compressive sampling. In a clinical setup, it might not be preferred to change the detector locations during a full scan, because of limited view or movement of the tissue under consideration. For this reason, the use of TV and TGV methods are promising. In this work, reconstructions were made using uniform compressive sampling. Future efforts could lie in finding the best sampling strategy for various data sets and regularisers.

Finally, research could be focussed on the reconstruction of anisotropic structures such as blood vessels and the difference between this reconstruction in 2D and 3D. It is unclear if the similarity between backprojection artefacts and vascular structures is preventing better reconstructions in the 2D case and if this problem is absent in the 3D case. For this reason, it is useful to develop or learn a dictionary that allows sparse reconstruction of vascular structures in 3D.

## Acknowledgements

The authors thank Maura Dantuma for the acquisition of the experimental data. Marinus J. Lagerwerf acknowledges financial support from the Netherlands Organization for Scientific Research (NWO), project 639.073.506. The authors acknowledge Stichting Achmea Gezondheidszorg for funding in project Z620, and the Pioneers in Healthcare Innovation (PIHC) fund 2014 for funding in project RAPACT.

## References

- [1] S. R. Arridge, P. Beard, M. M. Betcke, B. Cox, N. Huynh, F. Lucka, O. Ogunlade, and E. Zhang. Accelerated high-resolution photoacoustic tomography via compressed sensing. Phys. Med. Biol., 61(24):8908–8940, dec 2016.
- [2] S. R. Arridge, M. M. Betcke, B. T. Cox, F. Lucka, and B. E. Treeby. On the Adjoint Operator in Photoacoustic Tomography. Inverse Probl., 32(11):115012, nov 2016.
- [3] M. Benning, C. Brune, M. Burger, and J. Müller. Higher-order TV methods - Enhancement via Bregman iteration. J. Sci. Comput., 54(2-3):269–310, 2013.
- [4] S. Boyd and L. Vandenberghe. Convex Optimization. Cambridge University Press, 7th editio edition, 2004.
- [5] K. Bredies, K. Kunisch, and T. Pock. Total Generalized Variation. SIAM J. Imaging Sci., 3(3):492–526, 2010.
- [6] E. J. Candès and D. J. Donoho. Curvelet: A surprising effective non-adaptive representation for objects with edges. Department of Statistics. Technical report, 1999.
- [7] D. Chamberland, Y. Jiang, and X. Wang. Optical imaging: new tools for arthritis. Integr. Biol., 2(10):496–509, 2010.
- [8] A. Chambolle and T. Pock. A first-order primal-dual algorithm for convex problems with applications to imaging. J. Math. Imaging Vis., 40(1):120–145, 2011.
- [9] X. L. Deán-Ben, S. Gottschalk, B. Mc Larney, S. Shoham, and D. Razansky. Advanced optoacoustic methods for multiscale imaging of in vivo dynamics. Chem. Soc. Rev., 46:2158–2198, 2017.
- [10] J. Duan, W. Lu, C. Tench, I. Gottlob, F. Proudlock, N. N. Samani, and L. Bai. Denoising optical coherence tomography using second order total generalized variation decomposition. Biomed. Signal Process. Control, 24:120–127, 2016.
- [11] S. Esedoglu and S. J. Osher. Decomposition of Images by the Anisotropic Rudin – Osher – Fatemi Model. In Commun. Pure AppliedMathematics, volume LVII, pages 1609–1626. Wiley Periodicals, Inc., 2004.
- [12] Z. Guo, C. Li, L. Song, and L. V. Wang. Compressed sensing in photoacoustic tomography in vivo. J. Biomed. Opt., 15(2):021311, 2010.
- [13] E. Haber, L. Horesh, and L. Tenorio. Numerical methods for experimental design of large-scale linear ill-posed inverse problems. Inverse Probl., 24(5):2–17, 2008.
- [14] M. Haltmeier, T. Schuster, and O. Scherzer. Filtered backprojection for thermoacoustic computed tomography in spherical geometry. Math. Methods Appl. Sci., 28(16):1919–1937, 2005.
- [15] K. Hammernik, T. Pock, and R. Nuster. Variational photoacoustic image reconstruction with spatially resolved projection data. In A. A. Oraevsky and L. V. Wang, editors, Photons Plus Ultrasound Imaging Sens. 2017, Proc. SPIE, volume 10064, page 100643I, mar 2017.
- [16] P. C. Hansen. The L-Curve and its Use in the Numerical Treatment of Inverse Problems. Comput. Inverse Probl. Electrocardiology, ed. P. Johnston, Adv. Comput. Bioeng., 4:119–142, 2000.
- [17] M. Heijblom, D. Piras, M. Brinkhuis, J. C. G. van Hespen, F. M. van den Engh, M. van der Schaaf, J. M. Klaase, T. G. van Leeuwen, W. Steenbergen, and S. Manohar. Photoacoustic image patterns of breast carcinoma and comparisons with Magnetic Resonance Imaging and vascular stained histopathology. Sci. Rep., 5(May):11778, 2015.
- [18] Y. Hristova, P. Kuchment, and L. Nguyen. Reconstruction and time reversal in thermoacoustic tomography in acoustically homogeneous and inhomogeneous media. Inverse Probl., 24(5):55006, 2008.
- [19] N. Kingsbury. The dual-tree complex wavelet transform: a new efficient tool for image restoration and enhancement. Signal Process. Conf. (EUSIPCO 1998),, pages 2–5, 1998.
- [20] F. Knoll, K. Bredies, T. Pock, and R. Stollberger. Second order total generalized variation (TGV) for MRI. Magn. Reson. Med., 65(2):480–491, 2011.
- [21] W. Krauze, P. Makowski, M. Kujawińska, and A. Kuś. Generalized total variation iterative constraint strategy in limited angle optical diffraction tomography. Opt. Express, 24(5):4924, 2016.
- [22] R. Kruger, P. Liu, Y. Fang, and A. Robert. Photoacoustic ultrasound (PAUS)—Reconstruction tomography. Med. Phys., 22(10):1605, 1995.
- [23] Y. Lin and L. Huang. Ultrasound waveform tomography with the second-order total-generalized-variation regularization. Proc. SPIE 9783, Med. Imaging 2016 Phys. Med. Imaging, 9783(1):1–7, 2016.
- [24] J. Meng, L. V. Wang, L. Ying, D. Liang, and L. Song. Compressed-sensing photoacoustic computed tomography in vivo with partially known support. Opt. Express, 20(15):16510, jul 2012.
- [25] S. Niu, Y. Gao, Z. Bian, J. Huang, W. Chen, G. Yu, Z. Liang, and J. Ma. Sparse-view x-ray CT reconstruction via total generalized variation regularization. Phys. Med. Biol., 59(12):2997–3017, 2014.
- [26] T. Pock and A. Chambolle. Diagonal preconditioning for first order primal-dual algorithms in convex optimization. Proc. IEEE Int. Conf. Comput. Vis., pages 1762–1769, 2011.
- [27] J. Provost and F. Lesage. The application of compressed sensing for photo-acoustic tomography. IEEE Trans. Med. Imaging, 28(4):585–594, 2009.
- [28] A. Rosenthal, V. Ntziachristos, and D. Razansky. Acoustic Inversion in Optoacoustic Tomography: A Review. Curr. Med. Imaging Rev., 9:318–336, 2013.
- [29] B. R. Rubinstein, A. M. Bruckstein, and M. Elad. Dictionaries for Sparse Representation Modeling. 98(6):1045–1057, 2010.
- [30] L. I. Rudin, S. Osher, and E. Fatemi. Nonlinear total variation based noise removal algorithms. Phys. D Nonlinear Phenom., 60(1-4):259–268, 1992.
- [31] E. Y. Sidky, J. H. Jørgensen, and X. Pan. Convex optimization problem prototyping for image reconstruction in computed tomography with the Chambolle-Pock algorithm. Phys. Med. Biol., 57(10):3065–91, 2012.
- [32] J. Staal, M. D. Abràmoff, M. Niemeijer, M. A. Viergever, and B. Van Ginneken. Ridge-based vessel segmentation in color images of the retina. IEEE Trans. Med. Imaging, 23(4):501–509, 2004.
- [33] M. Toi, Y. Asao, Y. Matsumoto, H. Sekiguchi, A. Yoshikawa, M. Takada, M. Kataoka, T. Endo, N. Kawaguchi-Sakita, M. Kawashima, E. Fakhrejahani, S. Kanao, I. Yamaga, Y. Nakayama, M. Tokiwa, M. Torii, T. Yagi, T. Sakurai, K. Togashi, and T. Shiina. Visualization of tumor-related blood vessels in human breast by photoacoustic imaging system with a hemispherical detector array. Sci. Rep., 7(December 2016):41970, 2017.
- [34] B. E. Treeby and B. T. Cox. k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields. J. Biomed. Opt., 15(2):021314, 2010.
- [35] T. Valkonen, K. Bredies, and F. Knoll. Total Generalized Variation in Diffusion Tensor Imaging. SIAM J. Imaging Sci., 6(1):487–525, 2013.
- [36] P. van Es, S. K. Biswas, H. J. Bernelot Moens, W. Steenbergen, and S. Manohar. Initial results of finger imaging using photoacoustic computed tomography. J. Biomed. Opt., 19(6):060501–1–4, 2014.
- [37] P. Van Es, R. Vlieg, S. Biswas, E. Hondebrink, J. Van Hespen, H. Moens, W. Steenbergen, and S. Manohar. Coregistered photoacoustic and ultrasound tomography of healthy and inflamed human interphalangeal joints. Prog. Biomed. Opt. Imaging - Proc. SPIE, 9539:14–16, 2015.
- [38] Y. Wang, D. Xing, Y. Zeng, and Q. Chen. Photoacoustic imaging with deconvolution algorithm. Phys. Med. Biol., 49(14):3117–3124, 2004.
- [39] G. R. Willemink. Image reconstruction in a passive element enriched photoacoustic tomography setup. PhD thesis, University of Twente, Enschede, 2010.
- [40] M. Xu and L. V. Wang. Universal back-projection algorithm for photoacoustic computed tomography. Phys. Rev. E - Stat. Nonlinear, Soft Matter Phys., 71(1):1–7, 2005.
- [41] J. Yang, H. Yu, M. Jiang, and G. Wang. High Order Total Variation Minimization for Interior Tomography. Inverse Probl., 26(3):350131–3501329, 2010.
- [42] Y. Zhang, Y. Wang, and C. Zhang. Total variation based gradient descent algorithm for sparse-view photoacoustic image reconstruction. Ultrasonics, 52(8):1046–1055, 2012.
- [43] Y. Zhou, J. Yao, and L. V. Wang. Tutorial on photoacoustic tomography. J. Biomed. Opt., 21(6):061007, 2016.
- [44] M. Zhu and T. Chan. An efficient primal-dual hybrid gradient algorithm for total variation image restoration. UCLA CAM Rep., (1):1–29, 2008.