Analyzing Reconstruction Artifacts from Arbitrary Incomplete X-ray CT Data

This article provides a mathematical analysis of singular (nonsmooth) artifacts added to reconstructions by filtered backprojection (FBP) type algorithms for X-ray CT with arbitrary incomplete data. We prove that these singular 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. Object-dependent artifacts are generated by singularities of the object being scanned and these artifacts can extend along lines. They generalize the streak artifacts observed in limited-angle tomography. Object-independent artifacts, on the other hand, are essentially independent of the object and take one of two forms: streaks on lines if the boundary of the data set is not smooth at a point and curved artifacts if the boundary is smooth locally. We prove that these streak and curve artifacts are the only singular artifacts that can occur for FBP in the continuous case. 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 that arise in new practical applications. 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 on lines that have no relation to the object.


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 [22,27]; 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 [5,6] 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, artifacts 1 , such as the limited-angle streaks in Figure 2 can occur. Therefore, reconstruction quality suffers considerably, and this complicates the proper interpretation of images.
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 [24,40]. Subsequently, artifacts were extensively studied in the context of limited-angle tomography, e.g., [12,22]. 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 [1,13,14,32,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 [35,37,44], 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,10,15,43,47].
To the best of our knowledge, the mathematical literature up until now, e.g., [12,13,22,31,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 f is the density of the object to be reconstructed, then each CT measurement is modeled by a line integral of f over a line in the data set. As we will fully describe in Section 2.1, we parametrize lines by (ϕ, p) ∈ [0, 2π] × R, and the CT measurement of f over the line L(ϕ, p) is denoted Rf (ϕ, p). With complete data, where Rf (ϕ, p) is given over all (ϕ, p) ∈ [0, 2π] × R, accurate reconstructions can be produced by the FBP algorithm. In incomplete data CT problems, the data are taken over lines L(ϕ, p) for (ϕ, p) in a proper subset, A, of [0, 2π] × R 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, bd(A), and we show that there are two types of artifacts: object-dependent and objectindependent artifacts. The object-dependent artifacts are caused by singularities of the object being scanned. In this case, artifacts can spread all along a line L(ϕ 0 , p 0 ) (i.e., a streak) if (ϕ 0 , p 0 ) ∈ bd(A) 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 A is smooth near a point (ϕ 0 , p 0 ) ∈ bd(A), then we prove that artifacts can appear in the reconstruction along curves generated by bd(A) near (ϕ 0 , p 0 ), and they occur whether the object being scanned has singularities or not. We also prove if bd(A) is not smooth (see Definition 3.2) at a point (ϕ 0 , p 0 ) then, independent of the object, an artifact line can be generated all along L(ϕ 0 , p 0 ).
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.

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 L 2 (D), the set of absolutely square integrable functions supported in the open unit disk D = x ∈ R 2 : x < 1 . 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.
θ(ϕ) = (cos(ϕ), sin(ϕ)) and θ ⊥ (ϕ) = (− sin(ϕ), cos(ϕ)), so θ(ϕ) is the unit vector in the direction of ϕ and θ ⊥ (ϕ) is the unit vector π/2 radians counterclockwise from θ(ϕ). The line perpendicular to θ(ϕ) and containing pθ(ϕ) is denoted This line is parameterized by t → pθ(ϕ) + tθ ⊥ (ϕ). Note that L(ϕ, p) = L(ϕ + π, −p). We define the X-ray transform or Radon line transform of f ∈ L 2 (D) to be For functions g on [0, 2π] × R, 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 (x, ξdx) with the associated ordered pair of vectors (x, ξ). 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 x 0 ∈ R 2 and ξ 0 ∈ R 2 \ 0. A cutoff function at x 0 is any C ∞ -function of compact support that is nonzero at x 0 . Let f be a locally integrable function. We say f is smooth at x 0 in direction ξ 0 if there is a cutoff function ψ at x 0 and an open cone V containing ξ 0 such that the Fourier transform F(ψf )(ξ) is rapidly decaying at infinity. 2 If f is not smooth at x 0 in direction ξ 0 , then we say f has a singularity at x 0 in direction ξ 0 . The wavefront set of f is the set WF(f ) of all such (x 0 , ξ 0 ).
For example, if A is a subset of the plane with a smooth boundary, and f is equal to 1 on A and 0 off of A, then WF(f ) is the set of all points (x, ξ) where x is in the boundary of A and ξ is normal to the boundary at x.
We parametrize lines on [0, 2π] × R, rather than S 1 × R, because, in practice, the data space consists of points (ϕ, p) ∈ [0, 2π] × R and the sinogram is represented as a picture in the (ϕ, p)plane, not on a cylinder. We will consider only distributions on [0, 2π] × R that are restrictions of distributions for (ϕ, p) ∈ R 2 that are 2π-periodic in ϕ. The Radon transform and our other operators on [0, 2π] × R are all 2π-periodic in ϕ, and so this is no real restriction. For this reason, we can use the same definition of wavefront set as for R 2 , and we just apply it to the 2π-periodic extensions.
Our next definition allows us to express efficiently the correspondence between singularities of the object and of its Radon transform.
and the set of singularities of f normal to L(ϕ, p) is Let f be a locally integrable function on R 2 . We say f is smooth normal to the line It is important to understand each set introduced in Definition 2.2: N * (L(ϕ, p)) \ 0 is the set of all (x, ξ) such that x ∈ L(ϕ, p) and the vector ξ is normal to L(ϕ, p) at x. Therefore, WF L(ϕ,p) (f ) is the set of wavefront directions of f above points x ∈ L(ϕ, p) that are normal to this line.
If g is a locally integrable function on [0, 2π]×R, then WF (ϕ,p) (g) is the set of wavefront directions with base point (ϕ, p). The sets defined in Definition 2.2 have an important relationship that we will exploit starting in the next section.

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 [0, 2π] × R by setting it to zero off of the set A (cutoff region) over which data are taken. Therefore, the incomplete CT data can be modeled as where 1 A is the characteristic function of A. 3 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 R * is defined by (2.4). Our next assumption collects the conditions we will impose on the cutoff region A. There, we will use the notation int(A), bd(A), and ext(A) to denote the interior of A, the boundary of A, and the exterior of A, respectively. 3 The characteristic function of a set A is the function that is equal to one on the set and zero outside of A. 5 Assumption 3.1. Let A be a proper subset of [0, 2π] × R (i.e., A [0, 2π] × R) with a nontrivial interior and assume A is symmetric in the following sense: In addition, assume that A is the smallest closed set containing int(A), i.e. A = cl(int(A)).
We now justify using this assumption. Since A is proper, data over A are incomplete. Being symmetric means that the part of A in [0, π] × R determines the data set. We assume that A is the smallest closed set containing int(A) in order to exclude degenerate cases, such as when A includes an isolated curve.
Our next definition gives us the language to describe the geometry of bd(A).
-If this tangent line is vertical (i.e., of the form ϕ = ϕ 0 ), then we say the boundary has infinite slope at (ϕ 0 , p 0 ). -If this tangent line is not vertical, then bd(A) is defined near (ϕ 0 , p 0 ) by a smooth function p = p(ϕ). In this case, the slope of the boundary at (ϕ 0 , p 0 ) will be the slope of this tangent line, which is given by p (ϕ 0 ). • We say that bd(A) is not smooth at (ϕ 0 , p 0 ) if it is not a smooth curve at (ϕ 0 , p 0 ).
-We say that bd(A) has a corner at (ϕ 0 , p 0 ) if the curve bd(A) is continuous at (ϕ 0 , p 0 ), is smooth at all other points in some neighborhood of (ϕ 0 , p 0 ), and has one-sided slopes at (ϕ 0 , p 0 ) but those slopes are different.

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 L(ϕ, p) for (ϕ, p) ∈ bd(A). This material is either known (e.g., [40]) or it follows from [40].

Definition 3.3 (Artifacts and visible singularities)
. Every singularity of L A f that is not a singularity of f is called an artifact (i.e., any singularity in WF(L A f ) \ WF(f )). A streak artifact is a line of artifacts. Every singularity of L A f that is a singularity of f is called a visible singularity (from data on A) (i.e., any singularity in WF(L A f ) ∩ WF(f )). Other singularities of f are called invisible (from data on A). 4

Theorem 3.4 (Localization of artifacts to lines from bd(A)). Let
This theorem follows from the microlocal analysis in [40,42]. More generally, every singularity of L A f must come from a singularity of 1 A Rf , and if the singularity does not come from Rf , it must come from a singularity of 1 A , and these singularities are on bd(A).
Our next theorem gives an analysis of singularities in L A f corresponding to lines not in bd(A). It shows that visible singularities are along lines L(ϕ, p) when (ϕ, p) ∈ int(A) and invisible singularities are on lines for (ϕ, p) ∈ ext(A). 5 4 Invisible singularities are smoothed by LA and reconstruction of those singularities is in general is exponentially ill-posed in Sobolev scale. 5 Visible singularities of f can appear on lines L(ϕ0, p0) for (ϕ0, p0) ∈ bd(A), but these can cause artifacts all along the line L(ϕ0, p0) 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 f ∈ L 2 (D) and let
In this case, L A f is smooth normal to L(ϕ, p). That is, all singularities of f normal to L(ϕ, p) are invisible from data on A.
(iii) Now, let x ∈ D and assume that all lines through x are parameterized by points in int(A) (i.e., ∀ϕ ∈ [0, 2π], (ϕ, x · θ(ϕ)) ∈ int(A)). Then all singularities of f at x are visible in L A f : Therefore, artifacts occur only on lines in bd(A).
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.

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 bd(A) and, in some cases, singularities of the object f . Our next two theorems show that artifacts arise in two ways: • Object-independent artifacts: those caused essentially by the geometry of bd(A).
• Object-dependent artifacts: those caused by singularities of the object f that are normal to L(ϕ, p) for some (ϕ, p) ∈ bd(A). As a first generalization of the limited-angle characterizations in [12,22] we consider the case where bd(A) 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. I. Assume bd(A) has finite slope at (ϕ 0 , p 0 ) and let I be a neighborhood of ϕ 0 such that bd(A) is given by a smooth curve p = p(ϕ) near (ϕ 0 , p 0 ). Let A. An object-independent artifact can occur along the curve I ϕ → x b (ϕ).
B. If f has no singularity normal to L(ϕ 0 , p 0 ) then Rf is zero in a neighborhood of (ϕ 0 , p 0 ), then the x b (ϕ) artifact curve will not be visible near x b (ϕ 0 ). II. In all cases with smooth boundary at (ϕ 0 , p 0 ), if f has a singularity normal to L(ϕ 0 , p 0 ), then an object-dependent streak artifact can occur along L(ϕ 0 , p 0 ).
A. If f is smooth conormal to L(ϕ 0 , p 0 ) and bd(A) is vertical at (ϕ 0 , p 0 ), then L A f is smooth normal to L(ϕ 0 , p 0 ). B. If bd(A) is not vertical then L A f is smooth normal to L(ϕ 0 , p 0 ) except possibly at x b (ϕ 0 ). Theorem 3.6 part II. 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 objectdependent streak artifacts might not occur, even if there is a singularity of f normal to this line. This is why we state that artifacts can occur in parts I.A. and II. 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 bd(A) is smooth at with finite slope at (ϕ 0 , p 0 ). Let I be a neighborhood of ϕ 0 and p : I → R describes bd(A) near (ϕ 0 , p 0 ). Note that If the slope of bd(A) at (ϕ 0 , p 0 ) is small enough, i.e., holds, then the artifact curve ϕ → x b (ϕ) will be inside the unit disk, D, near (ϕ 0 , p 0 ). If not, this curve will be outside, D, at least at ϕ 0 . If bd(A) is smooth and vertical at (ϕ 0 , p 0 ) (infinite slope), then there will be no object-independent artifact on the line L(ϕ 0 , p 0 ) because x b (ϕ 0 ) is a "point at infinity" since "p (ϕ 0 )" is infinite. This is proven using the arguments in the proof of Theorem 3.6 B.
When bd(A) is not smooth at (ϕ 0 , p 0 ), the next theorem shows there can be a streak artifact in L A f along L(ϕ 0 , p 0 ) that is independent of the object f .
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 [12,22,31], but we are unaware of a general reference to Theorem 3.6, part II. 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 x b curve artifact in Theorem 3.6 part I. or to the corner artifacts from Theorem 3.8.

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 ϕ ∈ [0, π] on the horizontal axis and p ∈ [− √ 2, √ 2] on the vertical. Except for the center picture in Figure 3(A), the reconstruction is displayed on [−1, 1] 2 . 4.1. Limited-angle tomography. First, we analyze limited-angle tomography, a classical problem in which part II. of Theorem 3.6 applies. In this case bd(A) consists of vertical lines given by ϕ = ϕ 0 (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 L(ϕ 0 , p 0 ) in the data set with (ϕ 0 , p 0 ) ∈ bd(A) and which is normal to a singularity of f , then a streak artifact can be generated along L(ϕ 0 , p 0 ) in the reconstruction L A f . Therefore, Theorem 3.6 generalizes the results of [12,22] 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 bd(A) is smooth and the artifact curve ϕ → x b (ϕ) is not defined.  The boundary consists of the vertical lines ϕ = 4π/9 and ϕ = 5π/9. The artifact lines are exactly the lines with ϕ = 4π/9 or 5π/9 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 ϕ = 4π/9 and ϕ = 5π/9. 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 A which is defined as follows. It will be cut in the middle so that the left-most cutoff boundary occurs at ϕ = a := 4 9 π; the right-most boundary is constructed as ϕ = b := 5 9 π for p ≤ 0 and (4.1) for p > 0 such that the two parts join differentiably at (ϕ, p) = (0, 0). The steepness of the curved part of the right-most boundary is governed by the constant c (as seen in the two sinograms in Figure 3). According to the condition (3.7), the curved part of bd(A) is the only part that can potentially cause object-independent artifacts in D, since the other parts are vertical. In Figure 3, we consider two data sets A with smooth boundary, the top where the x b artifact curve is outside the unit disk and the bottom where it meets the object. Figure 3(A) provides a reconstruction with data set defined by c = 1.3 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 ϕ = 4π/9 and the lower circle with ϕ = 5π/9, 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 ϕ > 5π/9 so the corresponding artifact line has this larger angle, as seen in the reconstruction. The center reconstruction in Figure 3(A) shows the x b artifact curve, but it is far enough from D that it is not visible in the reconstruction on the far right. Figure 3(B) provides a reconstruction with data set defined by c = 0.65 in (4.1). In this case, the object-dependent artifacts are similar to those in Figure 3(A), but the lines for (ϕ, p) defined by (4.1) are different because bd(A) is different. The highlighted part of the boundary of A defined by (4.1) indicates the boundary points that create the part of the x b (ϕ) artifact curve that is in the reconstruction region. The highlighted curve in the right-hand reconstruction of Figure 3(B) is this part of the x b (ϕ) curve. Note that this curve is calculated using the formula (3.5) for x b (ϕ) 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 c > 0, the x b curve changes direction at x b (1/2 + 5π/9).
Let (ϕ 0 , p 0 ) be the coordinates of the circled point in the upper right of the sinogram in Figure  3(B). This circled point is on the boundary of supp(Rf ) so L(ϕ 0 , p 0 ) is tangent to the skull and an object-dependent artifact is visible along L(ϕ 0 , p 0 ) in the reconstruction. The x b curve ends at x b (ϕ 0 ) (as justified by Theorem 3.5 part (ii)) and so the x b curve seems to blend into this line. If supp(f ) were larger and the dotted part of the magenta curve on the sinogram were in supp(Rf ), the x b 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 (iii) implies that if the ROI is convex, then all singularities of f 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.
The boundary of the sinogram are given by horizontal lines p = ±p 0 where p 0 = 0.4 in Figure 4(A) and p 0 = 0.8 in Figure 4(B). The points given by (3.5) are x b (ϕ) = p 0 θ(ϕ) + 0θ ⊥ (ϕ) since p = 0, i.e., a circle of radius p 0 . The ring around the ROI in Figure 4(A) is exactly the curve ϕ → x b (ϕ) 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 (ϕ, p) ∈ bd(A).
More generally, if the ROI is smooth and strictly convex then the artifact curve x b traces the boundary of the ROI. The proof is an exercise using the parametrization in (ϕ, p) of tangent lines to this boundary.
Analyzing the reconstruction in Figure 4(B), one sees the top and bottom of an x b artifact circle which is analogous to the one in Figure 4(A). However, the artifact circle does not extend outside the object (as represented by the dotted magenta curve) because Rf is zero near the lines corresponding to p = ±0.8, ϕ < 0.7121 and ϕ > 2.429 (this is indicated by the dotted part of bd(A) in the sinogram). Therefore, the reconstruction is smooth there as explained by Theorem 3.6 part (ii).
One also sees object-dependent artifacts described by Theorem 3.6 part II. in Figure 4(B). These occur along lines L(ϕ 0 , p 0 ) where p 0 = ±0.8 that are tangent to the two boundaries of the skull. The angles for these lines are ϕ 0 = 0.7121 and ϕ 0 = 2.429 (these four points (ϕ 0 , p 0 ) are circled on the sinogram), and they are where the support of Rf intersects the line p = ±0.8; these are lines at the boundary of the data set at which Rf has wavefront set (so f has wavefront set normal to the line L(ϕ, p)). The same phenomenon happens on the line with p = ±0.8 that is tangent to the inside of the skull, which is what causes the second set of four visible line artifacts indicated in yellow. 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 p = p 0 = ±0.35 for φ ∈ 7 18 π, 11 18 π are displayed in solid magenta line. As in the ROI case, along these boundaries, we have p = 0 and thus circular arcs of radius p 0 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.

The general case. The reconstruction in
The circular arc between those lines corresponds to the top and bottom parts of bd(A) 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 ϕ = 7π 18 and at ϕ = 11π 18 as predicted by Theorem 3.6II., 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 (ϕ 0 , p 0 ) ∈ bd(A), and they fall into two categories.
• Streak artifacts along the line L(ϕ 0 , p 0 ): -Object-dependent streaks when bd(A) is smooth at (ϕ 0 , p 0 ) and a singularity of f is normal to L(ϕ 0 , p 0 ). -Object-independent streaks when bd(A) is nonsmooth at (ϕ 0 , p 0 ). • Artifacts on curves: -They are object-independent, and they are generated by the map ϕ → x b (ϕ) from parts of bd(A) that are smooth and of small slope.

Strength of added artifacts
Using the Sobolev continuity of Rf , one can measure the strength in Sobolev scale of added artifacts in several useful cases. First, we define the Sobolev norm [38,45]. We state it for distributions, so it applies to functions f ∈ L 2 (D).
Definition 5.1 (Sobolev wavefront set [38]). For s ∈ R, the Sobolev space H s (R n ) is the set of all distributions with locally square integrable Fourier transform and with finite Sobolev norm: Let f have locally square integrable Fourier transform and let x 0 ∈ R n and ξ 0 ∈ R n \ 0. We say f is in H s at x 0 in direction ξ 0 if there is a cutoff function ψ at x 0 and an open cone V containing ξ 0 such that the localized and microlocalized Sobolev seminorm On the other hand, if (5.2) does not hold for any cutoff at x 0 , ψ or conic neighborhood V of ξ 0 , then (x 0 , ξ 0 ) ∈ WF s (f ), the Sobolev wavefront set of f of order s.
An exercise using the definitions shows that WF(f ) = ∪ s∈R WF s (f ) [11]. The Sobolev wavefront set can be defined for periodic distributions on [0, 2π] × R by considering the periodic extension to (ϕ, p) ∈ [0, 2π] × R as discussed for C ∞ wavefront set in Section 2.2.
Note that this norm on distributions on [0, 2π]×R is not the typical H 0,s 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,  (i) Assume bd(A) is smooth and not vertical at (ϕ 0 , p 0 ). Let x b = x b (ϕ 0 ) be given by (3.5) and let ω = 0. Then, L A f is in H s for s < 0 at ξ 0 = (x b , ωθ(ϕ 0 )) and ξ 0 ∈ WF 0 (L A f ). Thus, there are singularities above x b in the 0-order wavefront set of L A f . (ii) Now, assume bd(A) has a corner at (ϕ 0 , p 0 ). Then for each (x, ξ) ∈ N * (L(ϕ 0 , p 0 )) \ 0, (x, ξ) ∈ WF 1 (L A f ) and, except for two points on L(ϕ 0 , p 0 ), L A f is in H s for s < 1 at (x, ξ). (iii) Now, assume (ϕ 0 , p 0 ) ∈ int(A). Then, for every (x, ξ) ∈ N * (L(ϕ 0 , p 0 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 (i) of this theorem, if bd(A) has a vertical tangent at (ϕ 0 , p 0 ), then, under the smoothness assumption on f , there are no added artifacts in L A f normal to L(ϕ 0 , p 0 ) (see Theorem 3.6 part II.). Part (i) of this theorem is a more precise version of the last assertion of Theorem 3.8. In cases (i) and (ii), bd(A) will cause specific singularities in specific locations on L(ϕ 0 , p 0 ). The two more singular points in part (ii) will be specified in equation (A.16). If one part of bd(A) is vertical at (ϕ 0 , p 0 ), then there is only one more singular point.
This theorem will be proven in Section A.4 of the appendix.

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 A ⊂ [0, 2π]×R to all of [0, 2π]×R by padding it with zeros off of A. This hard truncation can create discontinuities along bd(A) and that explains the artifacts. These jumps are stronger singularities than those of Rf One obvious way to get rid of the jump discontinuities of 1 A is to replace 1 A by a smooth function on [0, 2π] × R, ψ, that is equal to zero off of A and equal one on most of int(A) and smoothly transitions to zero near bd(A). We also assume ψ is symmetric in the sense ψ(ϕ, p) = ψ(ϕ + π, −p) for all (ϕ, p).
This gives the forward operator and the reconstruction operator Because ψ is a smooth function, R ψ is a standard Fourier integral operator and so L ψ is a standard pseudodifferential operator. This allows us to show that L ψ does not add artifacts. Therefore, L ψ does not add artifacts to the reconstruction. Let x ∈ D, ϕ ∈ [0, 2π], and ω = 0 and assume ψ(ϕ, x · θ(ϕ)) = 0. Then, Proof. We provide an outline using arguments from [14]. First, L ψ is a standard pseudodifferential operator (ΨDO) because it is the composition of R * , the ΨDO Λψ, and R, and because R satisfies the global Bolker assumption (see (A.7)). When the top-order symbol of L ψ is nonzero, the operator L ψ is elliptic. The symbol is calculated more generally in Theorems 6.1 and 6.2 in [14] or [39]. For reference, the symbol of L ψ is σ (L ψ ) (x, ωθ(ϕ)) = 4πψ ϕ, x · θ(ϕ) , and this follows from the symbol calculation in [14, (6.3) and (A.11)] with ξ = ωθ(ϕ), with weights µ = 1 and ν = 1, P = ξ , 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 [12,22] and more general tomography problems [13,14,46]. More sophisticated methods are discussed in [5,6] for the synchrotron problem that is described in Section 7.  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. Figure 7 shows tomographic data of a chalk sample (Figure 7(A) and 7(B)) that was acquired at a synchrotron experiment [5,6] (see [25] for related work) and a zoom of the corresponding reconstruction (Figure 7(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 [12-14, 22, 31]).

Application: a synchrotron experiment
Taking a closer look at the attenuation sinogram, Figures 7(A)-7(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 9(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 7(A)-7(B)). The authors of [5] then use a smooth cutoff at bd(A) that essentially eliminates the streaks. The resulting reconstruction is shown in Figure 9(B).

Discussion
We first make observations about our results for L A and then discuss generalizations.   For the detector position along the SW-NE diagonal, the beam is occluded completely by the metal bars and no signal is measured at the detector. In between one position is shown at which the outermost parts of the beam are occluded by two of the metal bars, leading to a projection with truncation from both sides.

Acquisition set-up
The motivating case for the present study is in situ X-ray micro-tomography imaging of fluid flow through porous chalk in which the goal is to recover oil from the North Sea underground. In situ X-ray tomography data was obtained for a cylindrical porous chalk sample of diameter 0.6 mm using beamline BL20XU of the SPring-8 Synchrotron Radiation Facility, Japan using a monochromatic (28 keV) parallel-beam scan configuration. Fluid is forced through the sample by a percolation cell, seen in Figure   8.1. Observations. Theorem 3.8 is valid as long as WF (ϕ 0 ,p 0 ) (1 A ) = {(ϕ 0 , p 0 )} × R 2 \ 0 , but the analogous theorem for Sobolev singularities, Theorem 5.2(ii), assumes that A has a corner at (ϕ 0 , p 0 ). If A has a weaker singularity at (ϕ 0 , p 0 ), then an analogous theorem would hold but one would need to factor in the Sobolev strength of the wavefront of 1 A above (ϕ 0 , p 0 ). 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., [2,26,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 [35,37,44] 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 [0, 2π]× R, 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 L A = R * (Λ (1 A R)), 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, L = − ∂ 2 ∂p 2 , in Lambda CT [9] satisfies this condition, and the only difference comes in our Sobolev Continuity Theorem 5.2. Since L is order two, the operator R * LR is of order 1 and the smoothness in Sobolev scale of the reconstructions would be one degree lower than for L A .
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 R 2 because they all have the same canonical relation, given by (A.4), and the proofs would be done as for L A 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 [13,14]. These generalizations are the subject of ongoing work. In incomplete data problems for R, the artifacts are either along curves ϕ → x b (ϕ) or they are streaks along the lines corresponding to points on bd(A). 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 [10,15].

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 x 0 ∈ R 2 . Let u and v be locally integrable functions or distributions.
(i) Let U be an open neighborhood of x 0 . Assume that u and v are equal on U , then WF x 0 (u) = WF x 0 (v). (ii) Let ψ be a smooth, compactly supported function, then WF x 0 (ψu) ⊂ WF x 0 (u). If ψ is nonzero at x 0 then WF x 0 (u) = WF x 0 (ψu). Our next definition will be useful to describe how wavefront set propagates under R and R * .
Our next proposition is the main technical theorem of the article. It provides the wavefront correspondences for R and R * which we will use in our proofs. Proposition A.3 (Microlocal correspondence of singularities [42]). The X-ray transform, R, is an elliptic Fourier integral operator with canonical relation Let f ∈ L 2 (D) and let g be a locally integrable function on [0, 2π] × R that is symmetric by (A.1). Let x 0 ∈ R 2 , ϕ 0 ∈ [0, 2π], and let p, α, and ω be real numbers with ω = 0.
Proof. The facts about R are directly from [40,Theorem 3.1] or [42,Theorem A.2], and they use the calculus of the FIO R [16,17] (see also [39]). Note that the crucial point is that R is an elliptic Fourier integral operator that satisfies the global Bolker Assumption: the natural projection (A.7) Π L : C → T * (Y ) is an injective immersion, so (A.5) holds for R. Note that we are using the identification (A.3) in asserting that (A.5) is an equality. The proofs for R * are parallel to those for R except they involve the canonical relation for R * , C t , rather than C. Our next remark will be used in ellipticity proofs that follow.
Remark A.5. The operator Λ is elliptic in all cotangent directions except dϕ because the symbol of Λ is |τ | where τ is the Fourier variable dual to p. However, the dϕ direction will not affect our proofs. This is true because, for any function f ∈ L 2 (D), the covector (ϕ, p, ωdϕ) is not in WF(Rf ) because WF(Rf ) = C • WF(f ) (use the definition of composition and (A.4)). So, for each f ∈ L 2 (D), WF(ΛRf ) = WF(Rf ). Because C t • {(ϕ, p, αdϕ)} = ∅ by (A.4), even if (ϕ, ωdϕ) ∈ WF(1 A Rf ), that covector will not affect the calculation of C t • WF(Λ1 A Rf ). Therefore, Λ is elliptic on all cotangent directions that are preserved when composed with C t , and these are all the directions we need in our proofs. Our theorems will be valid for any 2π-periodic ΨDO on [0, 2π] × R that is invariant under the symmetry condition (A.1) and satisfies this ellipticity condition (although the Sobolev results will depend on the order of the operator).
If (ϕ 0 , p 0 ) ∈ int(A) then there is a neighborhood U of (ϕ 0 , p 0 ) that is contained in A. By Lemma A.1 part (i), since Rf and 1 A Rf agree in U , WF (ϕ 0 ,p 0 ) (Rf ) = WF (ϕ 0 ,p 0 ) (1 A Rf ). The proof is finished using (A.6) and (A.9). The proof of part (ii) follows the same argument except that 1 A Rf is zero (so smooth) in a neighborhood of (ϕ 0 , p 0 ).
The last part of the theorem is proven using part (i) and the fact that every line through x is parameterized by points in int(A).
A.3. Proof of Theorems 3.6 and 3.8. In the proofs of Theorems 3.6 and 3.8, we use Proposition A.3 to analyze how multiplication by 1 A adds singularities to the data Rf and then to the reconstruction, L A f . We first make observations that will be useful in the proofs of both theorems.