Photoacoustic tomography in a rectangular reflecting cavity

In conventional PAT, time domain measurements of the photoacoustically generated acoustic waves are made by an array of ultrasound detectors positioned on a measurement surface. The inverse problem of PAT/TAT consists in finding the initial pressure distribution from the measurements. Several idealizing assumptions are typically made to simplify this problem. First, the speed of sound is frequently assumed to be known and constant throughout the whole space, i.e. the presence of any physical boundaries is neglected. Second, it is assumed that the measurements are done at all points lying on some imaginary observation surface surrounding the object of interest. Third, the detectors are considered small and non-reflecting (acoustically transparent). Under these assumptions the acoustic waves can be considered as propagating in free space ($\mathbb {R} ^{n}$, n = 2 or 3) and the acoustic pressure p = p(t, x) is a solution to the following (free-space) initial value problem (IVP):

$$\begin{equation} \left( \frac{\partial ^{2}}{\partial t^{2}}-c_{0}^{2}\Delta \right) p =0,\quad x\in \mathbb {R}^{n},\quad t\in [0,\infty ), \end{equation} \tag{ 1 }$$

$$\begin{equation} p|_{t=0} =f(x),\quad \left. \frac{\partial p}{\partial t}\right|_{t=0}=0 \end{equation} \tag{ 2 }$$

where f(x) is the initial acoustic pressure distribution. Given time series measurements g of the acoustic pressure at each point of the observation surface S surrounding the region of interest Ω (within which f(x) is supported)

$$\begin{equation*} g=p(t,x)\quad x\in S,\quad t\in [0,T], \end{equation*}$$

it is possible to reconstruct f(x) exactly. Indeed, since S is not a physical boundary, for a sufficiently large time T, the pressure p in Ω will vanish. (In 3D, due to the Huygens principle, the pressure will vanish after T = diam(Ω)/c₀. In 2D, the pressure does not completely vanish in finite time; however, it will decay fast enough that it can be approximately set to 0 at a sufficiently large value of T ⩾ diam(Ω)/c₀.)

Now one can solve in Ω the following initial/boundary value problem (IBVP):

$$\begin{equation} \left( \frac{\partial ^{2}}{\partial t^{2}}-c_{0}^{2}\Delta \right) p =0,\quad t\in [0,T],\quad x\in \Omega , \end{equation} \tag{ 3 }$$

$$\begin{equation} p|_{t=T} =0,\quad \left. \frac{\partial p}{\partial t}\right|_{t=T}=0,\quad p|_{x\in S}=g(t,x), \end{equation} \tag{ 4 }$$

backwards in time (from t = T to t = 0), thus obtaining f(x) = p(0, x). This technique is called time reversal. It is theoretically exact and works for an arbitrary closed surface S; many of the existing reconstruction methods are based on this idea. For instance, one can solve IBVP (3), (4) directly using numerical techniques [10, 13, 30]. Eigenfunction decomposition techniques [1, 19] and some of the known inversion formulas [20] yield reconstructions that are theoretically equivalent to those obtained by time reversal. In each case the crucial underlying assumption is that the values of the pressure are recorded until the acoustic wave vanishes within Ω.

The methods of the previous section can be used so long as the detectors (or anything else such as the tank walls) do not reflect the acoustic waves. When an array of detectors is used, instead of a single scanned detector, the situation will change as the array itself will reflect the wave. In the case of a single planar array, the traditional reconstruction techniques still apply, since the reflected waves will propagate away from the detector and will not affect the measurements. However, the images reconstructed from measurements made using a single finite-sized planar array configuration will contain 'partial view' or 'limited data' artifacts because the acoustic waves traveling parallel (or almost parallel) to the surface are not measured. In order to improve the reconstruction one could either reflect those waves back on the detector array using passive reflectors [5], or add new (perpendicular) detector arrays. If the region of interest is surrounded by the arrays and/or reflectors, a reverberant cavity is formed. (Note that the free surface of a liquid also reflects waves). Wave propagation in a reverberant cavity is no longer represented by the IVP (1), and a different mathematical model is needed.

The proper model should take into account that the domain Ω is surrounded by a reflecting boundary ∂Ω, and the corresponding boundary conditions should be imposed on the wave equation (as opposed to the free-space propagation discussed in the previous section). The measurement surface S is now a subset of the boundary ∂Ω.

When the boundary is treated as sound-hard the acoustic pressure p(t, x) is a solution to the following IBVP:

$$\begin{equation} \left( \frac{\partial ^{2}}{\partial t^{2}}-c_{0}^{2}\Delta \right) p =0,\quad x\in \Omega \end{equation} \tag{ 5 }$$

$$\begin{equation} p(0,x) =f(x),\quad \frac{\partial p(0,x)}{\partial t}=0, \end{equation} \tag{ 6 }$$

$$\begin{equation} \frac{\partial }{\partial n}p(t,x) =0,\quad x\in \partial \Omega. \end{equation} \tag{ 7 }$$

Now the measurements can be written as

$$\begin{equation*} g=p(t,x)\quad x\in S\subseteq \partial \Omega ,\quad t\in [0,T]. \end{equation*}$$

Other boundary conditions might be appropriate in some circumstances, such as for instance if the cavity were built as a tank and one part of the boundary was a water–air interface where p ≈ 0, i.e. Dirichlet conditions are imposed on this part of the boundary. Also, it is not necessary to measure the pressure on the whole boundary ∂Ω, and so S is only a subset of ∂Ω. (Taking this to an extreme, Cox et al [6] proposed using a single measurement point when the cavity is ray-chaotic. This scheme, however, is unlikely to yield a stable reconstruction.)

A distinct property of this model is the preservation of the acoustic energy trapped within the cavity. Since the model assumes that there is no absorption of the acoustic energy by the medium, and the (Neumann or Dirichlet) boundary conditions correspond to a complete reflection of waves, the oscillations within Ω will (theoretically) continue forever. In practice, of course, this is not the case and the waves will soon decrease to the extent that further measurements will become impossible due to the low signal-to-noise ratio. While we are not explicitly modeling the absorption, we have to assume that the measurement time T is bounded. It will be chosen to correspond, roughly, to several bounces of the acoustic wave between the cavity walls.

The preservation of acoustic energy within the reverberating cavity makes it impossible to solve the inverse problem by the time reversal techniques mentioned in section 1.1. Indeed, in order to obtain the accurate reconstruction of f one has to accurately prescribe conditions p(T, x) and ∂p/∂t(T, x) to initialize the time reversal. However, there is no way to measure these data within the object. Moreover, these values are of the same order of magnitude as f, and if one simply replaces them by a zero (as we can safely do in the free-space case), the induced error will also be of the same magnitude as f we seek to reconstruct. Therefore, almost none of the known reconstruction algorithms are applicable here (with the exception of [2, 5, 6, 29]). Below we propose inversion techniques suitable for PAT within the reverberant cavity.

We assume that the acoustic pressure p(t, x), x = (x₁, x₂) solves the IBVP (5)–(7) in 2D, in a square Ω = [0, 1] × [0, 1]. Without loss of generality, the speed of sound c₀ will be assumed to equal 1. The measurements are represented by the pair of time series g = (g₁(t, x₂), g₂(t, x₁)) corresponding to the pressure values on the two adjacent sides of the square Ω (figure 1):

$$\begin{equation*} g_{1}(t,x_{2})= p(t,x)|_{x_{1}=0},\quad g_{2} (t,x_{1})= p(t,x)|_{x_{2}=0},\quad t\in [0,T]. \end{equation*}$$

Zoom In Zoom Out Reset image size

Let us denote by $\mathcal {W}$ the linear operator transforming the initial condition f(x) into the boundary data g:

$$\begin{equation*} \mathbf {g}=\mathcal {W}(f). \end{equation*}$$

Our ultimate goal is to reconstruct the initial pressure f(x) = p(0, x) from g, i.e. to invert $\mathcal {W}$.

The simple shape of the domain Ω allows us to utilize the eigenfunctions

$$\begin{equation*} \varphi _{k,l}(x)=\cos (\pi kx_{1})\cos (\pi lx_{2}) \end{equation*}$$

of the Neumann Laplacian in a square

$$\begin{equation*} (\Delta +\omega _{k,l}^{2})\varphi _{k,l}=0,\quad x\in \Omega ,\quad \frac{\partial \varphi _{k,l}}{\partial n}=0,\quad x\in \partial \Omega, \end{equation*}$$

with the eigenfrequencies ω_{k, l}:

$$\begin{equation*} \omega _{k,l}=\pi \sqrt{k^{2}+l^{2}},\quad k,l=0,1,2,\ldots. \end{equation*}$$

Using the eigenpairs (φ_{k, l}, ω_{k, l}), the solution p(t, x) of the forward IBVP (5)–(7) can be represented by the following series:

$$\begin{equation} p(t,x)=\sum _{k=0}^{\infty }\sum _{l=0}^{\infty }f_{k,l}\cos (\pi kx_{1})\cos (\pi lx_{2})\cos (\omega _{k,l}t),\quad x\in \Omega ,\quad t\in [0,T], \end{equation} \tag{ 8 }$$

where numbers f_{k, l} are the coefficients of the 2D Fourier cosine series of the initial pressure f(x):

$$\begin{equation} f_{k,l} =\frac{1}{||\varphi _{k,l}||_{2}^{2}}\int \nolimits _{\Omega }f(x)\varphi _{k,l}(x) \,{\rm d}x,\quad k,l=0,1,2,\ldots , \end{equation} \tag{ 9 }$$

$$\begin{equation} ||\varphi _{k,l}||_{2}^{2} =\int \nolimits _{\Omega }\varphi _{k,l} ^{2}(x) \,{\rm d}x=\left\lbrace \begin{array}{@{}@{}cc@{}@{}} \frac{1}{4}, & k\ne 0 \,{\rm and }\, l\ne 0,\\ 1, & k=0 \,{\rm and }\, l=0,\\ \frac{1}{2}, & {\rm otherwise.} \end{array} \right. \end{equation} \tag{ 10 }$$

For future use we will denote the double-indexed sequence of the Fourier coefficients by $F=\lbrace f_{k,l}\rbrace _{k,l=0}^{\infty },$ and will use the notation $\mathcal {F}_{{series}}$ to refer to the transformation defined by equations (9) and (10), so that

$$\begin{equation*} F=\mathcal {F}_{{series}}f. \end{equation*}$$

Obviously, in order to obtain f(x) it is enough to reconstruct coefficients F; function f is then computed by inverting $\mathcal {F} _{{series}},$ i.e. by summing the Fourier cosine series:

$$\begin{equation*} f(x)=\mathcal {F}_{{series}}^{-1}F=\sum _{k=0}^{\infty }\sum _{l=0}^{\infty }f_{k,l}\cos (\pi kx_{1})\cos (\pi lx_{2}). \end{equation*}$$

Let us first analyze the connection between coefficients f_{k, l} and the data measured on one side of the square (say, g₂):

$$\begin{equation} g_{2}(t,x_{1})\equiv p(t,x)|_{x_{2}=0}=\sum _{k=0}^{\infty }\sum _{l=0}^{\infty }f_{k,l}\cos (\pi kx_{1})\cos (\omega _{k,l}t),\quad t\in [0,T]. \end{equation} \tag{ 11 }$$

Due to the orthogonality of the cosine functions (in variable x₁) the above equation splits:

$$\begin{equation} \sum _{l=0}^{\infty }f_{k,l}\cos (\omega _{k,l}t) =g_{2,k}(t),\quad t\in [0,T],\quad k=0,1,2,\ldots , \end{equation} \tag{ 12 }$$

$$\begin{equation} g_{2,0}(t) \equiv \int \nolimits _{0}^{1}g_{2}(x_{1},t)\,{\rm d}x_{1},\quad \! g_{2,k}(t)\equiv 2\int \nolimits _{0}^{1}g_{2}(x_{2},t)\cos (\pi kx_{1})\,{\rm d}x_{1},\!\quad k=1,2,3,\ldots . \end{equation} \tag{ 13 }$$

It follows from equation (12) that each of the k functions g_{2, k}(t) is a combination of cosine functions in t with known frequencies; the values f_{k, l} can be viewed as generalized Fourier coefficients. More precisely, equation (12) can be re-written as the inverse Fourier transform $\mathcal {F}^{-1}$ of a sequence of Dirac delta-functions:

$$\begin{equation*} g_{2,k}(t)=\sqrt{2\pi }\mathcal {F}^{-1}\left[ \sum _{l=0}^{\infty }f_{k,l} \frac{\delta (\xi -\omega _{k,l})+\delta (\xi +\omega _{k,l})}{2}\right], \end{equation*}$$

where the forward and inverse Fourier transforms of some function h(t) are defined as follows:

$$\begin{eqnarray*} \hat{h}(\xi ) & \equiv ( \mathcal {F}h) (\xi )\equiv \frac{1}{\sqrt{2\pi }}\int \nolimits _{\mathbb {R}}h(t)\,{\rm e}^{-{\rm i}\xi t}\,{\rm d}t,\\ h(t) & \equiv ( \mathcal {F}^{-1}\hat{h}) (t)\equiv \frac{1}{\sqrt{2\pi }}\int \nolimits _{\mathbb {R}}\hat{h}(\xi )\,{\rm e}^{{\rm i}\xi t}\,{\rm d}\xi . \end{eqnarray*}$$

One can conclude that the Fourier transform $\hat{g}_{2,k}(\xi )$ of g_{2, k}(t) is the following distribution:

$$\begin{equation} \hat{g}_{2,k}(\xi )=\sqrt{2\pi }\sum _{l=0}^{\infty }f_{k,l}\frac{\delta (\xi -\omega _{k,l})+\delta (\xi +\omega _{k,l})}{2}. \end{equation} \tag{ 14 }$$

It is tempting to try to recover coefficients f_{k, l} by applying the Fourier transform $\mathcal {F}$ to each g_{2, k}(t) and thus obtaining $\hat{g}_{2,k}(\xi )$. However, functions g_{2, k}(t) do not vanish at infinity, and they (in general) are not periodic on any finite interval [0, T]. Therefore, such a computation would be quite inaccurate, at best.

A technique well known in digital signal processing as a way to Fourier transform long time series is to multiply the signal by a window function. For convenience, let us extend g_{2, k}(t) as an even function to the interval [ − T, T] (formulas (11) and (12) will remain valid under such extension). Consider now an even, infinitely smooth function η(t) vanishing with all its derivatives at −1 and 1, and its scaled version η_T(t) = η(t/T). Since the product η_T(t)g_{2, k}(t) is periodic on [ − T, T], its Fourier transform

$$\begin{equation} \widehat{\eta _{T}g_{2,k}}(\xi )\equiv \frac{1}{\sqrt{2\pi }}\int \nolimits _{-T} ^{T}\eta _{T}(t)g_{2,k}(t)\,{\rm e}^{-{\rm i}\xi t}\,{\rm d}t, \end{equation} \tag{ 15 }$$

can be easily computed (for example, by applying the composite trapezoid rule to the discretized signal). Now, by the convolution theorem

$$\begin{eqnarray*} &&\widehat{\eta _{T}g_{2,k}}(\xi )=\frac{1}{\sqrt{2\pi }}\widehat{\eta _{T}} (\xi )\ast \widehat{g_{2,k}}(\xi )=\frac{1}{\sqrt{2\pi }}\int \nolimits _{\mathbb {R} }\widehat{\eta _{T}}(\zeta -\xi )\widehat{g_{2,k}}(\zeta )\,{\rm d}\zeta , \end{eqnarray*}$$

and, taking into account the equality $\widehat{\eta _{T}}(\xi )=T\widehat{\eta }(T\xi )$ and equation (14) one obtains

$$\begin{eqnarray*} \widehat{\eta _{T}g_{2,k}}(\xi )&=&\sum _{m=0}^{\infty }f_{k,m}\frac{\widehat{\eta _{T}}(\xi -\omega _{k,m})+\widehat{\eta _{T}}(\xi +\omega _{k,m})}{2}\\ &=&T\sum _{m=0}^{\infty }f_{k,m}\frac{\widehat{\eta }(T(\xi -\omega _{k,m}))+\widehat{\eta }(T(\xi +\omega _{k,m}))}{2}. \end{eqnarray*}$$

Function $\widehat{\eta _{T}}(\xi )$ is infinitely smooth (as a Fourier transform of a finitely supported function), and vanishes at infinity faster than any rational function in ξ (since η(t) is infinitely smooth). One may view the problem of reconstruction of $\hat{g}_{2,k}(\xi )$ from $\widehat{\eta _{T}}(\xi )\ast \widehat{g_{2,k}}(\xi )$ as a deconvolution problem. Due to the smoothness of the convolution kernel $\widehat{\eta _{T} }(\xi )$, such a deconvolution is an extremely ill-posed problem.

Fortunately, the problem at hand can be solved without resorting to standard deconvolution techniques. Since we know a priori that the $\hat{g}_{2,k}(\xi )$ is a combination of delta-functions (with unknown coefficients f_{k, l}) whose positions in the frequency space are known (see equation (14)), the problem becomes much simpler. Let us consider the values of $\widehat{\eta _{T}g_{2,k}}$ only at the frequencies ω_{k, l}. Then, for each fixed k we obtain the following system of linear equations with respect to unknown f_{k, m}:

$$\begin{eqnarray} &&T\sum _{m=0}^{\infty }f_{k,m}\frac{\widehat{\eta }(T(\omega _{k,l}-\omega _{k,m}))+\widehat{\eta }(T(\omega _{k,l}+\omega _{k,m}))}{2}=\widehat{\eta _{T}g_{2,k}}(\omega _{k,l}),\quad l=0,1,2,\ldots . \nonumber\\ \end{eqnarray} \tag{ 16 }$$

Further, by separating the 'diagonal' terms this system can be re-written in the form

$$\begin{eqnarray} &&f_{k,l}+\frac{\widehat{\eta }(2T\omega _{k,l})}{\widehat{\eta }(0)}f_{k,l} +\sum _{\underset{m\ne l}{m=0}}^{\infty }\frac{\widehat{\eta }(T(\omega _{k,l}-\omega _{k,m}))+\widehat{\eta }(T(\omega _{k,l}+\omega _{k,m} ))}{\widehat{\eta }(0)}f_{k,m}=\frac{2}{T}\frac{\widehat{\eta _{T}g_{2,k} }(\omega _{k,l})}{\widehat{\eta }(0)}\nonumber\\ \end{eqnarray} \tag{ 17 }$$

for l = 0, 1, 2, ..., k = 0, 1, 2, ..., except the case when k = l = 0 which assumes a simple form

$$\begin{equation} f_{0,0}+\sum _{\underset{m\ne l}{m=0}}^{\infty }\frac{\widehat{\eta } (T\omega _{0,m})}{\widehat{\eta }(0)}f_{0,m}=\frac{1}{T}\frac{\widehat{\eta _{T}g_{2,0}}(\omega _{0,0})}{\widehat{\eta }(0)}. \end{equation} \tag{ 18 }$$

Let us discuss the properties of these systems (one system for every value of k). Due to the fast decrease of $\widehat{\eta _{T}}(\xi )$ at infinity, as T goes to infinity, values of $\widehat{\eta }(T(\omega _{k,l}-\omega _{k,m}))$ and $\widehat{\eta }(T(\omega _{k,l}+\omega _{k,m}))$ converge to zero, and the system becomes diagonal. This suggests that for sufficiently large values of T one can try to approximately solve the problem by the formula

$$\begin{equation} f_{k,l} \approx \frac{2}{T}\frac{\widehat{\eta _{T}g_{2,k}} (\omega _{k,l})}{\widehat{\eta }(0)+\widehat{\eta }(2T\omega _{k,l})},\quad k,l=0,1,2,\ldots , \end{equation} \tag{ 19 }$$

$$\begin{equation} f_{0,0} \approx \frac{1}{T}\frac{\widehat{\eta _{T}g_{2,0}} (\omega _{0,0})}{\widehat{\eta }(0)} \end{equation} \tag{ 20 }$$

and thus to reconstruct p₀(x) from the measurements made on one side of the square Ω.

Examples of such approximate reconstructions can be found in [7]; it is clear that the images are distorted by noticeable artifacts. A closer look at the functions $\widehat{\eta }(T(\omega _{k,l}-\omega _{k,m}))$ reveals the reason for these distortions. The convergence of these functions to zero depends on the values of the difference ω_{k, l} − ω_{k, m} which is not uniform with respect to k and l. In particular, for large values of k and l = 0 and m = 1, this difference can be approximately estimated as follows:

$$\begin{equation*} \omega _{k,1}-\omega _{k,0}=\pi ( \sqrt{k^{2}+1}-\sqrt{k^{2}}) =\pi k( \sqrt{1+1/k^{2}}-1t) \approx \frac{\pi }{2k}, \end{equation*}$$

so that the value of T(ω_{k, l} − ω_{k, m}) is not large for large values of k. This implies that for any (large) value of T there are some values of k, l and m for which off-diagonal terms are not small, and approximation (19) is not accurate. Moreover, one may suspect (and numerical simulations confirm this) that the equations with large values of k and small values of l are a source of instability, so that attempts to solve system (17) numerically (for large k) lead to a significant amplification of the noise present in the data.

As shown below, one can eliminate this instability by using not all equations (17) but only those for l = k, k + 1, k + 2, ... . The missing information can be obtained from the data g₁(x₂, t) measured on the second side of the square corresponding to x₁ = 0. Define functions g_{1, l}(t) by the equations

$$\begin{equation*} g_{1,l}(t)\equiv 2\int \nolimits _{0}^{1}g_{1}(x_{2},t)\cos (\pi lx_{2})\,{\rm d}x_{2}. \end{equation*}$$

The derivation similar to the one in the previous section yields the following system of equations:

$$\begin{eqnarray} &&f_{k,l}+\widehat{\eta }(2T\omega _{k,l})f_{k,l}+\sum _{\underset{m\ne k}{m=0} }^{\infty }\frac{\widehat{\eta }(T(\omega _{k,l}-\omega _{m,l}))+\widehat{\eta }(T(\omega _{k,l}+\omega _{m,l}))}{\widehat{\eta }(0)}f_{k,m}=\frac{2}{T} \frac{\widehat{\eta _{T}g_{1,l}}(\omega _{k,l})}{\widehat{\eta }(0)}, \nonumber\\ \end{eqnarray} \tag{ 21 }$$

k, l = 0, 1, 2, ... (with a simpler formula for the case k = l = 0 which we will omit). We, however, are going to use only half of these equations, namely those corresponding to k = l + 1, l + 2, l + 3, ... in combination with a half of the equations (17), arriving at the following set of equations:

$$\begin{equation} \left\lbrace \begin{array}{@{}@{}ll@{}} f_{0,0}+\sum _{\underset{m\ne l}{m=0}}^{\infty }\frac{\widehat{\eta } (T(\omega _{0,m}))}{\widehat{\eta }(0)}f_{0,m}=\frac{1}{T}\frac{\widehat{\eta _{T}g_{2,0}}(\omega _{0,0})}{\widehat{\eta }(0)},\quad l=k=0,\\ f_{k,l}+\widehat{\eta }(2T\omega _{k,l})f_{k,l}+\sum _{\underset{m\ne l}{m=0} }^{\infty }\frac{\widehat{\eta }(T(\omega _{k,l}-\omega _{k,m}))+\widehat{\eta }(T(\omega _{k,l}+\omega _{k,m}))}{\widehat{\eta }(0)}f_{k,m}=\frac{2}{T} \frac{\widehat{\eta _{T}g_{2,k}}(\omega _{k,l})}{\widehat{\eta }(0)},\quad l\ge k,\\ f_{k,l}+\widehat{\eta }(2T\omega _{k,l})f_{k,l}+\sum _{\underset{m\ne k}{m=0} }^{\infty }\frac{\widehat{\eta }(T(\omega _{k,l}-\omega _{m,l}))+\widehat{\eta }(T(\omega _{k,l}+\omega _{m,l}))}{\widehat{\eta }(0)}f_{k,m}=\frac{2}{T} \frac{\widehat{\eta _{T}g_{1,l}}(\omega _{k,l})}{\widehat{\eta }(0)},\quad k>l. \end{array} \right. \end{equation} \tag{ 22 }$$

Since for very large values of T factors $\widehat{\eta }(T(\cdots ))$ vanish (while $\widehat{\eta }(0)$ is a constant), one can compute a crude approximation f⁽⁰⁾(x) by defining the coefficients F⁽⁰⁾ as follows:

$$\begin{equation} f_{k,l}^{(0)} =\left\lbrace \begin{array}{@{}@{}ll@{}} \displaystyle\frac{1}{T}\frac{\widehat{\eta _{T}g_{2,0}}(\omega _{0,0})}{\widehat{\eta } (0)}, & k=l=0,\\ \displaystyle\frac{2}{T}\frac{\widehat{\eta _{T}g_{2,k}}(\omega _{k,l})}{\widehat{\eta } (0)}, & l\ge k\\ \displaystyle\frac{2}{T}\frac{\widehat{\eta _{T}g_{1,l}}(\omega _{k,l})}{\widehat{\eta } (0)}, & k>l \end{array} \right. , \end{equation} \tag{ 23 }$$

$$\begin{equation} F^{(0)} =\big\lbrace f_{k,l}^{(0)}\big\rbrace _{k,l=0}^{\infty }, \end{equation} \tag{ 24 }$$

and by summing the Fourier series

$$\begin{equation} f^{(0)}(x)=\big[ \mathcal {F}_{{series}}^{-1}F^{(0)}\big] (x). \end{equation} \tag{ 25 }$$

Formulas (23) and (24) define a linear operator $\mathcal {R}$ that maps boundary values g into the set of Fourier coefficients F, so that

$$\begin{equation*} f^{(0)}=\mathcal {F}_{{series}}^{-1}F^{(0)}=\mathcal {F} _{{series}}^{-1}\mathcal {R}\mathbf {g=}\mathcal {F} _{{series}}^{-1}\mathcal {RW}f. \end{equation*}$$

As our numerical experiments show, approximation f⁽⁰⁾ yields qualitatively correct images even for moderate values of T. Our goal, however, is to obtain a quantitatively accurate, theoretically exact reconstruction. Forming and decomposing a matrix corresponding to the system (22) (after proper truncation of the spatial harmonics) is computationally prohibitive, since the number of unknowns in the 3D high-resolution images can reach hundreds of millions. Instead, we will solve this system iteratively.

Let us define an operator $\mathcal {A}$ that transforms a 2D sequence of numbers $H=\lbrace h_{k,m}\rbrace _{k,m=0}^{\infty }$ into a 2D sequence of numbers $E=\lbrace e_{k,m}\rbrace _{k,m=0}^{\infty }$ by the formulas

$$\begin{eqnarray} &&e_{0,0}=\sum _{\underset{m\ne l}{m=0}}^{\infty }\frac{\widehat{\eta } (T(\omega _{0,m}))}{\widehat{\eta }(0)}h_{0,m},\quad l=k=0,\nonumber\\ &&e_{k,l}=\frac{\widehat{\eta }(2T\omega _{k,l})}{\widehat{\eta }(0)}h_{k,l} +\sum _{\underset{m\ne l}{m=0}}^{\infty }\frac{\widehat{\eta }(T(\omega _{k,l}-\omega _{k,m}))+\widehat{\eta }(T(\omega _{k,l}+\omega _{k,m} ))}{\widehat{\eta }(0)}h_{k,m},\quad l\ge k,\nonumber\\ &&e_{l,k}=\frac{\widehat{\eta }(2T\omega _{k,l})}{\widehat{\eta }(0)}h_{k,l} +\sum _{\underset{m\ne k}{m=0}}^{\infty }\frac{\widehat{\eta }(T(\omega _{k,l}-\omega _{m,l}))+\widehat{\eta }(T(\omega _{k,l}+\omega _{m,l} ))}{\widehat{\eta }(0)}h_{k,m},\quad k>l, \end{eqnarray} \tag{ 26 }$$

or, in the operator notation,

$$\begin{equation*} E=\mathcal {A}H. \end{equation*}$$

Then equations (22) can be re-written as

$$\begin{equation} (\mathcal {I+A)}F=F^{(0)}, \end{equation} \tag{ 27 }$$

where F is the set of the Fourier coefficients of the sought f(x), and F⁽⁰⁾ are the Fourier coefficients of the initial, crude approximation f⁽⁰⁾(x). If the spectral radius of operator $\mathcal {A}$ is less than 1, the unique solution of equation (27) can be found in the form of a converging Neumann series

$$\begin{equation} F=\sum \limits _{K=0}^{\infty }(-\mathcal {A})^{K}F^{(0)}. \end{equation} \tag{ 28 }$$

Function f(x) is then found from F by summing the Fourier series (see (25)). Alternatively, this solution (F) can be represented as the limit of iterations defined by the recurrence relation

$$\begin{equation} F^{(K)}=F^{(0)}-\mathcal {A}F^{(K-1)},\quad K=1,2,3,\ldots . \end{equation} \tag{ 29 }$$

In section 2.4, we analyze operator $\mathcal {A}$ and find a sufficient condition for the convergence of the Neumann series (28).

In practical computations we are not evaluating series (28) directly, since this is still expensive. By summing the Fourier series (i.e. by applying the $\mathcal {F}_{{series}}^{-1}$ operator), iterations (29) can be replaced by the equivalent relation

$$\begin{equation*} f^{(K)}(x)=f^{(0)}(x)-\mathcal {F}_{{series}}^{-1}\mathcal {AF} _{{series}}f^{(K-1)},\quad K=1,2,3,\ldots \end{equation*}$$

or

$$\begin{equation} f^{(K)}(x)=\mathcal {F}_{{series}}^{-1}\mathcal {R}\mathbf {g} -\mathcal {F}_{{series}}^{-1}\mathcal {AF}_{{series} }f^{(K-1)},\quad K=1,2,3,\ldots . \end{equation} \tag{ 30 }$$

On the other hand, it follows from (27) that

$$\begin{equation*} (\mathcal {I+A)F}_{{series}}f=\mathcal {RW}f, \end{equation*}$$

and, since (27) holds for any initial condition f(x), the following operator identity holds:

$$\begin{equation*} \mathcal {F}_{{series}}^{-1}(\mathcal {I+A)F} _{{series}}=\mathcal {F}_{{series}}^{-1} \mathcal {RW}. \end{equation*}$$

Thus,

$$\begin{equation*} -\mathcal {F}_{{series}}^{-1}\mathcal {AF}_{{series} }=\mathcal {I-F}_{{series}}^{-1}\mathcal {RW}, \end{equation*}$$

and, therefore, recurrence relation (30) is equivalent to

$$\begin{equation*} f^{(K)}=\mathcal {F}_{{series}}^{-1}\mathcal {R}\mathbf {g+} f^{(K-1)}-\mathcal {F}_{{series}}^{-1}\mathcal {RW}f^{(K-1)},\quad K=1,2,3,\ldots \end{equation*}$$

or to

$$\begin{equation} f^{(K)} =f^{(K-1)}+\mathcal {F}_{{series}}^{-1}\mathcal {R(} \mathbf {g}-\mathcal {W}f^{(K-1)}),\quad K=1,2,3,\ldots , \end{equation} \tag{ 31 }$$

$$\begin{equation} f^{(0)} =\mathcal {F}_{{series}}^{-1}\mathcal {R}\mathbf {g.} \end{equation} \tag{ 32 }$$

Our computational algorithm is based on equations (31) and (32); these iterations are theoretically equivalent to those defined by equation (29) and have the same convergence properties. The computational advantage in using (31) and (32) lies in the possibility of easily implementing all the necessary operations using fast transforms (fast Fourier transforms (FFTs)). The implementation details are discussed in section 3. In the next section, we analyze the convergence of such iterations, and form (29) is more convenient for this purpose.

Let us consider the || · ||_∞ norm defined on the space of double-indexed infinite sequences $H=\lbrace h_{k,l}\rbrace _{k,l=0}^{\infty }$:

$$\begin{equation*} ||H||_{\infty }=\sup _{k,l=0.1,2,3, \ldots }|h_{k,l}|, \end{equation*}$$

and a space $\mathbb {L}$ of such sequences with bounded || · ||_∞ norm. Further, for a linear operator $\mathcal {B}:\mathbb {L\rightarrow L}$ define the induced ∞-norm:

$$\begin{equation*} ||\mathcal {B}||_{\infty }=\sup _{||H||_{\infty }\ne 0}\frac{||\mathcal {B} (H)||_{\infty }}{||H||_{\infty }}. \end{equation*}$$

Our goal is to estimate the induced ∞-norm of the operator $\mathcal {A}$ defined by equations (26) in the previous sections. These equations contain factors in the form $\widehat{\eta }(T(\omega _{k,l}-\omega _{k,m}))$. Note that for a sufficiently smooth cut-off function η(t) its Fourier transform declines fast at infinity:

$$\begin{equation} |\widehat{\eta }(\xi )|\le \frac{B}{1+|\xi |^{2}}. \end{equation} \tag{ 33 }$$

Let us find a lower bound on |ω_{k, l} − ω_{k, m}|. Consider first the case l ⩾ k:

$$\begin{equation*} |\omega _{k,l}-\omega _{k,m}|=\pi |\sqrt{k^{2}+l^{2}}-\sqrt{k^{2} +m^{2}}|=\frac{\pi |l-m|(l+m)}{\sqrt{k^{2}+l^{2}}+\sqrt{k^{2}+m^{2} }}. \end{equation*}$$

If, in addition, m ⩾ k, the above equation yields

$$\begin{equation*} |\omega _{k,l}-\omega _{k,m}|\ge \pi \frac{|l-m|(l+m)}{\sqrt{2}(l+m)}=\frac{\pi |l-m|}{\sqrt{2}}. \end{equation*}$$

If m < k (and l ⩾ k),

$$\begin{equation*} |\omega _{k,l}-\omega _{k,m}|\ge \frac{\pi |l-m|(l+m)}{\sqrt{2}(l+k)}\ge \frac{\pi |l-m|(l+m)}{2\sqrt{2}l}\ge \frac{\pi |l-m|}{2\sqrt{2}}, \end{equation*}$$

so that the uniform bound holds if l ⩾ k for all values of m:

$$\begin{equation*} |\omega _{k,l}-\omega _{k,m}|\ge \frac{\pi |l-m|}{2\sqrt{2}}. \end{equation*}$$

Now one can bound $|\widehat{\eta }(T(\omega _{k,l}-\omega _{k,m}))|$, still under assumption l ⩾ k:

$$\begin{eqnarray} |\widehat{\eta }(T(\omega _{k,l}-\omega _{k,m}))| &\le & \frac{B}{1+T^{2}(\omega _{k,l}-\omega _{k,m})^{2}}\le \frac{B}{1+\pi ^{2}T^{2}(l-m)^{2}/8}\nonumber \\ &\le & \frac{B}{\pi ^{2}T^{2}(l-m)^{2}/8},{\rm if }\;m\ne l. \end{eqnarray} \tag{ 34 }$$

Similarly, if l > 0 or k > 0,

$$\begin{equation} |\widehat{\eta }(T(\omega _{k,l}+\omega _{k,m}))|\le \frac{B}{1+T^{2}(\omega _{k,l}+\omega _{k,m})^{2}}\le \frac{B}{1+T^{2}\pi ^{2}(1+m^{2})}\le \frac{B}{T^{2}\pi ^{2}(1+m^{2})} \end{equation} \tag{ 35 }$$

and

$$\begin{equation} |\widehat{\eta }(2T\omega _{k,l})|\le \frac{B}{T^{2}\pi ^{2}}. \end{equation} \tag{ 36 }$$

Let us now estimate the ∞-norm of the vector $G=\mathcal {A}H$ resulting from applying operator $\mathcal {A}$ to a sequence H whose ∞-norm equals 1, so that

$$\begin{equation} |h_{k,l}|\le 1,\quad l,k=0,1,2,\ldots {\rm .} \end{equation} \tag{ 37 }$$

First, we use the second equation in (26) and bound e_{k, l} with l ⩾ k (excluding the case k = l = 0). Taking into account inequalities (34)–(37), one obtains

$$\begin{eqnarray*} |e_{k,l}| & \le \frac{B}{T^{2}\pi ^{2}}\left( 1+\sum _{\underset{m\ne l}{m=0}}^{\infty }\frac{1}{(1+m^{2})}+\frac{8}{(l-m)^{2}}\right) \\ & \le \frac{B}{T^{2}\pi ^{2}}\left( 1+\sum _{m=0}^{\infty }\frac{1}{(1+m^{2} )}+16\sum _{m=1}^{\infty }\frac{1}{m^{2}}\right) \\ & \le \frac{B}{T^{2}\pi ^{2}}\left( 2+17\sum _{m=1}^{\infty }\frac{1}{m^{2} }\right) . \end{eqnarray*}$$

The value of the latter series is well known (see [11], equation 0.233)

$$\begin{equation*} \sum _{m=1}^{\infty }\frac{1}{m^{2}}=\frac{\pi }{6}, \end{equation*}$$

resulting in the estimate

$$\begin{equation} |e_{k,l}|\le \frac{B}{T^{2}\pi ^{2}}\frac{12+17\pi ^{2}}{6}. \end{equation} \tag{ 38 }$$

The above computations can be replicated with very minor changes to show that (38) holds also in the cases k = l = 0 and k > l (corresponding to the first and third lines of (26)) which proves the following.

Lemma 1. Operator $\mathcal {A}$ is bounded in the induced ∞-norm by

$$\begin{equation*} ||\mathcal {A}||_{\infty }\le \frac{B}{T^{2}\pi ^{2}}\frac{12+17\pi ^{2}}{6}. \end{equation*}$$

The lemma implies that if the acquisition time T is large enough, e.g. if

$$\begin{equation} T>\frac{1}{\pi }\sqrt{\frac{B(12+17\pi ^{2})}{6}}, \end{equation} \tag{ 39 }$$

then the operator $\mathcal {A}$ is a contraction mapping in the induced ∞-norm, i.e. $||\mathcal {A}||_{\infty }<1$. In this case, the standard theory of contraction mappings applies and one arrives at the following.

Theorem 2. For a sufficiently large acquisition time T (satisfying (39)) the Neumann series (28) converges in the ∞-norm, implying convergence of iterations (31) and (32) in the following sense:

$$\begin{equation*} \sup _{k,l=0.1,2,3, \ldots }\big|f_{k,l}^{(K)}-f_{k,l}\big|\underset{K\rightarrow \infty }{\rightarrow }0. \end{equation*}$$

Remark. One may want to have an explicit estimate on the constant B in equations (33) and (39). Such an estimate would depend on function η(t). Although we assumed, for simplicity, that η(t) is infinitely smooth, the proof of the above theorem requires only relatively slow decay of $|\widehat{\eta }(\xi )|$. Therefore, a less smooth function can be used, whose Fourier transform is easily computed and analyzed. For example, if one chooses η(t) to be equal to cos ²(πt/2) for |t| < 1 and 0 otherwise, then

$$\begin{equation*} \widehat{\eta }(\xi )=\frac{\pi ^{2}}{\sqrt{2\pi }}\frac{\sin (\xi )}{\xi (\pi +\xi )(\pi -\xi )}. \end{equation*}$$

Clearly, for ξ ⩾ π + 1,

$$\begin{equation*} \left|\frac{\sin (\xi )}{\xi (\pi +\xi )(\pi -\xi )}\right|\le \frac{1}{\xi (\pi +\xi )}\le \frac{1}{\xi ^{2}+1}, \end{equation*}$$

so that

$$\begin{equation*} |\widehat{\eta }(\xi )|\le \frac{\pi ^{2}}{\sqrt{2\pi }}\frac{1}{\xi ^{2}+1}. \end{equation*}$$

Due to the evenness of $|\widehat{\eta }(\xi )|$ the same inequality holds for ξ ⩽ −(π + 1). In the interval |ξ| ⩽ π + 1 this inequality also holds, as can be easily verified numerically. Therefore, for the present choice of η(t), equations (33) and (39) hold with $B=\frac{\pi ^{2} }{\sqrt{2\pi }}.$

In the 3D case, the problem and the proposed algorithm are very similar to their 2D counterparts. Since the 3D case is both more important for the practical use and more challenging computationally, we provide in this section a brief outline of the derivation and convergence analysis of our technique.

The acoustic pressure p(t, x), x = (x₁, x₂, x₃) solves IBVP (5)–(7) in 3D, in a cube Ω = [0, 1] × [0, 1] × [0, 1]. As before, the speed of sound c₀ will be assumed to equal 1. The measurements are represented by the triple of time series g = (g₁(t, x₂), g₂(t, x), g₃(t, x)) corresponding to the pressure values on the three adjacent sides of the cube Ω:

$$\begin{eqnarray*} &&g_{1}(t,x_{2},x_{3})=\left. p(t,x)\right|_{x_{1}=0},\ g_{2} (t,x_{1},x_{3})=\left. p(t,x)\right|_{x_{2}=0},\\ &&g_{3}(t,x_{1} ,x_{2})=\left. p(t,x)\right|_{x_{3}=0},\ t\in [0,T]. \end{eqnarray*}$$

The solution of IBVP in the cube can be represented using the Neumann Laplacian eigenfunctions

$$\begin{equation*} \varphi _{k,l,n}(x)=\cos (\pi kx_{1})\cos (\pi lx_{2})\cos (\pi nx_{3}) \end{equation*}$$

with eigenfrequencies

$$\begin{equation*} \omega _{k,l,n}=\pi \sqrt{k^{2}+l^{2}+n^{2}},\quad k,l,n=0,1,2,\ldots {\rm .} \end{equation*}$$

It has the following form:

$$\begin{equation*} p(t,x)=\sum _{k=0}^{\infty }\sum _{l=0}^{\infty }\sum _{n=0}^{\infty }f_{k,l,n} \cos (\pi kx_{1})\cos (\pi lx_{2})\cos (\pi nx_{3})\cos (\omega _{k,l,n}t),\quad x\in \Omega ,\quad t\in [0,T]. \end{equation*}$$

Coefficients $F=\lbrace f_{k,l,n}\rbrace _{k,l,n=0}^{\infty }$ are related to f(x) through the Fourier cosine series, $F=\mathcal {F}_{{series}}f$. They are computed as follows:

$$\begin{equation} f_{k,l,n}=\frac{1}{||\varphi _{k,l,n}||_{2}^{2}}\int_{\Omega } f(x)\varphi _{k,l,n}(x)\, {\rm d}x,\quad k,l,n=0,1,2,\ldots {\rm .} \end{equation} \tag{ 40 }$$

Conversely, f is obtained from F by summing the standard 3D cosine Fourier series, $f=\mathcal {F}_{{series}}^{-1}F$.

Consider the side of the cube corresponding to x₂ = 0. Then

$$\begin{eqnarray} g_{2}(t,x_{1},x_{3})&\equiv & p(t,x)|_{x_{2}=0}\nonumber\\ &=&\sum _{k=0}^{\infty }\sum _{n=0}^{\infty }\cos (\pi kx_{1})\cos (\pi nx_{3})\left[ \sum _{l=0}^{\infty }f_{k,l,n}\cos (\omega _{k,l,n}t)\right],\quad t\in [0,T]. \end{eqnarray} \tag{ 41 }$$

Due to the orthogonality of the products of cosines, the functions in brackets (denote them g_{2, k, n}) can be easily found from g₂ by the 2D cosine Fourier series:

$$\begin{equation} g_{2,k,n}(t)\equiv \sum _{l=0}^{\infty }f_{k,l,n}\cos (\omega _{k,l,n}t),\quad t\in [0,T],\quad k,n=0,1,2,\ldots {\rm .} \end{equation} \tag{ 42 }$$

By finding the values of the Fourier transform $\widehat{g_{2,k,n}\eta _{T}}$ of the product g_{2, k, n}(t)η_T(t) at the frequencies ω_{k, l, n} one obtains the following infinite system of equations:

$$\begin{eqnarray*} &&f_{k,l,n}+\frac{\widehat{\eta }(2T\omega _{k,l,n})}{\widehat{\eta }(0)} f_{k,l,n}+\sum _{\underset{m\ne l}{m=0}}^{\infty }\frac{\widehat{\eta } (T(\omega _{k,l,n}-\omega _{k,m,n}))+\widehat{\eta }(T(\omega _{k,l,n} +\omega _{k,m,n}))}{\widehat{\eta }(0)}f_{k,m,n}\nonumber\\ &&=\frac{2}{T}\frac{\widehat{\eta _{T}g_{2,k,n}}(\omega _{k,l,n})}{\widehat{\eta }(0)}, \end{eqnarray*}$$

for k, l, n = 0, 1, 2, ..., except the case when k = l = n = 0 which yields a simpler equation

$$\begin{equation*} f_{0,0,0}+\sum _{\underset{m\ne l}{m=0}}^{\infty }\frac{\widehat{\eta } (T\omega _{0,m,0})}{\widehat{\eta }(0)}f_{0,m,0}=\frac{1}{T}\frac{\widehat{\eta _{T}g_{2,0,0}}(\omega _{0,0})}{\widehat{\eta }(0)}. \end{equation*}$$

As in the 2D case, for a stable reconstruction we need to use equations from three sides (equations for the remaining two sides are obtained in a similar fashion). We thus arrive at the system

$$\begin{equation*} \left\lbrace \begin{array}{@{}llll@{}} f_{0,0,0}+\left[ \sum _{\underset{m\ne l}{m=0}}^{\infty }\frac{\widehat{\eta }(T\omega _{0,m,0})}{\widehat{\eta }(0)}f_{0,m,0}\right] =f_{0,0,0}^{(0)},\quad l=k=n=0, \\ f_{k,l,n}+\left[\! \frac{\widehat{\eta }(2T\omega _{k,l,n})}{\widehat{\eta } (0)}f_{k,l,n}+\sum _{\underset{m\ne l}{m=0}}^{\infty }\frac{\widehat{\eta }(T(\omega _{k,l,n}-\omega _{m,l,n}))+\widehat{\eta }(T(\omega _{k,l,n} +\omega _{m,l,n}))}{\widehat{\eta }(0)}f_{m,l,n}\!\right] =f_{k,l,n} ^{(0)},\,k\ge l,\,k\ge n, \\ f_{k,l,n}+\left[\! \frac{\widehat{\eta }(2T\omega _{k,l,n})}{\widehat{\eta } (0)}f_{k,l,n}+\sum _{\underset{m\ne l}{m=0}}^{\infty }\frac{\widehat{\eta }(T(\omega _{k,l,n}-\omega _{k,m,n}))+\widehat{\eta }(T(\omega _{k,l,n} +\omega _{k,m,n}))}{\widehat{\eta }(0)}f_{k,m,n}\!\right] =f_{k,l,n} ^{(0)},\,l>k,\,l\ge n, \\ f_{k,l,n}+\left[\! \frac{\widehat{\eta }(2T\omega _{k,l,n})}{\widehat{\eta } (0)}f_{k,l,n}+\sum _{\underset{m\ne l}{m=0}}^{\infty }\frac{\widehat{\eta }(T(\omega _{k,l,n}-\omega _{k,l,m}))+\widehat{\eta }(T(\omega _{k,l,n} +\omega _{k,l,m}))}{\widehat{\eta }(0)}f_{k,l,m}\!\right] =f_{k,l,n} ^{(0)},\,n>l,\,n>k, \end{array} \right. \end{equation*}$$

where the right-hand side corresponds to the Fourier coefficients of the crude first approximation f⁽⁰⁾(x):

$$\begin{equation} f_{k,l,n}^{(0)}=\left\lbrace \begin{array}{@{}ccc@{}} \frac{1}{T}\frac{\widehat{\eta _{T}g_{2,0,0}}(\omega _{0,0})}{\widehat{\eta } (0)}, & l=k=n=0,\\ \frac{2\widehat{\eta _{T}g_{1,l,n}}(\omega _{k,l,n})}{T\widehat{\eta }(0)}, & \,k\ge l,\,k\ge n,\\ \frac{2\widehat{\eta _{T}g_{2,k,n}}(\omega _{k,l,n})}{T\widehat{\eta }(0)}, & l>k,\,l\ge n,\\ \frac{2\widehat{\eta _{T}g_{3,k,l}}(\omega _{k,l,n})}{T\widehat{\eta }(0)}, & n>l,\,n>k, \end{array} \right. \end{equation} \tag{ 43 }$$

and where the terms in the brackets define the operator $\mathcal {A}$. Now the system is in the form (27) and its solution can be found in the form of the Neumann series (29) if the series converges. Under the latter condition, the solution to the original inverse problem f(x) can be obtained by the iterative process (31), (32), where the operator $\mathcal {R}$ is now defined by (43) and the operator $\mathcal {W}$ maps solutions of the IBVP (5)–(7) in 3D corresponding to the initial condition f(x) into the values of p(t, x) at the three sides of the cube Ω.

The convergence analysis is also very similar to that in the 2D case. For example, the difference |ω_{k, l, n} − ω_{k, m, n}| can be estimated as follows:

$$\begin{equation*} |\omega _{k,l,n}-\omega _{k,m,n}|=\pi |\sqrt{k^{2}+l^{2}+n^{2}} -\sqrt{k^{2}+m^{2}+n^{2}}|=\pi \frac{|l-m|(l+m)}{\sqrt{k^{2} +l^{2}+n^{2}}+\sqrt{k^{2}+m^{2}+n^{2}}}. \end{equation*}$$

If l ⩾ k and l ⩾ n and m ⩾ l, then

$$\begin{equation*} |\omega _{k,l,n}-\omega _{k,m,n}|\ge \pi \frac{|l-m|(l+m)}{\sqrt{3l^{2}} +\sqrt{3m^{2}}}\ge \frac{\pi |l-m|}{\sqrt{3}}. \end{equation*}$$

If l ⩾ k and l ⩾ n and m < l, then

$$\begin{equation*} |\omega _{k,l,n}-\omega _{k,m,n}|\ge \pi \frac{|l-m|(l+m)}{\sqrt{3l^{2}} +\sqrt{3l^{2}}}\ge \pi \frac{|l-m|}{2\sqrt{3}}. \end{equation*}$$

The latter estimate therefore holds for all values of m, if l ⩾ k and l ⩾ n. Proceeding the same way as in the 2D case one obtains the following estimate on the ∞-norm of the operator $\mathcal {A}$.

Lemma 3 (3D case). Operator $\mathcal {A}$ is bounded in the induced ∞-norm by

$$\begin{equation*} ||\mathcal {A}||_{\infty }\le \frac{B}{T^{2}\pi ^{2}}\frac{12+25\pi ^{2}}{6}. \end{equation*}$$

This immediately leads to a sufficient condition for the convergence of the 3D algorithm if the acquisition time T satisfies the condition

$$\begin{equation} T>\frac{1}{\pi }\sqrt{\frac{B(12+25\pi ^{2})}{6}}. \end{equation} \tag{ 44 }$$

Theorem 4 (3D case). For a sufficiently large acquisition time T (satisfying (44)) the Neumann series (28) converges in the ∞-norm, implying the convergence of iterations (31) and (32) in the following sense:

$$\begin{equation*} \sup _{k,l,n=0.1,2,3,\ldots }\big|f_{k,l,n}^{(K)}-f_{k,l,n}\big|\underset{K\rightarrow \infty }{\rightarrow }0. \end{equation*}$$

We remind the reader that the reconstruction algorithm is defined by equations (31) and (32). $\mathcal {W}$ is one of the operators in (31); it maps the initial condition p(0, x) = f(x) (here f(x) is an arbitrary function) into values of p(t, x) at the cube faces. In our implementation, computation for each face is done separately. For example, in order to find $g_{2}(t,x_{1},x_{3})\equiv p(t,x)|_{x_{2}=0}$, the following steps are done.

• (i)  
  Expand f(x) into the 3D Fourier cosine series, obtain coefficients f_{k, l, n} (equation (40)), k, l, n = 0, 1, 2, ..., N.
• (ii)  
  On a uniform grid in t, compute functions g_{2, k, n}(t), k, n = 0, 1, 2, ..., N, (equation (42)).
• (iii)  
  For each value of t, find g₂(t, x₁, x₃) by summing 2D cosine Fourier series (41).

Let us obtain a crude estimate of the number of operations involved in such a computation. For simplicity, in all simulations we kept the grid step in time and in space the same. Therefore, for moderate measurement times T, the number of time nodes is $\mathcal {O}(N)$. Step 1 is done by the 3D cosine FFT, it requires $\mathcal {O}(N^{3}\log N)$ flops.

Step 2 is slightly more complicated. Evaluation of (42) can be viewed as computing on the uniform grid the 1D Fourier transform of a function given on a non-uniform grid (frequencies ω_{k, l, n} are not equispaced!). This can be done fast by applying the 1D Nonequispaced FFT (NSFFT, see, for example [8]). This algorithm is not exact (unlike the regular FFT). However, any needed accuracy can be obtained, at the expense of increased constant factor implicit in the $\mathcal {O}(N\log N)$ estimate of the complexity of this technique. In our simulations, the accuracy of the NSFFT was of the order of 10⁻⁵. Since there is N × N of such transforms, the complexity of step 2 is $\mathcal {O} (N^{3}\log N)$ flops.

Step 3 is implemented by summing the 2D cosine series (using the corresponding type of FFT of size N × N, $\mathcal {O}(N^{2}\log N)$ flops each) for each of $\mathcal {O}(N)$ time nodes. This step is clearly done in $\mathcal {O}(N^{3}\log N)$ flops.

To summarize, the expense of finding g₂(t, x₁, x₃) by our method is $\mathcal {O}(N^{3}\log N)$ flops. The data for the other faces of the cube are computed in a similar fashion (step 1 is exactly the same and is performed only once).

We also used this algorithm to compute the simulated measurements g for our experiments. In some of the simulations, a significant level of noise was added to g to model the measurement errors and to check the stability of the algorithm to such errors.

The iterations given by equation (31) can be viewed as repeated applications of the forward operator $\mathcal {W}$ and approximate inverse, $\mathcal {F} _{{series}}^{-1}\mathcal {R}$, in alternating order. The summation of the Fourier series $\mathcal {F}_{{series}}^{-1}$ is done by the 3D cosine FFT algorithm; it takes $\mathcal {O}(N^{3}\log N)$ flops. Operator $\mathcal {R}$ is defined by equations (43); it consists in expanding the data for each face and for each time sample in the 2D cosine Fourier series (i.e. finding g_{1, l, n}(t), g_{2, k, n}(t) and g_{3, k, l}(t)), in multiplying these functions by the scaled cut-off function η_T, computing the 1D Fourier transform of the products at frequencies ω_{k, l, n} and multiplying the results by some constant factors. Expanding the data in the Fourier series takes $\mathcal {O}(N)\times \mathcal {O}(N^{2}\log N)=\mathcal {O}(N^{3}\log N)$ flops. The only remaining non-trivial operation here is computing the 1D Fourier transforms $\widehat{\eta _{T}g_{...}}$ of products η_Tg_... at non-equally spaced frequencies. In our simulation, we simply applied the regular FFT to compute $\widehat{\eta _{T}g_{...}}$ on a uniform grid, and interpolated using third degree Lagrange polynomials, to obtain the values at the frequencies ω_{k, l, n}. The computational expense of this step is $\mathcal {O} (N^{3}\log N)$ flops. A more sophisticated approach to finding these values is to use the proper version of the NSFFT [8]. This would yield more accurate results at the price of longer computing time, although the complexity would still remain $\mathcal {O}(N^{3}\log N)$ flops for this step. We found, however, that the simple polynomial interpolation described above is sufficiently accurate for the purposes of the present work.

To summarize, one iteration of the proposed method requires $\mathcal {O} (N^{3}\log N)$ flops, i.e. the algorithm is indeed asymptotically fast. The total number of iterations needed to attain the convergence was from 1 to 4 in our simulations, so, it is fair to say that the whole technique is implemented in $\mathcal {O}(N^{3}\log N)$ flops.

Performance of the proposed reconstruction algorithm was tested in numerical simulations. The preliminary results in 2D were reported in [7]. Below we present 3D simulations.

3.3.1. Measurements done on three faces

As a numerical phantom representing the initial pressure f(x) we used a linear combination of several slightly smoothed characteristic functions of balls of various radii, similar to the phantom utilized in [22]. The centers of the balls were chosen to lie on pair-wise intersections of the planes x_j = 0.25, j = 1, 2, 3. The cross-section of the phantom by the plane x₃ = 0.25 is shown in figures 2(a) and 5(a). (If needed, additional cross-sections can be found in [22].) The measurement time T in the two simulations described below was equal to 2, which corresponds to the sound wave covering the shortest distance between the parallel faces twice. For comparison, the measurement time in the standard problem of TAT/PAT in 3D free space for a cube of this size would equal $\sqrt{3}$.

Zoom In Zoom Out Reset image size

The size of the discretization grid was 401 × 401 × 401, i.e. the reconstructed images contained about 64 million of unknowns. For simplicity the time step was chosen to equal the spatial step (which was reasonable since the model speed of sound was equal to 1). Therefore, 800 time steps of the forward problem were simulated. In both simulations the data were 'measured' on three mutually adjacent faces of the cubical domain.

Slices by the plane x₃ = 0.25 of the 3D images reconstructed from accurate data are shown in figure 2, parts (b) and (c). Part (b) corresponds to the initial approximation f⁽⁰⁾; part (c) represents the second iteration f⁽²⁾. The gray scale image in part (b) looks quite close to the slice of the phantom shown in part (a). However, let us look at figure 3 which shows the profiles of the functions along line x₂ = x₃ = 0.25. On can see that the first approximation f⁽⁰⁾ (shown by the dashed line) is not quite accurate. On the other hand, the profile corresponding to the second iteration f⁽²⁾ (represented by the gray solid line) practically coincides with that of f (black line).

Zoom In Zoom Out Reset image size

In order to test the sensitivity of our reconstruction technique to the noise always contained in real data, we added a strong noise component to the simulated data. The noise was modeled by a normally distributed random variable; it was scaled so that the noise intensity (as measured by the standard L₂ norm) coincided with the intensity of the unperturbed signal. Figure 4 demonstrates visually the high level of noise in the data; the black line shows exact signal g₁(t, 0.5, 0.5), while the gray line represents the same signal with added noise.

Zoom In Zoom Out Reset image size

It should be noted that most tomography modalities tend to amplify noise, and in the presence of 100% noise either the reconstructed image would be overwhelmed by noise-related artifacts, or significant blurring would occur due to the low-frequency filtration needed to regularize the reconstruction. This does not happen here, since the singularities of the solution p(t, x) do not get smoothed on their way to the detectors due to the properties of the wave equation. (Such low sensitivity to noise was observed previously [17, 22] in other inverse problems related to PAT/TAT.)

Figure 5, parts (b) and (c), contains images reconstructed from the noisy data (the initial approximation, f⁽⁰⁾, and the first iteration, f⁽¹⁾, respectively). The noise is almost unnoticeable in the grayscale images. In addition to the low noise sensitivity of the method, this phenomenon is partially explained by the nature of the phantom (the volume of the support of f is actually quite small compared to the volume of the cube).

Zoom In Zoom Out Reset image size

Figure 6 presents the profiles of the cross-sections by the line x₂ = x₃ = 0.25. One can see that a noticeable noise is present in the reconstructions, and that the first iteration (gray line) does represent a considerable improvement comparing to the initial approximation f⁽⁰⁾. In our simulation, in the presence of the strong noise the consecutive approximations f^(j), j = 2, 3, 4, ... did not show any improvement over f⁽¹⁾; moreover, a very slow growth in the level of noise was observed.

Zoom In Zoom Out Reset image size

The computation time in the above simulations was about 9 min per iteration (excluding the input/output time) on a 2.6 GHz processor of a desktop computer. The code was not highly optimized, and it ran as a single-thread (no parallelization). Most of the elapsed time (7 min) was spent computing the forward problem (i.e. $\mathcal {W}f^{(K-1)}$). This is due to a large constant factor implicit in the operation count of the NUFFT used on this step.

The computing time required by our algorithm is of the same order of magnitude as that of the fast algorithm for the free-space problem with a cube surface acquisition [19] (estimated 4 to 5 min on the 401 × 401 × 401 grid). Somewhat faster reconstruction time (about 1 min for a slightly larger grid) was reported in [21] for a fast algorithm for the free-space problem involving integrating linear detectors. (The latter two problems are somewhat simpler computationally; in particular, the algorithms are non-iterative.) It should be noted that the time reversal algorithm for the free-space problem implemented using finite differences (with complexity $\mathcal {O}(N^{4})$ flops) would take an estimated time of 3.5 h on a grid of our size, and methods based on a straightforward discretization of explicit inversion formulas (whose complexity equals $\mathcal {O}(N^{5})$ flops) would take several days to complete one reconstruction.

3.3.2. Measurements done from all six faces