Exact Series Reconstruction in Photoacoustic Tomography with Circular Integrating Detectors

Exact Series Reconstruction in Photoacoustic Tomography with Circular Integrating Detectors

Gerhard Zangerl Department of Mathematics University of Innsbruck, Technikerstraße 21a, 6020 Innsbruck (Gerhard.Zangerl@uibk.ac.at)    Otmar Scherzer Department of Mathematics University of Innsbruck, Technikerstraße 21a, 6020 Innsbruck, and Radon Institute of Computational and Applied Mathematics, Altenberger Strasße 69, 4040 Linz, Austria (Otmar.Scherzer@uibk.ac.at)    Markus Haltmeier Department of Mathematics, University of Innsbruck, Technikerstraße 21a, 6020 Innsbruck (Markus.Haltmeier@uibk.ac.at)
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.

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 non-ionizing 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 signal-to-noise 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 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.

Figure 1: Scanning geometry: A stack of circles centered on the boundary of is rotated around the -axis.

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 3.6).

The outline of this article is as follows. In 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.

2 Photoacoustic 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, 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, 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 (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.

Let and define and , by (2.4), (2.5). Then satisfies the axial symmetric wave equation

(2.7)
(2.8)
(2.9)

Moreover remains bounded as .

Proof.

Equations (2.8), (2.9) and the boundedness as immediately follow from (2.2), (2.3) 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 (2.7). ∎

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

Remark 2.2.

Proposition 2.1 is the basis of the following two-stage procedure which reconstructs the initial pressure in (2.1)-(2.3) from the data :

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

  2. 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 (back-propagation) 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., .

Proposition 2.3.

Let and define and by (2.5), (2.6). Then the relation

(2.11)

holds whenever .

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 .

Substituting in (2.12) and putting afterwards, leads to

(2.13)

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

Solving the last equation for shows (2.11). ∎

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

  1. The data are Fourier and cosine transformed, yielding to .

  2. According to (2.11), is mapped to .

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

Figure 2: Cross section of experimental buildup when (left) and (right). In both case object has to be supported in the gray disc.

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

In the following we fix , , define , by (2.5), (2.6), and let denote the zeros of the function .

Our first stable inversion formula is as follows:

Theorem 3.1 (The D’Hospital trick).

Assume . Then

(3.1)

for any .

Proof.

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

(3.2)

According to (2.11) we have

Applying the rule of D’Hospital gives

(3.3)

Inserting (3.3) in (3.2) and using the Fourier inversion formula shows (3.1). ∎

Remark 3.2.

[Stabilty of (3.1)] From the asymptotic approximation (see [1])

of the -th order Bessel function it follows that

Moreover the summands in (3.1) take the asymptotic form

Consequently, the parts that do not depend on the data are bounded, and (3.1) can be implemented in stable way.

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.

Let denote the unique solution of (2.1)-(2.3) and let and denote the Fourier coefficients of and

with respect to . Then

(3.4)

with denoting the -th order second kind Hankel function.


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

Theorem 3.4.

Assume . Then

(3.5)

Here is extended by for .

Proof.

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

(3.6)

Inserting (3.6) in (3.2) and using the Fourier inversion formula shows (3.5). ∎

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.

The derivation of (3.4) is based on the following Green function expansion in cylindrical coordinates [30]

(3.7)

with ,

and , denoting the -th order modified Bessel functions of first and second kind, respectively. Here its is assumed that .

Interchanging the roles of and implies that for the Green function expansion (3.7) holds with

Similar to the proof of (3.4) in [30] this leads to a formula for reconstructing in the case , however again an unstable one with in the denominator.

Figure 3: The first denominators in (3.10) for (left) and (right).
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 point-like 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 back-projection 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)
  • Finally, truncating the Fourier Bessel Series and approximating the inverse Fourier transform with the trapezoidal rule leads to discrete version of (3.1):

    (4.3)

    with in formula (4.3)

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 pre-computed 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 back-projection formula have to be applied times. Hence the total number of FLOPS is estimated as

(4.5)

Note that three dimensional back-projection type formulas which use point measurement data have complexity .

Figure 4: Left: Cross section of five absorbing spheres ( versus ). Right: The measurement data with Gaussian noise added ( versus ).
Figure 5: Reconstruction with (3.1) from simulated (left) and noisy data (right).
Figure 6: Reconstruction with (3.5) from simulated (left) and noisy data (right).

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 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, (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 S10505-N20. 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 quasi-reversibility 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 SPIE-OSA, 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 model-based 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 extreme-narrow 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 MHz-feasibility 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 X-ray 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. McGraw-Hill, 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]<