Fast algorithms enabling optimization and deep learning for photoacoustic tomography in a circular detection geometry

TL;DR

The paper tackles the computational bottleneck of reconstructing initial pressure in photoacoustic/thermoacoustic tomography with circular data by developing Hankel-free, FFT-based algorithms for the forward operator $\mathcal{A}$, its adjoint $\mathcal{A}^{\ast}$, and an accurate fast inverse $\mathcal{A}^{-1}$, all in $2$D with complexity $O(n^{2}\log n)$. It derives these operators in a way that avoids explicit Hankel function evaluations, enabling stable and scalable use inside iterative variational schemes and model-based deep learning, such as learned primal dual (LPD). The methods are validated through extensive numerical simulations across complete, partial, and incomplete data, showing that fast operator evaluation substantially speeds up NNLS, total variation regularization, and LPD reconstructions, with quantitative improvements in recovery quality and computational efficiency. The authors also release open-source Python/Torch implementations to promote adoption and further research.

Abstract

The inverse source problem arising in photoacoustic tomography and in several other coupled-physics modalities is frequently solved by iterative algorithms. Such algorithms are based on the minimization of a certain cost functional. In addition, novel deep learning techniques are currently being investigated to further improve such optimization approaches. All such methods require multiple applications of the operator defining the forward problem, and of its adjoint. In this paper, we present new asymptotically fast algorithms for numerical evaluation of the forward and adjoint operators, applicable in the circular acquisition geometry. For an $(n \times n)$ image, our algorithms compute these operators in $\mathcal{O}(n^2 \log n)$ floating point operations. We demonstrate the performance of our algorithms in numerical simulations, where they are used as an integral part of several iterative image reconstruction techniques: classic variational methods, such as non-negative least squares and total variation regularized least squares, as well as deep learning methods, such as learned primal dual. A Python implementation of our algorithms and computational examples is available to the general public.

Figures (8)

• Figure 1: Illustration of photoacoustic system MSOT inVision by iThera Medical with the data acquired by a ring of transducers.
• Figure 2: Reconstruction error in the fast forward computation
• Figure 3: Two samples (a) and (b) from the training data with each ground truth $f$ (left) and computed $g=\mathcal{A}f + \delta$ (right). The dotted white lines indicate the location of the transducers.
• Figure 4: 360 deg acquisition, accurate data
• Figure 5: 360 deg acquisition, 30% noise in the data