Photoacoustic Tomography with Direction Dependent Data: An Exact Series Reconstruction Approach

TL;DR

This paper develops exact series reconstruction formulas for photoacoustic tomography with direction-dependent data, accounting for detectors with directivity and frequency dependence, and demonstrates their robustness through numerical results.

Contribution

It introduces a novel exact frequency domain reconstruction method for data combining acoustic pressure and its normal derivative, specifically for spherical detection geometries.

Findings

Reconstruction formulas are robust and accurate.

Measurement model significantly influences the recovered initial pressure.

Numerical results validate the theoretical formulas.

Abstract

Photoacoustic image reconstruction often assumes that the restriction of the acoustic pressure on the detection surface is given. However, commonly used detectors often have a certain directivity and frequency dependence, in which case the measured data are more accurately described as a linear combination of the acoustic pressure and its normal derivative on the detection surface. In this paper, we consider the inverse source problem for data that are a combination of an acoustic pressure of the wave equation and its normal derivative For the special case of a spherical detection geometry we derive exact frequency domain reconstruction formulas. We present numerical results showing the robustness and validity of the derived formulas. Moreover, we compare several different combinations of the pressure and its normal derivative showing that used measurement model significantly affects the recovered initial pressure.

Equations

$${\partial }_t^{2}p(\mathbf{x},t)-{\Delta }_{\mathbf{x}}p(\mathbf{x},t)$$

$$p(\mathbf{x},0)$$

$${\partial }_{t}p(\mathbf{x},0)$$

$$m_{c_1,c_2}(\mathbf{x},t)=c_{1}p(\mathbf{x},t)+c_{2}n(\mathbf{x})\text{vbox}\text{scalebox}0.6\bullet \nabla p(\mathbf{x},t)\text{ for }(\mathbf{x},t)\in \partial \Omega \times [0,T],$$

$$\displaystyle g(\boldsymbol{\theta},t):=c_{1}p(R\boldsymbol{\theta},t)+c_{2}\Bigl{[}{\boldsymbol{\theta}}\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.6}{$\scriptscriptstyle\bullet$}}}}}{\nabla}p(\mathbf{x},t)\bigr{]}_{\mathbf{x}=R\boldsymbol{\theta}}\quad\text{ for }(\boldsymbol{\theta},t)\in\mathbb{S}^{n-1}\times[0,T]\,.$$

$$\hat{f}(\mathbf{\xi })$$

$$\mathbf{C}\left\{g\right\}(\lambda )$$

$$\mathbf{S}\left\{g\right\}(\lambda )$$

$$\forall f\in L^2({\mathbb{R}}^n):f(\rho \mathbf{\theta })=\sum _{\ell =0}^{\infty }\sum _{k=0}^{N(n,\ell )}{f}_{\ell ,k}(\rho )Y_{\ell ,k}(\mathbf{\theta }),$$

$$\forall (t,\mathbf{\theta })\in (0,\infty )\times {\mathbb{S}}^{n-1}:g(\mathbf{\theta },t)=\sum _{\ell =0}^{\infty }\sum _{k=0}^{N(n,\ell )}{g}_{\ell ,k}(t)Y_{\ell ,k}(\mathbf{\theta })$$

$$\forall (\lambda ,\mathbf{\omega })\in (0,\infty )\times {\mathbb{S}}^{n-1}:\hat{f}(\lambda \mathbf{\omega })=\sum _{\ell =0}^{\infty }\sum _{k=0}^{N(n,\ell )}{\hat{f}}_{\ell ,k}(\lambda )Y_{\ell ,k}(\mathbf{\omega }).$$

$$\begin{gathered} \forall \lambda >0:\mathbf{C}\{g_{\ell ,k}\}(\lambda )=(2R\pi {)}^{-\frac{n}{2}}\pi 2^{-1}{i}^{\ell } \\ \times \left[(R{c}_1+c_2\ell )J_{\ell +\frac{n-2}{2}}(R\lambda )-c_{2}R\lambda J_{\ell +\frac{n}{2}}(R\lambda )\right]{\hat{f}}_{\ell ,k}(\lambda ){\lambda }^{\frac{n}{2}}. \end{gathered}$$

$$g(\mathbf{\theta },t)$$

$$=(2\pi {)}^{-n}{\int }_{{\mathbb{S}}^{n-1}}{\int }_0^{\infty }[c_1+c_2(i\mathbf{\theta }\text{vbox}\text{scalebox}0.6\bullet \lambda \mathbf{\omega })]\cos (t\lambda )e^{i\lambda R\mathbf{\theta }\text{vbox}\text{scalebox}0.6\bullet \mathbf{\omega }}\hat{f}(\lambda \mathbf{\omega }){\lambda }^{n-1}\mathrm{d}\lambda \mathrm{d}S(\mathbf{\omega })$$

$$=(2\pi {)}^{-n}{\int }_{{\mathbb{S}}^{n-1}}{\int }_0^{\infty }\left(c_1{e}^{i\lambda R\mathbf{\theta }\text{vbox}\text{scalebox}0.6\bullet \mathbf{\omega }}+c_2\lambda \left({\partial }_{R\lambda }{e}^{i\lambda R\mathbf{\theta }\text{vbox}\text{scalebox}0.6\bullet \mathbf{\omega }}\right)\right)\cos (t\lambda )\hat{f}(\lambda \mathbf{\omega }){\lambda }^{n-1}\mathrm{d}\lambda \mathrm{d}S(\mathbf{\omega })$$

$$=(2\pi {)}^{-n}\sum _{\ell =0}^{\infty }\sum _{k=0}^{N(n,\ell )}{\int }_0^{\infty }\cos (t\lambda ){\hat{f}}_{\ell ,k}(\lambda ){\lambda }^{n-1}\left(c_1{A}_{\ell ,k}(\mathbf{\theta })+c_2{B}_{\ell ,k}(\mathbf{\theta })\right)\mathrm{d}\lambda ,$$

$$A_{\ell ,k}(\mathbf{\theta })$$

$$B_{\ell ,k}(\mathbf{\theta })$$

$$A_{\ell ,k}(\mathbf{\theta })=(2\pi {)}^{\frac{n}{2}}{i}^{\ell }(R\lambda {)}^{\frac{2-n}{2}}{J}_{\ell +\frac{n-2}{2}}(R\lambda )Y_{\ell ,k}(\mathbf{\theta }).$$

$$B_{\ell ,k}(\mathbf{\theta })$$

$$=\lambda (2\pi {)}^{\frac{n}{2}}{i}^{\ell }{Y}_{\ell ,k}(\mathbf{\theta }){\partial }_{R\lambda }\left[(R\lambda {)}^{-(n-2)/2}{J}_{\ell +\frac{n-2}{2}}(R\lambda )\right]$$

$$=\lambda (2\pi {)}^{\frac{n}{2}}{i}^{\ell }{Y}_{\ell ,k}(\mathbf{\theta }){\partial }_{R\lambda }\left[(R\lambda {)}^{\ell }(R\lambda {)}^{-(\ell +(n-2)/2)}{J}_{\ell +\frac{n-2}{2}}(R\lambda )\right]$$

$$=\lambda (2\pi {)}^{\frac{n}{2}}{i}^{\ell }{Y}_{\ell ,k}(\mathbf{\theta })\left[\ell (R\lambda {)}^{-n/2}{J}_{\ell +\frac{n-2}{2}}(R\lambda )-(R\lambda {)}^{-(n-2)/2}{J}_{\ell +n/2}(R\lambda )\right]$$

$$=(2\pi {)}^{\frac{n}{2}}{i}^{\ell }(R\lambda {)}^{\frac{2-n}{2}}\left[R^{-1}\ell J_{\ell +\frac{n-2}{2}}(R\lambda )-\lambda J_{\ell +n/2}(R\lambda )\right]Y_{\ell ,k}(\mathbf{\theta }).$$

$$\begin{gathered} g(\mathbf{\theta },t)=R^{-\frac{n}{2}}(2\pi {)}^{-\frac{n}{2}}\sum _{\ell =0}^{\infty }\sum _{k=0}^{N(n,\ell )}{i}^{\ell }{\int }_0^{\infty }\cos (t\lambda ){\hat{f}}_{\ell ,k}(\lambda )(\lambda {)}^{\frac{n}{2}} \\ \times \left((R{c}_1+c_2\ell )J_{\ell +\frac{n-2}{2}}(R\lambda )-c_{2}R\lambda J_{\ell +n/2}(R\lambda )\right)\mathrm{d}\lambda Y_{\ell ,k}(\mathbf{\theta }). \end{gathered}$$

$$\begin{gathered} g_{\ell ,k}(t)=R^{-\frac{n}{2}}(2\pi {)}^{-\frac{n}{2}}{i}^{\ell }{\int }_0^{\infty }\cos (t\lambda ){\hat{f}}_{\ell ,k}(\lambda ){\lambda }^{\frac{n}{2}} \\ \times \left((R{c}_1+c_2\ell )J_{\ell +\frac{n-2}{2}}(R\lambda )-c_{2}R\lambda J_{\ell +n/2}(R\lambda )\right)\mathrm{d}\lambda . \end{gathered}$$

$$\begin{gathered} f(\mathbf{x})=\frac{2}{\pi }\sum _{\ell =0}^{\infty }\sum _{k=0}^{N(n,\ell )}{\left(\frac{R}{|\mathbf{x}|}\right)}^{\frac{n-2}{2}} \\ \times {\int }_0^{\infty }\frac{\mathbf{C}\{g_{\ell ,k}\}(\lambda )J_{\ell +\frac{n-2}{2}}(\lambda |\mathbf{x}|)\mathrm{d}\lambda }{(c_1+c_2{R}^{-1}\ell )J_{\ell +\frac{n-2}{2}}(R\lambda )-c_2\lambda J_{\ell +\frac{n}{2}}(R\lambda )}{Y}_{\ell ,k}\left(\frac{\mathbf{x}}{|\mathbf{x}|}\right). \end{gathered}$$

$${\hat{f}}_{\ell ,k}(\lambda )$$

$$={\int }_{{\mathbb{S}}^{n-1}}{\int }_0^{\infty }f(\rho \mathbf{\theta }){\rho }^{n-1}\left({\int }_{{\mathbb{S}}^{n-1}}{e}^{-i\lambda \rho \mathbf{\theta }\text{vbox}\text{scalebox}0.6\bullet \mathbf{\omega }}(Y_{\ell ,k}(\mathbf{\omega }){)}^{\ast }\mathrm{d}S(\mathbf{\omega })\right)\mathrm{d}\rho \mathrm{d}S(\mathbf{\theta })$$

$$={\int }_{{\mathbb{S}}^{n-1}}{\int }_0^{\infty }f(\rho \mathbf{\theta }){\rho }^{n-1}{\left({\int }_{{\mathbb{S}}^{n-1}}{e}^{i\lambda \rho \mathbf{\theta }\text{vbox}\text{scalebox}0.6\bullet \mathbf{\omega }}{Y}_{\ell ,k}(\mathbf{\omega })\mathrm{d}S(\mathbf{\omega })\right)}^{\ast }\mathrm{d}\rho \mathrm{d}S(\mathbf{\theta })$$

$$=(2\pi {)}^{\frac{n}{2}}(i^{\ell }{)}^{\ast }{\int }_{{\mathbb{S}}^{n-1}}{\int }_0^{\infty }f(\rho \mathbf{\theta }){\rho }^{\frac{n}{2}}{\lambda }^{\frac{2-n}{2}}{J}_{\ell +\frac{n-2}{2}}(\lambda \rho )(Y_{\ell ,k}(\mathbf{\omega }){)}^{\ast }\mathrm{d}\rho \mathrm{d}S(\mathbf{\theta })$$

$${\hat{f}}_{\ell ,k}(\lambda )=(2\pi {)}^{\frac{n}{2}}(i^{\ell }{)}^{\ast }{\lambda }^{\frac{2-n}{2}}{\int }_0^{\infty }{f}_{\ell ,k}(\rho ){\rho }^{\frac{n}{2}}{J}_{\ell +\frac{n-2}{2}}(\lambda \rho )\mathrm{d}\rho .$$

$${\int }_0^{\infty }{f}_{\ell ,k}(\rho ){\rho }^{\frac{n}{2}}{J}_{\ell +\frac{n-2}{2}}(\lambda \rho )\mathrm{d}\rho =\frac{2}{\pi }\frac{R^{\frac{n}{2}}\mathbf{C}\{g_{\ell ,k}\}(\lambda )}{(R{c}_1+c_2\ell )\lambda J_{\ell +\frac{n-2}{2}}(R\lambda )-c_{2}R{\lambda }^2{J}_{\ell +\frac{n}{2}}(R\lambda )}.$$

$$f_{\ell ,k}(\rho )=R^{\frac{n}{2}}{\rho }^{-\frac{n-2}{2}}\frac{2}{\pi }{\int }_0^{\infty }\frac{\mathbf{C}\{g_{\ell ,k}\}(\lambda )J_{\ell +\frac{n-2}{2}}(\lambda \rho )\mathrm{d}\lambda }{(R{c}_1+c_2\ell )J_{\ell +\frac{n-2}{2}}(R\lambda )-c_{2}R\lambda J_{\ell +\frac{n}{2}}(R\lambda )}.$$

$$\begin{gathered} f(\mathbf{x})={\left(\frac{R}{|\mathbf{x}|}\right)}^{\frac{n-2}{2}}\frac{4}{c_2\pi }\sum _{\ell =0}^{\infty }\sum _{k=0}^{N(n,\ell )} \\ \left(\sum _{j=1}^{\infty }\frac{\mathbf{C}\{g_{\ell ,k}\}\left(\frac{w_{j,\ell }}{R}\right)}{w_{j,\ell }^2{J}_{\ell +\frac{n}{2}}(w_{j,\ell }{)}^3}{J}_{\ell +\frac{n-2}{2}}\left(\frac{w_{j,\ell }|\mathbf{x}|}{R}\right)\right)Y_{\ell ,k}\left(\frac{\mathbf{x}}{|\mathbf{x}|}\right). \end{gathered}$$

$$\begin{gathered} f(\mathbf{x})={\left(\frac{R}{|\mathbf{x}|}\right)}^{\frac{n-2}{2}}\frac{4}{c_1{R}^2\pi }\sum _{\ell =0}^{\infty }\sum _{k=0}^{N(n,\ell )} \\ \left(\sum _{j=1}^{\infty }\frac{\mathbf{S}\left\{t{g}_{\ell ,k}(t)\right\}\left(\frac{w_{j,\ell }}{R}\right)}{w_{j,\ell }{J}_{\ell +\frac{n}{2}}(w_{j,\ell }{)}^3}{J}_{\ell +\frac{n-2}{2}}\left(\frac{w_{j,\ell }|\mathbf{x}|}{R}\right)\right)Y_{\ell ,k}\left(\frac{\mathbf{x}}{|\mathbf{x}|}\right). \end{gathered}$$

$$\begin{gathered} f_{\ell ,k}(\rho ){\rho }^{\frac{n-2}{2}}=\sum _{j=1}^{\infty }\frac{2{\int }_0^{\infty }{f}_{\ell ,k}(r)J_{\ell +\frac{n-2}{2}}\left(\frac{w_{j,\ell }r}{R}\right)r^{\frac{n}{2}}\mathrm{d}r}{R^2{J}_{\ell +\frac{n}{2}}(w_{j,\ell }{)}^2}{J}_{\ell +\frac{n-2}{2}}\left(\frac{w_{j,\ell }\rho }{R}\right) \\ =(2\pi {)}^{-\frac{n}{2}}(i^{-\ell }{)}^{\ast }\sum _{j=1}^{\infty }\frac{2R^{-\frac{n+2}{2}}}{J_{\ell +\frac{n}{2}}(w_{j,\ell }{)}^2}{\hat{f}}_{\ell ,k}\left(\frac{w_{j,\ell }}{R}\right)w_{j,\ell }^{\frac{n-2}{2}}{J}_{\ell +\frac{n-2}{2}}\left(\frac{w_{j,\ell }\rho }{R}\right), \end{gathered}$$

$$\begin{gathered} \mathbf{C}\{g_{\ell ,k}\}\left(\frac{w}{R}\right)=(2R\pi {)}^{-\frac{n}{2}}\pi 2^{-1}{i}^{\ell }{\hat{f}}_{\ell ,k}\left(\frac{w}{R}\right)R^{-\frac{n}{2}}{w}^{\frac{n}{2}} \\ \times \left[(R{c}_1+c_2\ell )J_{\ell +\frac{n-2}{2}}(w)-c_{2}w{J}_{\ell +\frac{n}{2}}(w)\right]. \end{gathered}$$

Full text

Keywords: Photoacoustic tomography; spherical geometry; image reconstruction; wave equation; series inversion; reconstruction formula; Neumann data.

AMS subject classifications: 44A12, 65R32, 35L05, 92C55.

1 Introduction

Photoacoustic Tomography (PAT) is a hybrid imaging technique that combines high optical contrast with good ultrasonic resolution. It is based on the generation of an ultrasound wave inside an object of interest by pulses optical illumination. The initial pressure distribution of the induced sound wave encodes the optical absorption properties of the object, which are of great interest in medical diagnostics. PAT has proven to be very promising for medical applications like functional brain imaging of small animals, early cancer diagnostics, and imaging of vasculature [2, 29].

In a typical PAT setup, the induced acoustic waves are recorded by several point-like detectors outside the support of the investigates sample. Typically, the detectors are located on a surface $S$ that fully or partially encloses the volume $\Omega$ in which the object of interest is contained; compare Figure 1.1. In most reconstruction approaches, the measured data are identified with the restriction of the acoustic wave to the surface $\partial\Omega$ or sampled values of it, possible convolved in time with the detector impulse response function. Typical piezoelectric sensors, which are commonly used in PAT, however, have a directional dependence and are most sensitive in the normal direction to the sensor surface. Moreover, at the resonance frequency they are more sensitive than at lower frequencies. As noted in [30], such measurement data are rather modeled by the combination of the pressure field and its normal derivative, than the pressure alone. This is also suggested by visual comparison with real data [27, 26]. Systematic theoretical and experimental studies on PAT sensor modeling are interesting lines of future research.

Let $f$ be the smooth compactly supported initial pressure, and $p\colon\mathbb{R}^{n}\times(0,\infty)\to\mathbb{R}$ the induced pressure wave that satisfies

[TABLE]

In order to incorporate directional dependence, we model the data of a detector located at $\mathbf{x}\in\partial\Omega$ by

[TABLE]

where {n(\mathbf{x})}\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.6}{\displaystyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.6}{\textstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.6}{\scriptstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.6}{\scriptscriptstyle\bullet}}}}}{\nabla p(\mathbf{x},t)} is the normal derivative of the pressure, $n(\mathbf{x})$ the outwards pointing normal of the surface at $\mathbf{x}$, $T$ the final measurement time, and $c_{1},c_{2}\in\mathbb{R}$ are constants. The goal is to recover the initial data $f$ from data in (1.2). To the best of our knowledge, this paper is the first work investigating PAT with direction dependent data of the form (1.2) allowing arbitrary values of $c_{1},c_{2}$.

The standard image reconstruction problem in PAT corresponds to the special case $c_{1}\neq 0$ and $c_{2}=0$ in our data model (1.2). Many methods for reconstructing the initial pressure distribution have been derived in the recent years in various situations. This includes, for example, different detection geometries with variable or constant speed of sound, or limited view situations. Theoretical questions concerning uniqueness and stability of the inverse source problem have also been investigated [14, 16, 18, 3, 28]. A practically important case assumes a constant speed of sound. In this case, several exact reconstruction formulas have been developed [19, 20, 13, 12, 24, 8, 23, 25, 11, 1, 7, 31, 17, 15]. Among these formulas so-called series expansion formulas provide very fast and accurate reconstructions. They give an expansion of the initial pressure in terms of eigenfunctions of the Laplacian [11, 1, 31, 17, 15].

The case where measurements are modeled by the normal derivative of the pressure ($c_{1}=0$ and $c_{2}\neq 0$ in (1.2)) is studied much less. It has been considered in [6, 9], where an explicit inversion formula of the backprojection type is derived for the case that the detection surface is a sphere in three spatial dimensions. Besides that, we are not aware of any results for the case $c_{2}\neq 0$. We are not aware of any results when $c_{1},c_{2}$ are both non-vanishing.

In this paper, we provide series expansion reconstruction formulas for the direction dependent data model (1.2) allowing arbitrary values of $c_{1},c_{2}\in\mathbb{R}$, for the case that the measurement surface is a sphere. Our approach is based on expansions in spherical harmonics and an explicit formula relating the spherical harmonics coefficients of the direction dependent data, as a function of time, and the Fourier coefficients of the initial pressure distribution $f$, as a function of distance to the origin (see Section 2). By using Fourier Bessel series in the radial variable, in Section 3 we derive inversion formulas that allows for a stable implementation. Numerical implementation and results are presented in Section 4. The paper ends with a conclusion and outlook presented in Section 5

2 Preliminaries

In PAT, sound propagation is commonly described by an acoustic pressure wave $p\colon\mathbb{R}^{n}\times[0,\infty)\to\mathbb{R}$ that satisfies the initial value problem(1). In this work, we assume that the initial pressure distribution satisfies $f\in\mathcal{C}_{0}^{\infty}\left(\Omega\right)$. In PAT with direction dependent data we model measurement data by (1.2), where $p$ is the solution of (1). The aim is to recover the initial pressure distribution $f$ from such data. To the best of our knowledge, this is a new inverse problem for the wave equation that has not been considered so far. For the practical application, the cases $n=2$ and $n=3$ are of relevance. In particular, the case $n=2$ appears from direction dependent measurements with integrating line detectors [4, 15].

In particular, we study the case where $\Omega=B_{R}(0)$ is the ball of radius $R$ centered at the origin. In this situation, we can write the measurement data in the form

[TABLE]

We write $\mathbf{M}_{c_{1},c_{2}}$ for the operator that takes the initial data $f$ in (1) to the corresponding boundary data. The constants $c_{1},c_{2}\in\mathbb{R}$ are weights describing the contribution of the pressure and its normal derivative, respectively, to the measures data. The operator $\mathbf{M}_{1,0}$ is the standard PAT forward operator, and the operator $\mathbf{M}_{0,1}$ corresponds to the case where the data only consists of the normal derivative. Clearly, we have $\mathbf{M}_{c_{1},c_{2}}=c_{1}\mathbf{M}_{1,0}+c_{2}\mathbf{M}_{0,1}$.

For the derivation of the inversion formulas we use $T=\infty$. Because of the strict Huygens principle, in the case of odd spatial dimension, we have $p(\mathbf{x},t)=0$ for $t>2R$, which yields exact inversion for any measurement time $T\geq 2R$. For even spatial dimension, this is not the case and, strictly taken, our inversion formulas are exact only for $T=\infty$. Deriving exact formulas using data over a finite time interval $[0,T]$ in even dimension is an interesting open issue. For the standard PAT operator $\mathbf{M}_{1,0}$ such a formula has been derived in [7, Theorem 1.4].

2.1 Notation

Our approach is based on Fourier methods that lead to a relation between measurement data and the initial pressure distribution in the frequency domain. We denote the involved transforms by

[TABLE]

which are the Fourier, cosine and sine transforms, respectively.

Moreover, we denote by $Y_{\ell,k}(\boldsymbol{\theta})$ the spherical harmonics [22], which form a complete orthonormal system in $L^{2}(\mathbb{S}^{n-1})$. In particular,

[TABLE]

where $N(n,\ell)=(2\ell+n-2)(n+\ell-3)!/(\ell!(n-2)!)$ for $\ell\in\mathbb{N}$ and $N(n,0):=1$. We write $J_{\nu}$ for the $\nu$-th order Bessel function of the first kind and denote by $w_{j,\ell}$ the $j$-th positive root of $J_{\ell+(n-2)/2}$ .

2.2 Relations between transform coefficients

For the following we use the spherical harmonics expansions of the measurement data $g=\mathbf{M}_{c_{1},c_{2}}f$ in (2.1) and the Fourier transform $\hat{f}$ of the initial pressure,

[TABLE]

The following Lemma is the key to our results.

Lemma 1**.**

Let $f\in\mathcal{C}^{\infty}_{0}(\mathbb{R}^{n})$ be the initial pressure distribution in (1) and let $g=\mathbf{M}_{c_{1},c_{2}}f$ be the corresponding measurement data given by (2.1). Then,

[TABLE]

Proof.

Let $p$ be the solution of (1). Then p(\mathbf{x},t)=(2\pi)^{-n}\int_{\mathbb{R}^{n}}\cos(t|\boldsymbol{\xi}|)e^{i{\mathbf{x}}\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.6}{\displaystyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.6}{\textstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.6}{\scriptstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.6}{\scriptscriptstyle\bullet}}}}}{\boldsymbol{\xi}}}\hat{f}(\boldsymbol{\xi})\mathrm{d}\boldsymbol{\xi} and therefore {\boldsymbol{\theta}}\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.6}{\displaystyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.6}{\textstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.6}{\scriptstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.6}{\scriptscriptstyle\bullet}}}}}{\nabla p}(\mathbf{x},t)=(2\pi)^{-n}\int_{\mathbb{R}^{n}}({i}{\boldsymbol{\theta}}\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.6}{\displaystyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.6}{\textstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.6}{\scriptstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.6}{\scriptscriptstyle\bullet}}}}}{\boldsymbol{\xi}})\cos(t|\boldsymbol{\xi}|)e^{{i}{\mathbf{x}}\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.6}{\displaystyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.6}{\textstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.6}{\scriptstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.6}{\scriptscriptstyle\bullet}}}}}{\boldsymbol{\xi}}}\hat{f}(\boldsymbol{\xi})\mathrm{d}\boldsymbol{\xi} for all $(\mathbf{x},t)\in\mathbb{R}^{n}\times\mathbb{R}$. Changing to spherical coordinates, $\lambda\boldsymbol{\omega}\leftarrow\boldsymbol{\xi}$, and using the spherical harmonics expansion (2.2) we can write the direction dependent measurements as

[TABLE]

where we used the abbreviations

[TABLE]

According to the Fuck-Hecke theorem [22]

[TABLE]

Moreover, using the identity $\partial_{R\lambda}[(R\lambda)^{-\nu}J_{\nu}(R\lambda)]=-(R\lambda)^{-\nu}J_{\nu+1}(R\lambda)$ (see, for example, [10]), we further obtain

[TABLE]

The last two displayed equations imply $c_{1}A_{\ell,k}(\boldsymbol{\theta})+c_{2}B_{\ell,k}(\boldsymbol{\theta})=(2\pi)^{\frac{n}{2}}i^{\ell}\,(R\lambda)^{\frac{2-n}{2}}$ $Y_{\ell,k}(\boldsymbol{\theta})((c_{1}+R^{-1}c_{2}\ell)J_{\ell+\frac{n-2}{2}}(R\lambda)-c_{2}\lambda\,J_{\ell+n/2}(R\lambda))$ and therefore

[TABLE]

Comparing expansions (2.2) and (2.6) we conclude that

[TABLE]

This becomes (2.4) after applying the inversion formula $g=\frac{2}{\pi}\mathbf{C}\mathbf{C}{g}$ for the cosine transform; see [10, Equation (7.28)]. ∎

As a first application of Lemma 1, we obtain a range condition for the classical PAT forward operator $\mathbf{M}_{1,0}$ that maps the initial data $f$ to the solution of (1) on the boundary. More precisely, we have the following result.

Corollary 1** (Range condition).**

Let $f\in\mathcal{C}_{0}^{\infty}(B_{R}(0))$. Then $\mathbf{C}\{p_{\ell,k}\}(R^{-1}w_{j,\ell})=0$ where $(p_{\ell,k})_{\ell,k}$ are the spherical harmonics coefficients of $\mathbf{M}_{1,0}f$ and $w_{j,\ell}$ is the $j$-th positive zero of $J_{\ell+(n-2)/2}$.

Proof.

For $c_{1}=1$ and $c_{2}=0$, equation (2.4) reads $\mathbf{C}\{p_{\ell,k}\}(\lambda)=(2R\pi)^{-{\frac{n}{2}}}\pi 2^{-1}i^{\ell}\hat{f}_{\ell,k}(\lambda)$ $\lambda^{n/2}RJ_{\ell+(n-2)/2}(R\lambda)$. Because $\hat{f}_{\ell,k}(\lambda)$ is continuous, this implies $\mathbf{C}(p_{\ell,k})(R^{-1}w_{j,\ell})=0$ and concludes the proof. ∎

As another application of Lemma 1, we derive the following inversion formula.

Corollary 2**.**

Let $f\in\mathcal{C}^{\infty}_{0}\left(\mathbb{R}^{n}\right)$, let $c_{1},c_{2}\in\mathbb{R}$, and consider the measurement data $g:=\mathbf{M}_{c_{1},c_{2}}f$ as in (2.1). Then

[TABLE]

Proof.

Using (2.2), (2.5) and writing $(\,\cdot\,)^{*}$ for the complex conjugation we have

[TABLE]

and therefore

[TABLE]

Together with Lemma 1 this gives

[TABLE]

The left hand side in (2.9) is recognized as the Hankel transform of order $\ell+(n-2)/2$ of of the function $\rho\mapsto f_{\ell,k}(\rho)\rho^{(n-2)/2}$. Hence, by applying the inverse Hankel transform, we obtain

[TABLE]

Together with (2.3), this gives the desired result. ∎

Corollary 2 gives an exact inversion formula for reconstructing any smooth compactly supported function $f\in\mathcal{C}^{\infty}_{0}\left(\mathbb{R}^{n}\right)$ from data $\mathbf{M}_{c_{1},c_{2}}f$. In order to actually implement the formula one requires a proper discretization of the integral. In real situations only noisy data $g^{\delta}\in L^{2}\left(S_{R}\times(0,\infty)\right)$ are available and the cosine transform $\mathbf{C}\{g^{\delta}_{\ell,k}\}(\lambda)$ will not vanish at the roots of the denominator $(c_{1}+c_{2}\ell)J_{\ell+\frac{n-2}{2}}(\lambda)-c_{2}\lambda J_{\ell+n/2}(\lambda)$ that appears in formula (2.7). This means that Corollary 2 behaves unstable close to the roots and cannot be directly used to reconstruct the initial pressure density $f$. In this following section, we avoid the zeros of the denominator by using a Fourier Bessel series expansion similar as in [15, 31].

3 Stable series inversion formulas

In practice, reconstructing the initial pressure by the inversion formula stated in Corollary 2 is unstable due to the zeros of $(c_{1}+c_{2}\ell)J_{\ell+\frac{n-2}{2}}(\lambda)-c_{2}\lambda J_{\ell+n/2}(\lambda)$ in the denominator of (2.7). To avoid this issue, in this section we derive a series inversion formula using a Fourier Bessel series, similar to [15, 31]. We derive different inversion formulas for the cases $c_{2}\neq 0$ and $c_{2}=0$, respectively. For the following recall that $w_{j,\ell}$ denotes the $j$-th positive zero of $J_{\ell+(n-2)/2}$.

The following theorem is the main result of this paper.

Theorem 1** (Explicit series inversion formula).**

Let $f\in\mathcal{C}_{0}^{\infty}(B_{R}(0))$, let $c_{1},c_{2}\in\mathbb{R}$, and consider the measurement data $g:=\mathbf{M}_{c_{1},c_{2}}f$ as in (2.1).

1. (a)

If $c_{2}\neq 0$, then

[TABLE] 2. (b)

If $c_{2}=0$, then

[TABLE]

Proof.

Because $f_{\ell,k}$ is of class $C^{\infty}$ and is compactly supported, we can expand $f_{\ell,k}(\rho)\rho^{{\frac{n-2}{2}}}$ into a Fourier Bessel series [10],

[TABLE]

where for the second line, we used (2.8). From Lemma 1, we have

[TABLE]

• $\blacksquare$

For $c_{2}\neq 0$ we have $(Rc_{1}+c_{2}\ell)J_{\ell+\frac{n-2}{2}}(w_{j,\ell})-c_{2}w_{j,\ell}J_{\ell+\frac{n}{2}}(w_{j,\ell})=-c_{2}w_{j,\ell}J_{\ell+\frac{n}{2}}(w_{j,\ell})\neq 0$. We can therefore evaluate (3.4) at $w=w_{j,\ell}$ and solve for $\hat{f}_{\ell,k}\left(\frac{w}{R}\right)$. Together with (3.3), we obtain

[TABLE]

• $\blacksquare$

For the case $c_{2}=0$, equation (3.4) gives

[TABLE]

For the zeros $w=w_{j,\ell}$ of $J_{\ell+\frac{n-2}{2}}(w)$ this is an indeterminate form, which can be evaluated with L’Hospital’s rule. We have

[TABLE]

Using this and the identity $\partial_{w}\mathbf{C}\{g_{\ell,k}\}(wR^{-1})=-R^{-1}\mathbf{S}\{tg_{\ell,k}(t)\}(wR^{-1})$, and applying L’Hospital’s rule give

[TABLE]

Therefore equation (3.3) yields

[TABLE]

By combing (3.5) and (3.6) with the spherical harmonics expansion (2.3), we obtain the desired inversion formulas (3.1) and (3.2). ∎

Note that the inversion formula (3.1) for $\mathbf{M}_{c_{1},c_{2}}$ might be slightly surprising because it holds for any $c_{2}\neq 0$ regardless of the value of $c_{1}$. However, this independence is an immediate consequence of the decomposition $\mathbf{M}_{c_{1},c_{2}}f=c_{1}\mathbf{M}_{1,0}+c_{2}\mathbf{M}_{0,1}$, the range condition in Proposition 1 for the operator $\mathbf{M}_{1,0}$, and the particular form of the right hand side in (3.1). Equation (3.2), on the other hand, is a new series expansion formula for the standard PAT forward operator $\mathbf{M}_{1,0}$. In the special case $n=2$, it becomes the series inversion formula [15, Theorem 3.1].

4 Numerical experiments

In this section, we present reconstruction results with the inversion formulas in Theorem 1 using data $g=\mathbf{M}_{c_{1},c_{2}}f$ for different combinations of $c_{1},c_{2}$. We consider the case of $n=2$ spatial dimensions and use $R=1$. The 2D case arises in applications where the acoustic pressure is measured by integrating line detectors [4, 26].

4.1 Discretization and data simulation

In order to compute the pressure field and its gradient restricted to the unit circle, we use a discrete wave propagation model on a 2D equidistant grid. To numerically compute the solution at the discretization points we use the k-space method [5, 21] implemented as described in [14]. Using the k-space method, we obtain the values of the acoustic pressure field and compute its gradient by symmetric differences. The pressure and its normal derivative on the circle are obtained by linear interpolation. For the presented numerical results, the initial pressure (4.1) is given on $280\times 280$ Cartesian grid points in $[-1,1]^{2}$ and the data $\mathbf{M}_{c_{1},c_{2}}f$ are simulated for $\mathtt{{N_{\theta}}}=300$ equidistant sensor locations on unit circle. At each sensor we use $\mathtt{N}_{t}=1600$ equidistant time samples in the measurement interval $[0,6]$.

The top left image in Figure 4.1 shows the initial pressure distribution $f$ used our simulations. The top right image shows the classical PAT data $\mathbf{M}_{1,0}f+\xi_{1,0}$, and the bottom right image shows the normal derivative $\mathbf{M}_{0,1}f+\xi_{0,1}$, which has been re-scaled such that $\mathbf{M}_{1,0}f$ and $\mathbf{M}_{0,1}f$ have the same $\ell^{2}$-norm. The bottom left image shows the mixed data $\mathbf{M}_{1,1}f+\xi_{1,1}$. In all cases we have added Gaussian white noise $\xi_{c_{1},c_{2}}$ with a standard deviation of $50\text{\,}\mathrm{\char 37\relax}$ of the $\ell^{2}$-norm $\lVert\mathbf{M}_{c_{1},c_{2}}f\rVert_{2}:=(\sum_{\mathtt{\mathbf{m}},\mathtt{\mathbf{n}}}\lvert\mathbf{M}_{c_{1},c_{2}}f[\mathtt{\mathbf{m}},\mathtt{\mathbf{n}}]\rvert^{2}/(\mathtt{{N_{\theta}}}\mathtt{N}_{t}))^{1/2}$, resulting in a relative $\ell^{2}$-data error $\lVert\xi_{c_{1},c_{2}}\rVert_{2}/\lVert\mathbf{M}_{c_{1},c_{2}}f\rVert_{2}$ of $45\text{\,}\mathrm{\char 37\relax}$ in all cases. Note that we simulated the data until $T=6$ but only show data until $t=2.2$ in Figure 4.1, because after $T=2$ the data $\mathbf{M}_{c_{1},c_{2}}f$ are smooth and monotonically decreasing.

4.2 Implementation of the inversion formulas

In order to reconstruct the initial pressure distribution, we implement discrete versions of the reconstruction formulas (3.1) and (3.2), which for $n=2$ and $R=c_{1}=c_{2}=1$ are given by

[TABLE]

Here $g_{k}(t)=\frac{1}{\sqrt{2\pi}}\int_{0}^{2\pi}g(\theta(\varphi),t)e^{-ik\varphi}\mathrm{d}\varphi$ are the Fourier coefficients in the angular variable, and $\mathbf{x}=\rho\theta(\varphi)$ with $\theta(\varphi)=(\cos\varphi,\sin\varphi)$ and $\varphi\in[0,2\pi)$. Moreover, $\omega_{j,k}$ are the positive roots of $J_{k}$. Formulas (4.1) and (4.2) are implemented following [11] where the inversion formula (4.2) for standard PAT data are considered.

Below we briefly describe the discretization of inversion formula (4.1). The numerical reconstruction uses discrete samples $\mathtt{g}\left[\mathtt{\mathbf{m}},\mathtt{\mathbf{n}}\right]=g\left(\theta_{\mathtt{\mathbf{m}}},t_{\mathtt{\mathbf{n}}}\right)$ for $(\mathtt{\mathbf{m}},\mathtt{\mathbf{n}})\in\left\{1,\dots,\mathtt{{N_{\theta}}}\right\}\times\left\{1,\dots,\mathtt{N}_{t}\right\}$ where $\theta_{\mathtt{\mathbf{m}}}=2\pi\left(\mathtt{\mathbf{m}}-1\right)/\mathtt{{N_{\theta}}}$ and $\mathtt{\mathbf{n}}=T\left(\mathtt{\mathbf{n}}-1\right)/\mathtt{N}_{t}$. First, the Fourier coefficients $\mathtt{g}_{k}$ in the angular variable are computed by the FFT algorithm. Next, the cosine transform is approximated by $\mathtt{C}\{\mathtt{g}_{k}\}[j]:=\frac{T}{\mathtt{N}_{t}}\sum_{\mathtt{\mathbf{n}}=1}^{\mathtt{N}_{t}}\mathtt{g}_{k}[\mathtt{\mathbf{n}}]\cos(\omega_{j,k}t_{\mathtt{\mathbf{n}}})$. This is implemented by matrix-vector multiplication where the entries $\cos(\omega_{j,k}t_{\mathtt{\mathbf{n}}})$ are pre-computed and stored. Finally, we evaluate (4.1) by truncating both sums. Formula (4.2) is discretized in an analogous manner. This results in the discrete versions of (4.1), (4.2)

[TABLE]

In both formulas, the inner sums are evaluated by matrix-vector multiplications where the required matrix elements are pre-computed and stored. The outer sums are evaluated with the inverse FFT algorithms. Application of (4.3) (or (4.4)) with a standard Matlab implementation on a desktop PC with $16\text{\,}\mathrm{GB}$ RAM and $3.40\text{\,}\mathrm{MHz}$ eight-core processor with $\mathtt{{N_{\theta}}}=300$, $\mathtt{N}_{r}=180$ and $\mathtt{N}_{t}=1200$ takes about $1.2\text{\,}\sec$. Note that $\omega_{j,k}$ are the positive roots of $J_{k}$ and therefore $J_{k+1}(\omega_{j,k})\neq 0$ which implies that neither (4.3) nor (4.4) suffer from a division by zero problem.

4.3 Reconstruction results

Figure 4.2 shows reconstruction results for the simulated data without noise using both formulas (4.3), (4.4) applied to all three data cases $\mathbf{M}_{1,0}f$ (top), $\mathbf{M}_{1,1}f$ (middle) and $\mathbf{M}_{0,1}f$ (bottom). The left column show the reconstruction with $c_{2}=0$ formula and the right column the results with the $c_{2}\neq 0$ formula. Note that up to discretization error and truncation error at time $T=6$, Theorem 1 shows that (4.3) is exact for $\mathbf{M}_{1,0}f$, and that (4.3) is exact for $\mathbf{M}_{1,1}f$ and $\mathbf{M}_{0,1}f$. The numerical results shown in Figure 4.2 support these theoretical findings. Reconstruction results for noisy are shown in Figure 4.3, which again support exactness of inversion formula and further shows their stability with respect to Gaussian noise. A quantitative error analysis is shown in Figure 4.4, where the relative reconstruction error for cases is shown in dependence of the noise level of the data.

4.4 Discussion

First, note that the theoretically exact inversion formulas (3.1) and (3.2) use data for all times whereas the discrete counterparts only use data up to finite time. Nevertheless, the discrete formula (4.3) applied to $\mathbf{M}_{1,0}f$ and (4.4) applied to $\mathbf{M}_{0,1}f$ give visually satisfactory results. This is consistent with our previous observations [4, 7, 15]. However, the relative reconstruction error is quite large, even for data without noise. We mainly address this due to the required truncation of the data. In future work we will therefore address this issue and aim deriving exact inversion formula which only use data until a finite time $T$.x

Second, up to discretization and truncation errors, due to the range condition in Corollary 1, the inversion formula (4.3) applied to classical PAT data $\mathbf{M}_{1,0}f$ should yield the zero image. From Figures 4.2 (top right) and 4.3 (top right) we see that this is clearly not the case. Again, this is mainly due to the data truncation, which we have verified (result not shows) by varying $T$. The introduced violation of the range condition is also the reason that the inversion formula applied to $\mathbf{M}_{1,1}f$ yields worse results then the result for $\mathbf{M}_{0,1}f$. Finally, applying the inversion formula (3.2) to data $\mathbf{M}_{1,1}f$ results in amplified boundary structures. Visually the amplification is quite appealing which may be the reason that commonly data are modeled by $\mathbf{M}_{1,0}f$ instead of the more general model $\mathbf{M}_{c_{1},c_{2}}f$. Quantitative analyzing the effect of $\mathbf{M}_{0,1}f$ component to (3.2) is an interesting open issue.

5 Conclusion

We investigated PAT with the direction dependent data model (2.1), which uses linear combinations $\mathbf{M}_{c_{1},c_{2}}f=c_{1}\mathbf{M}_{1,0}f+c_{2}\mathbf{M}_{0,1}f$ of the acoustic pressure and its normal derivative. We developed an exact and stable reconstruction formula for the special case of spherical detection geometry. Numerical results show the validity of the proposed approach. Investigating such data for more general detection geometries in PAT is an interesting line of future research. Moreover, deriving an explicit inversion formula that only used data for $t\leq T<\infty$ is an important future aspect. Finally, in future work we will test our formulas on experimental data and investigate which values of $c_{1}$, $c_{2}$ actually are the best to accurately model experimental PAT data, for example using small piezoelectric sensors or piezoelectric line sensors.

Acknowledgments

The work of G.Z. and M.H. has been supported by the Austrian Science Fund (FWF), project P 30747-N32. The work of S.M. has been supported by the National Research Foundation of Korea grant funded by the Korea government (MSIP) (NRF-2018R1D1A3B07041149).