Detectable singularities from dynamic Radon data

Detectable singularities from dynamic Radon data

B. N. Hahn111E-mail: Hahn@num.uni-sb.de   Department of Applied Mathematics, Saarland University, Saarbrücken, Germany E. T. Quinto222E-mail: Todd.Quinto@tufts.edu  Department of Applied Mathematics, Saarland University, Saarbrücken, Germany
Abstract

In this paper, we use microlocal analysis to understand what X-ray tomographic data acquisition does to singularities of an object which changes during the measuring process. Depending on the motion model, we study which singularities are detected by the measured data. In particular, this analysis shows that, due to the dynamic behavior, not all singularities might be detected, even if the radiation source performs a complete turn around the object. Thus, they cannot be expected to be (stably) visible in any reconstruction. On the other hand, singularities could be added (or masked) as well. To understand this precisely, we provide a characterization of visible and added singularities by analyzing the microlocal properties of the forward and reconstruction operators. We illustrate the characterization using numerical examples.

1 Introduction

The data collection in X-ray computerized tomography takes a certain amount of time due to the time-dependent rotation of the radiation source around the specimen. A crucial assumption in the classical mathematical theory (including modeling, analysis and derivation of reconstruction algorithms) is that the investigated object does not change during this time period. However, this assumption is violated in many applications, e.g. in medical imaging due to internal organ motion. In this case, the measured data suffer from inconsistencies. Especially, the application of standard reconstruction techniques leads to motion artifacts in the resulting images [39, 40].

Analytic reconstruction methods to compensate for these inconsistencies have been developed for special types of motion, including affine deformations, see e.g. [3, 36, 5]. An inversion formula for the dynamic forward operator in case of affine motion has been stated in [15], which also serves as basis for suitable reconstruction methods. For general, non-affine deformations, no inversion formula is known so far. Besides iterative methods, e.g. [2, 19], approximate inversion formulas that accurately reconstruct singularities exist for fan beam and parallel beam data in the plane [23] and for cone beam data in space [24]. They are based on the observation, that operators of the form

(1)

with forward operator , specially designed pseudodifferential operator and backprojection operator (which is, typically, related to the formal dual to ), are known to reconstruct singularities of the object. In addition, methods developed in the general context of dynamic inverse problems have been successfully applied in computerized tomography [38, 16].

Nevertheless, there can still arise artifacts in the reconstructions, even if the motion is known and the compensation method is exact, as e.g. [15]. On the other hand, the dynamic behavior of the object can lead to a limited data problem even if the radiation source rotates completely around the object. This means that some singularities of the object might not be visible in the reconstruction.

To guarantee reliable diagnostics in practice, it is essential to study these limitations carefully. Therefore, our aim is to analyse which singularities are detected by the measured data in the dynamic case and to characterize which of them can be reliably reconstructed or whether they create additional artifacts in the reconstruction process.

In this research, we understand the motion problem using generalized Radon transforms and microlocal analysis. The mathematical model of X-ray tomography with stationary specimen is integration along straight lines [28]. If the object moves during the data acquisition, the measured data can be interpreted as data for a (static) reference object where the integration now takes place along curves rather than straight lines [23, 2, 15]. Microlocal analysis is the rigorous theory of singularities and the study of how Fourier Integral Operators (FIO) transform them. Guillemin [12] was first to make the connection between microlocal analysis and Radon transforms (see also [14, 13]) when he showed that many generalized Radon transforms, , are FIO. He showed that, under the Bolker Assumption (Def. 2.9) and an extra smoothness assumption related to our definition of smoothly periodic (see Sect. 4.1), is an elliptic pseudodifferential operator (DO). This means that images all singularities of functions and does not add artifacts. This theorem was exploited in [1] to show that a broad range of Radon transforms on surfaces in can be “inverted” modulo lower order terms. Greenleaf and Uhlmann [11] and others developed the microlocal analysis of generalized Radon transforms that occur in X-ray CT [33, 26], cone beam CT [6, 21, 24], seismics [4], sonar [34], radar [31], and other applications in tomography.

Microlocal analysis has begun to be used in motion compensated CT. In [22], Katsevich proved that under certain completeness conditions on the motion model, the reconstruction operator in (1) detects all singularities of the object. This is related to theorems of Beylkin [1] showing that operators like are elliptic pseudodifferential operators. In [7] uniqueness is proven for a broad range of Radon transforms on curves. The cone beam CT case is more subtle since artifacts can be added to backprojection reconstructions, even with stationary objects [11, 6]. Katsevich characterized the added artifacts for this case and developed reconstruction algorithms to, at least locally, decrease the effect of those added artifacts. He uses this information to develop motion estimation algorithms when the motion model is not known [24].

Motivated by large field of view electron microscopy, the article [35] presents the microlocal analysis of general curvilinear Radon transform as well as local reconstruction methods. Analyzing added artifacts for X-ray tomography without motion has been done in [20, 8, 29], and generalizations to other types of tomography have been done in [9, 30].

In this article, we consider general motion models with less restrictive completeness assumptions. To develop our characterization of detectable and added singularities, we describe in Section 2 the mathematical model for the dynamic case as generalized Radon transform. We also present the mathematical bases of our work, including microlocal analysis. In Section 3, we assume the model is exact and study which object singularities are encoded in the measured data. In Section 4 we consider the reconstruction operator in the case of smoothly periodic motion, so the object is in the same state at the end of the scan as the start. Based on these results, in Section 5 we analyze the case when limited data arise and characterize visible and added singularities in reconstruction methods of type filtered backprojection. Our theoretical results are evaluated on numerical examples in Section 6. The more intricate proofs are in the appendix and we show in A.5 that our theorems are true even if the weights are arbitrary on the Radon transforms.

2 Mathematical basis

We use the following notation for function spaces. The space of all smooth (i.e., ) functions of compact support is denoted . A distribution is an element of the dual space with the weak-* topology and pointwise convergence (i.e., in if for every , in ). Further, will denote the set of smooth functions on ; its dual space, is the set of distributions that have compact support. See [37] for a description of the topologies and properties of these spaces.

A data set in computerized tomography can be interpreted as a function (or distribution) with domain . In the static case, the data are -periodic in the first variable, but this does not necessarily hold in the dynamic case since the object does not necessarily return to its initial state at the end of the scanning.

Generally, smooth functions (and hence distributions) are defined on open sets because derivatives will then be well defined. With this in mind, we make the following definition.

Definition 2.1

Let be a function with domain , where is an open subset of .

  1. We call smoothly periodic if extends to a smooth function on that is -periodic in the first variable.

  2. In the non-periodic case, we call smooth if, for some , extends to a smooth function on .

If is smoothly periodic, then can be viewed as a smooth function on the unit circle by identifying and . We define as the set of all smoothly periodic compactly supported functions on , and is its dual space with the weak-* topology. The set of smoothly periodic functions on , , and its dual space are defined in a similar way. Including the condition of -periodicity in these definitions will simplify the mapping properties of the dynamic forward operator and its dual (see Sect. 4.1).

In general, the object does not return to its initial state at the end of the scanning, i.e. its motion is not -periodic. For this case, we will state our theorems and definitions using open domains with for some . Finally, distributions can be restricted to open subsets and microlocal properties that hold on the larger set (e.g., smoothness) hold on the smaller set. So, our theorems are also true when mapping to distributions on (i.e., when the data are restricted to ) when is open.

In computerized tomography with stationary specimen, the given data correspond to integrals along straight lines of the distribution describing the x-ray attenuation coefficients of the investigated object. Hence, the mathematical model in the 2D parallel scanning geometry is given by the Radon line transform

(2)

with , and the -distribution. For fixed source and detector position , the integration takes place over the line

(3)

The data acquisition in computerized tomography is time-dependent, since the rotation of the radiation source around the object takes a certain amount of time. The source rotation is the only time-dependent part of the scanning procedure since, in modern CT scanners, detector panels are used such that all detector points record simultaneously for each fixed source position. Concerning the mathematical model, this means that a time instance can be uniquely identified with a source position and vice versa. In terms of the Radon transform, the source position is given by the angle , and there is the unique relation to a time instance via

with being the rotation angle of the radiation source. Therefore, throughout the paper, we interpret also as a time instance, and as time interval.

2.1 Mathematical model for moving objects in computerized tomography

We now derive the mathematical model for the case when the investigated object changes during the measuring process. A dynamic object is described by a time-dependent function . In the application of computerized tomography, for a fixed time corresponds to the x-ray attenuation coefficient of the specimen at this particular time instance.

The dynamic behavior of the object is considered to be due to particles which change position in a fixed coordinate system of . This physical interpretation of object movement is now incorporated into a mathematical model.

Let denote the state of the object at the initial time. We call a reference function. Please note that is a distribution since . Further, let be a motion model describing the dynamic behavior of the specimen, where and denotes which particle is at position at the time instance (in other words, is the location at time of the particle that is at at time ). For fixed , we write

(4)

to simplify the notation. Using this motion model and the reference function , the state of the object at time instance is given by

(5)
Remark 2.2

In the model (5), each particle keeps its initial intensity over time. However, this means that the mass of the object may no longer be conserved. If the density varies due to the deformation, this can be taken into account by the mathematical model

(6)

In both cases, the respective Fourier Integral Operators describing the dynamic setting have the same phase function and hence the same canonical relation. Thus, our results provided in this paper hold equivalently for the mass preserving model (6), see also A.5.

Throughout the paper, we make the following assumptions on the motion model , which we justify by physical properties of moving objects and imaging systems.

Hypothesis 2.3

Let and let . Assume . Then, is called a motion model and a motion function if there is an such that

  1. extends smoothly to (so is smooth by Def. 2.1).

  2. For each is a diffeomorphism.

A motion model is smoothly periodic if it satisfies these conditions for some and is smoothly periodic.

Remark 2.4
  1. Note that if the motion model is smoothly periodic and satisfies (1) and (2) of this hypothesis for some , then it does for any because is -periodic in in this case.

  2. Under these hypotheses, the trajectory of a fixed particle, which is the map

    is a smooth curve.
    In practical applications of computerized tomography, only discrete data are measured. Thus, the object’s motion is ascertained for finitely discrete time instances only, which justifies this (theoretical) assumption of smooth trajectories.

  3. Hypothesis (2) ensures that two particles cannot move into the same position, and no particle gets lost (or added). The relocation is smooth because is a smooth function.

With the mathematical model of a dynamic object (5), the operator of the dynamic setting is given by

(7)

Using the change of coordinates , we obtain the representation

(8)

Thus, integrates the respective intensity-corrected reference function along the curve

(9)

So, for each , . Because is a diffeomorphism, each is a smooth simple unbounded curve, and for each , the curves for cover the plane and they are mutually disjoint (they foliate the plane).

2.2 Microlocal analysis and Fourier integral operators

In this section we will outline the basic microlocal principles used in the article. We refer to [17, 41, 42, 18, 25] for more details.

The key to understanding singularities and wavefront sets is the relation between smoothness and the Fourier transform: a distribution is smooth if and only if its Fourier transform is rapidly decreasing at infinity. However, to make the definition invariant on manifolds (such as with and identified), we need to define the wavefront set as a set in the cotangent bundle [41]. So, we will introduce some notation.

Let and . Now let be a smooth scalar function of variables including and let be a function with codomain , then we define

where is the cotangent space at ,

and the other derivatives (using ) and differentials (using ) are defined in a similar way; for example, .

Definition 2.5

Let and let . Then is smooth at in direction if there is a cutoff function at , (i.e., ) and an open cone containing such that is rapidly decreasing at infinity for all .

On the other hand, if is not smooth at in direction , then , the wavefront set of .

We now define the fundamental class of operators on which our analysis is based: Fourier integral operators. Note that we define them only for the special case we use. For other applications, one would use the definition for general spaces in [42, Chapter VI.2] or [17].

Definition 2.6 (Fourier Integral Operator (FIO))

Let . Now let be a smooth function on , then is an amplitude of order if it satisfies the following condition. For each compact subset in and , there exists a positive constant such that

(10)

for , and all and all .

The real-valued function is called a phase function if is positive homogeneous of degree in and both and are nonzero for all . The phase function is called non-degenerate if on the zero-set

(11)

one has that . In this case, the operator defined for by

(12)

is a Fourier Integral Operator (FIO) of order . The canonical relation for is

(13)

Note that since the phase function is non-degenerate, the sets and are smooth manifolds. Because of the conditions on and , and is continuous in both cases [42]. If the amplitude and phase function are smoothly periodic, then the conditions in this definitions are valid on where and are identified. In this case, is periodic in for all .

To state the theorems that form the key to our proofs, we need the following definitions. Let and be sets and let , , and . Then,

(14)

We will use these rules for sets of cotangent vectors to calculate wavefront sets.

Theorem 2.7 ([17, Theorem 4.2.1])

Let be an FIO with canonical relation . Then the formal dual operator, to is an FIO with canonical relation .

FIO transform wavefront sets in precise ways, and our next theorem, a special case of the Hörmander Sato Lemma, is a key to our analysis.

Theorem 2.8 ([17, Theorems 2.5.7 and 2.5.14])

Let be an FIO with canonical relation . Let . Then .

To understand the more subtle properties of FIO, we investigate the mapping properties of the canonical relation . Let and be the natural projections. Then we have the following diagram:

(15)

First, note that if and then

(16)

These statements are proven using the definitions of composition and the projections.

Definition 2.9

Let be an FIO with canonical relation . Then, satisfies the Bolker Assumption if the projection is an injective immersion.

Recall that an immersion is a smooth map with injective differential. Victor Guillemin [12, 14] named this assumption after Ethan Bolker who gave a similar assumption for finite Radon transforms.

Definition 2.10

The FIO in (12) is elliptic of order if its amplitude, , is of order and satisfies, for each compact set there are constants and such that for all and , .

Now, we apply these ideas to dynamic tomography.

3 Microlocal analysis of the dynamic forward operator

In this section, we study the microlocal properties of the forward operator in dynamic computerized tomography. We show that it is an FIO and provide conditions under which it fulfills the Bolker Assumption. Theorem 3.6 gives the relationship between singularities of and those of which is then analysed in more detail, especially with respect to the importance of the Bolker Assumption. Our theorems are true for more general FIO, but the proofs are easier in our special case.

We now introduce some notation and describe its geometric meaning. Here is a motion model that satisfies Hyp. 2.3 and let be as in that hypothesis. For define

(17)

Then, the integration curve in (9) can be written

Now, define

(18)

Our next lemma gives the geometric meaning of this covector.

Lemma 3.1

Let and let be a point on the integration curve . The vector is normal the curve at , and therefore the covector is conormal to this curve at .

Proof: The curve is defined by the equation . Therefore the gradient in of at each , which is , is normal to this curve at . So, its dual covector, which is , is conormal to at (i.e., in the conormal space of above ).  

3.1 The canonical relation of

We first prove that the forward operator (8) for the dynamic setting is an elliptic FIO.

Theorem 3.2

Under Hypothesis 2.3, the operator is an elliptic FIO of order with phase function

(19)

and amplitude

(20)

which is elliptic of order zero.

The proof is given in the appendix A.1.

Since is an FIO, we can determine its canonical relation using Definition 2.6, eq. (13).

Lemma 3.3

Under Hypothesis 2.3, the canonical relation of is

(21)

where is as given in Hypothesis 2.3.

If the motion model is smoothly periodic in then the condition on in (21) is replaced by and is still a smooth manifold without boundary when is identified with the unit circle, .

Proof: According to Definition 2.6, (13), the canonical relation of is given by

where . Using the representation of the phase function (19) along with (17), , and thus if . The representation of then follows from the representation of the differentials and , as noted in the proof of Theorem 3.2. 

In the following theorem, we find conditions on the motion model under which satisfies the Bolker Assumption.

Theorem 3.4

Assume the motion model satisfies Hypothesis 2.3.

  1. If, for each , the map

    (22)

    is one-to-one, then is injective.

  2. If the motion model fulfills the condition

    (23)

    for all , then the projection is an immersion.

Thus, under these two conditions, satisfies the Bolker Assumption (Definition 2.9).

If the motion is smoothly periodic, then can be replaced by in this theorem.

To illustrate the geometric meaning of condition (22) for the motion model, we assume there exist two points and with and for some . The first equality implies that the two points are on the same integration curve, i.e. the data at angle cannot distinguish between them. The second equality means, if the angle of view changes infinitesimally, also the new curve cannot distinguish the two points because they both stay on the same curve (at least infinitesimally). An example for a motion model not satisfying (22), is any dynamic behavior, where two particles, which are on the same integration curve for a time instance , are rotated with the same speed and in the same direction as the radiation source

Condition (23), also referred to as an immersion condition, is equivalent to the condition

The property means that, at least infinitesimally at , the line normal to the curve
at is stationary at , i.e. the curves near are infinitesimally rigid at (these statements are justified in a related case in [35, Remarks 2 and 5]).

We should remark that the conditions in Theorem 3.4 are essentially equivalent to the conditions of Theorem 2.1 in [22] for the fan beam case. There is an additional assumption in his theorem that ensures that all singularities are visible.

Proof of Theorem 3.4 On the set , we introduce global coordinates by the map

(24)

In these coordinates, the projection is given by

(25)

Using the representation (25) of , one sees that is injective if for each , the map in (22) is injective.

The map is an immersion if its differential has constant rank , and in coordinates

Thus, condition (23) is equivalent to for all , and thus is equivalent to, being an immersion. 

The importance of this Bolker Assumption for the detection of object singularities in dynamic Radon data is discussed in the next section.

3.2 Visible Singularities

Now, we classify singularities of functions that appear in the data, both algebraically and geometrically.

Theorem 3.5

Assume the motion model, , satisfies Hypothesis 2.3. Let . Then,

(26)

Now assume, in addition, that satisfies the Bolker Assumption. Then,

(27)

We will prove this theorem in the appendix, §A.2.

The explicit correspondence between object and data singularities is given in the following corollary.

Corollary 3.6

Let , and let be a motion model satisfying Hypothesis 2.3. Let be an open subset of and let , , .

If then there is an such that

where is the integration curve given by (9) and is given by (18).

Now assume in addition satisfies the Bolker Assumption. For ,

(28)

Furthermore, if such a point exists, then it is unique.

The proof follows immediately from Theorem 3.5 and the expression (21) for the canonical relation . In particular, the first statement follows from (26), and the equivalence (28) follows from (27).

For define

(29)

Corollary 3.6 justifies our next definition.

Definition 3.7

Let and let be a motion model satisfying Hypothesis 2.3. Assume the associated Radon transform, , satisfies the Bolker Assumption. Let and let . Then, is a visible singularity from data above if has the representation

(30)

for some and .

We call

(31)

the set of all possible visible singularities from above . Covectors in

will be called invisible singularities from .

Using (16), it follows that

(32)

Corollary 3.6 justifies the definition: if the motion model satisfies Hypothesis 2.3 and satisfies the Bolker Assumption, then a singularity causes a singularity from the data above the open set (i.e., in ) if and only if it is in . The singularities of that are in are smoothed by . Note that the singularities of in are problematic because we will show they are in directions that can be added singularities or that can be visible or masked by added singularities.

We can now describe the geometric meaning of the visible singularities.

Corollary 3.8

Let the motion model fulfill the Bolker Assumption. The dynamic operator detects a singularity of at a point in direction if and only if there is an integration curve passing through with conormal to the curve at (i.e., the curve has tangent line at this point that is normal to ).

Proof: Let . Corollary 3.6 shows that, under the Bolker Assumption, a singularity of at is visible if and only if for some . Furthermore, Lemma 3.1 establishes that for each and each , the covector is conormal to at . Thus a singularity of at is visible if and only if is conormal to at . 

Remark 3.9

In general, each data singularity at a point in data space, , stems from an object singularity with direction , where is perpendicular to the curve at . However, in case the Bolker Assumption is not fulfilled by the motion model, two object singularities could cancel in the data and thus, not lead to a corresponding data singularity.
In contrast, under the Bolker Assumption, every singularity in the data comes from a singularity in the object. Note that Example 3.11 shows that not all singularities of the object necessarily show up in the data.

Another way to understand visible singularities is the following. if there is some and , such that , where is the map

(33)

for (see (30)). If this map is not injective, the object singularity can cause two different data singularities, resulting in redundant data, as illustrated by our next example.

Example 3.10

Let the dynamic behavior of be given by the rotation with rotation matrix

This describes an object which rotates in the opposite direction as the radiation source with the same rotational speed. In particular, it holds for , so this is a smoothly periodic motion model. Since is a unitary matrix for all , it is

By a calculation using it’s definition, , and the map

is one-to-one since the matrix is nonsingular. Thus, the dynamic operator satisfies the Bolker Assumption, and .

Now, let with . Since it holds that

as well as

this one singularity in object space causes two singularities

This is according to the fact that the projection is not injective due to the motion introduced data redundancy.

If the map in (33) is surjective for all then all singularities and all directions are gathered in the measured data, and we speak of complete data. In the static case, this corresponds to the fact that the radiation source rotates around the complete circle (e.g., [32]). If is not surjective, when the point is only probed by data from a limited angular range. The following example illustrates that the dynamic behavior of the object can lead to incomplete data, even if the full angular range is covered by the source.

Example 3.11

We consider the rotational movement with

In this setting, the object rotates in the same direction as the radiation source with half of its rotation speed. In particular, this is a non-periodic motion model. It is

One shows the injectivity condition, (22), is fulfilled in the same way as in Example 3.10. Computing the derivatives, we obtain . So, the Bolker Assumption holds.

Now, assume with . According to Theorem 3.5, a corresponding singularity is seen in the data if there exists an angle with , or . Since for all , an angle with the required property does not exist. Hence, the singularity cannot be seen in the data.

4 The dynamic reconstruction operator for smoothly periodic motion

In this section, we prove the main theorem for smoothly periodic motion. Basically, under our assumptions, the reconstruction operator is well behaved and reconstructs all singularities of the object without introducing new artifacts. First, we define the backprojection operator.

4.1 Backprojection for Smoothly Periodic Motion

In general, we denote the backprojection operator by and define it as

(34)

Note that, for smoothly periodic motion, this backprojection operator is the formal dual, , to for . A generalization to arbitrary weights is explained in section A.5.

Proposition 4.1

If the motion model is smoothly periodic, then the backprojection operator, , can be composed with for and, if is a pseudodifferential operator, then the reconstruction operator

is defined and continuous on domain .

Proof: The proof will now be outlined. First, we show when ,