Fast reconstruction approaches for photoacoustic tomography with smoothing Sobolev/Matérn priors

Jaakko Kultima; Ronny Ramlau; Teemu Sahlström; Tanja Tarvainen

doi:10.3934/ipi.2026015

Article Contents

2026, Volume 25: 51-66. Doi: 10.3934/ipi.2026015

This volume Previous Article A direct sampling method for inverse elastic scattering in cavity structure Next Article Factorization method for near-field inverse scattering problems in elastodynamics

Fast reconstruction approaches for photoacoustic tomography with smoothing Sobolev/Matérn priors

1.
Johann Radon Institute for Computational and Applied Mathematics (RICAM), Linz, Altenbergerstraße 69, A-4040 Linz, Austria
2.
Industrial Mathematics Institute, Johannes Kepler University, Linz, Altenbergerstraße 69, A-4040 Linz, Austria
3.
Department of Technical Physics, University of Eastern Finland, P.O. Box 1627, 70211 Kuopio, Finland

^*Corresponding author: Jaakko Kultima
^*Corresponding author: Jaakko Kultima

Received: July 3, 2025

Revised: December 19, 2025

Early access: January 26, 2026

Published online: January 26, 2026

Abstract / Introduction Full Text(HTML) Figure(6) / Table(2) Related Papers Cited by

Abstract

In photoacoustic tomography (PAT), the computation of the initial pressure distribution within an object from its time-dependent boundary measurements over time is considered. This problem can be approached from two well-established points of view: deterministically using regularization methods, or stochastically using the Bayesian framework. Both approaches frequently require the solution of a variational problem. In this paper, we elaborate the connection between these approaches by establishing the equivalence between a smoothing Matérn class of covariance operators and Sobolev embedding operator $ E_s: H^s \hookrightarrow L^2 $. We further discuss the use of a Wavelet-based implementation of the adjoint operator $ E_s^* $, which also allows for efficient evaluations for certain Matérn covariance operators, leading to efficient implementations both in terms of computational effort as well as memory requirements. The proposed methods are validated with reconstructions for the photoacoustic problem.

Keywords:

Mathematics Subject Classification: Primary: 92C55, 65M12; Secondary: 35R30, 62P10.

Citation:

Full Text(HTML)

Figure 1. The condition numbers of the acoustic forward operator $ K $ vs. the number of sensors at the boundary as described above

Download: Full-size image PowerPoint slide

Figure 2. Simulated initial pressure distribution $ p_0 $, consisting of a vein-like pattern within a square-shaped domain with non-constant background. The phantom was modified from the high-resolution Fundus (HRF) image database published in [7]. Sensors are evenly distributed along one and two sides of the target. Cross-sections along the diagonal white line of the reconstructions are plotted in Figure 6

Download: Full-size image PowerPoint slide

Figure 3. Estimated $ p_0 $ in the two-sided sensor geometry. Top row, reference reconstructions using total variation, $ \ell^1 $ for wavelet coefficients, and standard $ L^2 $- Tikhonov on the top row. Bottom row, reconstructions with $ s = 1 $ (Matérn covariance $ \Gamma_{0, 0} $ in $ \mathbb{{R}}^2 $), $ s = 3/2 $ (corresponding to Matérn covariance $ \Gamma_{1/2, 1} $ in $ \mathbb{{R}}^2 $), and $ s = 2 $ (Matérn covariance $ \Gamma_{1, \sqrt{2}} $ in $ \mathbb{{R}}^2 $)

Download: Full-size image PowerPoint slide

Figure 4. Estimated $ p_0 $ in the one-sided sensor geometry. Top row, reference reconstructions using total variation, $ \ell^1 $ for wavelet coefficients, and standard $ L^2 $- Tikhonov on the top row. Bottom row, reconstructions with $ s = 1 $ (Matérn covariance $ \Gamma_{0, 0} $ in $ \mathbb{{R}}^2 $), $ s = 3/2 $ (corresponding to Matérn covariance $ \Gamma_{1/2, 1} $ in $ \mathbb{{R}}^2 $), and $ s = 2 $ (Matérn covariance $ \Gamma_{1, \sqrt{2}} $ in $ \mathbb{{R}}^2 $)

Download: Full-size image PowerPoint slide

Figure 5. Estimated $ p_0 $ in the two-sided sensor geometry, with Sobolev smoothness $ s = 1 $, and noise-levels 1%, 5%, and 10%

Download: Full-size image PowerPoint slide

Figure 6. Diagonal cross sections of the reconstructions in two-sided sensor geometry. The orange plot is the underlying ground truth, and blue line is the reconstructed value

Download: Full-size image PowerPoint slide

Table 1. The computational efficiency as measured by the number of evaluations of the wave propagation using pseudo spectral method. Computational cost is dominated by these evaluations; convergence of the chosen method is critical. Algorithmic runtimes and memory usage are comparable for all examples

	$ \ell^1 $ (FISTA)	TV (Barzilai-Borwein)	Tikhonov
Number of iterations	200	150	25
Number of simulations	401	301	51

| Show Table

DownLoad: CSV

Table 2. Reconstruction quality at 5% noise level measured in terms of relative error (RE), peak signal-to-noise ratio (PSNR), and structural similarity index measure (SSIM) for TV, $ \ell^1 $, and $ H^s $-Tikhonov with $ s = 0, 1, 3/2, 2 $

		TV	$ \ell^1 $	$ s = 0 $	$ s = 1 $	$ s = 3/2 $	$ s = 2 $
	RE(%)	48.70	39.88	40.12	38.95	40.31	44.96
One side	PSNR	13.38	15.11	15.06	15.32	15.02	14.07
	SSIM	0.57	0.61	0.58	0.62	0.64	0.58
	RE(%)	19.31	15.00	12.16	12.42	17.32	33.93
Two side	PSNR	21.41	23.61	25.43	25.25	22.35	16.51
	SSIM	0.74	0.72	0.70	0.71	0.72	0.64

| Show Table

DownLoad: CSV

Related Papers

Cited by

References

[1]	M. Agranovsky and P. Kuchment, Uniqueness of reconstruction and an inversion procedure for thermoacoustic and photoacoustic tomography with variable sound speed, Inv Probl, 23 (2007), 2089-2102. doi: 10.1088/0266-5611/23/5/016.
[2]	N Aronszajn and K. T. Smith, Theory of bessel potentials, Ⅰ, Annales de l'institut Fourier, 11 (1961), 385-475. doi: 10.5802/aif.116.
[3]	S. Arridge, P. Beard, 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 (2016), 8908-8940. doi: 10.1088/1361-6560/61/24/8908.
[4]	P. Beard, Biomedical photoacoustic imaging, Interface Focus, 1 (2011), 602-631. doi: 10.1098/rsfs.2011.0028.
[5]	M. Bergounioux, X. Bonnefond, T. Haberkorn and Y. Privat, An optimal control problem in photoacoustic tomography, Mathematical Models and Methods in Applied Sciences, 24 (2014), 2525-2548. doi: 10.1142/S0218202514500286.
[6]	Y. E. Boink, M. J. Lagerwerf, W. Steenbergen, S. A. van Gils, S. Manohar and C. Brune, A framework for directional and higher-order reconstruction in photoacoustic tomography, Physics in Medicine & Biology, 63 (2018), 045018.
[7]	A. Budai, R. Bock, A. Maier, J. Hornegger and G. Michelson, Robust vessel segmentation in fundus images, International Journal of Biomedical Imaging, 2013 (2013), 154860. doi: 10.1155/2013/154860.
[8]	A. Buehler, A. Rosenthal, T. Jetzfellner, A. Dima, D. Razansky and V. Ntziachristos, Model-based optoacoustic inversions with incomplete projection data, Medical Physics, 38 (2011), 1694-1704. doi: 10.1118/1.3556916.
[9]	B. Cox, J. G. Laufer, S. R. Arridge and P. C. Beard, Quantitative spectroscopic photoacoustic imgaging: A review, J Biomed Opt, 17 (2012), 061202.
[10]	B. T. Cox and B. E. Treeby, Artifact trapping during time reversal photoacoustic imaging for acoustically heteregeneous media, IEEE Trans Med Imag, 29 (2010), 387-396. doi: 10.1109/TMI.2009.2032358.
[11]	H. W. Engl and R. Ramlau, Regularization of Inverse Problems, In B. Engquist, editor, Encyclopedia of Applied and Computational Mathematics. Springer, 2015.
[12]	H. W. Engl, M. Hanke and A. Neubauer, Regularization of Inverse Problems, Math. Appl., 375, Kluwer Academic Publishers, Dordrecht, Boston, London, 2000.
[13]	D. Finch, M. Haltmeier and Ra kesh, Inversion of spherical means and the wave equation in even dimensions, SIAM Journal on Applied Mathematics, 68 (2007), 392-412. doi: 10.1137/070682137.
[14]	D. Finch, S. K. Patch and Ra kesh, Determining a function from its mean values over a family of spheres, SIAM J Math Anal, 35 (2004), 1213-1240. doi: 10.1137/S0036141002417814.
[15]	J. Frikel and M. Haltmeier, Efficient regularization with wavelet sparsity constraints in photoacoustic tomography, Inverse Problems, 34 (2018), 024006, 28pp. doi: 10.1088/1361-6420/aaa0ac.
[16]	M. Haltmeier, T. Berer, S. Moon and P. Burgholzer, Compressed sensing and sparsity in photoacoustic tomography, Journal of Optics, 18 (2016), 114004. doi: 10.1088/2040-8978/18/11/114004.
[17]	M. Haltmeier, T. Schuster and O. Scherzer, Filtered backprojection for thermoacoustic computed tomography in spherical geometry, Mathematical Methods in the Applied Sciences, 28 (2005), 1919-1937. doi: 10.1002/mma.648.
[18]	A. Hauptmann and B. Cox, Deep learning in photoacoustic tomography: Current approaches and future directions, Journal of Biomedical Optics, 25 (2020), 112903-112903.
[19]	Y. Hristova, P. Kuchment and L. Nguyen, Reconstruction and time reversal in thermoacoustic tomography in acoustically homogeneous and inhomogeneous media, Inv Probl, 24 (2008), 055006, 25pp. doi: 10.1088/0266-5611/24/5/055006.
[20]	S. Hubmer, E. Sherina and R. Ramlau, Characterizations of adjoint sobolev embedding operators with applications in inverse problems, Electronic Transactions on Numerical Analysis, 59 (2023), 116-144. doi: 10.1553/etna_vol59s116.
[21]	A. Javaherian and S. Holman, A multi-grid iterative method for photoacoustic tomography, IEEE Trans Med Imag, 36 (2017), 696-706. doi: 10.1109/TMI.2016.2625272.
[22]	J. Kaipio and E. Somersalo, Statistical and Computational Inverse Problems, Appl. Math. Sci., 160, Springer-Verlag, New York, 2005.
[23]	P. Kuchment and L. Kunyansky, Mathematics of photoacoustic and thermoacoustic tomography, In Handbook of Mathematical Methods in Imaging, Springer New York, 19 (2011), 817-865. doi: 10.1007/978-0-387-92920-0_19.
[24]	L. A. Kunyansky, Explicit inversion formulae for the spherical mean Radon transform, Inv Probl, 23 (2007), 373-383. doi: 10.1088/0266-5611/23/1/021.
[25]	L. A. Kunyansky, A series solution and a fast algorithm for the inversion of the spherical mean Radon transform, Inv Probl, 23 (2007), S11-S20. doi: 10.1088/0266-5611/23/6/S02.
[26]	C. Li and L. V. Wang, Photoacoustic tomography and sensing in biomedicine, Phys Med Biol, 54 (2009), R59-R97. doi: 10.1088/0031-9155/54/19/R01.
[27]	A. K. Louis, Inverse und Schlecht Gestellte Probleme, Teubner Studienbücher Mathematik, B. G. Teubner, Stuttgart, 1989.
[28]	Y. Meyer, Wavelets and Operators, Cambridge Studies in Advanced Mathematics, 37, Cambridge University Press, Cambridge, 1992.
[29]	W. Naetar and O. Scherzer, Quantitative photoacoustic tomography with piecewise constant material parameters, SIAM J Imag Sci, 7 (2014), 1755-1774. doi: 10.1137/140959705.
[30]	J. Poudel, Y. Lou and M. A. Anastasio, A survey of computational frameworks for solving the acoustic inverse problem in three-dimensional photoacoustic computed tomography, Physics in Medicine & Biology, 64 (2019), 14TR01.
[31]	J. Qian, P. Stefanov, G. Uhlmann and H. Zhao, An efficient neumann series–based algorithm for thermoacoustic and photoacoustic tomography with variable sound speed, SIAM Journal on Imaging Sciences, 4 (2011), 850-883. doi: 10.1137/100817280.
[32]	R. Ramlau and G. Teschke, Regularization of Sobolev embedding operators and applications to medical imaging and meteorological data. Part Ⅰ: Regularization of sobolev embedding operators, Sampling Theory in Signal and Image Processing, 3 (2004), 175-195.
[33]	R. Ramlau and G. Teschke, Regularization of Sobolev embedding operators and applications to medical imaging and meteorological data. Part Ⅱ: Regularization incorporating noise with applications in medical imaging and meteorological data, Sampling Theory in Signal and Image Processing, 3 (2004), 205-226.
[34]	C. E. Rasmussen and C. K. I. Williams, Gaussian Processes for Machine Learning, Adapt. Comput. Mach. Learn, MIT Press, Cambridge, MA, 2006.
[35]	T. Sahlström, A. Pulkkinen, J. Tick, J. Leskinen and T. Tarvainen, Modeling of errors due to uncertainties in ultrasound sensor locations in photoacoustic tomography, IEEE Transactions on Medical Imaging, 39 (2020), 2140-2150. doi: 10.1109/TMI.2020.2966297.
[36]	O. Scherzer, M. Grasmair, H. Grossauer, M. Haltmeier and F. Lenzen, Variational Methods in Imaging, Appl. Math. Sci., 167, Springer, New York, 2009.
[37]	P. Stefanov and G. Uhlmann, Thermoacoustic tomography with variable sound speed, Inverse Problems, 25 (2009), 075011, 16pp. doi: 10.1088/0266-5611/25/7/075011.
[38]	M. Suhonen, A. Pulkkinen and T. Tarvainen, One-step approach for estimating optical parameters in spectral quantitative photoacoustic tomography, Biomed Opt Express, (2023).
[39]	T. Tarvainen and B. Cox, Quantitative photoacoustic tomography: Modeling and inverse problems, Journal of Biomedical Optics, 29 (2024), S11509-S11509.
[40]	J. Tick, A. Pulkkinen, F. Lucka, R. Ellwood, B. T. Cox, J. P. Kaipio, S. R. Arridge and T. Tarvainen, Three dimensional photoacoustic tomography in Bayesian framework, J Acoust Soc Am, 144 (2018), 2061-2071. doi: 10.1121/1.5057109.
[41]	J. Tick, A. Pulkkinen and T. Tarvainen, Image reconstruction with uncertainty quantification in photoacoustic tomography, J Acoust Soc Am, 139 (2016), 1951-1961. doi: 10.1121/1.4945990.
[42]	J. Tick, A. Pulkkinen and T. Tarvainen, Bayesian approach to image reconstruction in photoacoustic tomography, In A.A. Oraevsky and L.V. Wang, editors, Photons Plus Ultrasound: Imaging and Sensing 2017, Proc of SPIE, 10064 (2017), 100643M.
[43]	B. E. Treeby and B. T. Cox, k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields, J Biomed Opt, 15 (2010), 021314.
[44]	B. E. Treeby, E. Z. Zhang and B. T. Cox, Photoacoustic tomography in absorbing acoustic media using time reversal, Inv Probl, 26 (2010), 115003, 20pp. doi: 10.1088/0266-5611/26/11/115003.
[45]	K. Wang and M. A. Anastasio, Photoacoustic and thermoacoustic tomography: Image formation principles, In Handbook of Mathematical Methods in Imaging, Springer New York, New York, NY, (2011), 781-815. doi: 10.1007/978-0-387-92920-0_18.
[46]	K. Wang, R. Su, A. A. Oraevsky and M. A Anastasio, Investigation of iterative image reconstruction in three-dimensional optoacoustic tomography, Phys Med Biol, 57 (2012), 5399-5423. doi: 10.1088/0031-9155/57/17/5399.
[47]	M. Xu and L. V. Wang, Universal back-projection algorithm for photoacoustic computed tomography, Phys Rev E, 71 (2005), 016706. doi: 10.1103/PhysRevE.71.016706.
[48]	Y. Xu and L. V. Wang, Time reversal and its application to tomography with diffracting sources, Phys Rev Lett, 92 (2004), 033902. doi: 10.1103/PhysRevLett.92.033902.
[49]	Y. Xu, L. V. Wang, G. Ambartsoumian and P. Kuchment, Reconstructions in limited-view thermoacoustic tomography, Med Phys, 31 (2004), 724-733. doi: 10.1118/1.1644531.