Exact Series Reconstruction in Photoacoustic Tomography with Circular Integrating Detectors
Abstract
A method for photoacoustic tomography is presented that uses circular integrals of the acoustic wave for the reconstruction of a threedimensional image. Image reconstruction is a twostep process: In the first step data from a stack of circular integrating detectors are used to reconstruct the circular projection of the source distribution. In the second step the inverse circular Radon transform is applied. In this article we establish inversion formulas for the first step, which involves an inverse problem for the axially symmetric wave equation. Numerical results are presented that show the validity and robustness of the resulting algorithm.
Keywords. Radon transform; Photoacoustic tomography; photoacoustic microscopy; Hankel transform; image reconstruction; integrating detectors; axially symmetric; wave equation;
AMS classifications. 44A12, 65R32, 35L05, 92C55.
1 Introduction
The principle of Photoacoustic tomography (PAT), also called Thermoacoustic tomography, is based on the excitation of high bandwidth acoustic waves with pulsed nonionizing electromagnetic energy inside tissue [9, 17, 23, 25, 28]. PAT presents a hybrid imaging technique that combines the advantages of optical (high contrast) and ultrasound imaging (high resolution). It has proven great potential for important medical applications including cancer diagnostics [16, 21] and imaging of vasculature [7, 15].
The common approach uses small conventional piezoelectric transducers that approximate point detectors to measure the acoustic pressure [28]. Reconstruction algorithms which are based on the point detector assumption yield images with a spatial resolution that is essentially limited by the size of the used piezoelectric transducers [27]. The size of the detector can in principle be reduced, but only at the cost of also reducing the signaltonoise ratio. In order to overcome this limitation large size planar or linear integrating detectors have been proposed in [3, 12]. Line shaped detectors integrate the acoustic pressure over its length and can be implemented by a MachZehnder [22] or a FabryPerot interferometer [10], which naturally integrate the acoustic pressure over the length of a laser beam.
A drawback of linear integrating detectors is that because of attenuation parts of the detector which are distant from the object may be less influenced by the pressure wave. Moreover, linear integrating detectors do not provide a compact experimental buildup. In [32] so called circular integrating detectors where introduced, which integrate the acoustic pressure over circles. A circular integrating detector can, similar to a linear integrating detector, be implemented by an interferometer where the laser beam is guided along a circle in an optical fiber. It is free of any aperture effect and can provide a uniformly and high resolution throughout the image area. Since it is possible to fabricate optical fibers out of materials which have nearly the same acoustical density like the fluid in which they are contained [10] no shadowing effects due to the circular integrating detectors are expected. The use of circular integrating detectors has been proposed independently in [31]. However, their study was limited to two spatial dimensions, where the circular shaped detector can be used as virtual point detector.
In [32] we derived an inversion formula based on the expansion of the involved functions in special basisfunctions. This formula, is hard to implement directly due to possible division by a zero (see Subsection 2.2). In this article we prove two novel exact series solutions that allow for stable implementation. As a byproduct we obtain a novel inversion formula for PAT using pointlike detectors on a cylindrical recording surface (see Remark 3.6).
2 Photoacoustic Tomography with Circular Integrating Detectors
In PAT an acoustic pressure wave , inside an object, is generated by a pulse of nonionizing electromagnetic radiation. In the case of spatially constant sound speed, the induced acoustic pressure satisfies the initial value problem [25, 28]
(2.1)  
(2.2)  
(2.3) 
PAT is concerned with recovering the initial pressure from measurements of outside the support of .
In [32] it is proposed to measure the acoustic signals with a stack of parallel circles which is rotated around a single axis, see Figure 1. In such a situation, threedimensional imaging involves the inversion of the classical circular Radon transform and the inversion of a reduced wave equation, as outlined in the following.
Throughout this article it is assumed that is smooth and supported in the cylinder , where is a fixed positive number. Let denote the unique solution of (2.1)(2.3) and, for , define
(2.4)  
(2.5) 
where
The stack of circular integrating detectors measures
(2.6) 
with denoting the fixed radius of the detectors.
The goal is to recover the unknown initial data from measured data .
2.1 Two Stage Reconstruction
Reconstruction with circular integrating detectors is based on the following reduction to the axial symmetric wave equation:
Proposition 2.1.
Proof.
Remark 2.2.
Proposition 2.1 is the basis of the following twostage procedure which reconstructs the initial pressure in (2.1)(2.3) from the data :

For (fixed position of the stack of circles) determine the initial pressure of (2.7)(2.9) using data .
Repeating this procedure for every , one obtains a family of functions , , corresponding to averages over circles centered on .

Next one recognizes that for fixed , the function
is the circular mean transform of with centers on a circle. For the circular mean transform stable analytic inversion formulas have been discovered recently [2, 8, 11, 18, 19]. Exemplarily, one of the formulas in [8, Theorem 1.1] reads
(2.10) where denotes the coordinates in the plane .
The key task for reconstruction is to derive stable and fast algorithms to reconstruct the initial data in (2.7)(2.9) from measurement data . A possible reconstruction method is based of time reversal (backpropagation) similar to [4, 5, 14]. However, the degeneration of at and the open detector set may cause difficulties in such procedures. The inversion approach in this paper is based on analytic inversion formulas for reconstructing .
2.2 Exact Inversion Formula
In the following denote by
the sine, cosine, Fourier, and Hankel transform, respectively. (Here is the zero order Bessel function [1].) When the above transforms are applied to functions depending on several variables then the transformed variable is added as subscript, e.g., .
Proof.
The zero order Bessel function satisfies on . Hence the chain rule implies that for every . Separation of variables shows that the functions
solve (2.7), (2.9). Employing the initial condition (2.8) and the inversion formulas for the Fourier and Hankel transforms it follows that the unique bounded solution of (2.7)(2.9) is given by
(2.12) 
with .
Proposition 2.3 implies that can be reconstructed from data as follows:

The data are Fourier and cosine transformed, yielding to .

According to (2.11), is mapped to .

Finally, application of the inverse Fourier and Hankel transforms yields
Remark 2.4 (Instability of (2.11)).
Inversion formula (2.11) is not defined when equals . From the proof of the above theorem it is clear that for exact data
(2.14) 
with denoting the zeros of . In practice, however, only noisy (approximately measured) data are available. In general,
It is therefore difficult to stably evaluate the quotient in (2.11) in practice.
Because the acoustic pressure is measured outside of the investigated object, only the following two situations occur in practical applications, see Figure 2:

The stack of circles is strictly outside the object. In this case and .

The object is enclosed in the stack of circles. In this case .
For the case we will provide two stable formulas based on expansions in bases of special functions. The case turns out to be harder, we currently do not have a stable alterative to (2.11). In this case the function is not supported in the interval and thus it cannot be expanded into a Fourier Bessel series which is crucial in the proves of theorems 3.1 and 3.4. However, in the limiting case a stable reconstruction formula is obtained, compare with Remark 3.6 below.
3 Stable Inversion Formulas
Our first stable inversion formula is as follows:
Theorem 3.1 (The D’Hospital trick).
Assume . Then
(3.1) 
for any .
Proof.
Remark 3.2.
In the sequel we derive an additional inversion formula that circumvents the division by zero problem. In fact, our formula will be a consequence of the following result derived in [30].
Proposition 3.3.
Now the second stable inversion formula can be stated as follows:
Theorem 3.4.
Assume . Then
(3.5) 
Here is extended by for .
Proof.
Equation (3.6) is quite similar to (2.11). However, in the denominator in (3.6) the zero order second kind Hankel function appears (instead of the the zero order Bessel function) which cannot be zero for a finite argument [1]. Moreover, the asymptotic expansion of the Bessel and the second kind Hankel function show that the summands in (3.5) that do not depend on the data remain bounded as .
Remark 3.5.
Remark 3.6.
Suppose that and that is supported in . For let denote the zeros of . Then one can expand in a Fourier Bessel series (see [26])
(3.8) 
According to (2.11) we have
(3.9) 
If we assume that is a integer multiple of , i.e., , then for any . Therefore, inserting (3.9) in (3.8) yields
(3.10) 
In general, (3.10) is a again sensible to noise when gets close to a zero of .
In the limiting case and for not to large, the denominators are well bounded from below (see right image in Figure 3). In this case, truncating (3.10) leads to a stable inversion formula. In particular, for one obtains
(3.11) 
for any . Together with (2.10) this provides a novel inversion formula for PAT using pointlike detectors on a cylindrical recording surface.
4 Numerical Experiments
In practice one deals with discrete measurement data
where is as in (2.6), and where , and are discrete samples of the angle, height and time, respectively. Here represents the finite height of the stack of circular integrating detectors (see Figure 1) and is such that for and .
In this section we outline how to implement (3.1) and (3.5) in order to find an approximation
with . Having calculated such an approximation, one can reconstruct a discrete approximation to by applying the filtered backprojection algorithm of [8] for fixed , see Remark 2.2.
A numerical reconstruction method based on (3.1) is as follows:

The discrete Fourier transform (with respect to the first component) of the data
(4.1) with , is considered as an approximation to .

The sine transform , evaluated at
is approximated by the trapezoidal rule, leading to
(4.2)
A numerical reconstruction method using (3.1) is be obtained in an analogous manner. In this case one replaces (4.3) by
(4.4) 
which is the discrete analogue of (3.5).
To give a rough estimate of the computational complexity for the previous calculations let us assume and that the values of the sine function and the Bessel function are precomputed and stored in lookup tables. Then the evaluation (4.1) needs floating point operations (FLOPS) whereas (4.2) and (4.3) require FLOPS. The filtered back projection formula (2.10) also requires FLOPS, see [8]. For three dimensional reconstruction (4.1), (4.2), (4.3) and the filtered backprojection formula have to be applied times. Hence the total number of FLOPS is estimated as
(4.5) 
Note that three dimensional backprojection type formulas which use point measurement data have complexity .
In the following numerical experiments we take , and . The synthetic initial data is assumed to be a superposition of radially symmetric objects around centers , i.e.,
The acoustic pressure generated by a single radially symmetric object at position and time is given by (see [13])
(4.6) 
By the superposition principle the total pressure is
The measurement data , see (2.4), (2.6), were generated by evaluating of (4.6) followed by numerical integration over . The radius of the circular integrating detectors is chosen to be . In this case the stack of circular integrating detectors fully encloses the synthetic initial data , see right image in in Figure 2.
Figure 4 shows a vertical cross section of the initial pressure and the measurement data where Gaussian noise with a variance of of the maximal data valued is added. The stack of circular integrating detectors is centered to the left of the objects. The presented discrete implementation measurements in time and in space. In both reconstructions the value was chosen to be
The reconstructions of with (4.3) from exact and noisy data are depicted in Figure 5 from formula (4.3) and with formula (4.4) in Figure 6. In the reconstructed images one notices some blurred boundaries which are limited data artifacts [20, 24, 29] arising from the finite height of the stack of circular integrating detectors. Moreover the images reconstructed with (4.4) are less sensitive to noise.
5 Conclusion
In this article a novel experimental buildup for PAT using circular integrating detectors was proposed. For collecting measurement data a fiberbased MachZehnder or FabryPerot interferometer can be used as an circular integrating detector. We showed that the 3D imaging problem reduces to a series of 2D problems. This decomposition can be used to reduce the operation count of derived reconstruction algorithms. We derived two stable exact reconstruction formulas, (3.1) and (3.5), for the case that the object is contained in the stack of detecting circles. In the case where the object is outside the detecting circles, a stable reconstruction formula is obtained for the limiting case . As a byproduct, this leads to a novel reconstruction formula (3.11) for PAT using point detectors on a cylindrical surface.
Acknowledgement
This work has been supported by the Austrian Science Foundation (FWF) within the framework of the NFN “Photoacoustic Imaging in Biology and Medicine”, Project S10505N20. Moreover, the work of M. Haltmeier has been supported by the technology transfer office of the University Innsbruck (transIT).
References
 [1] M. Abramowitz and I.A. Stegun. Handbook of Mathematical Functions. Dover, New York, 1972.
 [2] M. L. Agranovsky, K. Kuchment, and E. T. Quinto. Range descriptions for the spherical mean Radon transform. J. Funct. Anal., 248(2):344–386, 2007.
 [3] P. Burgholzer, C. Hofer, G. Paltauf, M. Haltmeier, and O. Scherzer. Thermoacoustic tomography with integrating area and line detectors. IEEE Trans. Ultrason., Ferroeletr., Freq. Control, 52(9):1577–1583, 2005.
 [4] P. Burgholzer, G. J. Matt, M. Haltmeier, and G. Paltauf. Exact and approximate imaging methods for photoacoustic tomography using an arbitrary detection surface. Phys. Rev. E, 75(4):046706, 2007.
 [5] C. Clason and M. Klibanov. The quasireversibility method for thermoacoustic tomography in a heterogeneous medium. SIAM J. Sci. Comp., 2007. accepted.
 [6] C. Depeursinge, editor. Novel Optical Instrumentation for Biomedical Applications III, volume 6631 of Proceedings of SPIEOSA, 2007.
 [7] R. O. Esenaliev, I. V. Larina, K. V. Larin, D. J. Deyo, M. Motamedi, and D. S. Prough. Optoacoustic technique for noninvasive monitoring of blood oxygenation: a feasibility study. App. Opt., 41(22):4722–4731, 2002.
 [8] D. Finch, M. Haltmeier, and Rakesh. Inversion of spherical means and the wave equation in even dimensions. SIAM J. Appl. Math., 68(2):392–412, 2007.
 [9] D. Finch and Rakesh. The spherical mean value operator with centers on a sphere. Inverse Probl., 23(6):37–49, 2007.
 [10] H. Grün, M. Haltmeier, G. Paltauf, and P. Burgholzer. Photoacoustic tomography using a fiber based Fabry–Perot interferometer as an integrating line detector and image reconstruction by modelbased time reversal method. In [6], 2007.
 [11] M. Haltmeier, O. Scherzer, P. Burgholzer, R. Nuster, and G. Paltauf. Thermoacoustic tomography & the circular Radon transform: Exact inversion formula. Math. Models Methods Appl. Sci., 17(4):635–655, 2007.
 [12] M. Haltmeier, O. Scherzer, P. Burgholzer, and G. Paltauf. Thermoacoustic imaging with large planar receivers. Inverse Probl., 20(5):1663–1673, 2004.
 [13] M. Haltmeier, T. Schuster, and O. Scherzer. Filtered backprojection for thermoacoustic computed tomography in spherical geometry. Math. Methods Appl. Sci., 28(16):1919–1937, 2005.
 [14] Y. Hristova, P. Kuchment, and L. Nguyen. Reconstruction and time reversal in thermoacoustic tomography in acoustically homogeneous and inhomogeneous media. Inverse Problems, 24(5):055006 (25pp), 2008.
 [15] R. G. M. Kolkman, E. Hondebrink, W. Steenbergen, and F. F. M. De Mul. In vivo photoacoustic imaging of blood vessels using an extremenarrow aperture sensor. IEEE J. Sel. Topics Quantum Electron., 9(2):343–346, 2003.
 [16] R. A. Kruger, K. D. Miller, H. E. Reynolds, W. L. Kiser, D. R. Reinecke, and G. A. Kruger. Breast cancer in vivo: contrast enhancement with thermoacoustic CT at 434 MHzfeasibility study. Radiology, 216(1):279–283, 2000.
 [17] P. Kuchment and L. A. Kunyansky. Mathematics of thermoacoustic and photoacoustic tomography. European J. Appl. Math., 19:191–224, 2008.
 [18] L. A. Kunyansky. Explicit inversion formulae for the spherical mean Radon transform. Inverse Probl., 23(1):373–383, 2007.
 [19] L. A. Kunyansky. A series solution and a fast algorithm for the inversion of the spherical mean radon transform. Inverse Probl., 23(6):S11–S20, 2007.
 [20] A.K. Louis and E.T. Quinto. Local tomographic methods in sonar. In Surveys on solution methods for inverse problems, pages 147–154. Springer, Vienna, 2000.
 [21] S. Manohar, A. Kharine, J. C. G. van Hespen, W. Steenbergen, and T. G. van Leeuwen. The twente photoacoustic mammoscope: system overview and performance. Physics in Medicine and Biology, 50(11):2543–2557, 2005.
 [22] G. Paltauf, R. Nuster, M. Haltmeier, and P. Burgholzer. Experimental evaluation of reconstruction algorithms for limited view photoacoustic tomography with line detectors. Inverse Probl., 23(6):81–94, 2007.
 [23] S. K. Patch and O. Scherzer. Special section on photo and thermoacoustic imaging. Inverse Probl., 23:S1–S122, 2007.
 [24] E. T. Quinto. Singularities of the Xray transform and limited data tomography in and . SIAM Journal on Mathematical Analysis, 24(5):1215–1225, 1993.
 [25] O. Scherzer, M. Grasmair, H. Grossauer, M. Haltmeier, and F. Lenzen. Variational Methods in Imaging, volume 167 of Applied Mathematical Sciences. Springer, New York, 2008.
 [26] I. N. Sneddon. The Use of Integral Transforms. McGrawHill, New York, 1972.
 [27] M. Xu and L. V. Wang. Analytic explanation of spatial resolution related to bandwidth and detector aperture size in thermoacoustic or photoacoustic reconstruction. Phys. Rev. E, 67(5):0566051–05660515 (electronic), 2003.
 [28] M. Xu and L. V. Wang. Photoacoustic imaging in biomedicine. Rev. Sci. Instruments, 77(4):041101, 2006.
 [29]<