Photoacoustic Tomography in a Rectangular Reflecting Cavity

Almost all known image reconstruction algorithms for photoacoustic and thermoacoustic tomography assume that the acoustic waves leave the region of interest after a finite time. This assumption is reasonable if the reflections from the detectors and surrounding surfaces can be neglected or filtered out (for example, by time-gating). However, when the object is surrounded by acoustically hard detector arrays, and/or by additional acoustic mirrors, the acoustic waves will undergo multiple reflections. (In the absence of absorption they would bounce around in such a reverberant cavity forever). This disallows the use of the existing free-space reconstruction techniques. This paper proposes a fast iterative reconstruction algorithm for measurements made at the walls of a rectangular reverberant cavity. We prove the convergence of the iterations under a certain sufficient condition, and demonstrate the effectiveness and efficiency of the algorithm in numerical simulations.


Introduction
Photoacoustic tomography (PAT) and the closely related modality thermoacoustic tomography (TAT) [3,14,15,26,28] are based on the photoacoustic effect, in which an acoustic wave is generated by the absorption of an electromagnetic (EM) pulse. The distinction between PAT and TAT is that the former uses visible or near-infrared light, while the latter uses energy in the microwave region. There are several endogenous chromophores which absorb in these wavelength ranges, the most important of which are oxy-and deoxy-hemoglobin, and externally administered, molecularlytargeted, chromophores can be used as contrast agents. These emerging modalities are therefore attracting considerable interest for molecular and functional imaging applications in pre-clinical and clinical imaging.
When EM pulses at these wavelengths are sent into biological tissue, the EM energy will be absorbed and then thermalised. When this happens sufficiently rapidly, there will be simultaneous, small, increases in the temperature and pressure in the regions where the energy was absorbed. As tissue is an elastic medium, the pressure increase propagates as an acoustic (ultrasonic) pulse, and can be detected as a time series at the boundary of the tissue. The image of the initial pressure distribution can subsequently be reconstructed from these measurements by solving a certain inverse problem. Since acoustic waves propagate through soft tissues with little absorption or scattering, PAT and TAT yield high-resolution images related to the EM properties of the tissue. Such images cannot be obtained by purely optical or electrical (or electromagnetic) techniques such as, e.g. diffuse optical tomography or electrical impedance tomography.
In order to facilitate the acoustic measurements the object is usually submerged in water or a gel with acoustic properties close to those of the tissues. One method of measuring the acoustic pressure is to use a small single-element detector to record the time-dependent acoustic pressure at a point. In order to obtain enough data for the reconstruction, the detector is scanned across an imaginary acquisition surface surrounding the region of interest; the EM pulse must be re-emitted for each new location of the detector. The wave propagation within such an acquisition scheme can be modelled by the free-space wave equation, as the reflection of the waves from the single detector does not affect the measurement at that detector. The corresponding inverse problem has been studied extensively (see, for example reviews [16,27,28] and references therein), and various reconstruction techniques have been developed to solve it. These techniques include time reversal [4,13,30], series expansion methods [1,12,19,21,24,25], and several explicit inversion formulas [9,10,18,20,23,31] (the references above are by no means exhaustive).
In order to significantly speed up the data acquisition, the measurements may be made by an array of detectors. (Both piezoelectric and optically-addressed arrays have been used). Such an array typically forms a solid surface, a side effect of which is that the sound waves are reflected back into the region of interest. In the case of a single planar array the free space reconstruction methods still apply, since the reflected waves will propagate away from the detector and will not affect the measurements. However, it is not possible to obtain exact reconstructions from data measured on a finite section of a plane. To overcome this limitation, more advanced measuring systems will have the object surrounded by the detector arrays (or perhaps intentionally installed acoustic mirrors, e.g. [5]). One such system, based on optically-addressed Fabry-Perot detection arrays, is currently being investigated. The waves within the reverberant cavity formed by the detectors and passive reflectors will undergo multiple reflections. Under the idealized assumption of negligible acoustic absorption used by most classical methods, the oscillations within such a cavity will continue forever. All the standard reconstruction techniques assume that the signal either has a final duration (due to Huygens' principle in 3D) or dies out sufficiently fast (in 2D or in the presence of nonuniform speed of sound). As a result, these techniques are not applicable in the presence of such multiple reflections.
In the present paper we develop an algorithm for reconstructing the initial pressure distribution within a rectangular reverberant cavity from acoustic pressure time series measurements made on at least three mutually adjacent walls. The algorithm finds first a crude initial approximation that can be further refined iteratively. The derivation of the algorithm is presented in Section 2, together with the convergence analysis for the iterations (Section 2.4). The algorithm was tested in numerical simulations (Section 3) that showed both the high quality of the reconstruction and very low noise sensitivity of the proposed technique.
High resolution images in 3D utilized in modern biomedical practice are described by hundreds millions of unknowns. In order to be useful in practice, a reconstruction algorithm must be fast. The proposed technique is asymptotically fast; it requires O(N 3 log N ) floating point operations (flops) per iteration on an N ×N ×N computational grid. In our simulations the reconstruction time was comparable with that of known fast algorithms for various free-space problem (e.g. [19,21]).
1 Formulation of the problem 1

.1 Photoacoustic Tomography in Free Space
In conventional photoacoustic tomography, 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 (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): 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) 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 0 . 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 0 .) Now one can solve in Ω the following initial/boundary value problem (IBVP) 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 the ones 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 Ω.

Photoacoustic Tomography in a Reflective Cavity
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 travelling 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 initial boundary value problem: Now the measurements can be written as 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, this 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 be also 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 photoacoustic tomography within the reverberant cavity.

Derivation of the reconstruction algorithm
In this section we develop a fast and accurate reconstruction algorithm for photoacoustic tomography within a rectangular reverberant cavity. On one hand, such an acquisition geometry is one of the simplest from the mathematical standpoint, allowing one to design a simple and fast Fourier  Figure 1: The 2D reverberant cavity with two passive reflecting walls, and two reflecting detector arrays on which the data g 1 (t, x 2 ) and g 2 (t, x 1 ) are measured. based reconstruction technique. On the other hand, photoacoustic scanners with this particular configuration are quite conceivable -indeed a design based on optically-addressed Fabry-Perot ultrasound arrays [32] is currently under development -and algorithms will be needed to process the data.
The reconstruction technique we propose can be used in both 2D and 3D settings, with only minor changes. The 3D case is, obviously, more interesting from the applied point of view, while the 2D case is somewhat easier to present. Thus, in the next section we provide a derivation for the 2D algorithm, with the caveat that the 3D case is quite similar. In Section 3 we test the 3D version of the algorithm in numerical simulations.

2D formulation and the series solution of the forward problem
We assume that the acoustic pressure p(t, x), x = (x 1 , x 2 ) solves the initial/boundary value problem (5)-(7) in 2D, in a square Ω = [0, 1] × [0, 1]. Without loss of generality the speed of sound c 0 will be assumed to equal 1. The measurements are represented by the pair of time series g = (g 1 (t, x 2 ), g 2 (t, x 1 )) corresponding to the pressure values on the two adjacent sides of the square Ω ( Figure 1): Let us denote by W the linear operator transforming the initial condition f (x) into the boundary data g: Our ultimate goal is to reconstruct the initial pressure f (x) = p(0, x) from g, i.e. to invert W . The simple shape of domain Ω allows us to utilize the eigenfunctions of the Neumann Laplacian in a square with the eigenfrequencies ω k,l : ω k,l = π k 2 + l 2 , k, l = 0, 1, 2, ...
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: where numbers f k,l are the coefficients of the two-dimensional Fourier cosine series of the initial pressure f (x): For future use we will denote the double-indexed sequence of the Fourier coefficients by F = {f k,l } ∞ k,l=0 , and will use the notation F series to refer to the transformation defined by equations (9) and (10), so that F = F series f.
Obviously, in order to obtain f (x) it is enough to reconstruct coefficients F ; function f is then computed by inverting F series , i.e. by summing the Fourier cosine series:

What can be found from the data measured on one side?
Let us first analyze the connection between coefficients f k,l and the data measured on one side of the square (say, g 2 ): Due to the orthogonality of the cosine functions (in variable x 1 ) the above equation splits: 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 F −1 of a sequence of Dirac delta-functions: where the forward and inverse Fourier transforms of some function h(t) are defined as followŝ One can conclude that the Fourier transformĝ 2,k (ξ) of g 2,k (t) is the following distribution It is tempting to try to recover coefficients f k,l by applying the Fourier transform F to each g 2,k (t) and thus obtainingĝ 2,k (ξ). 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 η can be easily computed (for example, by applying the composite trapezoid rule to the discretized signal). Now, by the convolution theorem and, taking into account the equality η T (ξ) = T η(T ξ) and equation (14) one obtains Function η T (ξ) 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ĝ 2,k (ξ) from η T (ξ) * g 2,k (ξ) as a deconvolution problem. Due to the smoothness of the convolution kernel η T (ξ) 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ĝ 2,k (ξ) 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 η 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 : Further, by separating the 'diagonal' terms this system can be re-written in the form for l = 0, 1, 2, ..., k = 0, 1, 2..., except the case when k = l = 0 which assumes a simple form Let us discuss the properties of these systems (one system for every value of k). Due to the fast decrease of η T (ξ) at infinity, as T goes to infinity, values of η(T (ω k,l − ω k,m )) and η(T (ω k,l + ω 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 and thus to reconstruct p 0 (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 η(T (ω k,l − ω 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 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.

Reconstruction using from data measured on two adjacent sides
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 1 (x 2 , t) measured on the second side of the square corresponding to x 1 = 0. Define functions g 1,l (t) by the equations The derivation similar to the one in the previous section yields the following system of equations: 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: , k > l.
Since for very large values of T factors η(T (...)) vanish (while η(0) is a constant), one can compute a crude approximation f (0) (x) by defining the coefficients F (0) as follows and by summing the Fourier series Formulas (23), (24) define a linear operator R that maps boundary values g into the set of Fourier coefficients F , so that f (0) = F −1 series F (0) = F −1 series Rg =F −1 series RWf. As our numerical experiments show, approximation f (0) 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 A that transforms a two-dimensional sequence of numbers H = {h k,m } ∞ k,m=0 into a two-dimensional sequence of numbers E = {e k,m } ∞ k,m=0 by the formulas or, in the operator notation, E = AH.
Then equations (22) can be re-written as where F is the set of the Fourier coefficients of the sought f (x), and F (0) are the Fourier coefficients of the initial, crude approximation f (0) (x). If the spectral radius of operator A is less than 1, the unique solution of equation (27) can be found in the form of a converging Neumann series 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 In section 2.4 we analyze operator 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 F −1 series operator) iterations (29) can be replaced by the equivalent relation On the other hand, it follows from (27) that and, since (27) holds for any initial condition f (x), the following operator identity holds Thus, −F −1 series AF series = I − F −1 series RW, and, therefore, recurrence relation (30) is equivalent to Our computational algorithm is based on equations (31), (32); these iterations are theoretically equivalent to those defined by equation (29) and have the same convergence properties. The computational advantage in using (31), (32) lies in the possibility of easily implementing all the necessary operations using fast transforms (FFTs). The implementation details are discussed in Section 3. In the next section we analyze the convergence of such iterations, and the form (29) is more convenient for this purpose.

Convergence analysis
Let us consider the || · || ∞ norm defined on the space of double-indexed infinite sequences H = {h k,l } ∞ k,l=0 : ||H|| ∞ = sup k,l=0.1,2,3... |h k,l |, and a space L of such sequences with bounded ||·|| ∞ norm. Further, for a linear operator B : L → L define the induced ∞-norm: Our goal is to estimate the induced ∞-norm of the operator A defined by equations (26) in the previous sections. These equations contain factors in the form η(T (ω k,l − ω k,m )). Note that for a sufficiently smooth cut-off function η(t) its Fourier transform declines fast at infinity: Let us find a lower bound on |ω k,l − ω k,m |. Consider first the case l ≥ k: If, in addition m ≥ k, the above equation yields If m < k (and l ≥ k): so that the uniform bound holds if l ≥ k for all values of m: Now one can bound | η(T (ω k,l − ω k,m ))|, still under assumption l ≥ k: Similarly, if l > 0 or k > 0, and Let us now estimate the ∞-norm of the vector G = AH resulting from applying operator A to a sequence H whose ∞-norm equals 1, so that |h k,l | ≤ 1, l, k = 0, 1, 2, ... .
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 (33)-(36) one obtains The value of the latter series is well known (see [11], equation 0.233) resulting in the estimate The above computations can be replicated with very minor changes to show that (37) 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 A is bounded in the induced ∞-norm by The Lemma implies that if the acquisition time T is large enough, e.g. if then operator A is a contraction mapping in the induced ∞-norm, i.e. ||A|| ∞ < 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 (38)) the Neumann series (28) converges in the ∞-norm, implying convergence of iterations (31), (32) in the following sense

3D version of the method
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 1 , x 2 , x 3 ) solves IBVP (5)- (7) in 3D, in a cube Ω = [0, 1] × [0, 1] × [0, 1]. As before, the speed of sound c 0 will be assumed to equal 1. The measurements are represented by the triple of time series g = (g 1 (t, x 2 ), g 2 (t, x), g 3 (t, x)) corresponding to the pressure values on the three adjacent sides of the cube Ω: The solution of IBVP in the cube can be represented using the Neumann Laplacian eigenfunctions ϕ k,l,n (x) = cos(πkx 1 ) cos(πlx 2 ) cos(πnx 3 ) with eigenfrequencies ω k,l,n = π k 2 + l 2 + n 2 , k, l, n = 0, 1, 2, ... .
Conversely, f is obtained from F by summing the standard 3D cosine Fourier series, f = F −1 series F . Consider side of the cube corresponding to x 2 = 0. Then f k,l,n cos(ω k,l,n t) , t ∈ [0, T ]. (40) 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 2 by the 2D cosine Fourier series.
By finding the values of the Fourier transform g 2,k,n η 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 f k,l,n + η(2T ω k,l,n ) η(0) f k,l,n + ∞ m=0 m =l η(T (ω k,l,n − ω k,m,n )) + η(T (ω k,l,n + ω k,m,n )) η(0) f k,m,n = 2 T η T g 2,k,n (ω k,l,n ) η(0) , for k, l, n = 0, 1, 2, ..., except the case when k = l = n = 0 which yields a simpler 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 obtain in a similar fashion). We thus arrive at the system where the right hand side corresponds to the Fourier coefficients of the crude first approximation , l > k, l ≥ n, , n > l, n > k, and where the terms in the brackets define the operator 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 R is now defined by (42) and the operator 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: |ω k,l,n − ω k,m,n | = π k 2 + l 2 + n 2 − k 2 + m 2 + n 2 = π |l − m|(l + m) √ k 2 + l 2 + n 2 + √ k 2 + m 2 + n 2 .
If l ≥ k and l ≥ n and m ≥ l, then If l ≥ k and l ≥ n and m < l, then 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 A.
This immediately leads to a sufficient condition for the convergence of the 3D algorithm if the acquisition time T satisfies the condition

Numerical implementation and simulations
In this section the proposed reconstruction algorithm is implemented numerically and tested in numerical simulations. Since real measurements are made in 3D, we will concentrate on the 3D version of the method. The 2D algorithm is very similar. In 3D tomography reconstructions, numerical complexity becomes a major issue, since in a fine 3D mesh the number of unknowns can easily exceed hundreds of millions. Thus, it is highly desirable to have an asymptotically fast reconstruction algorithm reconstructing images on a N × N ×N computational grid in (preferably) O(N 3 ) or (at most) O(N 3 log N ) floating point operations (flops). Our technique achieves the latter operation count by utilizing the various Fast Fourier Transforms (FFTs) on all computationally intensive steps of the algorithm, as described below.

Forward problem
We remind the reader that the reconstruction algorithm is defined by equations (31), (32). 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 ) ≡ p(t, x)| x 2 =0 the following steps are done: 1. Expand f (x) into the 3D Fourier cosine series, obtain coefficients f k,l,n (equation (39)), k, l, n = 0, 1, 2, ...N .
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 O(N ). Step 1 is done by the 3D cosine FFT, it requires O(N 3 log N ) flops.
Step 2 is slightly more complicated. Evaluation of (41) 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 Fast Fourier Transform (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 O(N log N ) estimate of the complexity of this technique. In our simulations the accuracy of NSFFT was of order of 10 −5 . Since there is N × N of such transforms, the complexity of Step 2 is O(N 3 log N ) flops.
Step 3 is implemented by summing the 2D cosine series (using the corresponding type of To summarize the expense of finding g 2 (t, x 1 , x 3 ) by our method is 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.

Approximate inversion
The iterations given by equation (31) can be viewed as repeated applications of the forward operator W and approximate inverse, F −1 series R, in alternating order. The summation of the Fourier series F −1 series is done by the 3D cosine FFT algorithm; it takes O(N 3 log N ) flops. Operator R is defined by equations (42); 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 O(N ) × O(N 2 log N ) = O(N 3 log N ) flops. The only remaining non-trivial operation here is computing the 1D Fourier transforms η T g ... of products η T g ... at non-equally spaced frequencies. In our simulation we simply applied the regular FFT to compute η 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 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 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 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 O(N 3 log N ) flops.

Simulations
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.

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 3 = 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 √ 3. 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 3 = 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 (0) ; parts (c) represents the second iteration f (2) . 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 the Figure ( 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 modelled by a normally distributed random variable; it was scaled so that the noise intensity (as measured by the standard L 2 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 1 (t, 0.5, 0.5), while the gray line represents the same signal with added noise.
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 (0) , and the first iteration, f (1) , 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).  Figure 6 presents the profiles of the cross sections by the line x 2 = x 3 = 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 (0) . 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 (1) ; moreover, a very slow growth in the level of noise was observed.
The computation time in the above simulations was about 9 minutes 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. Wf (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 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 O(N 4 ) flops) would take an estimated 3.5 hours on a grid of our size, and methods based on a straightforward discretization of explicit inversion formulas (whose complexity equals O(N 5 ) flops) would take several days to complete one reconstruction.

Measurements done from all six faces
We have also conducted simulations aimed at estimating the impact of additional measurements (done on additional faces of the cube). Clearly, if the data were measured on the other three faces of the cube (corresponding to planes x j = 1, j = 1, 2, 3), an algorithm very similar to described above could be used for reconstruction. When the data are measured on all six faces, the image is obtained simply by averaging the results obtained from the two sets of three-face measurements. Our simulations demonstrate that with the additional data the measurement time T can be decreased. With measurement time T = 1 and a six-face acquisition the results of reconstruction were quite similar to those corresponding to the three-face acquisition and T = 2, as described at the beginning of this section.