Artifacts from Arbitrary Incomplete CT Data

Analyzing Reconstruction Artifacts from Arbitrary Incomplete X-ray CT Data

Leise Borg Jürgen Frikel Jakob Sauer Jørgensen  and  Eric Todd Quinto

This article provides a mathematical analysis of artifacts from arbitrary incomplete X-ray computed tomography (CT) data when using the classical filtered backprojection algorithm. We prove that artifacts arise from points at the boundary of the data set. Our results show that, depending on the geometry of this boundary, two types of artifacts can arise: object-dependent and object-independent artifacts. The object-dependent artifacts are generated by singularities of the object being scanned and these artifacts can extend all along lines. They generalize the streak artifacts observed in limited-angle tomography. The article also characterizes two new types of artifacts that are essentially independent of the object; they occur along lines if the boundary of the data set is not smooth at a point and along curves if the boundary is smooth locally. In addition to the geometric description of artifacts, the article provides characterizations of their strength in Sobolev scale in certain cases. The results of this article apply to the well-known incomplete data problems, including limited-angle and region-of-interest tomography, as well as to unconventional X-ray CT imaging setups. Reconstructions from simulated and real data are analyzed to illustrate our theorems, including the reconstruction that motivated this work–a synchrotron data set in which artifacts appear along lines that have no relation to the object.

11footnotetext: Department of Computer Science, University of Copenhagen ( Department of Computer Science and Mathematics, OTH Regensburg, ( Her work was partially funded by: Innovation Fund Denmark and Maersk Oil and Gas 33footnotetext: School of Mathematics, University of Manchester, Manchester, M13 9PL, United Kingdom ( HC Ørsted Postdoc programme, co-funded by Marie Curie Actions at the Technical University of Denmark (Frikel); ERC Advanced Grant No. 291405 and EPSRC Grant EP/P02226X/1 44footnotetext: Department of Mathematics, Tufts University, Medford, MA 02155 USA, ( partial support from NSF grants DMS 1311558 and 1712207, as well as support from the Otto Mønsteds Fond, DTU, and Tufts University Faculty Research Awards Committee

Keywords: X-ray tomography, incomplete data tomography, limited angle tomography, region of interest tomography, reconstruction artifact, wavefront set, microlocal analysis, Fourier integral operators

2010 AMS subject clasifications: 44A12, 92C55, 35S30, 58J40

1. Introduction

Over the past decades computed tomography (CT) has established itself as a standard imaging technique in many areas, including materials science and medical imaging. Its principle is based on collecting numerous X-ray measurements of the object along all possible directions (lines) and then reconstructing the interior of the object by using an appropriate mathematical algorithm. In classical tomographic imaging setups, this procedure works very well because the data can be collected all around the object (i.e., the data are complete) and standard reconstruction algorithms, such as filtered backprojection (FBP), provide accurate reconstructions [29, 36]. However, in many recent applications of CT, some data are not available, and this leads to incomplete (or limited) data sets. The reasons for data incompleteness might be health-related (e.g., to decrease dose) or practical (e.g., when the scanner cannot image all of the object as in digital breast tomosynthesis).

Classical incomplete data problems have been studied from the beginning of tomography, including limited-angle tomography, where the data can be collected only from certain view-angles [27, 22]; interior or region-of-interest (ROI) tomography, where the X-ray measurements are available only over lines intersecting a certain subregion of the object [9]; or exterior tomography, where only measurements over all lines outside a certain subregion are available [28, 41].

In addition, new scanning methods generate novel data sets, such as the synchrotron experiment [6, 5] in Section 7 that motivated this research. That reconstruction, in Figure 1, includes dramatic streaks that are independent of the object and were not described in the mathematical theory at that time but are explained by our main theorems. A thorough practical investigation of this particular problem was recently presented in [5].

Regardless of the type of data incompleteness, in most practical CT problems a variant of FBP is used on the incomplete data to produce reconstructions [36]. It is well-known that incomplete data reconstruction problems that do not incorporate a priori information (as is the case in all FBP type reconstructions) are severely ill-posed. Consequently, certain image features cannot be reconstructed reliably [40] and, in general, artifacts111Artifacts are image features (singularities), e.g. streaks or points, that are added by the algorithm and that are not part of the original object., such as the limited-angle streaks in Figure 2 can occur. Therefore, reconstruction quality suffers considerably, and this complicates the proper interpretation of images.

Figure 1. Left: A small part of the sinogram of the chalk sample analyzed in Section 7. Notice the rough boundary. Right: Small central section of a reconstruction of the chalk. Monochromatic parallel beam data were taken of the entire cross section of the chalk over views covering degrees, and there were detector elements with a mm field of view, providing micrometer resolution of the sample. Data [51] obtained, with thanks from the Japan Synchrotron Radiation Research Institute from beam time on beamline BL20XU of SPring-8 (Proposal 2015A1147). For more details, see Section 7 and [5, ©IOP Publishing. Reproduced by permission of IOP Publishing. All rights reserved].

1.1. Related research in the mathematical literature

Our work is based on microlocal analysis, a deep theory that describes how singularities are transformed by Fourier integral operators, such as the X-ray transform. In the 1990s microlocal analysis was used to characterize visible singularities from X-ray CT data and their strength was given in [40, 24]. Subsequently, artifacts were extensively studied in the context of limited-angle tomography, e.g., [22, 12]. The strength of added artifacts in limited-angle tomography was analyzed in [31]. Similar characterizations of artifacts in limited-angle type reconstructions have also been derived for the generalized Radon line and hyperplane transforms as well as for other Radon transforms (such as circular and spherical Radon transform), see [13, 14, 32, 1, 33].

Metal in objects can corrupt CT data and create dramatic streak artifacts [3]. This can be dealt with as an incomplete data problem by excluding data over lines through the metal. Recently, this problem has been mathematically modeled in a sophisticated way using microlocal analysis in [37, 44, 35], all of which are in the spirit of our article. A related problem is studied in [34], where the authors develop a streak reduction method for quantitative susceptibility mapping (see also [7]). Moreover, microlocal analysis has been used to characterize properties of related integral transforms in pure and applied settings [4, 15, 47, 10, 43].

To the best of our knowledge, the mathematical literature up until now, e.g., [22, 31, 12, 13, 33] used microlocal and functional analysis to explain streak artifacts along lines that are generated by singularities of the object, and they were for specific problems, primarily limited-angle tomography. Even classical setups, such as region-of-interest tomography, had not yet been thoroughly explored microlocally, not to mention nonstandard imaging setups such as the one presented in Figure 1.

1.2. Basic mathematical setup and our results

In this article, we present a unified approach to characterize reconstruction artifacts for arbitrary incomplete X-ray CT data that are caused by the choice of data set. We not only consider all of the above mentioned classical incomplete data problems but also emerging imaging situations with incomplete data. To this end, we employ tools from microlocal analysis and the calculus of Fourier integral operators.

If is the density of the object to be reconstructed, then each CT measurement is modeled by a line integral of over a line in the data set. As we will fully describe in Section 2.1, we parametrize lines by , and the CT measurement of over the line is denoted . With complete data, where is given over all , accurate reconstructions can be produced by the FBP algorithm. In incomplete data CT problems, the data are taken over lines for in a proper subset, , of and, even though FBP is designed for complete data, it is still one of the preferred reconstruction methods in practice. As a result, incomplete data CT reconstructions usually suffer from artifacts.

We prove that incomplete data artifacts arise from points at the boundary or “edge” of the data set, , and we show that there are two types of artifacts: object-dependent and object-independent artifacts. The object-dependent artifacts are caused by singularities of the object being scanned. In this case, artifacts can spread all along a line (i.e., a streak) if and if there is a singularity of the object along the line (such as a jump or object boundary tangent to the line)—this singularity of the object “generates” the artifact. The streak artifacts observed in limited-angle tomography are special cases of this type of artifact.

In addition, we characterize two new phenomena that are (essentially) independent of the object being scanned. If the boundary of is smooth near a point , then we prove that artifacts can appear in the reconstruction along curves generated by near , and they occur whether the object being scanned has singularities or not. We also prove if is not smooth (see Definition 3.2) at a point then, independent of the object, an artifact line can be generated all along .

We will illustrate our results using classical problems including limited-angle tomography and region-of-interest tomography, as well as problems with novel data sets, including the synchrotron data set in Figure 1. In addition to the geometric characterization of artifacts, we also provide a characterization of their strength in Sobolev scale, and we review a general method to suppress the artifacts.

1.3. Organization of the article

In Section 2, we provide notation and some of the basic ideas about distributions and wavefront sets. In Section 3 we give our main theoretical results, and in Section 4, we apply them to classical examples to explain added artifacts. In Section 5, we describe the strength of added artifacts in Sobolev scale. Then, in Section 6, we describe a straightforward way to decrease the added artifacts. We provide more details of the synchrotron experiment in Section 7 and observations and generalizations in Section 8. Finally, in the appendix, we give some technical theorems and then prove the main theorems.

2. Mathematical basis

Much of our theory can be made rigorous for distributions of compact support (see [11, 45] for an overview of distributions), but we will consider only functions in , the set of absolutely square integrable functions supported in the open unit disk . This setup is realistic in practice, and our theorems are simpler in this case than for general distributions. Remark A.4 provides perspective on this.

2.1. Notation

For , we define


so is the unit vector in the direction of and is the unit vector radians counterclockwise from . The line perpendicular to and containing is denoted


This line is parameterized by . Note that .

We define the X-ray transform or Radon line transform of to be


For functions on , the dual Radon transform or backprojection operator is defined


2.2. Wavefront sets

In this section, we define some important concepts needed to describe singularities in general. Sources, such as [11], provide introductions to microlocal analysis. Generally cotangent spaces are used to describe microlocal ideas, but they would complicate this exposition, so we will identify a covector with the associated ordered pair of vectors . The book chapter [23] provides some basic microlocal ideas and a more elementary exposition adapted for tomography.

The wavefront set is a deep concept that makes singularities of functions precise, and we take the definition from [11].

Definition 2.1 (Wavefront set).

Let and . A cutoff function at is any -function of compact support that is nonzero at . Let be a locally integrable function. We say is smooth at in direction if there is a cutoff function at and an open cone containing such that the Fourier transform is rapidly decaying at infinity.222That is, for every , there is a constant such that for all .

If is not smooth at in direction , then we say has a singularity at in direction . The wavefront set of is the set of all such .

For example, if is a subset of the plane with a smooth boundary, and is equal to on and off of , then is the set of all points where is in the boundary of and is normal to the boundary at .

We parametrize lines on , rather than , because, in practice, the data space consists of points and the sinogram is represented as a picture in the -plane, not on a cylinder. We will consider only distributions on that are restrictions of distributions for that are -periodic in . The Radon transform and our other operators on are all -periodic in , and so this is no real restriction. For this reason, we can use the same definition of wavefront set as for , and we just apply it to the -periodic extensions.

Our next definition allows us to express efficiently the correspondence between singularities of the object and of its Radon transform.

Definition 2.2.

Let . The 6 normal space of the line is


and the set of singularities of normal to is


Let be a locally integrable function on . We say is smooth normal to the line if .

For , we let

Now, let be a locally integrable function on and . We define


It is important to understand each set introduced in Definition 2.2: is the set of all such that and the vector is normal to at . Therefore, is the set of wavefront directions of above points that are normal to this line.

If is a locally integrable function on , then is the set of wavefront directions with base point . The sets defined in Definition 2.2 have an important relationship that we will exploit starting in the next section.

3. Main results

In contrast to limited-angle characterizations in [12], our main results characterize arbitrary incomplete data reconstructions, including classical problems such as region-of-interest tomography, exterior tomography, and limited-angle tomography. Our results are formulated in terms of the wavefront set (Definition 2.1).

In many applications, reconstructions from incomplete CT data are calculated by the filtered backprojection algorithm (FBP), which is designed for complete data (see [36] for a practical discussion of FBP in practice). In this case, the incomplete data is often implicitly extended by the algorithm to a complete data set on by setting it to zero off of the set (cutoff region) over which data are taken. Therefore, the incomplete CT data can be modeled as


where is the characteristic function of .333The characteristic function of a set is the function that is equal to one on the set and zero outside of . Thus, using the FBP algorithm to calculate a reconstruction from such data gives rise to the reconstruction operator:


where is the standard FBP filter (see e.g., [29, Theorem 2.5] and [30, §5.1.1] for numerical implementations) and is defined by (2.4).

Our next assumption collects the conditions we will impose on the cutoff region . There, we will use the notation , , and to denote the interior of , the boundary of , and the exterior of , respectively.

Assumption 3.1.

Let be a proper subset of (i.e., ) with a nontrivial interior and assume is symmetric in the following sense:


In addition, assume that is the smallest closed set containing , i.e. .

We now justify using this assumption. Since is proper, data over are incomplete. Being symmetric means that the part of in determines the data set. We assume that is the smallest closed set containing in order to exclude degenerate cases, such as when includes an isolated curve.

Our next definition gives us the language to describe the geometry of .

Definition 3.2 (Smoothness of ).

Let and let .

  • We say that is smooth near if this boundary is a curve near . In this case, there is a unique tangent line in -space to at .

    • If this tangent line is vertical (i.e., of the form ), then we say the boundary has infinite slope at .

    • If this tangent line is not vertical, then is defined near by a smooth function . In this case, the slope of the boundary at will be the slope of this tangent line, which is given by .

  • We say that is not smooth at if it is not a smooth curve at .

    • We say that has a corner at if the curve is continuous at , is smooth at all other points in some neighborhood of , and has one-sided slopes at but those slopes are different.

3.1. Basic theorems on artifacts and visible singularities

In this section we begin to characterize singularities and artifacts, and we show that artifacts occur only on lines for . This material is either known (e.g., [40]) or it follows from [40].

Definition 3.3 (Artifacts and visible singularities).

Every singularity of that is not a singularity of is called an artifact (i.e., any singularity in ). A streak artifact is a line of artifacts.

Every singularity of that is a singularity of is called a visible singularity (from data on ) (i.e., any singularity in ). Other singularities of are called invisible (from data on ).444Invisible singularities are smoothed by and reconstruction of those singularities is in general is exponentially ill-posed in Sobolev scale.

Theorem 3.4 (Localization of artifacts to lines from ).

Let satisfy Assumption 3.1. If , then there is a , such that and is normal to . That is, if a singularity of at is an artifact, then it, necessarily, comes from .

This theorem follows from the microlocal analysis in [40, 42]. More generally, every singularity of must come from a singularity of , and if the singularity does not come from , it must come from a singularity of , and these singularities are on .

Our next theorem gives an analysis of singularities in corresponding to lines not in . It shows that visible singularities are along lines when and invisible singularities are on lines for .555Visible singularities of can appear on lines for , but these can cause artifacts all along the line as noted in Theorem 3.6 and seen in the reconstructions in Section 4. In addition, artifacts along such lines can mask visible singularities on those lines.

Theorem 3.5 (Visible and invisible singularities in the reconstruction).

Let and let satisfy Assumption 3.1.

  1. If then . That is, all singularities of normal to are visible from data on .

  2. If , then . In this case, is smooth normal to . That is, all singularities of normal to are invisible from data on .

  3. Now, let and assume that all lines through are parameterized by points in (i.e., , ). Then all singularities of at are visible in :


Therefore, artifacts occur only on lines in .

This theorem follows directly from [40, Theorem 3.1] or [42], and the proof will be sketched in Section A.2 since it is not exactly given in the literature.

3.2. Characterization of artifacts

We now characterize artifacts in FBP reconstructions from incomplete data given by (3.2). In particular, we show that the nature of artifacts depends on the smoothness and geometry of and, in some cases, singularities of the object . Our next two theorems show that artifacts arise in two ways:

  • Object-independent artifacts: those caused essentially by the geometry of .

  • Object-dependent artifacts: those caused by singularities of the object that are normal to for some .

As a first generalization of the limited-angle characterizations in [22, 12] we consider the case where is locally smooth. This theorem describes the range of artifacts that can occur in this case, and it guarantees they will occur in certain cases.

Theorem 3.6 (Artifacts for locally smooth boundary).

Let and satisfy Assumption 3.1. Let and assume that is smooth near .

  1. Assume has finite slope at and let be a neighborhood of such that is given by a smooth curve near . Let

    1. An object-independent artifact can occur along the curve .

    2. If has no singularity normal to then


      That is, is smooth normal to , except perhaps at the single point .

      1. If, in addition, , then equality holds in (3.6) and the artifact curve will be visible near .

      2. If is zero in a neighborhood of , then the artifact curve will not be visible near .

  2. In all cases with smooth boundary at , if has a singularity normal to , then an object-dependent streak artifact can occur along .

    1. If is smooth conormal to and is vertical at , then is smooth normal to .

    2. If is not vertical then is smooth normal to except possibly at .

Theorem 3.6 part 2 gives a necessary condition under which there can be an object-dependent artifact. This is observed in practice (see Figure 2). However, Remark A.6 shows such object-dependent streak artifacts might not occur, even if there is a singularity of normal to this line. This is why we state that artifacts can occur in parts 11a and 2 of Theorem 3.6.

The proof is given in Section A.3 and we will present reconstructions illustrating this theorem in Section 4.2.

Remark 3.7.

Assume is smooth at with finite slope at . Let I be a neighborhood of and describes near . Note that

If the slope of at is small enough, i.e.,


holds, then the artifact curve will be inside the unit disk, , near . If not, this curve will be outside, , at least at .

If is smooth and vertical at (infinite slope), then there will be no object-independent artifact on the line because is a “point at infinity” since “” is infinite. This is proven using the arguments in the proof of Theorem 3.61b

When is not smooth at , the next theorem shows there can be a streak artifact in along that is independent of the object .

Theorem 3.8 (Artifacts for nonsmooth boundary).

Let and satisfy Assumption 3.1. Let and assume that is not smooth at . Then, can have a streak artifact on independent of .

If is smooth normal to , and has a corner at (see Definition 3.2), then does have a streak artifact all along :

This theorem is proven in Section A.3.

The streak artifacts (artifacts along lines) in this theorem are object-independent, and they are illustrated in a simple case in Section 4.4 and on real data in Section 7.

In Theorem 5.2, we will describe the strength of the artifacts in Sobolev scale in specific cases of Theorems 3.6 and 3.8.

Object-dependent streak artifacts along lines were analyzed for limited-angle tomography in articles such as [22, 12, 31], but we are unaware of a general reference to Theorem 3.6, part 2 for general incomplete data problems, although this is not surprising given the limited-angle result. We are not aware of a previous reference in the mathematical literature to the curve artifact in Theorem 3.6 part 1 or to the corner artifacts from Theorem 3.8.

4. Numerical illustrations of our theoretical results

We now consider a range of well-known incomplete data problems as well as unconventional ones to show how the theoretical results in Section 3 are reflected in practice. All sinograms have on the horizontal axis and on the vertical. Except for the center picture in Figure 2(A), the reconstruction is displayed on .

4.1. Limited-angle tomography

First, we analyze limited-angle tomography, a classical problem in which part 2 of Theorem 3.6 applies. In this case consists of vertical lines given by (infinite slope). Taking a closer look at the statement of Theorem 3.6 and the results of [12, 14] one can observe that, locally, they describe the same phenomena, namely: whenever there is a line in the data set with and which is normal to a singularity of , then a streak artifact can be generated along in the reconstruction . Therefore, Theorem 3.6 generalizes the results of [22, 12] as it also apply to cutoffs with non-vertical slope.

It is important to note that, with limited-angle data, there are no object-independent artifacts since is smooth and the artifact curve is not defined.

Figure 2. Left: Limited-angle data (cutoff with infinite slope). Center: FBP reconstruction. Right: Reconstruction highlighting object-dependent artifact lines tangent to skull corresponding to the four circled points in the sinogram.

Figure 2 illustrates limited-angle tomography. The boundary consists of the vertical lines and . The artifact lines are exactly the lines with or that are tangent to boundaries in the object (i.e., wavefront directions are normal to the line). The four circled points on the sinogram correspond to the object-dependent artifact lines at the boundary of the skull. The corresponding lines are tangent to the skull and have angle and . One can also observe artifact lines tangent to the inside of the skull with these same angles.

4.2. Smooth boundary with finite slope

We now consider the general case in Theorem 3.6 by analyzing the artifacts for a specific set which is defined as follows. It will be cut in the middle so that the left-most cutoff boundary occurs at ; the right-most boundary is constructed as for and


for such that the two parts join differentiably at . The steepness of the curved part of the right-most boundary is governed by the constant (as seen in the two sinograms in Figure 3).

According to the condition (3.7), the curved part of is the only part that can potentially cause object-independent artifacts in , since the other parts are vertical. In Figure 3, we consider two data sets with smooth boundary, the top where the artifact curve is outside the unit disk and the bottom where it meets the object.

(A) Left: Boundary of with large slope (). Center: FBP reconstruction over the larger region to show that the artifact is outside of the region displayed in the right frame. Right: Reconstruction highlighting object-dependent artifact lines tangent to the skull corresponding to the four circled points in the sinogram.
(B) Smooth boundary with small slope ( causing a prominent curve of artifacts in . Left: Sinogram with part of boundary that could causes artifacts in the reconstruction region highlighted in magenta. The solid part of the curve indicates the artifacts that are realized in the reconstruction. The dotted curve at the right end of the sinogram indicates potential artifacts that are not realized because the corresponding part of is outside (see Theorem 3.52). Center: FBP reconstruction. Right: The magenta curve in the reconstruction is the curve of artifacts and the yellow artifact lines are object-dependent artifacts similar to those in Figure 2(A).
Figure 3. Illustration of artifacts with smooth boundary given by (4.1). The artifact curve is outside the reconstruction region in the top figure and the curve meets the object in the bottom picture.

Figure 2(A) provides a reconstruction with data set defined by in (4.1). Many artifacts in the reconstruction region are the same as in Figure 2 because the boundaries of the cutoff regions are substantially the same: the artifacts corresponding to the circles with and the lower circle with , are the same limited-angle artifacts as in Figure 2 because those parts of the boundaries are the same. However, the upper right circled point in the sinogram has so the corresponding artifact line has this larger angle, as seen in the reconstruction. The center reconstruction in Figure 2(A) shows the artifact curve, but it is far enough from that it is not visible in the reconstruction on the far right.

Figure 2(B) provides a reconstruction with data set defined by in (4.1). In this case, the object-dependent artifacts are similar to those in Figure 2(A), but the lines for defined by (4.1) are different because is different. The highlighted part of the boundary of defined by (4.1) indicates the boundary points that create the part of the artifact curve that is in the reconstruction region. The highlighted curve in the right-hand reconstruction of Figure 2(B) is this part of the curve. Note that this curve is calculated using the formula (3.5) for rather than by tracing the physical curve on the reconstruction. That they are substantially the same shows the efficacy of our theory. A simple exercise shows that, for any , the curve changes direction at .

Let be the coordinates of the circled point in the upper right of the sinogram in Figure 2(B). This circled point is on the boundary of so is tangent to the skull and an object-dependent artifact is visible along in the reconstruction. The curve ends at (as justified by Theorem 3.5 part 2) and so the curve seems to blend into this line. If were larger and the dotted part of the magenta curve on the sinogram were in , the curve would be longer.

4.3. Region-of-interest (ROI) tomography

The ROI problem, also known as interior tomography, is a classical incomplete data tomography problem in which a part of the object (the ROI) is imaged using only data over lines that meet the ROI. ROI data are used when the detector width is not large enough to contain the complete object or when researchers would like a higher resolution scan of a small part of the object. We now demonstrate how our theoretical results predict precisely which reconstruction artifacts occur for two choices of detector width.

Note that Theorem 3.5 part 3 implies that if the ROI is convex, then all singularities of in the interior of the ROI are recovered. The reason is that, in this case, all wavefront directions at all points in the interior are normal to lines in the data set. This is observed in both reconstructions in Figure 4.

(A) Sinogram and ROI reconstruction within a disk of radius centered at the origin, .
(B) Sinogram and ROI reconstruction within a disk of radius centered at the origin, .
Figure 4. ROI reconstructions with small and large ROI. In each subfigure top left is the sinogram, top right the reconstruction and bottom the same with sinogram boundary and reconstruction artifacts plotted. Note that the part of the sinogram for is shown.

The boundary of the sinogram are given by horizontal lines where in Figure 3(A) and in Figure 3(B). The points given by (3.5) are since , i.e., a circle of radius . The ring around the ROI in Figure 3(A) is exactly the curve for this ROI. There are no object-dependent streak artifacts in this picture because there are no artifacts of the object tangent to lines for .

More generally, if the ROI is smooth and strictly convex then the artifact curve traces the boundary of the ROI. The proof is an exercise using the parametrization in of tangent lines to this boundary.

Analyzing the reconstruction in Figure 3(B), one sees the top and bottom of an artifact circle which is analogous to the one in Figure 3(A). However, the artifact circle does not extend outside the object (as represented by the dotted magenta curve) because is zero near the lines corresponding to , and (this is indicated by the dotted part of in the sinogram). Therefore, the reconstruction is smooth there as explained by Theorem 3.6 part 1(b)ii.

One also sees object-dependent artifacts described by Theorem 3.6 part 2 in Figure 3(B). These occur along lines where that are tangent to the two boundaries of the skull. The angles for these lines are and (these four points are circled on the sinogram), and they are where the support of intersects the line ; these are lines at the boundary of the data set at which has wavefront set (so has wavefront set normal to the line ). The same phenomenon happens on the line with that is tangent to the inside of the skull, which is what causes the second set of four visible line artifacts indicated in yellow.

4.4. The general case

The reconstruction in Figure 5 illustrates all our cases in one. In that figure, we consider a general incomplete data set with a rectangular region cut out of the sinogram leading to all considered types of artifacts. Now, we describe the resulting artifacts. In Figure 5 the horizontal sinogram boundaries at for are displayed in solid magenta line. As in the ROI case, along these boundaries, we have and thus circular arcs of radius for the given interval for are added in the reconstruction (as indicated by solid magenta). As predicted by Theorem 3.8, each of the four corners produce a line artifact as marked by the yellow solid lines in the right-hand reconstruction, and they align tangentially with the ends of the curved artifacts.

Figure 5. Left: The sinogram for a general incomplete data problem in which the cutoff region, , has a locally smooth boundary with zero and infinite slope as well as corners. The cutout from the sinogram is at and , . Center: FBP reconstruction. Right: Reconstruction with the circular curve artifacts highlighted in magenta and object-independent “corner” streak artifacts highlighted in yellow.

The circular arc between those lines corresponds to the top and bottom parts of as the data are, locally, constrained as in ROI CT (see Section 4.3).

In Figure 5, there are other object-dependent streaks corresponding to the vertical lines in the sinogram at and at as predicted by Theorem 3.62, but they are less pronounced and more difficult to see.

4.5. Summary

We have presented reconstructions that illustrate all of types of incomplete data and each of our theorems from Section 3. All artifacts arise because of points , and they fall into two categories.

  • Streak artifacts along the line :

    • Object-dependent streaks when is smooth at and a singularity of is normal to .

    • Object-independent streaks when is nonsmooth at .

  • Artifacts on curves:

    • They are object-independent, and they are generated by the map from parts of that are smooth and of small slope.

5. Strength of added artifacts

Using the Sobolev continuity of , one can measure the strength in Sobolev scale of added artifacts in several useful cases. First, we define the Sobolev norm [45, 38]. We state it for distributions, so it applies to functions .

Definition 5.1 (Sobolev wavefront set [38]).

For , the Sobolev space is the set of all distributions with locally square integrable Fourier transform and with finite Sobolev norm:


Let have locally square integrable Fourier transform and let and . We say is in at in direction if there is a cutoff function at and an open cone containing such that the localized and microlocalized Sobolev seminorm


On the other hand, if (5.2) does not hold for any cutoff at , or conic neighborhood of , then , the Sobolev wavefront set of of order .

An exercise using the definitions shows that [11].

The Sobolev wavefront set can be defined for periodic distributions on by considering the periodic extension to as discussed for wavefront set in Section 2.2.

Note that this norm on distributions on is not the typical norm used in elementary continuity proofs for the Radon transform (see e.g., [19, equation (2.11)]), but this is the appropriate norm for the continuity theorems for general Fourier integral operators [20, Theorem 4.3.1], [8, Corollary 4.4.5].

Our next theorem gives the strength in Sobolev scale of added singularities of under certain assumptions on . It uses the relation between microlocal Sobolev strength of and of , [40, Theorem 3.1] and of and , which is given in Proposition A.7.

Theorem 5.2.

Let and let satisfy Assumption 3.1. Let and assume and is smooth normal to , i.e., .

  1. Assume is smooth and not vertical at . Let be given by (3.5) and let . Then, is in for at and . Thus, there are singularities above in the -order wavefront set of .

  2. Now, assume has a corner at . Then for each , and, except for two points on , is in for at .

  3. Now, assume . Then, for every , if and only if .

This theorem provides estimates on smoothness for more general data sets than the limited-angle case, which was thoroughly considered in [22, 31]. In contrast to part 1 of this theorem, if has a vertical tangent at , then, under the smoothness assumption on , there are no added artifacts in normal to (see Theorem 3.6 part 2). Part 1 of this theorem is a more precise version of the last assertion of Theorem 3.8. In cases 1 and 2, will cause specific singularities in specific locations on . The two more singular points in part 2 will be specified in equation (A.16). If one part of is vertical at , then there is only one more singular point.

This theorem will be proven in Section A.4 of the appendix.

6. Artifact reduction

In this section, we briefly describe a simple method to suppress the added streak artifacts described in Theorems 3.6 and 3.8. As outlined in Section 3, the application of FBP to incomplete data extends the data from to all of by padding it with zeros off of . This hard truncation can create discontinuities along and that explains the artifacts. These jumps are stronger singularities than those of for since .

One obvious way to get rid of the jump discontinuities of is to replace by a smooth function on , , that is equal to zero off of and equal one on most of and smoothly transitions to zero near . We also assume is symmetric in the sense for all .

This gives the forward operator


and the reconstruction operator


Because is a smooth function, is a standard Fourier integral operator and so is a standard pseudodifferential operator. This allows us to show that does not add artifacts.

Theorem 6.1 (Artifact Reduction Theorem).

Let and satisfy Assumption 3.1. Then


Therefore, does not add artifacts to the reconstruction.

Let , , and and assume . Then,


We provide an outline using arguments from [14]. First, is a standard pseudodifferential operator (DO) because it is the composition of , the DO , and , and because satisfies the global Bolker assumption (see (A.7)). When the top-order symbol of is nonzero, the operator is elliptic. The symbol is calculated more generally in Theorems 6.1 and 6.2 in [14] or [39]. For reference, the symbol of is

and this follows from the symbol calculation in [14, (6.3) and (A.11)] with , with weights and , , and one uses that is symmetric. ∎

Many practitioners include a smooth cutoff function in incomplete data algorithms, but others do not. Theorem 6.1 shows the advantages of including such a cutoff, and it has been suggested in several settings, including limited-angle X-ray CT [22, 12] and more general tomography problems [14, 13, 46]. More sophisticated methods are discussed in [6, 5] for the synchrotron problem that is described in Section 7.

Figure 6. Left: Smoothed sinogram. Center: Smoothed reconstruction with suppressed artifacts. Right: Reconstruction using , with sharp cutoff.

Figure 6 illustrates the effects of our smoothing algorithm. At the end of Section 7, we will illustrate the efficiency of this simple artifact reduction method on real synchrotron data.

7. Application: a synchrotron experiment

Figure 7 shows tomographic data of a chalk sample (Figure 6(A) and 6(B)) that was acquired at a synchrotron experiment [6, 5] (see [25] for related work) and a zoom of the corresponding reconstruction (Figure 6(C)). As can be clearly observed, the reconstruction includes dramatic streaks that are independent of the object. These streaks motivated the research in this article since they were not explained by the mathematical theory at that time (such as in [22, 31, 12, 13, 14]).

(A) Segmented attenuation
sinogram (processed Radon
transform data) with the
rectangular zoom region in Fig. 6(B) highlighted.
(B) Zoom of segmented sinogram.
(C) Zoomed synchrotron reconstruction given in Figure 1.
Figure 7. Left and Center: The truncated attenuation sinogram (after processing to get Radon transform data) and an enlargement of a central section of . Right: Reconstruction given in Figure 1 [5, ©IOP Publishing. Reproduced by permission of IOP Publishing. All rights reserved].
Figure 8. Data acquisition setup for the synchrotron experiment [5, ©IOP Publishing. Reproduced by permission of IOP Publishing. All rights reserved].

Taking a closer look at the attenuation sinogram, Figures 6(A)-6(B), a staircasing is revealed with vertical and horizontal boundaries. This is a result of X-rays being blocked by four metal bars that help stabilize the percolation chamber (sample holder) as the sample is subjected to high pressure during data acquisition, see Figure 8. More details are given in [5].

Because the original reconstructions of this synchrotron data used a sharp cutoff, the original reconstructions suffer from severe streak artifacts, see Figure 8(A). These artifacts are exactly described by Theorem 3.8 in that each corner of the sinogram gives rise to a line artifact in the reconstruction (cf. Figures 6(A)-6(B)). The authors of [5] then use a smooth cutoff at that essentially eliminates the streaks. The resulting reconstruction is shown in Figure 8(B).

(A) Standard FBP reconstruction
(B) FBP reconstruction with artifact reduction (cf. Theorem 6.1).
Figure 9. Reconstructions from synchrotron data without smoothing (top) and with smoothing (bottom)[5, ©IOP Publishing. Reproduced by permission of IOP Publishing. All rights reserved].

8. Discussion

We first make observations about our results for and then discuss generalizations.

8.1. Observations

Theorem 3.8 is valid as long as , but the analogous theorem for Sobolev singularities, Theorem 5.22, assumes that has a corner at . If has a weaker singularity at , then an analogous theorem would hold but one would need to factor in the Sobolev strength of the wavefront of above .

Our artifact reduction method, which is motivated by Theorem 6.1, works well for the synchrotron data as was shown in Figure 9. The article [5] provides more elaborate artifact reduction methods that are even more successful. We would also like to point out that in other incomplete data tomography problems, this simple technique might not work as efficiently as in the presented problem. In general, the artifact reduction methods need to be designed for each particular incomplete data tomography problem, but Theorem 6.1 implies that smooth cutoffs are better than abrupt cutoffs.

There are other methods to deal with incomplete data. In other approaches, authors have completed incomplete data so that the completed data is in the range of the Radon transform with full data (e.g., [26, 2, 49]). In [34] and [7], the authors develop artifact reduction methods using quantitative susceptibility mapping. For metal artifacts, there is vast literature (see, e.g., [3]) for artifact reduction methods, and we believe that those methods might also be useful for certain other incomplete data tomography problems. Authors [37, 44, 35] have effectively used microlocal analysis to understand these related problems.

Our theory is developed based on the continuous case – we view the data as functions on , not just defined at discrete points. As shown in this article, our theory predicts and explains the artifacts and visible and invisible singularities. In practice, real data are discrete, and discretization may also introduce artifacts, such as undersampling streaks. Discretization in our synchrotron experiment could be a factor in the streaks in Figure 7. Furthermore, numerical experiments have finite resolution, and this can cause (and sometimes de-emphasize) artifacts. For all these reasons, further analysis is needed to shed light on the interplay between the discrete and the continuous theory for CT reconstructions from incomplete data.

8.2. Generalizations

Our theorems were proven for , but the results hold for more general filtering operators besides . One key to our proofs is that satisfies the ellipticity condition in Remark A.5, but many other operators satisfy this condition. For example, the operator, , in Lambda CT [9] satisfies this condition, and the only difference comes in our Sobolev Continuity Theorem 5.2. Since is order two, the operator is of order and the smoothness in Sobolev scale of the reconstructions would be one degree lower than for .

Our theorems hold for fan-beam data when the source curve is smooth and convex and the object is compactly supported inside . This is true because, in this case, the fan-beam parameterization of lines is diffeomorphic to the parallel beam parametrization we use and the microlocal theorems we use are invariant under diffeomorphisms.

Theorems 3.6 and 3.8 hold verbatim for generalized Radon transforms with smooth measures on lines in because they all have the same canonical relation, given by (A.4), and the proofs would be done as for but using the basic microlocal analysis in [39].

Analogous theorems hold for other Radon transforms including the hyperplane transform and the spherical transform of photoacoustic CT. The proofs would use our arguments here plus the proofs in [14, 13]. These generalizations are the subject of ongoing work. In incomplete data problems for , the artifacts are either along curves or they are streaks along the lines corresponding to points on . However, in these higher dimensional cases, the results will be more subtle because artifacts can spread on proper subsets of the surface over which data are taken, not necessarily the entire set (see [13, Remark 4.7]).

Analogous theorems should hold for cone beam CT, but this type of CT is more subtle because the reconstruction operator itself can add artifacts, even with complete data [15, 10].

Appendix A Proofs

Here, we provide some basic microlocal analysis and then use this to prove our theorems. In this section, we adapt the standard terminology of microlocal analysis and consider wavefront sets as subsets of cotangent spaces [50]. Elementary presentations of microlocal analysis for tomography are in [23] and [18, Section 2.2]. Standard references include [11, 48].

a.1. Building blocks

Our first lemma gives some basic facts about wavefront sets.

Lemma A.1.

Let . Let and be locally integrable functions or distributions.

  1. Let be an open neighborhood of . Assume that and are equal on , then .

  2. Let be a smooth, compactly supported function, then . If is nonzero at then .

  3. if and only if there is an open neighborhood of on which is a smooth function.

The analogous statements hold for functions on .666This lemma is valid on since we consider only distributions on that are the restriction to of distributions for that are -periodic in .

These basic properties are proven using the arguments in Section 8.1 of [21], in particular, Lemma 8.1.1, Definition 8.1.2, and Proposition 8.1.3.

Our next definition will be useful to describe how wavefront set propagates under and .

Definition A.2.

If then,

If , then

The function on will be called symmetric if


If , then and