# 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 three-dimensional image. Image reconstruction is a two-step 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.

Radon transform; Photoacoustic tomography; photoacoustic microscopy; Hankel transform; image reconstruction; integrating detectors; axially symmetric; wave equation;

44A12, 65R32, 35L05, 92C55.

## 1Introduction

The principle of Photoacoustic tomography (PAT), also called Thermoacoustic tomography, is based on the excitation of high bandwidth acoustic waves with pulsed non-ionizing electromagnetic energy inside tissue [9]. 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] and imaging of vasculature [7].

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 signal-to-noise ratio. In order to overcome this limitation large size planar or linear integrating detectors have been proposed in [3]. Line shaped detectors integrate the acoustic pressure over its length and can be implemented by a Mach-Zehnder [22] or a Fabry-Perot 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 basis-functions. 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 point-like detectors on a cylindrical recording surface (see Remark ?).

The outline of this article is as follows. In Section Section 2 we review PAT with circular integrating detectors and recall the main results of [32]. In the following Sections we will be concerned with derivation and implementation of the novel stable inversion formulas.

## 2Photoacoustic Tomography with Circular Integrating Detectors

In PAT an acoustic pressure wave , inside an object, is generated by a pulse of non-ionizing electromagnetic radiation. In the case of spatially constant sound speed, the induced acoustic pressure satisfies the initial value problem [25]

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 ?. In such a situation, three-dimensional 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 (Equation 1)-( ?) and, for , define

where

The stack of circular integrating detectors measures

with denoting the fixed radius of the detectors.

The goal is to recover the unknown initial data from measured data .

### 2.1Two Stage Reconstruction

Reconstruction with circular integrating detectors is based on the following reduction to the axial symmetric wave equation:

Equations ( ?), ( ?) and the boundedness as immediately follow from ( ?), ( ?) and the definitions of and . In the cylindrical coordinates , the Laplace operator is given by the well known expression

Integrating the equation with respect to yields ( ?).

Note that ( ?)-( ?) is uniquely solvable if we require that its solution remains bounded as .

The key task for reconstruction is to derive stable and fast algorithms to reconstruct the initial data in ( ?)-( ?) from measurement data . A possible reconstruction method is based of time reversal (back-propagation) similar to [4]. 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.2Exact 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., .

The zero order Bessel function satisfies on . Hence the chain rule implies that for every . Separation of variables shows that the functions

solve ( ?), ( ?). Employing the initial condition ( ?) and the inversion formulas for the Fourier and Hankel transforms it follows that the unique bounded solution of ( ?)-( ?) is given by

with .

Substituting in (Equation 4) and putting afterwards, leads to

The inversion formulas for the Fourier and Cosine transforms now imply that

Solving the last equation for shows ( ?).

Proposition ? implies that can be reconstructed from data as follows:

The data are Fourier and cosine transformed, yielding to .

According to ( ?), is mapped to .

Finally, application of the inverse Fourier and Hankel transforms yields

Because the acoustic pressure is measured outside of the investigated object, only the following two situations occur in practical applications, see Figure ?:

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 ( ?). 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 ? and ?. However, in the limiting case a stable reconstruction formula is obtained, compare with Remark ? below.

## 3Stable Inversion Formulas

In the following we fix , , define , by ( ?), (Equation 3), and let denote the zeros of the function .

Our first stable inversion formula is as follows:

The assumptions and imply that is compactly supported in . It can therefore be expanded in a Fourier Bessel series [26]

According to ( ?) we have

Applying the rule of D’Hospital gives

Inserting (Equation 7) in (Equation 6) and using the Fourier inversion formula shows ( ?).

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

Now the second stable inversion formula can be stated as follows:

We use again the Fourier-Bessel series (Equation 6) of proof of Theorem ?. Recalling the definitions of , and the Fourier coefficients , one notices that , . Therefore ( ?) for implies

Inserting (Equation 8) in (Equation 6) and using the Fourier inversion formula shows ( ?).

Equation (Equation 8) is quite similar to ( ?). However, in the denominator in (Equation 8) 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 ( ?) that do not depend on the data remain bounded as .

## 4Numerical Experiments

In practice one deals with discrete measurement data

where is as in (Equation 3), 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 ?) and is such that for and .

In this section we outline how to implement ( ?) and ( ?) in order to find an approximation

with . Having calculated such an approximation, one can reconstruct a discrete approximation to by applying the filtered back-projection algorithm of [8] for fixed , see Remark ?.

A numerical reconstruction method based on ( ?) is as follows:

The discrete Fourier transform (with respect to the first component) of the data

with , is considered as an approximation to .

The sine transform , evaluated at

is approximated by the trapezoidal rule, leading to

Finally, truncating the Fourier Bessel Series and approximating the inverse Fourier transform with the trapezoidal rule leads to discrete version of ( ?):

with in formula (Equation 11)

A numerical reconstruction method using ( ?) is be obtained in an analogous manner. In this case one replaces (Equation 11) by

which is the discrete analogue of ( ?).

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 pre-computed and stored in lookup tables. Then the evaluation (Equation 9) needs floating point operations (FLOPS) whereas (Equation 10) and (Equation 11) require FLOPS. The filtered back projection formula ( ?) also requires FLOPS, see [8]. For three dimensional reconstruction (Equation 9), (Equation 10), (Equation 11) and the filtered back-projection formula have to be applied times. Hence the total number of FLOPS is estimated as

Note that three dimensional back-projection 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])

By the superposition principle the total pressure is

The measurement data , see (Equation 2), (Equation 3), were generated by evaluating of (Equation 13) 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 ?.

Figure 2 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 (Equation 11) from exact and noisy data are depicted in Figure 3 from formula (Equation 11) and with formula (Equation 12) in Figure Figure 4. In the reconstructed images one notices some blurred boundaries which are limited data artifacts [20] arising from the finite height of the stack of circular integrating detectors. Moreover the images reconstructed with (Equation 12) are less sensitive to noise.

## 5Conclusion

In this article a novel experimental buildup for PAT using circular integrating detectors was proposed. For collecting measurement data a fiber-based Mach-Zehnder or Fabry-Perot 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, ( ?) and ( ?), 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 ( ?) 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 S10505-N20. Moreover, the work of M. Haltmeier has been supported by the technology transfer office of the University Innsbruck (transIT).

### References

*Handbook of Mathematical Functions*.

M. Abramowitz and I.A. Stegun. Dover, New York, 1972.**Range descriptions for the spherical mean Radon transform.**

M. L. Agranovsky, K. Kuchment, and E. T. Quinto.*J. Funct. Anal.*, 248(2):344–386, 2007.**Thermoacoustic tomography with integrating area and line detectors.**

P. Burgholzer, C. Hofer, G. Paltauf, M. Haltmeier, and O. Scherzer.*IEEE Trans. Ultrason., Ferroeletr., Freq. Control*, 52(9):1577–1583, 2005.**Exact and approximate imaging methods for photoacoustic tomography using an arbitrary detection surface.**

P. Burgholzer, G. J. Matt, M. Haltmeier, and G. Paltauf.*Phys. Rev. E*, 75(4):046706, 2007.**The quasi-reversibility method for thermoacoustic tomography in a heterogeneous medium.**

C. Clason and M. Klibanov.*SIAM J. Sci. Comp.*, 2007.*Novel Optical Instrumentation for Biomedical Applications III*, volume 6631 of*Proceedings of SPIE-OSA*, 2007.

C. Depeursinge, editor.**Optoacoustic technique for noninvasive monitoring of blood oxygenation: a feasibility study.**

R. O. Esenaliev, I. V. Larina, K. V. Larin, D. J. Deyo, M. Motamedi, and D. S. Prough.*App. Opt.*, 41(22):4722–4731, 2002.**Inversion of spherical means and the wave equation in even dimensions.**

D. Finch, M. Haltmeier, and Rakesh.*SIAM J. Appl. Math.*, 68(2):392–412, 2007.**The spherical mean value operator with centers on a sphere.**

D. Finch and Rakesh.*Inverse Probl.*, 23(6):37–49, 2007.**Photoacoustic tomography using a fiber based Fabry–Perot interferometer as an integrating line detector and image reconstruction by model-based time reversal method.**

H. Grün, M. Haltmeier, G. Paltauf, and P. Burgholzer. In*, 2007.***Thermoacoustic tomography & the circular Radon transform: Exact inversion formula.**

M. Haltmeier, O. Scherzer, P. Burgholzer, R. Nuster, and G. Paltauf.*Math. Models Methods Appl. Sci.*, 17(4):635–655, 2007.**Thermoacoustic imaging with large planar receivers.**

M. Haltmeier, O. Scherzer, P. Burgholzer, and G. Paltauf.*Inverse Probl.*, 20(5):1663–1673, 2004.**Filtered backprojection for thermoacoustic computed tomography in spherical geometry.**

M. Haltmeier, T. Schuster, and O. Scherzer.*Math. Methods Appl. Sci.*, 28(16):1919–1937, 2005.**Reconstruction and time reversal in thermoacoustic tomography in acoustically homogeneous and inhomogeneous media.**

Y. Hristova, P. Kuchment, and L. Nguyen.*Inverse Problems*, 24(5):055006 (25pp), 2008.**In vivo photoacoustic imaging of blood vessels using an extreme-narrow aperture sensor.**

R. G. M. Kolkman, E. Hondebrink, W. Steenbergen, and F. F. M. De Mul.*IEEE J. Sel. Topics Quantum Electron.*, 9(2):343–346, 2003.**Breast cancer in vivo: contrast enhancement with thermoacoustic CT at 434 MHz-feasibility study.**

R. A. Kruger, K. D. Miller, H. E. Reynolds, W. L. Kiser, D. R. Reinecke, and G. A. Kruger.*Radiology*, 216(1):279–283, 2000.**Mathematics of thermoacoustic and photoacoustic tomography.**

P. Kuchment and L. A. Kunyansky.*European J. Appl. Math.*, 19:191–224, 2008.**Explicit inversion formulae for the spherical mean Radon transform.**

L. A. Kunyansky.*Inverse Probl.*, 23(1):373–383, 2007.**A series solution and a fast algorithm for the inversion of the spherical mean radon transform.**

L. A. Kunyansky.*Inverse Probl.*, 23(6):S11–S20, 2007.**Local tomographic methods in sonar.**

A.K. Louis and E.T. Quinto. In*Surveys on solution methods for inverse problems*, pages 147–154. Springer, Vienna, 2000.**The twente photoacoustic mammoscope: system overview and performance.**

S. Manohar, A. Kharine, J. C. G. van Hespen, W. Steenbergen, and T. G. van Leeuwen.*Physics in Medicine and Biology*, 50(11):2543–2557, 2005.**Experimental evaluation of reconstruction algorithms for limited view photoacoustic tomography with line detectors.**

G. Paltauf, R. Nuster, M. Haltmeier, and P. Burgholzer.*Inverse Probl.*, 23(6):81–94, 2007.**Special section on photo- and thermoacoustic imaging.**

S. K. Patch and O. Scherzer.*Inverse Probl.*, 23:S1–S122, 2007.**Singularities of the X-ray transform and limited data tomography in and .**

E. T. Quinto.*SIAM Journal on Mathematical Analysis*, 24(5):1215–1225, 1993.*Variational Methods in Imaging*, volume 167 of*Applied Mathematical Sciences*.

O. Scherzer, M. Grasmair, H. Grossauer, M. Haltmeier, and F. Lenzen. Springer, New York, 2008.*The Use of Integral Transforms*.

I. N. Sneddon. McGraw-Hill, New York, 1972.**Analytic explanation of spatial resolution related to bandwidth and detector aperture size in thermoacoustic or photoacoustic reconstruction.**

M. Xu and L. V. Wang.*Phys. Rev. E*, 67(5):0566051–05660515 (electronic), 2003.**Photoacoustic imaging in biomedicine.**

M. Xu and L. V. Wang.*Rev. Sci. Instruments*, 77(4):041101, 2006.**Reconstructions in limited-view thermoacoustic tomography.**

Y. Xu, L. V. Wang, G. Ambartsoumian, and P. Kuchment.*Med. Phys.*, 31(4):724–733, 2004.**Exact frequency-domain reconstruction for thermoacoustic tomography–II: Cylindrical geometry.**

Y. Xu, M. Xu, and L. V. Wang.*IEEE Trans. Med. Imag.*, 21:829–833, 2002.**Ring-based ultrasonic virtual point detector with applications to photoacoustic tomography.**

X. Yang and L. V. Wang.*Applies Physics Letters*, 90(25):251103, 2007.**Circular integrating detectors in photo and thermoacoustic tomography.**

G. Zangerl, M. Haltmeier, and O. Scherzer.*Inverse Probl. Sci. Eng.*, 17(1):133–142, 2009.