Circular integrating detectors in photo and thermoacoustic tomography

Abstract

Thermoacoustic computed tomography provides a novel and promising technology in medical imaging. It is based on the generation of acoustic waves inside tissues by stimulation with electromagnetic waves. Acoustic waves are measured outside the stimulated sample and converted into a three-dimensional (3D) image with high contrast and high resolution. Since optical fibres can be used to guide a laser beam on an arbitrary curve, it is possible to realize an experimental buildup in such a way that the acoustic pressure is measured over circles. The aim of this article is to present the mathematical theory when the acoustic pressure is measured over a stack of circles.

Keywords:

• Radon transform
• photoacoustic tomography
• Hankel transform
• image reconstruction
• integrating detectors

AMS Subject Classifications:

• 44A12
• 65R32
• 35L05
• 92C55

1. Introduction

Thermoacoustic computed tomography (TACT) is an emerging technology for non-invasive medical imaging. Based on the excitation of high bandwidth acoustic waves with pulsed non-ionizing electromagnetic energy 1–4 it combines the high contrast of optical imaging with the high resolution of ultrasonic imaging.

The common approach has been to use small conventional piezoelectric transducers 3,4 that approximate point-detectors. A point-detector measures the mean values of the absorption density function over spherical surfaces with the point-detector in their centre 4,5. Reconstruction algorithms that are based on point-detector data yield images with a spatial resolution essentially determined by the size and bandwidth of the used piezoelectric transducers 6. A prerequisite for achieving a high-resolution image is therefore that the detector size in all three dimensions should not exceed the desired minimum resolution. Small detectors with high bandwidth and sensitivity are technically hardly realizable 7,8. Different approaches have been introduced in 9,10, where the detector size is larger (in one or two dimensions) than the resolution limit.

Line-shaped detectors 9,11,12 integrate the acoustic pressure over its length and can be implemented by a Mach–Zehnder interferometer (which is an optical device where the the thermoacoustic pressure is measured with a laser beam). One drawback of the approach with linear detectors is that the lines are not aligned to the observed object, i.e. the points on a line move away from the object. Another disadvantage is that measuring the integrated pressure over laser beams around the object does not provide a simple and compact experimental buildup.

In this article, we propose to use integrals of the pressure over a stack of parallel closed circles instead of infinite lines. Optical fibres allow it to guide laser beams along curved paths. Therefore, circular integrating detectors can be implemented by fibre-based interferometers as acoustic detector 13. Since it is possible to fabricate optical fibres out of materials, which have nearly the same acoustical density like the surrounding liquid in which they are contained 14, there are little refraction effects expected.

The outline of this article is as follows. In Section 2, we present the proposed setup for PAT with circular integrating detectors. It turns out that the three-dimensional imaging problem reduces to a sequence of two-dimensional ones (similar as in 12). In Section 3 we derive an analytic reconstruction formula.

2. Photoacoustic tomography with circular detectors

Photoacoustic tomography (PAT) is based on the excitation of acoustic pressure waves p(x, t) inside an object of interest, when exposed to non-ionizing electromagnetic radiation 1,2,4,15. It combines the advantages of purely optical tomography (high contrast) and pure ultrasound imaging 16,17. It has demonstrated great promise for important medical applications including cancer diagnostics 18–22 and imaging of vasculature 23–26. For recent experimental buildups of PAT scanners, the reader is referred to 4,11.

Mathematically, the thermoacoustic pressure is described by the initial value problem (IVP) (1) (2) (3) (We assume that the speed of sound is spatially constant and equal to one.) Suppose the initial pressure f has compact support in some region Ω ⊂ ℝ³. The aim of PAT is to determine the initial pressure from the measurements p(x, t) outside Ω.

The classical approach with point-detector measurements leads to problem of reconstruction p from p(x, t), for (x, t) ∈ S × ℝ⁺, where S ⊂ ∂Ω. A different approach 9,11,12 investigates so-called integrating line detectors, which measure the overall pressure (4) on a line. This approach is motivated by the superior measuring accuracy of integrating line detectors, which can be realized by a thin laser beam that is part of a Fabry–Perot or Mach–Zehnder interferometer 11,13. By integrating (1) over the family of lines parallel to a given direction (i.e. taking the Radon transform) one obtains a wave equation in two dimensions. In this situation, the inverse problem of recovering f decomposes to a sequence of two-dimensional (2D) problems: first, the Radon transform of f (in a given direction) is reconstructed from values of the 2D wave equation on a curve, while in the the second stage a three-dimensional image is reconstructed by inverting the Radon transform.

A drawback of integrating line detectors compared to ultrasonic transducers is that they are not aligned to the object since points on the line ‘move away’ from the object quickly. Moreover, taking measurements with laser beams around the object does not lead to a compact experimental buildup. The purpose of this article is to combine the advantages of both technologies by using compact circular integrating detectors.

2.1. Circular integrating detectors

We propose the following experimental setup for circular integrating detectors (Figure 1):

Figure 1. Parallel circles with centre on z passing through C.

Consider a stack of parallel circles, which have possibly different radii and whose centres are on a line pointing in z-axis, which intersects the circle with radius R centred at the origin of a reference coordinate system. Then suppose the stack is rotated around an object along this circle and that the object is surrounded by the stacks. Any circle is assumed to collect integrals of pressure wave field over its length. This is the analog of the experimental buildup in 12, Section 2] for circular integrating detectors.

For our further considerations, we introduce the following notations: Let V the volume surrounded by the stacks and let σ ∈ S¹ denote the angle between the x′ axis of the reference coordinate system and the centre of the stack, see the picture in Figure 1(a). Then D_σ denotes the cross-section of V with the plane spanned by σ and the z-axis, and C_σ the intersection curve of this pane and the circles on the stack. Sometimes it will be convenient to use the notation {σ} × C for the set of points (σ, y₀) on C_σ. Finally, (5) denote the parallel circles that are centred on the z-axis (Figure 1(b)). Here y ≔ (r, z) ∈ ℝ², where z and r are the height and radius of the circles, respectively. Note that for each point on the curve C_σ there is a circle c_σ,y passing through it and that this circle is unique up to the sign of the radial variable r. For , let G_ev(z, r) ≔ G(z, |r|) and define the spaces of smooth compactly supported functions which are even in the r variable.

In the following, we shall prove that circular detectors yield a similar two-step reconstruction procedure as for integrating line detectors 2, Theorem 2.1]. Before we can state the corresponding theorem, we make the following definitions.

Definition 2.1 (Circular Radon transform)

By (6) we denote the circular Radon transform, which maps a function f to integrals over circles in the planes ℝ² × {z}, which centres on the circle with radius R and centre 0 in this plane.

Due to the definition of c_σ,y, ℛf is an even function in the r variable. Circular integrating detectors are modelled as a device performing the circular Radon transform of the pressure wave field for restricted radii. This leads to the following definitions:

Definition 2.2 (Forward operator)

The forward operator in cylindrical coordinates is (7) Here p is a solution of IVP (1–3).

Definition 2.3

The operator (8) maps a function to G ∈ , where G(σ, ·) is the restriction to C × [0, ∞) of the solution of (9) (10) (11)

Initial value problem (9–11) has a unique solution that can be given in an explicit form, which is due to the connection of (9–11) to the 2D wave equation and, in particular, to the existence of transformations that interwine the second derivative with the Bessel operator r⁻¹∂_rr∂r. The operator we are dealing with is the Weyl transformation (12) in the r variable. The Weyl transform is linear and satisfies the interwining property 27 (13)

This leads to the following proposition:

PROPOSITION 2.1

For a fixed σ ∈ S¹, the unique solution of (9–11) is given by (14) with

Proof

From the interwining property (13), it follows that P is a solution of (9–11), if and only if 𝒲₀P is a solution of the IVP (15) (16) (17) The unique solution of (15–17) is given by 𝒜𝒲₀F_σ, so we have that 𝒲₀ P = 𝒜𝒲₀F_σ and the claim follows.▪

Theorem 2.1 is a generalization in 2, Theorem 2.1].

THEOREM 2.1

The forward operator (7) decomposes into the product of the linear operators 𝒬 and ℛ, i.e., (18)

Proof

Let p be the solution of (1–3), fix σ ∈ S¹ and define (This is the same expression as in (7) but in this case y is not restricted to C_σ.) From the definitions of 𝒬 and P, it follows that we it is sufficient to show that 𝒫_σ satisfies (9–11) with F = ℛf.

First note that initial data (10) and (11) immediately follow from (2) and (3). It remains to show that 𝒫_σ solves the modified 2D wave Equation (9). By inserting the well-known expression for the Laplacian in cylindrical coordinates in (1) and integrating over α, we obtain The last term cancels because c_σ,y is a closed curve; therefore 𝒫_σ satisfies (9).▪

The decomposition of the original problem into a sequence of lower dimensional problems is the basis of the two-stage procedure for three-dimensional image reconstruction with circular integrating detectors. Before going into details we explain how it works in principle.

Figure 2. Reconstruction of the initial pressure from its circular Radon transform.

1.	For each of the directions σ ∈ S1 solve (19) That is, reconstruct the initial data Fσ in (9–11) from the restriction of the solution to Cσ and note that (20) is nothing else than the circular integrals of the initial pressure, see Figure 2.
2.	Again from Figure 2, it is clear that for a fixed z = z0, the map (21) is nothing other than the circular Radon transform (6) of the initial pressure f. Therefore, apply the (recently discovered) inversion formula for the circular Radon transform 28, theorem 1.1], which is given by (22) to (21) where x′ denotes the coordinates in the plane .

A novel analytic reconstruction formula for reconstructing F_σ in (19) is derived in the following Section 3.

3. Analytic inversion formula

The first step in the reconstruction algorithm is to solve Equation (19). We outline a strategy that uses the Fourier, Hankel and Cosine transformation to derive a reconstruction in the case when the circular stack is a cylinder. In this case, C_σ is a line consisting of the points (r₀, z). In the following computations, we write F(y) instead F_σ(y) in order to keep our notation simple.

It is well known that after going into Equation (9) with the separation ansatz (23) one obtains the ordinary differential equation (24) for φ. Its solution is given by the Bessel function φ(r) = J₀(rv) of order 0 and therefore, we end up with the family of solutions (25)

Now consider the function P(y, t) = P(r, z, t), which satisfies (9–11) and try to find a such that the equation holds. Because of the initial condition (10), i.e. at t = 0, this leads to which is the inverse Hankel Fourier transform of the function 29 and therefore the function is given by Putting all things together, the desired solution of (9–11) is given by (26) Since our goal is to determine the initial pressure from the measurements G(z, t) = P(r₀, z, t), it is desirable to eliminate the dependence of dv in (26) because the function g does not depend on r. Remember that ω² = k² + v² and therefore (27) Plugging this into Equation (26) and setting r = r₀, we have (28) After this change of the integration variables, we recognize that (28) can be written in the form where (29) (30) denote the Fourier and Cosine transform in the z and t variable, respectively. Using that 𝒞⁻¹ = 𝒞, gives (31)

Equation (31) is an explicit formula for the operator 𝒬. Solving (31) for , provided the denominator is not zero, finally gives (32)

Summarizing the computations above gives the following theorem:

THEOREM 3.1

Let and let G be the restriction of the solution of (9–11). Then G and F are related by (32).

4. Conclusion

In this article, we proposed a novel experimental buildup for PAT using circular integrating detectors. For collecting measurement data, a fibre-based 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 an exact reconstruction formula for circular integrating detectors.

The reconstruction formula in Section 3 needs further analysis to handle the effect of the zeros of Bessel functions in denominator. Future work will also be concerned with the numerical implementation of (32), which will investigate the behaviour of the derived reconstruction algorithms for noisy measurement data and 3D reconstruction with circular integrating detectors. There is considerable hope for stable implementation like in 30, since inversion formulas in 30 for integrating line detectors are quite similar.

Acknowledgements

This work has been supported by the Austrian Science Fund (FWF), Project Y123-INF and project P18172-N02.