## Learned backprojection for sparse and limited view photoacoustic tomography

Johannes Schwab, Stephan Antholzer, Markus Haltmeier

Department of Mathematics, University of Innsbruck Technikerstraße 13, 6020 Innsbruck, Austria, johannes.schwab@uibk.ac.at

## Abstract

Filtered backprojection (FBP) is an efficient and popular class of tomographic image reconstruction methods. In photoacoustic tomography, these algorithms are based on theoretically exact analytic inversion formulas which results in accurate reconstructions. However, photoacoustic measurement data are often incomplete (limited detection view and sparse sampling), which results in artefacts in the images reconstructed with FBP. In addition to that, properties such as directivity of the acoustic detectors are not accounted for in standard FBP, which affects the reconstruction quality, too. To account for these issues, in this papers we propose to improve FBP algorithms based on machine learning techniques. In the proposed method, we include additional weight factors in the FBP, that are optimized on a set of incomplete data and the corresponding ground truth photoacoustic source. Numerical tests show that the learned FBP improves the reconstruction quality compared to the standard FBP.

Keywords: Photoacoustic tomography, image reconstruction, sparse data, limited view, machine learning, detector directivity, learned backprojection

## 1 Introduction

Photoacoustic tomography (PAT) is a promising imaging method for medical diagnosis based on the photoacoustic effect. The process of data acquisition for PAT can be described as follows. Illumination of the sample with a short laser pulse leads to heating of the tissue, and along with that to a

thermoelastic expansion. This expansion generates a pressure wave, which is measured outside of the tissue. In the case of full data, several efficient reconstruction methods have been developed. In particular, this includes filtered backprojection algorithms bases on analytic inversion formulas known for special geometries[19, 9, 10, 11, 16, 25].

In practice it is often not possible to collect measurement data of the resulting pressure wave on a measurement surface surrounding the whole object, and the data is known only on parts of such a surrounding surface. This problem of image reconstruction from incomplete data of this type, is called limited view problem. Furthermore, because for each spatial measurement an individual, often costly, detector is needed, the number of detectors placed on the detection area is limited, leading to spacially sparsly sampled data. In particular, in the reconstruction of the initial pressure distribution, the application of existing FBP algorithms to the limited view problem as well as the sparse data problem leads to specific artifacts. It is known, that the reconstruction of singularities not visible from the detector positions is unstable, whereas singularities contained within the convex hull can theoretically be stably recovered in appropriate spaces [15].

For PAT from limited view data several iterative reconstruction methods have been proposed to reduce limited view artifacts and improve the reconstruction quality [22, 7, 13]. In [21] the authors proposed using weight factors depending on the angle between the reconstruction- and detection point for an inversion formula, which is exact on continuously sampled data given on the whole boundary. This method tries to compensate the zero-extension of the data on parts of the boundary, where the data is not available. In [18] an iterative algorithm, updating these angle-depending weight factors has been proposed. Deep learning methods for artifact reduction in limited view PAT have been proposed in [24] and a framework to learn an extension of the data was introduced in [8]. In this paper we propose to learn weight factors from a big dataset of incomplete data with simulated directional detector sensitivity and the corresponding ground truth reconstructions. We compute appropriate weights via a gradient descent method in the framework of neural networks and the available software and implementations on GPUs.

## 2 Methods

## 2.1 Photoacoustic tomography (PAT)

In PAT the investigated sample is illuminated with a short laser pulse, which induces heating in the tissue which in turn causes a thermoelastic expansion.

This expansion generates a pressure wave, which satisfies the initial value problem

$$\partial _ { t } ^ { 2 } u ( x, t ) - c ^ { 2 } \Delta _ { x } u ( x, t ) & = 0, \\ u ( x, 0 ) & = f ( x ), \\ \partial _ { t } u ( x, 0 ) & = 0.$$

Here for d ∈ { 2 3 , } the spatial location is denoted by x ∈ R d , t ∈ R denotes the time, ∆ x the spatial Laplacian and c the speed of sound. We further assume, that the photoacoustic (PA) source satisfies supp f ⊂ Ω for a bounded domain Ω ⊂ R d .

Spatial dimension d = 3 is relevant for PAT, when point-wise measurements are taken on the detection surface ∂ Ω . In this paper we consider two spatial dimension, which is relevant for PAT using integrating line detectors [3, 4, 23, 20]. In the case of integrating line detectors, the PA source f corresponds to a projection image of the 3D source. By reconstructing projection images of sufficiently many directions, the 3D source can be reconstructed using the inverse Radon transform.

In two spatial dimensions the solution of (1) denoted by u f for given initial pressure f is given by the well known solution formula

$$u _ { f } ( x, t ) \coloneqq \frac { 1 } { 2 \pi } \frac { \partial } { \partial t } \int _ { 0 } ^ { t } \int _ { \partial B _ { 1 } ( 0 ) } \frac { r f ( x + r \omega ) } { \sqrt { t ^ { 2 } - r ^ { 2 } } } \text{d} \omega \text{ d} r,$$

where ∂B 1 (0) denotes the unit sphere. From data as in (2), the PA source can be reconstructed by exact inversion formulas. In this paper, we consider the following three issues affecting PAT image reconstruction. We denote points on detection curve by s ∈ ∂ Ω .

## glyph[squaresolid] Limited view PAT

For a measurement surface ∂ Ω containing the support of the source f , the problem of recovering f from measurements on parts of ∂ Ω is called limited view PAT. For limited view PAT, the given data is of the form

$$u _ { f } ( s _ { i }, t _ { j } ) \quad i = 1, \dots, N _ { s } ; \, j = 1, \dots N _ { t } \text{ with } s _ { i } \in \Gamma, \quad ( 3 )$$

where the detector locations for i = 1 , . . . N s are contained in a subset Γ ⊂ ∂ Ω , which is not surrounding the source f . (see Figure 2 A) )

## glyph[squaresolid] Sparse data PAT

The number of spatial measurements of the solution of (1) on the measurement curve ∂ Ω is limited, which is called sparse data problem. The data for the reconstruction of a acoustic source f is therefore given by

$$\{ u _ { f } ( s _ { i }, t _ { j } ) \ | \ i = 1, \dots, N _ { s } ; j = 1, \dots N _ { t } \} \text{ with } s _ { i } \in \Gamma, \quad ( 4 )$$

where the number of detectors N s is small, whereas the number of measurements in time N t is typically large enough. (see Figure 2 B) )

## glyph[squaresolid] Detector directivity

To model the directional sensitivity of the acoustic detectors we consider a weight function depending on the angle between the detection points s and source locations. In particular we assume, that the wavedata is given by [26, 27]

$$u _ { f } ( s, t ) \coloneqq \frac { 1 } { 2 \pi } \frac { \partial } { \partial t } \int _ { 0 } ^ { t } \int _ { \partial B _ { 1 } ( 0 ) } \frac { r f ( s + r \omega ) } { \sqrt { t ^ { 2 } - r ^ { 2 } } } \varphi ( \nu _ { s }, \omega ) \, \mathrm d \omega \, \mathrm d r, \quad ( 5 )$$

where ϕ : R 2 × R 2 → [0 , 1] is a function depending on the angle between source location and the outer normal on the detection surface at s . In our simulations with a circle as detection curve, we modelled the sensitivity by choosing the function ϕ as

$$\varphi ( s, \omega ) \coloneqq \begin{cases} \frac { \langle s - \omega, \omega \rangle ^ { 2 } } { \| s - \omega \| ^ { 2 } } = \cos ( \alpha ) ^ { 2 }, & \quad | \alpha | < \frac { \pi } { 2 }, \\ 0, & \quad | \alpha | \geq \frac { \pi } { 2 }, \end{cases}$$

where α denotes the angle between source and detection point.

## 2.2 Weighted universal backprojection (UBP)

In the case, where the data of a source f is given as a continuous function g : ∂ Ω × [0 , ∞ → ) R , defined on the whole boundary, the inverse problem of determining the source f is uniquely and stably solvable. Many efficient and robust reconstruction methods for recovering the PA source including filtered backprojection, Fourier methods or time reversal have been developed []. In particular an exact inversion formula for continuous sampled data on certain detection surfaces ∂ Ω , including circular geometries, is the universal backprojection (UBP) fomula. In 2D the UBP has first been derived in [4], and [25], for the 3D case.

Figure 1: Illustration of the directivity model used for our simulations.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00593v1-figure-a53e3e7a-3549-4972-bc92-77c0df713032.jpg)

For wave data g , given on a subset of the boundary Γ ⊂ Ω , the 2D universal backprojection (UBP) including a weight function w : R 2 × ∂ Ω → R , is defined by

$$\Psi ( w, g ) ( x ) = - \frac { 2 } { \pi } \int _ { \Gamma } w ( x, s ) ^ { 2 } \langle \nu _ { s }, x - s \rangle \int _ { \| x - s \| } ^ { \infty } \frac { \partial _ { t } ( t ^ { - 1 } g ( s, t ) ) } { \sqrt { t ^ { 2 } - \| x - s \| ^ { 2 } } } d t \, \mathrm d s. \quad ( 7 )$$

Here x denotes the reconstruction point, ν s the exterior normal vector to a point s ∈ Γ ⊂ ∂ Ω on the detection surface and 〈 · , · 〉 the inner product in R 2 . The weight function w in the backprojection depends on the the reconstruction point and the detector point. In the standard UBP the weight function is constant, w ≡ 1 .

In practice, for PA image reconstruction, the measurements of the pressure are only given for a finite number of detector locations, often not surrounding the whole object, and for a finite number of time samples. We denote by G ∈ R N t × N s the data for time samples { t k = kT/N t | k = 1 , . . . , N t } ⊂ [0 , T ] , with final measurement time T ∈ R , and detector locations { s k | k = 1 , . . . , N s } ⊂ Γ on the measurement curve. The image reconstruction problem of PAT can now be modelled as estimating the discretized PA source F ∈ R N x × N x corresponding to the initial pressure f in (1), from discrete data

$$G = A ( F ).$$

Here A : R N x × N x → R N t × N s the discretized operator mapping a source to the corresponding discrete wave data. Discetizing the weighted UBP formula

Ψ leads to a mapping P ( W , · ) : R N t × N s → R N x × N x with a weight vector W ∈ R N x × N x × N s . The case W ≡ 1 corresponds to the discretization of the ordinary UBP.

## 2.3 Learning the weights in the UBP

Deep learning and in particular deep convolutional neural networks (CNNs) have been very successfully applied to a great variety of pattern recognition and image processing tasks. Recently a lot of research has been done in solving inverse problems, incorporating deep learning techniques, including efficient and accurate image reconstruction methods in tomographic problems [2, 5, 12, 13, 14, 28, 17].

In the context of a FBP in limited data, the task of learning is to find optimal weights in the discretized UBP formula, with respect to a data set ( G , F i i ) M i =1 ⊂ R N s × N t × R N x × N x , where G , i i = 1 , . . . , M denote the incomplete wave data of the sources F i . The high dimensional weight vector of the UBP formula should minimize some error function

$$E ( \mathbb { W }, ( G _ { i }, F _ { i } ) _ { i = 1 } ^ { M } ) = \frac { 1 } { M } \sum _ { i = 1 } ^ { M } \mathcal { D } ( F _ { i }, P ( \mathbb { W }, G _ { i } ) ),$$

where D : R N x × N x × R N x × N x → [0 , ∞ ) is some distance measure. The implementation of the network was done in keras, running on top of tensorflow, where a first layer of the network does the filtering in time without any trainable weights and a second layer implements the backprojection with the adjustable weights W .

## 3 Results and discussion

In our numerical studies we consider the case where Ω = B 1 is the ball around the origin with radius one. We learn weights in the UBP for three different limited data scenarios, shown in Figure 2. To efficiently optimize the weights we use the deep learning framework keras [6] which is based on tensorflow [1]. Using this efficient implementation the weighted UBP formula can be combined with CNNs to further improve reconstructions.

For learning the weights in all three cases, we numerically compute direction dependent wave data (5) for M = 1800 data-source pairs, with discretized sources F i ∈ R 256 × 256 . These generated initial sources are variations of a Shepp-Logan phantom. In particular the positions and amplitudes of the

Figure 2: Considered measurement scenarios. The detectors { s k } N s k =1 ⊂ Γ are represented by dots and ∂ Ω by the solid line.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00593v1-figure-7a2fa179-44c1-4094-a9f1-1ef723009ec8.jpg)

ellipses in the Shepp-Logan are randomly selected. Additionally, the phantoms were rotated by a randomly chosen angle and some finer structures were added. Finally we performed elastic deformations on the phantoms generated this way. The time discretization for all three scenarios is given by 400 equidistant samples in the interval [0 , 3] . For the error function we choose the distance measure induced by ‖·‖ 2 , which gives the mean squared error

$$E ( \mathbb { W }, ( G _ { i }, F _ { i } ) _ { i = 1 } ^ { M } ) = \frac { 1 } { M } \sum _ { i = 1 } ^ { M } \| F _ { i } - P ( \mathbb { W }, G _ { i } ) \| _ { 2 }.$$

The network incorporates about two million free parameters and the minimization is done by the stochastic gradient descent algorithm with batch size one over 100 epochs.

## 3.1 Limited view

For the limited view case we assumed the detectors to lie on the half circle as shown in Figure 2 A) . The simulated photoacoustic data was computed on 100 detector locations. The results for a phantom randomly selected according to the model described above but not contained in the training set is shown in Figure 3.

As can be seen in Figure 3 the reconstructions with the weighted UBP leads slightly better results than the reconstruction with the unmodified

Figure 3: Reconstruction results for the limited view case in the reconstruction domain [ -1 1] , 2 . In the top row a simulated phantom not contained in the training set (left), the reconstruction using the UBP without trained weights (center) and the reconstruction using the weighted UBP (right) are shown. The bottom row displays the learned weights for measurements on a half circle (left) and the absolute differences from the plain UBP (center) and weighted UBP (right) to the ground truth.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00593v1-figure-b1939912-72b3-44e8-8ea9-87eedb87e0ea.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00593v1-figure-46c84bee-4cca-4598-9c45-a9a4d7f3940a.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00593v1-figure-0fbe3612-87ae-4176-a7e5-390abb3bc61b.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00593v1-figure-1b284c26-b33a-4e06-859e-244da5671e9e.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00593v1-figure-d6f76355-9949-415e-8373-11bb0b09433b.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00593v1-figure-861aca50-4ca4-4e04-8c00-9194c8d7f38a.jpg)

UBP. We also observe, that the learned weights follow a certain structure depending on the angle between the reconstruction and detection point.

## 3.2 Sparse sampling

In the sparse data scenario the measurement geometry was taken to be a circle with N s = 20 equidistant measurement positions, see Figure 2 B) . Reconstruction results in the sparse data case are presented in Figure 4.

Again, the reconstructions using the learned UBP are visually superior to the results obtained with the unweighted UBP formula.

## 3.3 Limited view and sparse sampling

In the limited view and sparse sampling case the detector locations were chosen to lie on a half circle as shown in 2 C) . The wave data was computed on N s = 20 detector positions. The reconstruction results for these case are

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00593v1-figure-6d2340d1-fb72-4423-97a0-ed3c137912ef.jpg)

-1

Figure 4: Reconstruction results for sparse data on the reconstruction domain [ -1 1] , 2 . In the top row a simulated phantom not contained in the training set (left), the reconstruction using the UBP without trained weights (center) and the reconstruction using the weighted UBP (right) are shown. The third row displays the learned weights for sparse measurements and the absolute differences from the plain UBP (center) and weighted UBP (right) to the ground truth.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00593v1-figure-2dcaa7ab-e771-4d00-a631-962e060c40d7.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00593v1-figure-87a770ac-9d5b-4c74-88f3-ceb5115af0ea.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00593v1-figure-c74b29a4-3aea-4d37-aa2c-4cac88f3a1f4.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00593v1-figure-04ca5045-1c1f-4f3b-8e35-f72346aa29c1.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00593v1-figure-56e86653-2e72-4d30-bd3e-546cb8b7fcf2.jpg)

displayed in Figure 5 and are once more superior to the reconstructions with the ordinary UBP.

## 3.4 Quantitative evaluation

To evaluate of the learned universal backprojection formula not just visually the relative squared glyph[lscript] 2 -error averaged over 200 reconstructions ( ˆ ) F i 200 i =1 of phantoms ( F i ) 200 i =1 not contained in the training set

$$\text{error} \left ( ( F _ { i } ) _ { i = 1 } ^ { 2 0 0 }, ( \hat { F } _ { i } ) _ { i = 1 } ^ { 2 0 0 } \right ) \coloneqq \frac { 1 } { 2 0 0 } \sum _ { i = 1 } ^ { 2 0 0 } \frac { \| \hat { F } _ { i } - F _ { i } \| _ { 2 } } { \| F _ { i } \| _ { 2 } },$$

was computed. The averaged reconstruction errors for all scenarios A) B) , , C) are displayed in the following table.

Table 1 shows that the reconstruction error is roughly halved by the use of adjusted weights in the UBP formula.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00593v1-figure-0126885e-cabc-49e5-97e2-364d62f0dcd8.jpg)

5

10

15

20

Figure 5: Reconstruction results for limited view and sparse data on the reconstruction domain [ -1 1] , 2 . In the top row a simulated phantom not contained in the training set, the reconstruction using the UBP without trained weights and the reconstruction using the weighted UBP are shown. The third row shows the learned weights for sparse measurements on the half circle and the absolute differences to the ground truth.

## 4 Conclusion

In this paper, we proposed learned weight factors to improve the reconstruction quality of filtered backprojection algorithms in PAT accounting for limited view data, sparse sampling and detector directivity. The optimization procedure of the weight factors is very flexible and can be done for arbitrary surfaces, including geometries, where no exact inversion formula is known. One only needs an appropriate training set consisting of PA sources and corresponding PA data. The presented results demonstrate, that for data given on a half circle, the learned FBP clearly improves the results compared the standard FBP.

We note that the learnable backprojection could also be used as first layer in a deep CNN as proposed in [24]. Further, it is possible to improve the reconstruction quality by learning temporal filters in the UBP, or to use correction weights to reduce the error on arbitrary convex and bounded domains.

|              |     A) |     B) |     C) |
|--------------|--------|--------|--------|
| UBP          | 0.2002 | 0.3461 | 0.3546 |
| weighted UBP | 0.0912 | 0.1806 | 0.1649 |

Table 1: Reconstruction errors for the different scenarios A) limited view, B) sparse sampling and C) sparse sampling and limited view.

## Acknowledgement

The research of S.A. and M.H. has been supported by the Austrian Science Fund (FWF), project P 30747-N32; the work of R.N. has been supported by the FWF, project P 28032. We acknowledge the support of NVIDIA Corporation with the donation of the Titan Xp GPU used for this research.

## References

- [1] Martín Abadi, Paul Barham, Jianmin Chen, Zhifeng Chen, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Geoffrey Irving, Michael Isard, et al. Tensorflow: a system for large-scale machine learning. In OSDI , volume 16, pages 265-283, 2016.
- [2] Stephan Antholzer, Markus Haltmeier, and Johannes Schwab. Deep learning for photoacoustic tomography from sparse data. Inverse Problems in Science and Engineering , pages 1-19, 2018.
- [3] Johannes Bauer-Marschallinger, Karoline Felbermayer, and Thomas Berer. All-optical photoacoustic projection imaging. Biomedical optics express , 8(9):3938-3951, 2017.
- [4] Peter Burgholzer, J Bauer-Marschallinger, H Grün, Markus Haltmeier, and G Paltauf. Temporal back-projection algorithms for photoacoustic tomography with integrating line detectors. Inverse Problems , 23(6):S65, 2007.
- [5] Hu Chen, Yi Zhang, Weihua Zhang, Peixi Liao, Ke Li, Jiliu Zhou, and Ge Wang. Low-dose ct via convolutional neural network. Biomedical optics express , 8(2):679-694, 2017.
- [6] François Chollet et al. Keras. https://github.com/fchollet/keras , 2015.

| [7]   | X Luís Dean-Ben, Andreas Buehler, Vasilis Ntziachristos, and Daniel Razansky. Accurate model-based reconstruction algorithm for three- dimensional optoacoustic tomography. IEEE Transactions on Medical Imaging , 31(10):1922-1928, 2012.                                       |
|-------|----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| [8]   | Florian Dreier, Sergiy Pereverzyev Jr, and Markus Haltmeier. Oper- ator learning approach for the limited view problem in photoacoustic tomography. Computational Methods in Applied Mathematics , 2017.                                                                         |
| [9]   | David Finch, Markus Haltmeier, and Rakesh. Inversion of spherical means and the wave equation in even dimensions. SIAM Journal on Applied Mathematics , 68(2):392-412, 2007.                                                                                                     |
| [10]  | David Finch and Sarah K Patch. Determining a function from its mean values over a family of spheres. SIAM journal on mathematical analysis , 35(5):1213-1240, 2004.                                                                                                              |
| [11]  | Markus Haltmeier. Universal inversion formulas for recovering a func- tion from spherical means. SIAM Journal on Mathematical Analysis , 46(1):214-232, 2014.                                                                                                                    |
| [12]  | Yoseob Han, Jawook Gu, and Jong Chul Ye. Deep learning inte- rior tomography for region-of-interest reconstruction. arXiv preprint arXiv:1712.10248 , 2017.                                                                                                                      |
| [13]  | Andreas Hauptmann, Felix Lucka, Marta Betcke, Nam Huynh, Jonas Adler, Ben Cox, Paul Beard, Sebastien Ourselin, and Simon Arridge. Model-based learning for accelerated, limited-view 3-d photoacoustic tomography. IEEE transactions on medical imaging , 37(6):1382-1393, 2018. |
| [14]  | Brendan Kelly, Thomas P Matthews, and Mark A Anastasio. Deep learning-guided image reconstruction from incomplete data. arXiv preprint arXiv:1709.00584 , 2017.                                                                                                                  |
| [15]  | Peter Kuchment and Leonid Kunyansky. Mathematics of photoacoustic and thermoacoustic tomography. In Handbook of Mathematical Methods in Imaging , pages 817-865. Springer, 2011.                                                                                                 |
| [16]  | Leonid A Kunyansky. Explicit inversion formulae for the spherical mean radon transform. Inverse problems , 23(1):373, 2007.                                                                                                                                                      |

| [17]   | Housen Li, Johannes Schwab, Stephan Antholzer, and Markus Halt- meier. Nett: Solving inverse problems with deep neural networks. arXiv preprint arXiv:1803.00092 , 2018.                                                                     |
|--------|----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| [18]   | Xueyan Liu, Dong Peng, Xibo Ma, Wei Guo, Zhenyu Liu, Dong Han, Xin Yang, and Jie Tian. Limited-view photoacoustic imaging based on an iterative adaptive weighted filtered backprojection approach. Applied optics , 52(15):3477-3483, 2013. |
| [19]   | Linh V Nguyen. A family of inversion formulas in thermoacoustic to- mography. arXiv preprint arXiv:0902.2579 , 2009.                                                                                                                         |
| [20]   | Guenther Paltauf, Petra Hartmair, Georgi Kovachev, and Robert Nuster. Piezoelectric line detector array for photoacoustic tomography. Photoacoustics , 8:28-36, 2017.                                                                        |
| [21]   | Guenther Paltauf, Robert Nuster, and Peter Burgholzer. Weight factors for limited angle photoacoustic tomography. Physics in Medicine & Biology , 54(11):3303, 2009.                                                                         |
| [22]   | Guenther Paltauf, Robert Nuster, Markus Haltmeier, and Peter Burgholzer. Experimental evaluation of reconstruction algorithms for limited view photoacoustic tomography with line detectors. Inverse Problems , 23(6):S81, 2007.             |
| [23]   | Guenther Paltauf, Robert Nuster, Markus Haltmeier, and Peter Burgholzer. Photoacoustic tomography using a mach-zehnder interfer- ometer as an acoustic line detector. Applied optics , 46(16):3352-3358, 2007.                               |
| [24]   | Johannes Schwab, Stephan Antholzer, Robert Nuster, and Markus Halt- meier. Real-time photoacoustic projection imaging using deep learning. arXiv preprint arXiv:1801.06693 , 2018.                                                           |
| [25]   | Minghua Xu and Lihong V Wang. Universal back-projection algo- rithm for photoacoustic computed tomography. Physical Review E , 71(1):016706, 2005.                                                                                           |
| [26]   | Yuan Xu and Lihong V Wang. Time reversal and its application to to- mography with diffracting sources. Physical review letters , 92(3):033902,                                                                                               |

- [27] Gerhard Zangerl, Sunghwan Moon, and Markus Haltmeier. Photoacoustic tomography with direction dependent data: An exact series reconstruction approach. Inverse Problems , 2019.
- [28] Hanming Zhang, Liang Li, Kai Qiao, Linyuan Wang, Bin Yan, Lei Li, and Guoen Hu. Image prediction for limited-angle tomography via deep learning with convolutional neural network. arXiv preprint arXiv:1607.08707 , 2016.
Johannes Schwab Stephan Antholzer Markus Haltmeier
Learned backprojection for sparse and limited view photoacoustic tomography
Filtered backprojection (FBP) is an efficient and popular class of
tomographic image reconstruction methods. In photoacoustic tomography, these
algorithms are based on theoretically exact analytic inversion formulas which
results in accurate reconstructions. However, photoacoustic measurement data
are often incomplete (limited detection view and sparse sampling), which
results in artefacts in the images reconstructed with FBP. In addition to that,
properties such as directivity of the acoustic detectors are not accounted for
in standard FBP, which affects the reconstruction quality, too. To account for
these issues, in this papers we propose to improve FBP algorithms based on
machine learning techniques. In the proposed method, we include additional
weight factors in the FBP, that are optimized on a set of incomplete data and
the corresponding ground truth photoacoustic source. Numerical tests show that
the learned FBP improves the reconstruction quality compared to the standard
FBP.