FULL CHARACTERIZATION OF RECONSTRUCTION ARTIFACTS FROM ARBITRARY INCOMPLETE X-RAY CT DATA
arXiv:1707.03055v3 [math.FA] 24 Jan 2018
¨ LEISE BORG1 , JURGEN FRIKEL2 , JAKOB SAUER JØRGENSEN3 , AND ERIC TODD QUINTO4 Abstract. This article provides a mathematical classification of artifacts from arbitrary incomplete X-ray tomography data when using the classical filtered backprojection algorithm. Using microlocal analysis, we prove that all artifacts arise from points at the boundary of the data set. Our results show that, depending on the geometry of the data set 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. This is a generalization of the streak artifacts observed in limited angle CT. The article also characterizes two new phenomena: the object-independent artifacts are caused only by the geometry of the data set boundary; they occur along lines if the boundary of the data set is not smooth and along curves if the boundary of the data set is smooth. In addition to the geometric description of artifacts, the article also provides characterizations of their strength in Sobolev scale in certain cases. Moreover, numerical reconstructions from simulated and real data are presented illustrating our theorems. This work is motivated by a reconstruction we present from a synchrotron data set in which artifacts along lines appeared that were independent of the object. The results of this article apply to a wide range of well-known incomplete data problems, including limited angle CT and region of interest tomography, as well as to unconventional x-ray CT imaging setups. Some of those problems are explicitly addressed in this article, theoretically and numerically.
1. Introduction Over the past decades Computed Tomography (CT) has established itself as a standard imaging technique in many areas such as materials sciences, medical imaging, and several other fields. Its principle is based on collecting numerous x-ray measurements of the object along all possible directions (lines) and then reconstruct 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 [24, 30]. However, in many recent applications of CT, the data cannot be collected all around the object leading to incomplete (or limited) data. The reasons for data incompleteness might be technical; for example, in the in situ experiment discussed in Section 7, the scanning equipment obscures some x-ray lines. Alternatively, incomplete data might be taken for health-related issues, such as for dose reduction, or for other limitations. Concrete examples of such imaging situations include limited angle tomography, where the data can be collected only from certain view-angles [22]; interior or 1
Department of Computer Science, University of Copenhagen (
[email protected]) Department of Computer Science and Mathematics, OTH Regensburg, (
[email protected]) 3 School of Mathematics, University of Manchester, Manchester, M13 9PL, United Kingdom (
[email protected]). Part of the work was carried out while the author was with the Department of Applied Mathematics and Computer Science, Technical University of Denmark. 4 Department of Mathematics, Tufts University, Medford, MA 02155 USA, (
[email protected]) 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 2
1
region of interest (ROI) tomography, where the x-ray measurements are available only over lines intersecting a certain subregion of the object [8, 19]; or exterior tomography, where only measurements over all lines outside a certain subregion are available [23, 35]. Another important example of incomplete data tomography is given by any tomographic imaging situation where the imaged object has metal inclusions that destroy parts of the data and lead to so-called metal artifacts. In all of above mentioned imaging situations, a (significant) portion of the data is missing or destroyed leading to instabilities during the reconstruction (see, for example, [3]). In particular, it is wellknown that certain image features cannot be reconstructed reliably [34] and that artifacts might be generated [10, 19]. As a consequence, the reconstruction quality suffers, and this complicates the proper interpretation of images. In this article we use microlocal analysis to understand artifacts and visible singularities in the general incomplete data problem for X-ray tomography in two dimensions. Throughout the article, by artifacts we mean features (streaks or points) that are added by the algorithm and that are not part of the original object. Visible singularities are singularities of the object that are visible in the reconstruction and invisible singularities are singularities of the object that are smoothed in the reconstruction and cannot be stably imaged. We study the generation of artifacts in FBP type reconstruction methods with arbitrary incomplete data, and now, we will describe our results in general. Starting in Section 2.1, we will make these ideas precise by associating each line L (representing an x-ray) to a point (ϕ0 , p0 ) ∈ [0, 2π] × R, such that L = L(ϕ0 , p0 ) and consider the (incomplete) data set, A, as a subset of [0, 2π] × R. We will define singularities rigorously and prove our theorems using microlocal analysis – the analysis of singularities of functions and what operators do to them. We prove that incomplete data artifacts arise from points at the boundary or “edge” of the data set, bd(A). Depending on the geometry of bd(A) we show that there are two types of artifacts: object-dependent and object-independent artifacts. The object-dependent artifacts are caused by singularities (e.g., density jumps) of the object being scanned. They can only be generated if the boundary curve is smooth near a point (ϕ0 , p0 ) ∈ bd(A) and if the slope of this curve at the point (ϕ0 , p0 ) fulfills certain restrictions. In this case, artifacts can be generated and spread all along a line L(ϕ0 , p0 ) (streak) that is defined by the boundary point (ϕ0 , p0 ) if this line is tangent some singularity of the scanned object (we say that this singularity generates the artifact). This is a generalization of the artifacts observed in limited angle CT. Furthermore, we characterize two new phenomena that are independent of the object being scanned. If the boundary of A is smooth near a point (ϕ0 , p0 ) ∈ bd(A), then we prove that artifacts can appear in the reconstruction along curves (related to the slope of the curve bd(A)), and they are (essentially) independent of the object being scanned. We also prove if bd(A) is not smooth (see Definition 3.2) at a point (ϕ0 , p0 ) then, completely independently of the object, an artifact line can be created all along L(ϕ0 , p0 ). We illustrate our results using classical and well studied problems such as limited angle tomography and region of interest tomography, as well as problems with novel data sets obtained by micro-CT. In addition to the geometric characterization of artifacts, we also provide a characterization of their strength in Sobolev scale. Characterization of artifacts and visible singularities has been extensively studied in the in the context of limited angle tomography and Lambda tomography, see [10, 19]. The strength of added artifacts in limited angle CT was analyzed in [26]. Similar characterizations 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, 11, 12, 27, 28]. Moreover, metal artifacts have been mathematically modeled and evaluated using microlocal analysis in [31, 37], both of which are in the spirit of this article. A related problem is studied in [29], where authors develop a streak reduction method for quantitative susceptibility mapping (see also [6]). 2
This research project was motivated by reconstructions from [4, 5] of chalk samples in a synchrotron experiment described in Section 7. The reconstructions include dramatic streaks that are independent of the object, as can be seen in Figure 1. These streaks could not be explained by the microlocal theory at that time (e.g., [10–12, 19, 26]) since they are independent of the object, and that theory explains only streaks that were generated by singularities of the object. A thorough practical investigation of this particular problem was recently presented in [4]. In this article, we present a unified approach to characterize artifacts, visible and invisible singularities for arbitrary incomplete data, including the synchrotron experiment (see Figure 1 and Section 7) and the classical limited data CT problems such as limited angle, region of interest and exterior tomography.
Figure 1. Left: A small part of the sinogram of a chalk sample. Notice the rough boundary. Right: Small central section of a reconstruction of the chalk. Monochromatic parallel beam data were taken over 1800 views covering 180 degrees, and there were 2048 × 2048 detector elements with a 0.5 mm field of view, providing micrometer resolution of the sample. Data [43] obtained, with thanks from the Japan Synchrotron Radiation Research Institute from beam time on beamline BL20XU c of SPring-8 (Proposal 2015A1147). For more details, see Section 7 and [4, IOP Publishing. Reproduced by permission of IOP Publishing. All rights reserved]. In section 2, we provide notation and some of the basic ideas about distributions and wavefront sets. In section 3 we give our main 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. 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 [38] for an overview of distributions), but we will consider only absolutely square integrable functions supported in the open unit disk D = x ∈ R2 : kxk < 1 because this is more realistic in practice and because our theorems are simpler in this case than for general distributions. Remark A.5 provides some more perspective on this point. We denote this set of functions as L2 (D). 3
A function f is locally integrable if it isR Lebesgue measurable and its absolute value, |f | is integrable on all compact sets, ∀K ⋐ R2 , K |f | < ∞, and f is locally square integrable if it is measurable and |f |2 is integrable on all compact sets. We use the analogous definitions on [0, 2π]×R. We use the standard convention that two locally integrable functions are equal if they are equal almost everywhere (a.e.). This means that they are equal as distributions. We say a function is smooth on an open set U if it is equal to a smooth function a.e. on U . 2.1. Notation. For (ϕ, p) ∈ [0, 2π] × R, we define θ(ϕ) = (cos(ϕ), sin(ϕ)) and
(2.1)
θ ⊥ (ϕ) = (− sin(ϕ), cos(ϕ)),
so θ(ϕ) is the unit vector in the direction of ϕ and θ ⊥ (ϕ) is the unit vector π/2 radians counterclockwise from θ(ϕ). We denote the line perpendicular to θ(ϕ) and containing pθ(ϕ) by (2.2) L(ϕ, p) = x ∈ R2 : x · θ(ϕ) = p .
This line is parameterized by t 7→ pθ(ϕ) + tθ ⊥ (ϕ). Note that L(ϕ, p) = L(ϕ + π, −p). We denote the integral of f ∈ L2 (D) over the line L(ϕ, p) by Z ∞ f (pθ(ϕ) + tθ ⊥ (ϕ)) dt. (2.3) Rf (ϕ, p) = −∞
This transform is called the X-ray transform or Radon line transform of f . The sinogram of data Rf is a picture representing the data on [0, 2π] × R where the color (gray scale) value of the pixel at (ϕ, p) represents the value of the data Rf (ϕ, p). For functions g on [0, 2π] × R, the dual Radon transform or backprojection operator is defined as Z 2π ∗ g(ϕ, x · θ(ϕ)) dϕ. (2.4) R g(x) = 0
2.2. Cotangent spaces and microlocal analysis. In this section, we define some important concepts needed to describe singularities in general. Sources, such as [40] and [20], provide an introduction to microlocal analysis in tomography. An introduction to differential geometry, including tangent spaces and cotangent spaces is in [42]. We introduce these ideas because this is the standard notation for singularities and so that singularities can be defined on manifolds, not just R2 . In Remark 2.1 we provide a more elementary way to view these concepts. space of R2 at x ∈ R2 , Tx (R2 ) is the vector space with basis the differential operators o n The tangent ∂ ∂ 2 ∂x1 , ∂x2 . They evaluate on smooth functions at x, and Tx (R ) consists of tangent vectors ~y = y1
∂ ∂ + y2 where (y1 , y2 ) ∈ R2 . ∂x1 ∂x2
The tangent bundle of R2 is the set of all tangent vectors above all points in R2 : ∂ ∂ 2 2 2 + y2 , (y1 , y2 ) ∈ R . T (R ) = (x, ~y ) : x ∈ R , ~y = y1 ∂x1 ∂x2
The cotangent space of R2 at x ∈ R2 , Tx∗ (R2 ) is the dual vector space to Tx (R2 ) which is the set of all linear functionals maps from on Tx (R2 ) to R). For i = 1, 2, let dxi be the linear (linear functional satisfying dxi ∂x∂ j is equal to one if i = j and zero otherwise. Then, {dx1 , dx2 } is a basis of Tx∗ (R2 ). Therefore, an arbitrary vector in Tx∗ (R2 ) is of the form
ξdx := ξ1 dx1 + ξ2 dx2 where ξ = (ξ1 , ξ2 ) ∈ R2 . 4
The cotangent bundle of R2 is the set of all covectors above all points in R2 : T ∗ (R2 ) = (x, ξdx) : x ∈ R2 , ξ ∈ R2 .
To define the cotangent space on [0, 2π] × R we need to make [0, 2π] × R a manifold. To do this, we identify 0 and 2π so that, for each p ∈ R, the point (0, p) and (2π, p) are considered the same. This is equivalent to identifying [0, 2π] with the unit circle S 1 so that [0, 2π] × R is identified with the cylinder S 1 × R. We will consider only functions on [0, 2π] × R that are 2π-periodic in ϕ; thus they are well defined under this identification. Tangent and cotangent bundles are defined in a similar way on other manifolds using coordinates [42]. For example, the cotangent bundle T ∗ ([0, 2π] × R) = (ϕ, p, adϕ + bdp) : (ϕ, p) ∈ [0, 2π] × R, (a, b) ∈ R2 and for (ϕ0 , p0 ) ∈ [0, 2π] × R, the cotangent space above (ϕ0 , p0 ) is ∗ T(ϕ ([0, 2π] × R) = (ϕ0 , p0 , adϕ + bdp) : (a, b) ∈ R2 . 0 ,p0 )
Remark 2.1. As noted at the start of the section, we introduce tangent and cotangent spaces so that we can define singularities and wavefront sets in an invariant way on manifolds (i.e., independent of coordinate changes). However, one can view tangent and cotangent vectors as just vectors, so the tangent vector ~y is identified with (y1 , y2 ) in the plane: y = y1 ~
∂ ∂ + y2 7→ (y1 , y2 ) ∈ R2 . ∂x1 ∂x2
Cotangent vectors ξdx can be identified with ξ = (ξ1 , ξ2 ) ∈ R2 rigorously and (2.5)
ξdx := ξ1 dx1 + ξ2 dx2 7→ (ξ1 , ξ2 ) ∈ R2 .
So, when we say ξdx is conormal to a line L, this means that the vector ξ is perpendicular (normal) to L. 2.3. Wavefront sets and R. The wavefront set is a deep concept that makes singularities of functions precise, and we take the definition from [40]. Definition 2.2 (Wavefront Set). Let (x0 , ξ0 dx) ∈ T ∗ (R2 ) \ 0 (so ξ0 ∈ R2 \ 0). A cutoff function at x0 is any C ∞ function of compact support that is nonzero at x0 . Let f be a locally integrable function. We say f is smooth at x0 in direction ξ0 if there is a cutoff function ψ at x0 and an open cone V containing ξ0 such that the Fourier transform F(ψf )(ξ) is rapidly decaying at infinity (decaying faster at infinity than any power of 1/ |ξ|) for ξ ∈ V . If f is not smooth at x0 in direction ξ0 , then we say (x0 , ξ0 ) is in the wavefront set of f : (x0 , ξ0 ) ∈ WF(f ).
As noted in Remark 2.1, cotangent vectors can be viewed (or identified) as vectors in R2 . This allows one to understand a point (x, ξdx) in the wavefront set, WF(f ) in the following way. The point x is where f is not smooth, and the vector ξ = (ξ1 , ξ2 ) points in a direction in which f is not smooth at x. The wavefront set is conic in ξ (if (x0 , ξ0 ) ∈ WF(f ), then so is (x0 , tξ0 ) for any t > 0) because the neighborhood V of ξ0 that determines microlocal smoothness is conic. We define the wavefront set as a subset of the cotangent bundle because it transforms under diffeomorphism as the cotangent bundle does. Note that we can define the wavefront set above (ϕ, p) ∈ [0, 2π] × R locally by choosing a cutoff function that is only nonzero in a small neighborhood of (ϕ, p). Our next definition allows us to express efficiently the correspondence between wavefront set of f and that of Rf . 5
Definition 2.3. Let (ϕ, p) ∈ [0, 2π] × R and let f be a locally integrable function. The conormal space, N ∗ (L(ϕ, p)) \ 0, of the line L(ϕ, p) is defined to be (2.6) N ∗ (L(ϕ, p)) \ 0 = (x, ωθ(ϕ)dx) : x ∈ L(ϕ, p), ω 6= 0 and
(2.7)
WFL(ϕ,p) (f ) = WF(f ) ∩ (N ∗ (L(ϕ, p)) \ 0) .
We will say f is smooth conormal to the line L(ϕ, p) if WFL(ϕ,p) (f ) = ∅. For x0 ∈ R2 , we let WFx0 (f ) = WF(f ) ∩ Tx∗0 (R2 ). Now, let g be a locally integrable function on [0, 2π] × R and (ϕ, p) ∈ [0, 2π] × R. We define
(2.8)
∗ WF(ϕ,p) (g) = WF(g) ∩ T(ϕ,p) ([0, 2π] × R).
It is important to understand each set introduced in the last definition: N ∗ (L(ϕ, p)) \ 0 is the set of all (x, ξdx) such that x ∈ L(ϕ, p) and the vector ξ is normal to L(ϕ, p) at x. So, WFL(ϕ,p) (f ) is the set of wavefront directions of f above points x ∈ L(ϕ, p) that are conormal to this line. If g is a locally integrable function on [0, 2π] × R, WF(ϕ,p) (g) is the set of wavefront directions with base point (ϕ, p). These sets have an important relationship that we will exploit starting in the next proposition. Proposition 2.4. Let f ∈ L2 (D) and let (ϕ0 , p0 ) ∈ [0, 2π] × R. Then,
(2.9)
WFL(ϕ0 ,p0 ) (f ) 6= ∅ if and only if WF(ϕ0 ,p0 ) (Rf ) 6= ∅.
That is, f is smooth conormal to the line L(ϕ0 , p0 ) if and only if Rf is a smooth function near (ϕ0 , p0 ). This proposition says that singularities of Rf above (ϕ0 , p0 ) are caused by singularities of f on L(ϕ0 , p0 ) that are conormal to the line. The proof will be given in the appendix at the end of section A.1, just after the proof of Proposition A.3. 3. Main Results In contrast to limited angle characterizations in [10], our main results characterize arbitrary incomplete data situations, including classical problems such as region of interest tomography, exterior tomography, and limited angle tomography. Our results are formulated with respect to the wavefront set (Definition 2.3) because it precisely describes singularities. We start with the basic setup. By definition, incomplete CT data are given over a closed, proper, nonempty subset A ⊂ [0, 2π] × R. We will sometimes refer to the region A over which data are given as the cutoff region. If A is a subset of [0, 2π]×R, then the interior of A, int(A), is the set of all points (ϕ0 , p0 ) ∈ A such that there is an open ball in (ϕ, p)-space centered at (ϕ0 , p0 ) that is contained in A. The boundary of A is the set of all points in (ϕ0 , p0 ) ∈ [0, 2π] × R such that every open ball centered at (ϕ0 , p0 ) contains points of A and points of the complement, ([0, 2π] × R) \ A, and ext(A) is everything else: ext(A) = [0, 2π] × R \ (int(A) ∪ bd(A)). In calculating these sets for points with ϕ = 0 and ϕ = 2π, we identify 0 and 2π. To avoid trivial counterexamples, we assume bd(int(A)) = bd(A). Because of the symmetry condition on lines L(ϕ, p) = L(ϕ + π, −p), we will assume the following symmetry condition on A: (3.1)
if (ϕ, p) ∈ A then (ϕ + π, −p) ∈ A.
In many situations (in particular in most of the commercial devices), reconstructions from incomplete CT data are obtained by using classical algorithms such as filtered backprojection (see [30] for a practical discussion of FBP in commercial scanners). In this case, the incomplete data is 6
implicitly extended to all of [0, 2π] × R by padding it with zeros off of A. Then, one can view the incomplete data as (3.2)
RA f (ϕ, p) = 1A (ϕ, p)Rf (ϕ, p),
where 1A is the characteristic function of A. Therefore, RA f (ϕ, p) gives the measured data on A and zero off of A. This would be represented in the sinogram by gray values corresponding to values of Rf (ϕ, p) for (ϕ, p) ∈ A and zero for (ϕ, p) ∈ / A. Using filtered backprojection (FBP) on this padded data is a standard reconstruction method that gives rise to the reconstruction operator: (3.3)
LA f = R∗ (ΛRA f ) = R∗ (Λ1A Rf )
where Λ is the standard FBP filter (see e.g., [24, Theorem 2.5] and [25, §5.1.1] for numerical implementations) and R∗ is defined by (2.4). Our first theorem clarifies that singularities of 1A and of Rf above (ϕ0 , p0 ) can create singularities and artifacts of LA f on N ∗ (L(ϕ0 , p0 )) \ 0 and not elsewhere. Theorem 3.1 (Localization of Singularities and Artifacts). Assume A ⊂ [0, 2π] × R is closed and satisfies the symmetry condition (3.1). Let f ∈ L2 (D) and (x0 , ξ0 dx) ∈ T ∗ (R2 ) \ 0. Choose ϕ0 ∈ [0, 2π] and ω0 6= 0 such that ξ0 = ω0 θ(ϕ0 ). Finally, let p0 = x0 · θ(ϕ0 ). (i) If (x0 , ξ0 dx) ∈ WF(LA f ) then there is a singularity of either Rf or of 1A above (φ0 , p0 ) that caused it. (ii) Singularities of Rf or of 1A above (ϕ0 , p0 ) only cause singularities or artifacts in LA f at points on L(ϕ0 , p0 ) and directions conormal to this line (and not at other points or other directions). The above theorem shows that all reconstructed singularities in LA f are generated either by singularities in the data Rf or by singularities that are generated by the multiplication of Rf with 1A (hard truncation of the data). In the latter case, the reconstructed singularities are in general not present in the original object which is why they are called artifacts. In what follows, our goal is to characterize those artifacts for general incomplete data sets A ⊂ [0, 2π] × R. Definition 3.2. Let A ⊂ [0, 2π] × R. Geometrically, we think of A being represented by part of a Cartesian coordinate system with φ on the horizontal axis and p on the vertical axis. Now, let (ϕ0 , p0 ) ∈ bd(A). We say that bd(A) is smooth near (ϕ0 , p0 ) if this boundary is a C ∞ curve near (ϕ0 , p0 ). In this case, there is a unique tangent line in (ϕ, p)-space to bd(A) at (ϕ0 , p0 ). If this tangent line is vertical (i.e., of the form ϕ = ϕ0 ), then we say the boundary has infinite slope at (ϕ0 , p0 ). If this tangent line is not vertical, then bd(A) is defined near (ϕ0 , p0 ) by a smooth function p = p(ϕ). In this case, the slope of the boundary at (ϕ0 , p0 ) will be the slope of this line, p′ (ϕ0 ). If bd(A) is not a smooth curve near (ϕ0 , p0 ), then we say bd(A) is not smooth at (ϕ0 , p0 ). We say that bd(A) has a corner at (ϕ0 , p0 ) if the curve bd(A) is continuous at (ϕ0 , p0 ), is smooth at other points sufficiently near (ϕ0 , p0 ), but has different one-sided slopes at (ϕ0 , p0 ). 3.1. Visible and invisible singularities. Our next theorem gives an analysis of singularities in LA f corresponding to lines not in bd(A). It shows that there are added streak artifacts along lines L(ϕ, p) only when (ϕ, p) ∈ bd(A).
Theorem 3.3 (Visible and Invisible Singularities in the Reconstruction). Let f ∈ L2 (D) and let A be a closed subset of [0, 2π] × R satisfying the symmetry condition (3.1). If (ϕ, p) ∈ / bd(A) then LA f does not have streak artifacts along L(ϕ, p), i.e., the only singularities of LA f at covectors conormal to those lines are singularities of f . (i) If (ϕ, p) ∈ int(A) then WFL(ϕ,p) (f ) = WFL(ϕ,p) (LA f ). (ii) If (ϕ, p) ∈ / A then WFL(ϕ,p) (LA f ) = ∅. 7
(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 LA f : (3.4)
WFx (f ) = WFx (LA f ).
Remark 3.4 (ROI CT). We now apply Theorem 3.3 to show that in region of interest (ROI) tomography, all singularities of f in the interior of the ROI are recovered. This is well-known (e.g., [34]), and we revisit the ROI problem in Section 4.5. Let K be a closed convex subset of D and let A be the set of all (ϕ, p) for lines meeting K (i.e., (ϕ, p) ∈ A if and only if L(ϕ, p) ∩ K 6= ∅). Let x ∈ int(K). Then, every line through x (i.e., every line L(ϕ, x · θ(ϕ)) for all ϕ ∈ [0, 2π]) is in int(A). So, the condition of Theorem 3.3 part (iii) holds and (3.4) holds so all singularities of f at x are visible in LA f . 3.2. Characterization of artifacts. We now state our main results. In what follows, we assume the cutoff region A has nontrivial interior so that 1A is not zero as a distribution. We let (ϕ0 , p0 ) be an arbitrary point in bd(A) and evaluate the behavior of the reconstruction depending on the geometric properties of bd(A) near (ϕ0 , p0 ). We show that artifacts are created in incomplete data FBP reconstructions depending on the local smoothness of bd(A) and, in some cases, singularities of the object f . Theorem 3.1 and Theorem 3.3 demonstrate that artifacts appear only on lines L(ϕ, p) for (ϕ, p) ∈ bd(A) (i.e., with wavefront conormal to such lines). Our next set of theorems show that artifacts arise in two ways: • Object-dependent Artifacts: artifacts caused by singularities of the object f that are conormal to L(ϕ, p) for some (ϕ, p) ∈ bd(A).
• Object-independent Artifacts: artifacts caused essentially only by the geometry of bd(A). The presentation of our results is organized as follows. I. Boundary of A is smooth at (ϕ0 , p0 ) ∈ bd(A):
A. bd(A) has infinite slope (i.e., vertical tangent line) or large slope (Theorem 3.5). In this case, the only artifacts that appear in the reconstruction region D are objectdependent and they will be along the line L(ϕ0 , p0 ). This case includes limited angle CT. B. bd(A) has small slope at (ϕ0 , p0 ) ∈ bd(A) (Theorem 3.7). This includes region of interest (ROI) CT. In this case, two types of artifacts can appear in the reconstruction region D: (i) Artifacts on L(ϕ0 , p0 ): The same object-dependent artifacts as in I.A. (ii) Artifacts on Curves: Smooth boundary at (ϕ0 , p0 ) with small slope. These artifacts are generated by bd(A) and are object-independent.
II. Boundary of A is not smooth at (ϕ0 , p0 ) ∈ bd(A) (Theorem 3.8): In this case, an artifact line can run along L(ϕ0 , p0 ). It is object-independent and can occur even if f is smooth. This case explains the artifacts in the reconstructions from the synchrotron data in Figure 1 and Section 7.
3.2.1. Locally Smooth Cutoff Boundary with Large or Infinite Slope. As a first generalization of the limited angle characterizations in [10] we consider the case where the incomplete data A ⊂ [0, 2π]×R has a smooth boundary bd(A) locally near (ϕ0 , p0 ) and a large slope. 8
Theorem 3.5 (Locally Smooth Boundary, Large Slope). Let f ∈ L2 (D) and A ⊂ [0, 2π] × R be closed and satisfy the symmetry condition (3.1). Moreover, let (ϕ0 , p0 ) ∈ bd(A) and assume that bd(A) is smooth near (ϕ0 , p0 ). Assume either that the boundary curve is given by p = p(ϕ) locally near (ϕ0 , p0 ) and it has a relatively large slope at (ϕ0 , p0 ), i.e. ′ q (3.5) p (ϕ0 ) > 1 − p20 , or bd(A) has a vertical tangent at (ϕ0 , p0 ) (see Definition 3.2). (i) If f has a singularity conormal to L(ϕ0 , p0 ) (i.e., for some x ∈ L(ϕ0 , p0 ), (x, λθ(ϕ0 )) ∈ WF(f )), then WFL(ϕ0 ,p0 ) (LA f ) ⊂ N ∗ (L(ϕ0 , p0 )) \ 0.
This means that a streak artifact can be generated all along L(ϕ0 , p0 ) by the reconstruction operator LA . (ii) If f has no singularity conormal to L(ϕ0 , p0 ) then LA f is smooth conormal to L(ϕ0 , p0 ). That is, WFL(ϕ0 ,p0) (f ) = ∅ implies that WFL(ϕ0 ,p0) (LA f ) is empty above points in D. Hence, there is no artifact along (i.e., conormal to) L(ϕ0 , p0 ) in D.
We explain the slope condition (3.5) in Section 4.2 and in particular Remark 4.1. Note that for large slope, the only artifacts that occur in D related to the boundary point (ϕ0 , p0 ) are objectdependent since they occur only if f has a singularity conormal to L(ϕ0 , p0 ). This theorem is given in the appendix as parts I.A. and I.B. of Theorem A.6, and its proof is given there. We will present reconstructions illustrating this case in Section 4.1 and we provide a quick illustration of the theorem in the following remark. Remark 3.6 (Limited Angle Tomography). A classical problem where Theorem 3.5 applies is given by limited angle tomography. In this case the boundary curve is a vertical line given by ϕ = ϕ0 (infinite slope). The results of Theorem 3.5 generalize the results of [10, 12] as they also apply to cutoffs with non-vertical slope. Taking a closer look at the statement of Theorem 3.5 and the results of [10,12] one can observe that, locally, they describe the same phenomena, namely: whenever there is a line L(ϕ0 , p0 ) in the data set with (ϕ0 , p0 ) ∈ bd(A) and which is conormal to a singularity of f , then a streak artifact can be generated in the reconstruction LA f . It is important to note that in this case, the only artifacts are object-dependent: a streak artifact is generated only if f has a singularity conormal to a line L(ϕ0 , p0 ) with (ϕ0 , p0 ) ∈ bd(A). In fact, a streak artifact will be generated all along L(ϕ0 , p0 ) if 1A Rf has wavefront set in all directions at (ϕ0 , p0 ). This shows that the generation of streak artifacts depends on the object that is imaged. This is a classical phenomenon known from limited angle tomography and Theorem 3.5 shows that it is still valid in this more general setting. 3.2.2. Locally Smooth Cutoff Boundary with Small Slope. Here, we present a second generalization of results in [10–12] by considering incomplete data sets A ⊂ [0, 2π] × R with a smooth boundary bd(A) with a small slope. In contrast to the large slope situation (as presented in Theorem 3.5), the next theorem shows that in this case there can be singularities in LA f that are generated independently of the object f . Theorem 3.7 (Locally Smooth Boundary, Small Slope). Let f ∈ L2 (D) and A ⊂ [0, 2π] × R be closed and satisfies the symmetry condition (3.1). Moreover, let (ϕ0 , p0 ) ∈ bd(A) and assume that bd(A) is smooth near (ϕ0 , p0 ) and is defined locally near (ϕ0 , p0 ) by p = p(ϕ). Assume the boundary curve has a relatively small slope at (ϕ0 , p0 ), i.e. ′ q p (ϕ0 ) ≤ 1 − p2 . (3.6) 0 9
Let I be a neighborhood of ϕ0 such that p = p(ϕ) describes bd(A) near (ϕ0 , p0 ). For ϕ ∈ I let
(3.7)
xb = xb (ϕ) = p(ϕ)θ(ϕ) + p′ (ϕ)θ ⊥ (ϕ).
(i) If f has a singularity conormal to L(ϕ0 , p0 ), then WFL(ϕ0 ,p0 ) (LA f ) ⊂ N ∗ (L(ϕ0 , p0 )) \ 0.
This means that a streak artifact can be generated all along L(ϕ0 , p0 ) by the reconstruction operator LA . (ii) If f has no singularity conormal to L(ϕ0 , p0 ) then (3.8) WFL(ϕ0 ,p0 ) (LA f ) ⊂ (xb , ωθ(ϕ0 )) : ω 6= 0 .
That is, LA f is smooth conormal to L(ϕ0 , p0 ), except perhaps at xb (ϕ0 ). • If, in addition, Rf (ϕ0 , p0 ) 6= 0, then equality holds in (3.8). • If Rf is zero near (ϕ0 , p0 ), then there will be no singularities conormal to L(ϕ0 , p0 ), WFL(ϕ0 ,p0) (LA f ) = ∅.
Theorem 3.7 is covered in the appendix as part I.A. of Theorem A.6 and its proof is given there. We will present reconstructions illustrating this case in Section 4.3. In both the large slope (Theorem 3.5) and small slope case (Theorem 3.7), there are object dependent artifacts. If (ϕ0 , p0 ) ∈ bd(A) and f has a singularity conormal to L(ϕ0 , p0 ), then a streak artifact can occur on L(ϕ0 , p0 ). The only difference between Theorem 3.5 and Theorem 3.7 is in item 3.7(ii). This item states that only when the boundary curve is smooth and has small slope, can there be additional singularities in the reconstruction in D. They are on the curve ϕ 7→ xb (ϕ) (see (3.7)), and they are objectindependent (essentially independent of f ). This phenomenon was not explained by the previous mathematical theory. We will analyze the points xb more carefully in Section 4.2 and explain the slope conditions 3.5 and 3.6. 3.2.3. Locally Nonsmooth Boundary. We now consider the case where the incomplete data A ⊂ [0, 2π] × R has nonsmooth boundary. The next theorem shows that in this situation there can be streak artifacts in LA f along L(ϕ0 , p0 ) that are independent of the object f . Theorem 3.8 (Nonsmooth Boundary). Let f ∈ L2 (D) and A ⊂ [0, 2π] × R be closed (with nonempty interior) that satisfy the symmetry condition (3.1). Moreover, let (ϕ0 , p0 ) ∈ bd(A) and assume that bd(A) is not smooth at (ϕ0 , p0 ). Then, WFL(ϕ0 ,p0 ) (LA f ) ⊂ N ∗ (L(ϕ0 , p0 )) \ 0.
Thus, LA f can have a streak artifact on L(ϕ0 , p0 ) independent of f . If f is smooth conormal to L(ϕ0 , p0 ), Rf (ϕ0 , p0 ) 6= 0 and bd(A) has a corner at (ϕ0 , p0 ), then LA f does have a streak artifact all along L(ϕ0 , p0 )–equality holds in (A.19).
This theorem is part II. of Theorem A.6 and its proof is in section A.5 of the appendix. These line artifacts 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 present a stronger version of the parts of Theorems 3.7 and 3.8 when Rf (ϕ0 , p0 ) 6= 0 and f is smooth conormal to L(ϕ0 , p0 ) because it gives the strength in Sobolev scale of the added singularities in these cases. We point out that we have taken care of all cases. Namely, either bd(A) is smooth near (ϕ0 , p0 ) and either it has large slope or is vertical at (ϕ0 , p0 ) and Theorem 3.5 applies, bd(A) has small slope at (ϕ0 , p0 ) and Theorem 3.7 applies, or bd(A) is not smooth at (ϕ0 , p0 ) ∈ bd(A) and Theorem 3.8 applies. 10
4. Numerical Illustrations of our Theoretical Results To show how our theoretical results are reflected in practice, we now analyze reconstructions of the Shepp Logan phantom from incomplete data. We consider a range of well-known incomplete data problems as well as unconventional ones in order to demonstrate the generality of our theoretical results. The set A we use in Sections 4.1, 4.2, and 4.3 is defined as follows. It will be cut in the middle so that the left-most cut-off boundary occurs at ϕ = a = 94 π; the right-most boundary is constructed as ϕ = b = 59 π for p ≤ 0 and p (4.1) p(ϕ) = c ϕ − b, ϕ > b
for p > 0 such that the two parts join smoothly at (ϕ, p) = (0, 0). The steepness of the curved part of the right-most boundary is governed by the constant c, in particular c , ϕ > b. (4.2) p′ (ϕ) = √ 2 ϕ−b According to Theorems 3.5 and 3.7, the only part of the data boundary that can potentially cause object-independent reconstruction artifacts is the curved part, since the other parts are vertical. The small slope condition (3.6) for potential existence of added artifacts can be evaluated for the parametrized curve by substituting expressions (4.1) and (4.2) into the condition for xb (ϕ) to be in D, (3.6). The resulting expression can simplified and recognized as a quadratic inequality in ϕ0 − b, and, for c ∈ [0, 1], the inequality has a real solution when # " √ √ 1 + 1 − c4 1 − 1 − c4 ,b + , for c ∈ [0, 1], (4.3) ϕ0 ∈ [ϕ1 , ϕ2 ] = b + 2c2 2c2 where ϕ1 and ϕ2 represent the end points of interval. If ϕ2 > π, then we set ϕ2 = π because the symmetry condition (3.1) determines A on [π, 2π] × R from A on [0, π] × R.
4.1. Smooth boundary with large slope. In this section, we consider the case c = ∞ (limited angle–vertical boundary curves) and the large slope case c = 1.3. Because c > 1 in both cases, there will be no artifact curves ϕ 7→ xb (ϕ) visible in D. The only artifacts will be of the object-dependent type, caused by singularities of f conormal to lines L(ϕ, p) for (ϕ, p) ∈ bd(A). In Figure 2(A), 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 set is conormal 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 the same angles. In Figure 2(B), most artifacts are as with Figure 2(A) because the boundaries of the cutoff regions are the same at those point. However, the upper part of the right boundary, which corresponds to ϕ > 5π/9 and p > 0, are different. This discussion is illustrated by the four circled points on the sinogram and the corresponding artifacts: 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(A) because 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. 4.2. The role of xb . We now analyze the points xb and the resulting curved artifacts. These points can occur when bd(A) is smooth and not vertical near (ϕ0 , p0 ). In this case, one can define bd(A) locally near (ϕ0 , p0 ) by a smooth function p = p(ϕ) for ϕ in some open neighborhood I of ϕ0 : I ∋ ϕ 7→ (ϕ, p(ϕ)) ∈ bd(A). 11
(A) 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.
(B) Left: Boundary of A with large positive slope (boundary defined by (4.1) for c = 1.3). Center: FBP reconstruction. Right: Reconstruction highlighting object-dependent artifact lines tangent to the skull corresponding to the four circled points in the sinogram.
Figure 2. Left: Incomplete data sinograms for which bd(A) has large or infinite slope. Note that the part of the sinogram for ϕ ∈ [0, π] is shown. Center and Right: Reconstructions from the sinograms on the left. All streak artifacts in the reconstructions are object-dependent: they are generated by points (ϕ0 , p0 ) ∈ bd(A) for which L(ϕ0 , p0 ) is tangent to some boundary of f (i.e., conormal to the associated singularity). A closer look at proofs of Theorems 3.5 and 3.7 (Theorem A.6 part I. and its proof, in particular, the derivation (A.22)) reveals that whenever bd(A) is smooth with finite slope at (ϕ, p(ϕ)), then an additional singularity can be generated at xb given by (3.7): (4.4)
I ∋ ϕ 7→ xb (ϕ) = p(ϕ)θ(ϕ) + p′ (ϕ)θ ⊥ (ϕ).
This means that a curve of artifacts ϕ 7→ xb (ϕ) can be generated for ϕ ∈ I. Note that for each ϕ ∈ I, xb (ϕ) is a point on the line L(ϕ, p(ϕ)) and any artifact at xb (ϕ) is object-independent since it can occur even if f is smooth conormal to L(ϕ, p(ϕ)) as noted in Theorem 3.7, part (ii). Figure 3 shows the sinogram with boundary given by (4.1) and c = 0.5. A specific point (ϕ0 , p0 ) in the curved part of bd(A) is circled and the conormal to the boundary at (ϕ0 , p0 ) is shown. This conormal is in WF(1A ) and, as shown in the derivation (A.22), it causes the artifact at xb (ϕ0 ). The 12
right-hand picture in Figure 3 shows the reconstruction including the corresponding point xb (ϕ0 ) and the line L(ϕ0 , p0 ) that it is on. In the case of large slope, the expression for xb , (4.4), and inequality (3.5) imply that xb 6∈ D. Therefore, this singularity is not seen the reconstruction in the region D. However, this point will be visible if one reconstructs a sufficiently large region. In the case of Theorem 3.7, inequality (3.6) holds, and the point xb (ϕ0 ) will be in the unit disk, D. As p′ (ϕ0 ) gets larger in magnitude, (4.4) shows that xb becomes farther from the origin and one might suggest that for vertical tangents (“p′ (ϕ0 ) = ∞”), xb given by (4.4) is a point at infinity. Remark 4.1 (Recap). To paraphrase this discussion, basically the curve artifacts ϕ 7→ xb (ϕ) occurs in the large slope case (Theorem 3.5 with non-vertical slope) as well as the small slope case (Theorem 3.7). However, one sees the curve artifacts only if one reconstructs a large enough region to see them. This explains the slope conditions (3.5) and (3.6): those conditions determine if the artifact point xb (ϕ0 ) is in the reconstruction region D or not.
Figure 3. Left: A data set, A, with smooth boundary. The set A is described at the start of this section and has constant c = 0.5. A specific boundary point (ϕ0 , p0 ) is circled and the perpendicular vector (−p′ (ϕ0 ), 1) is drawn. The associated covector η0 = (ϕ0 , p0 , −p′ (ϕ0 )dϕ + dp) is in WF(1A Rf ) because η0 is conormal to bd(A). The singularity at xb (ϕ0 ) comes from the singularity of 1A Rf at η0 as shown in the proof of Theorem A.6 and in particular the calculation (A.22). Right: The reconstruction showing the location of the corresponding artifact xb (ϕ0 ) and the line L(ϕ0 , p0 ).
4.3. Smooth boundary with small slope. In this section, we consider two cutoff regions A defined by (4.1) with c = 0.82 and c = 0.5. In both cases, there are ϕ, given by (4.3), for which the point xb (ϕ) is in D. We show the resulting reconstruction in Figure 4(A) for c = 0.82 and in Figure 4(B) for c = 0.5. For each of these two values of c, let ϕ1 and ϕ2 be the limit angles in (4.3). The sinograms on the left of Figure 4 show the points on the boundary for ϕ ∈ [ϕ1 , ϕ2 ] that potentially cause added artifacts as cyan curves. Since artifacts can only arise from points (ϕ0 , p0 ) in the support of the sinogram by Theorem 3.3 part (ii), we distinguish between the relevant part of the boundary 13
(solid curve) and the part outside the support of 1A Rf (dashed curve). The center figures are the reconstructions without added markings so the reader clearly sees the artifacts. We take the analysis one step further and construct an explicit parametrization of the added artifacts in the reconstructions on the right in Figure 4. We evaluate (4.4) and generate the curve ϕ 7→ xb (ϕ) which is shown in cyan on the right-hand reconstructions in Figure 4, again using solid and (resp.) dashed lines corresponding to the sinogram boundaries in and (resp.) out of supp(Rf ). In both cases the cyan artifact curves closely trace the visible added artifacts in their full extent. For c = 0.82, the artifact curve disappears at the points xb (ϕ) for which (ϕ0 , p0 ) ∈ bd(A) is outside the support of 1A Rf . This occurs because WFL(ϕ,p) (LA f ) = ∅ for such (ϕ0 , p0 ) by Theorem 3.3 part (ii). This is indicated on the right-hand reconstruction in Figure 4(A) by a dashed cyan curve that curves down. There is a nearby line artifact caused by the skull. Note that object-dependent artifacts similar to those in Figure 2 show up in the reconstruction, and they are most prominent at the boundary of the skull corresponding to the points in the sinogram at the intersection of the boundary of the support of Rf and bd(A). The two vertical artifact lines on the left side of the reconstruction in Figure 4(B) are explained when one extends the set A to [π, 2π] using p the symmetry condition (3.1) as shown in Figure 5. The curve (4.1) for c = 0.5 ends at (π, 0.5 4π/9) because the symmetry condition determines A on [π, 2π] × [−1, 1]. This forces bd(A) to include a vertical segment o n p ℓ = (π, p) : p > 0.5 4π/9 .
p The yellow circle at ϕ = π and repeated at ϕ = 2π represents the resulting corner at (π, 0.5 4π/9) between the square root boundary curve and ℓ. By Theorem 3.8, the vertical line artifact highlighted in yellow is an object-independent artifact that is caused by this corner. The vertical line artifact tangent to the left side of the skull is an object-dependent artifact. As can be seen from Figure 5, it is a line L(π, p0 ) where (π, p0 ) is in the vertical boundary, ℓ. Since L(π, p0 ) is tangent to the skull, and (π, p0 ) ∈ bd(A), Theorem 3.5(i) justifies that this line is an object-dependent artifact line.
4.4. Nonsmooth boundary. Theorem 3.8 explains that, if bd(A) is not smooth at (ϕ0 , p0 ), then an object-independent artifact can be generated all along L(ϕ0 , p0 ). In Figure 6, we show a simple cutoff region with two corners and we plot these lines and confirm that they agree precisely with the observed artifact lines. Since the slope of the cut-off is numerically large (±5), which exceeds the maximal attainable value of the right hand side of (3.5), we know from Theorem 3.5 that there cannot be additional curve artifacts added in D. This is in agreement with none being seen in the reconstruction. Note that object-dependent artifacts similar to those in Figure 2 show up in the reconstruction, and they are most prominent at the boundary of the skull corresponding to the points in the sinogram at the intersection of the boundary of supp(Rf ) and bd(A). 4.5. Region of interest (ROI) tomography. A classical incomplete-data tomography problem where Theorem 3.7 applies is the ROI problem, also known as interior tomography, in which the detector width is not large enough to contain the full projection, see Figure 7. Here we demonstrate how our theoretical results predict precisely which reconstruction artifacts occur for two choices of detector width. The boundary of the sinogram are given by horizontal lines p = ±p0 where p0 = 0.4 in Figure 7(A) and p0 = 0.8 in Figure 7(B). The points given by (4.4) are xb (ϕ) = p0 θ(ϕ) + 0θ ⊥ (ϕ) since p′ = 0, i.e., a circle of radius p0 . The ring around the ROI in Figure 7(A) is exactly the curve ϕ 7→ xb (ϕ) for this ROI. 14
(A) Smooth boundary with part of small slope (c = 0.82) causing curve of artifacts, ϕ 7→ xb (ϕ) in D.
(B) Smooth boundary with smaller slope (c = 0.5) causing a more prominent curve of artifacts in D.
Figure 4. Illustration of artifacts with smooth boundary and small slope and boundary given by (4.1). Left: Sinogram with part of boundary that potentially causes artifacts in cyan. The part that actually causes artifacts (where the sinogram is nonzero) is a solid curve, and the part outside supp(Rf ) is a dashed curve. Center: Reconstructions without added markings. Right: The cyan curve in the reconstruction is the curve of artifacts: ϕ 7→ xb (ϕ). The solid part indicates the artifacts that are realized in the picture. The dashed curve indicates potential artifacts that are not realized because the corresponding part of bd(A) is outside supp(Rf ) (see Theorem 3.3(ii)). More generally, if the ROI is smooth and strictly convex then an exercise using the correspondence between tangent lines and their parametrization in (ϕ, p) shows that the points xb in (4.4) at which these point artifacts might occur trace the boundary of the ROI. Analyzing the reconstruction in Figure 7(B), one sees the top and bottom of a circle analogous to the one in Figure 7(A). However, the sides of the circle are not visible in this larger ROI because Rf is zero near the lines corresponding to p = ±0.8, ϕ < 0.7121 and ϕ > 2.429, and therefore the reconstruction is smooth there (see the case when Rf = 0 in Theorem 3.7 part (ii)). One also sees object-dependent artifacts described by Theorem 3.7 part (i) in Figure 7(B). These occur along lines L(ϕ0 , p0 ) where p0 = ±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 points (ϕ0 , p0 ) 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 conormal to 15
Figure 5. Left: The sinogram from Figure 4(B) extended to ϕ ∈ [0, 2π] and showing the square root boundary curve in cyan and with the corner at ϕ = π circled. Right: the reconstruction in Figure 4(B) with one vertical artifact line highlighted.
Figure 6. The nonsmooth boundary case covered in Theorem 3.8. Left: Sino8 11 gram with two corners (circled), which are at (ϕ0 , p0 ) = ( 18 π, 0) and ( 18 π, 0.3). Center: Reconstruction from LA f . Somewhat faint artifact lines (compared to the object-dependent artifacts at the boundary of the skull) occur in the reconstruction corresponding to the two corners in the boundary of the sinogram. Right: reconstruction with the “corner” artifacts highlighted. 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. All features of the object inside the ROI are visible in our reconstruction as predicted in Remark 3.4. 4.6. The General case. The reconstruction in Figure 8 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 8 11 7 π, 18 π are displayed in solid the horizontal sinogram boundaries at p = p0 = ±0.35 for φ ∈ 18 cyan line. As in the ROI case, along these boundaries, p′ = 0 and thus circular arcs of radius p0 for the interval for ϕ are added in the reconstruction as indicated by solid cyan. 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. 16
(A) ROI data with p ∈ [−0.4, 0.4].
(B) ROI data with p ∈ [−0.8, 0.8].
Figure 7. 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 ϕ ∈ [0, π] is shown. 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.5). In Figure 8, there are other object-dependent streaks corresponding to the vertical lines in the 11π sinogram at ϕ = 7π 18 and at ϕ = 18 as predicted by Theorem 3.5(i), but they are less pronounced and more difficult to see. 4.7. Summary. We have presented reconstructions that illustrate all of types of incomplete data and each of our theorems from Section 3. To summarize, all artifacts arise because of points (ϕ0 , p0 ) ∈ bd(A), and they fall into two categories.
• Streak artifacts along the line L(ϕ0 , p0 ) when (ϕ0 , p0 ) ∈ bd(A): – Object-dependent streaks when bd(A) is smooth at (ϕ0 , p0 ) and a singularity of f is conormal to L(ϕ0 , p0 ). – Object-independent streaks when bd(A) is nonsmooth at (ϕ0 , p0 ). • Artifacts on curves: – They are object-independent, and they are generated by the map ϕ 7→ xb (ϕ) from parts of bd(A) that are smooth and of small slope. 17
Figure 8. Left: The sinogram for a general incomplete data problem in which the cutoff region, A, has a locally smooth boundary with zero and infinite slope as well 11π as corners. The cutout from the sinogram is at 7π 18 and 18 , p = ±0.35. Center: FBP reconstruction. Right: Reconstruction with the circular xb curve artifacts highlighted in cyan and object-independent “corner” line artifacts highlighted. Furthermore, by Theorem 3.1 any artifacts caused by (ϕ0 , p0 ) ∈ bd(A) show up only on the line L(ϕ, p). 5. Strength of Added Artifacts Using 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 [32,38]. We state it for distributions, but functions f ∈ L2 (D) are just distributions.
Definition 5.1. For s ∈ R, the Sobolev space Hs (Rn ) is the set of all distributions with locally square integrable Fourier transform and with finite Sobolev norm: 1/2 Z < ∞. |Ff (y)|2 (1 + kyk2 )s dy (5.1) kf ks := y∈Rn
Let f have locally square integrable Fourier transform and let x0 ∈ Rn and ξ0 ∈ Rn \ 0. We say f is in H s at x0 in direction ξ0 if there is a cutoff function ψ at x0 and an open cone V containing ξ0 such that the localized and microlocalized Sobolev seminorm 1/2 Z 2 2 s < ∞. |F (ψf ) (y)| (1 + kyk ) dy (5.2) kf ks,ψ,V := y∈V
On the other hand, if (5.2) does not hold for any cutoff at x0 , ψ or conic neighborhood V of ξ0 , then (x0 , ξ0 dx) ∈ WFs (f ), the Sobolev wavefront set of f of order s. If (x0 , ξ0 dx) ∈ / WF(f ), then an exercise using definitions shows that (x0 , ξ0 dx) ∈ / WFs (f ) for all s ∈ R. Note that WFs can be defined on [0, 2π] × R. Let g be a distribution on [0, 2π] × R. If ψ ∈ D([0, 2π] × R) has sufficiently small support, then ψg can be viewed as a function on an open subset of R2 and then one can use the definition of Hs norm on this open subset of R2 . Note that this norm on distributions on [0, 2π] × R is not the typical H0,s norm used in elementary continuity proofs for the Radon transform (see e.g., [16, equation (2.11)]), but this is the appropriate space for the continuity theorems for general Fourier integral operators [17, Theorem 4.3.1], [7, Corollary 4.4.5]. 18
Our next theorem gives the strength in Sobolev scale of added singularities of LA f under certain assumptions on f . It uses the relation between microlocal Sobolev strength of f and of Rf , [34, Theorem 3.1] and of g and R∗ g, which is given in the useful Theorem A.8. Theorem 5.2. Let f ∈ L2 (D) and let A be a closed symmetric subset of [0, 2π] × R. Let (ϕ0 , p0 ) ∈ bd(A) and assume Rf (ϕ0 , p0 ) 6= 0 and f is smooth conormal to L(ϕ0 , p0 ), i.e., WFL(ϕ0 ,p0 ) (f ) = ∅. (i) Assume bd(A) is smooth at (ϕ0 , p0 ). (a) Assume bd(A) is non-vertical near (ϕ0 , p0 ) and let α be the slope of bd(A) at (ϕ0 , p0 ) as a function of ϕ. Let xb = xb (ϕ0 , p0 , α) be given by (3.7). Then, LA f is in Hs for s < 0 at η0 = (xb , tθ(ϕ0 )dx) and η0 ∈ WF0 (LA f ) for t 6= 0. Thus, there are singularities above xb in the 0-order wavefront set of LA f . (b) If bd(A) has a vertical tangent at (ϕ0 , p0 ), then, under the smoothness assumption on f , there are no added artifacts in LA f conormal to L(ϕ0 , p0 ). (ii) Assume bd(A) comes to a corner at (ϕ0 , p0 ). Then for each (x, ξ) ∈ N ∗ (L(ϕ0 , p0 )) \ 0, (x, ξ) ∈ WF1 (LA f ) and, except for two points on L(ϕ0 , p0 ), LA f is in Hs for s < 1 at (x, ξ). (iii) Now, assume (ϕ0 , p0 ) ∈ int(A). Then, for every (x, ξ) ∈ N ∗ (L(ϕ0 , p0 )) \ 0, (x, ξ) ∈ WFs (f ) if and only if (x, ξ) ∈ WFs (LA f ). Part (i) of this theorem is a more precise version of the last assertion of Theorem 3.8. This theorem will be proven in section A.7 of the appendix. In cases (i) and (ii), bd(A) will cause specific singularities in specific locations on L(ϕ0 , p0 ). The two more singular points in part (ii) will be specified in equation (A.25). 6. Artifact Reduction In this section, we describe a simple method to suppress the added streak artifacts described in Theorems 3.5, 3.7 and 3.8. More sophisticated methods are discussed in [4, 5] for the synchrotron problem that is described in Section 7. As outlined in Section 3, the application of FBP to incomplete data implicitly extends the data from A ⊂ [0, 2π] × R to all of [0, 2π] × R by padding it with zeros off of A. Therefore, the incomplete data RA f is obtained from full data by hard truncation via RA f (ϕ, p) = 1A (ϕ, p)Rf (ϕ, p). As we showed in section 3, the lines L(ϕ0 , p0 ) corresponding to (ϕ0 , p0 ) ∈ bd(A) can cause point and streak artifacts in the reconstruction. The reason is that multiplying by 1A creates jump discontinuities in the data. These are stronger singularities than those of Rf since Rf ∈ H1/2 ([0, 2π] × R) since f ∈ L2 (D) = H0 (D). One natural way to get rid of the jump discontinuities of 1A is to replace 1A 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 (6.1)
Rψ f (ϕ, p) = ψ(ϕ, p)Rf (ϕ, p)
and the reconstruction operator (6.2)
Lψ f = R∗ (ΛRψ f ) = R∗ (ΛψRf ) .
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 Lψ does not add artifacts. Theorem 6.1 (Artifact Reduction Theorem). Let f ∈ L2 (D). Then
(6.3)
WF(Lψ f ) ⊂ WF(f ). 19
(A) Left: Sinogram. Right: Reconstruction without smoothing.
(B) Left: Smoothed sinogram. Right: Smoothed reconstruction with suppressed artifacts.
Figure 9. Sinograms and reconstructions without smoothing in Figure 9(A) and with smoothing in Figure 9(B) (cf. Theorem 6.1). Therefore, Lψ does not add artifacts to the reconstruction. Let x ∈ D, ϕ ∈ [0, 2π], and ω 6= 0 and assume ψ(ϕ, x · θ(ϕ)) 6= 0. Then, (6.4)
(x, ωθ(ϕ)dx) ∈ WF(Lψ f ) if and only if (x, ωθ(ϕ)dx) ∈ WF(f ).
Many practitioners include a smooth cutoff function in incomplete data algorithms, but others do not. This theorem shows the advantages of including such a cutoff, and it has been suggested in several settings, including X-ray CT [10, 12, 19] and more general tomography problems [11, 39]. Figure 9 illustrates the effects of our smoothing algorithm. Proof of Theorem 6.1. This follows from the proof of Theorems 6.1 and 6.2 in [12]. Those results are for the generalized hyperplane transform and our theorem is a special case: the line transform with standard weights. In [12], the cutoff function is only a function of the angular variable. However, the symbol calculation is the same as the proof of [12, (6.3)] except that the cutoff function ϕ(ω) in that proof is replaced by the cutoff ψ(ϕ, p). (see [33] for the general case). From that calculation, the symbol of Lψ as a pseudodifferential operator is σ (Lψ ) (x, ωθ(ϕ)) = 4πψ ϕ, x · θ(ϕ) . 20
To prove this, one applies the symbol calculation in [12, (6.3) and (A.11)] with ξ = ωθ(ϕ), with weights µ = 1 and ν = 1, P = kξk, and one uses that ψ is symmetric. Using this expression [12, (6.3)] for the symbol σ (Lψ ), one sees that σ (Lψ ) is elliptic whenever ψ ϕ, x · θ(ϕ) 6= 0 and this proves (6.4). We illustrate the efficiency of this simple artifact reduction method on real synchrotron data at the end of Section 7. In practice, the smooth cut-off function ψ can be constructed in different ways depending on the set A. If A is a product of sets in ϕ and in p, then ψ can be a product of a cutoff in ϕ and one in p. This is how we made the cutoff for Figure 9(B) using the function [10, equation (39)]. In general, this is not possible and we now describe one general construction. Let ε > 0 and let Aε be the set of points in A at least ε from the boundary bd(A). If ε is small enough, Aε contains most of A. Let h be a smooth function supported in the ε/2 ball about (0, 0) ∈ [0, 2π] × R and with integral 1 (e.g., if h is an approximate identity). Then the convolution ψ = 1Aε ∗ h will be a smooth function supported in A which is equal to one on A2ε . This construction can be done numerically by numerical convolution. 7. Application: A Synchrotron Experiment Figure 10 shows tomographic data of a chalk sample (Figure 10(A) and 10(B)) that was acquired at a synchrotron experiment [4, 5] and a zoom of the corresponding reconstruction (Figure 10(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 theory at that time (such as in [10–12, 19, 26]).
(A) Segmented attenuation (B) Zoom of segmented sinogram. (C) Zoomed synchrotron reconsinogram (processed Radon struction given in Figure 1. transform data) with the rectangular zoom region in Fig. 10(B) highlighted.
Figure 10. Left and Center: The truncated attenuation sinogram (after processing to get Radon transform data) and an enlargement of a central section of bd(A). c Right: Reconstruction given in Figure 1 [4, IOP Publishing. Reproduced by permission of IOP Publishing. All rights reserved]. Taking a closer look at the attenuation sinogram, Figures 10(A)-10(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 11. The scanner is fixed and the sample holder, 21
Reduction Reductionofofvariable-truncation variable-truncationartifacts artifactsduring duringininsitu situX-ray X-raytomography tomography nonosignal signal
metal metalbar bar position positionofofsample sample
truncated truncated projection projection full full projections projections sample sample
metal metalbar bar (a) (a)
(b) (b)
c Figure 11. Data acquisition setup for the synchrotron experiment [4, IOP PubFigure 2: Left: Side-view of the percolation cell. Right: Top-view of the setup, where lishing. Reproduced by permission of IOP Publishing. All rights reserved]. Figure 2: Left: Side-view of the percolation cell. Right: Top-view of the setup, where the thespecimen specimenisisplaced placedininthe thecenter centerbetween betweenfour fourmetal metalbars bars(not (nottotoscale). scale).Four Fourangular angular positions positionsfor forthe thedetector detector(the (theblack blackbar) bar)isisshown: shown:For Fortwo twoofofthem, them,the theprojections projectionsare are including the metal bars, rotates. As the holder rotates, the bars block different amounts of X-rays complete. For the detector position along the SW-NE diagonal, the beam is occluded complete. For the detector position along the SW-NE diagonal, the beam is occluded from the source, as the sample holder rotates. Discretization of the angular measurements results in completely by bars is measured the In between one bythe themetal metalHowever, barsand andno nosignal signal measuredatat thedetector. detector. between one 10(B) the sharp completely vertical boundaries. the sharpishorizontal boundaries in pInseen in Figure position is shown at which the outermost parts of the beam are occluded by two of the position is following shown at which theTo outermost of the beam are occluded by two of theto Radon are required for the reason: get fromparts the measured sinogram (photon counts) metal bars, leading totoaaprojection with truncation from both sides. metal bars, leading projection with truncation from both sides. data, one takes negative logarithms of the raw photon counts. Because the photon count goes slowly
to zero near the horizontal boundaries, the negative log (Radon data) blows up. Therefore, the researchers sharply truncate the data in p while there is still signal, getting the sharp horizontal 2. Data boundaries in Figure 10(B) and the sharp corners. In practice, Acquisition the reconstructions from the measured incomplete data are computed by using the set-up standard FBP algorithm (which is the core of many commercial tomographs). As a result, the The motivating case for from the present is in situ X-ray imaging ofanalysis, obtained reconstructions suffer severestudy streak artifacts, seemicro-tomography Figure 12(A). Microlocal fluid flow through porous chalk in which the goal is to recover oil from the North in particular, Theorem 3.8, explains that the streaks in the reconstruction are caused by the myriad underground. In 10(A)-10(B)), situ X-ray tomography data was obtained for a More cylindrical corners inSea the data (cf. Figures and each corner produces a streak. importantly, the Artifact Reduction Theorem, 6.1, provides simple methodBL20XU to improve theSPring-8 classical builtporous chalk sample of diameter 0.6 mma using beamline of the in reconstruction algorithms to reduce To achieve that, Theorem 6.1 suggests to Synchrotron Radiationand Facility, Japanartifacts. using a monochromatic (28 keV) parallel-beam smoothly scan truncate the data along bd(A) a preprocessing step. implemented this in procedure configuration. Fluid is forcedinthrough the sample by aWe percolation cell, seen and applied it to synchrotron data. ofThe resulting reconstruction is shown in Figure 12(B). Figure 2a,the by applying a pressure 50 bars imitating the underground conditions. The Figure 3: a nice plot Here, one goal can is observe that our artifact reduction works very well in this situation. In the zoomedto model the structural changes of the sample during the fluid flow and a series in imagesofone canareobserve that all of theover streaks are gone while all of changes the details seem to be scans acquired continuously the experiment. Structural are slow preserved. 2. Data to the acquisition time of each complete scan and any sample deformations compared We would like to point out that the authors of [4, 5] also mitigate the streaks using the methodwithin each scan can be neglected. The percolation cell is equipped with four metal ology of Section 6 and by using more subtle methods. Acquisition set-up
bars which can sustain pressures of 200 bar and temperatures of 100 C. The metal bars havemotivating a radius ofcase 1 mm andpresent are positioned square and at approx. distance The for the study is in asitu X-rayaround micro-tomography imaging of 8. Discussion of 15.6 mm from the sample. The number of detector pixels is 2048 2048, providing fluid flow through porous chalk in which the goal is to recover oil from the North in make each horizontal slice a field of our view (FOV) for of data approx. 0 obtained 5∗ (Λ mm(1in Rf diameter [24]. Sea underground. In situ X-ray tomography for))aand cylindrical We first some observations on results LA was = R thenThe we discuss A detector is positioned outside the percolation cell and 1800 projections are collected porous chalk sample of diameter 0.6 mm using beamline BL20XU of the SPring-8 generalizations. covering 0 toRadiation 180 degrees. As seen in Figure most projections arekeV) complete, some are Synchrotron Facility, Japan using 2b a monochromatic (28 parallel-beam ∗ fully occluded by the metal bars, while some projections are partially occluded. These scan configuration. Fluid is forced through the sample by a percolation cell, seen in \ 0, but 8.1. Observations. Theorem 3.8 is valid as long as WF(ϕ0 ,p0 ) (1A ) = T ([0, 2π] × R) partial projections are the focus of the present work. In addition, since the sample Figuretheorem 2a, by applying a pressure of 50 bars imitating underground The the analogous for Sobolev singularities, Theoremthe 5.2(ii), assumesconditions. that A has a iscorner at than the FOV, all projections truncated. to model the singularity structural changes sample the fluid flow and a series (ϕ0 , p0 ). goal Iflarger A ishas a weaker at (ϕare ), then anduring analogous theorem would hold but one 0 ,ofpslightly 0the
would need to factor in the Sobolev strength of the wavefront of 1A above (ϕ0 , p0 ). Our artifact reduction method, which is motivated by Theorem 6.1, works well for the synchrotron data as was shown in Figure 12. The article [4] provides more elaborate artifacts 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 22
(A) Standard FBP reconstruction
(B) FBP reconstruction with artifact reduction (cf. Theorem 6.1).
Figure 12. Reconstructions from synchrotron data without smoothing (top) and c with smoothing (bottom) [4, IOP Publishing. Reproduced by permission of IOP Publishing. All rights reserved]. 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, 21, 41]). In [29] and [6], 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 [31, 37] 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 23
may also introduce artifacts, such as undersampling streaks. Discretization in our synchrotron experiment could be a factor in the streaks in Figure 10. 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 LA = R∗ (Λ (1A R)), but the results hold for more general filtering operators besides Λ. One key to our proofs is that Λ satisfies the ellipticity condition A.7, but many other operators satisfy this condition. For example, the operator, L = ∂2 − ∂p 2 , in Lambda CT [8] satisfies this condition, and the only difference come in our Sobolev continuity statements in Section 5 since L is order two and Λ is order one. The operator R∗ LR is of order 1 and so the statements in Theorem 5.2 would involve wavefront sets that are one degree lower than in Theorem 5.2. Theorems 3.5, 3.7 and 3.8 hold verbatim for generalized Radon transforms on lines in R2 because they all have the same canonical relation, given by (A.4), and the proofs would be done as those in this article with the symbol calculations in [12]. 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 [12] for the hyperplane transform or those in [11] for the spherical transform of PAT. These generalizations are the subject of ongoing work. In incomplete data problems for R, the artifacts are either along curves ϕ 7→ xb (ϕ) 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, in these cases, artifacts can spread on proper subsets of the surface over which data are taken, not necessarily the entire set. One such example is discussed in [11, Remark 4.7]. Analogous theorems should hold for cone beam CT, but the problem is more subtle because the reconstruction operator itself can add artifacts, even with complete data [9, 13]. Acknowledgements We thank the Japan Synchrotron Radiation Research Institute for the allotment of beam time on beamline BL20XU of SPring-8 (Proposal 2015A1147) that provided the raw data described in section 7. Leise Borg received partial support from P 3 , a joint project between Maersk Oil and Gas, Innovation Fund Denmark, and University of Copenhagen. The work of J¨ urgen Frikel was partially supported by the HC Ørsted Postdoc programme, co-funded by Marie Curie Actions at the Technical University of Denmark. Jakob Sauer Jørgensen was supported in part by the project “High-Definition Tomography” which was funded by Advanced Grant No. 291405 from the European Research Council, and in part by EPSRC grant EP/P02226X/1. Todd Quinto thanks John Schotland and Guillaume Bal for a stimulating discussion about multiplying distributions that relate to Remark A.5 and the reasons we consider only functions in this article. He thanks the Technical University of Denmark and DTU Compute for a wonderful semester during which this research was being done and he is indebted to his colleagues there for stimulating, enjoyable discussions that influenced this work. He received partial support from NSF grants DMS 1311558 and DMS 1712207, as well as a fellowship from the Otto Mønsteds Fond, and support from the Tufts University Faculty Research Awards Committee. Appendix A Here, we provide some basic microlocal analysis and then use this to prove our theorems. 24
A.1. Building blocks. Our first lemma gives some basic facts about wavefront sets. Lemma A.1. Let x0 ∈ R2 Let u and v be locally integrable functions. (i) Let U be an open neighborhood of x0 . Assume that u and v are equal on U , then WFx0 (u) = WFx0 (v). (ii) Let ψ be a smooth, compactly supported function, then WFx0 (ψu) ⊂ WFx0 (u). If ψ is nonzero at x0 then WFx0 (u) = WFx0 (ψu). (iii) WFx0 (u) = ∅ if and only if there is an open neighborhood U of x0 on which u is a smooth function. The analogous statements hold for functions on [0, 2π] × R. Proof. To prove part (i), we use the fact that if the Fourier transform F(u) is rapidly decreasing in a conic open neighborhood, V , of ξ0 , then, for any smooth compactly supported function ψ, F(ψu) is rapidly decreasing in some (perhaps smaller) conic neighborhood V ′ of ξ0 [18, Lemma 8.1.1]. Let u and v be equal on the neighborhood U of x0 and let u be smooth in direction ξ0 at x0 . Therefore, there is a cutoff function ψ1 at x0 and an open conic neighborhood V1 of ξ0 such that F(ψ1 u) is rapidly decreasing in V1 . Let ψ2 be a cutoff function at x0 that is supported in U , then F(ψ1 ψ2 u) is rapidly decreasing in some (perhaps smaller) conic neighborhood V1′ of ξ0 . Since, ψ1 ψ2 u = ψ1 ψ2 v, we see that F ((ψ1 ψ2 ) v) is rapidly decreasing in V1′ and (x0 , ξdx) ∈ / WF(v). The other containment follows by symmetry – switching u and v in this proof. To prove part (ii) we assume ψ is a smooth function of compact support. By Lemma 8.1.1 and Definition 8.1.2 of [18], (A.1)
WFx0 (ψu) ⊂ WFx0 (u).
Now, assume ψ(ϕ0 , p0 ) 6= 0 and choose a neighborhood U of x0 such that |ψ| is bounded away from zero in U . Let ̺ be a smooth, compactly supported function that is equal to 1/ψ in U , then, on U , u = ̺ψu and so WFx0 (u) = WFx0 (̺ψu) by part (i) of this lemma and then WFx0 (̺ψu) ⊂ WFx0 (ψu) by containment (A.1). Part (iii) is true as the singular support is the projection of the wavefront set on the first coordinate [18, Proposition 8.1.3], so the wavefront set is empty above x0 if and only if x0 ∈ / sing supp(u) if and only if u is smooth in a neighborhood of x0 . Our next definition will be useful to describe how wavefront set propagates under R and R∗ . Definition A.2. If C ⊂ T ∗ ([0, 2π] × R) × T ∗ (R2 ) then,
C t = {(x, ξ, ϕ, p, η) : (ϕ, p, η, x, ξ) ∈ C} .
If B ⊂ T ∗ ([0, 2π] × R), then C t ◦ B = (x, ξ) ∈ T ∗ (R2 ) : (x, ξ, ϕ, p, η) ∈ C t for some (ϕ, p, η) ∈ B
Our next theorem is the main technical theorem of the article. Its proof uses the theory of Fourier integral operators (FIO). If F : E ′ (X) → D ′ (Y ) is a Fourier integral operator between manifolds X and Y with associated to the canonical relation C ⊂ T ∗ (Y ) × T ∗ (X), then for f ∈ E ′ (X), F maps singularities in a precise way given by the H¨ormander-Sato Lemma [18, Theorem 8.2.13] plus the fact that Fourier integral operators satisfy the hypotheses of that theorem: WF(F f ) ⊂ C ◦ WF(f ). If F is elliptic and satisfies the microlocal Bolker Assumption: the natural projection (A.2) then equality holds: (A.3)
πL : C → T ∗ (Y ) is an injective immersion, WF(F f ) = C ◦ WF(f ).
More information about FIO can be found in [17, 18, 32, 40]. 25
Guillemin [14,15] proved, under some geometric assumptions that generalized Radon transforms, R, are elliptic FIO. If, in addition, they satisfy the full Bolker Assumption (this microlocal assumption plus a proper support assumption), then R∗ R is an elliptic pseudodifferential operator and (A.3) holds. It is well-known [15, 20, 33, 36] that R satisfies these conditions, and we will show how R∗ essentially satisfies them, too. Proposition A.3. The X-ray transform, R, is an elliptic Fourier integral operator with canonical relation C = ϕ, x · θ(ϕ), ω(−x · θ(ϕ)dϕ + dp), x, ωθ(ϕ) (A.4) : (ϕ, p) ∈ [0, 2π] × R, x ∈ R2 , ω 6= 0
and WF(Rf ) = C ◦ WF(f ) for f ∈ L2 (D). The specific correspondence is as follows. Let x0 ∈ R2 and let ξ0 ∈ R2 \ 0. Now, let ϕ0 ∈ [0, 2π] and ω ∈ R \ 0 such that ξ0 = ωθ(ϕ0 ). Then
(A.5)
(x0 , ξ0 dx) ∈ WF(f ) if and only if
ϕ0 , x0 · θ(ϕ0 ), ω(−x0 · θ ⊥ (ϕ0 )dϕ + dp) ∈ WF(Rf ).
Let C t be the transpose of C given by (A.4). Let (ϕ, p) ∈ [0, 2π] × R, ω 6= 0, and α ∈ R, and let (A.6)
x0 = x0 (ϕ, p, α) = αθ ⊥ (ϕ) + pθ(ϕ).
Then (A.7) Furthermore, (A.8)
C t ◦ {(ϕ, p, ω(−αdϕ + dp))}
= C t ◦ {(ϕ + π, −p, (−ω)(−αdϕ + dp))} = (x0 , ωθ(ϕ)) . C t ◦ {(ϕ, p, αdϕ)} = ∅.
The dual transform, R∗ is an elliptic Fourier integral operator associated with C t . Let g be a locally integrable function on [0, 2π] × R and assume g satisfies the symmetry condition
(A.9)
Then, (A.10) and (A.11) Therefore, (A.12)
∀(ϕ, p) ∈ [0, 2π] × R, g(ϕ, p) = g(ϕ + π, −p). WF(R∗ g) = C t ◦ WF(g)
WFL(ϕ0 ,p0) (R∗ g) = C t ◦ WF(ϕ0 ,p0 ) (g) for (ϕ0 , p0 ) ∈ [0, 2π] × R. (ϕ0 , p0 , ω0 (−αdϕ + dp)) ∈ WF(g) if and only if x0 (ϕ0 , p0 , α), ω0 θ(ϕ0 ) ∈ WF(R∗ g).
Note that, because g satisfies the symmetry condition (A.9), then (A.13)
(ϕ0 , p0 , ω0 (−αdϕ + dp)) ∈ WF(g) if and only if (ϕ0 + π, −p0 , −ω0 (αdϕ + dp)) ∈ WF(g). 26
Also note that, if f ∈ L2 (D), then Rf and Λ1A Rf both satisfy the symmetry condition (A.9) and are at least locally integrable. Proof. The statements about R are directly from Theorem A.2 in [36]. Note that θ (respectively θ ⊥ ) in [36] is the same as θ(ϕ) (respectively: θ ⊥ (ϕ)) in this article. Also, note that the crucial points are that R is an elliptic Fourier integral operator and that R satisfies the full Bolker assumption so (A.3) holds for R. Then, (A.5) is proven using the definition of composition and (A.4). The proofs for R∗ follow from arguments in the proof of Theorem A.2 and Theorem 2.8 in [36] and they will be outlined. First, (A.14) C t = x0 (ϕ, p, α), ωθ(ϕ), ϕ, p, ω(−αdϕ + dp) : (ϕ, p) ∈ [0, 2π] × R, α ∈ R, ω 6= 0
by the Definition A.2 and the expression (A.4) for C. The calculations (A.7) and (A.8) follow from the definition of composition and the expression (A.14) for C t . From [36, Theorem A.2], we know that R∗ is an elliptic Fourier integral operator associated to C t . The projection πL : C t → T ∗ (R2 ) is two to one, and wavefront of g at covector (ϕ0 , p0 , ω0 (−αdϕ + dp)) can, in general, be cancelled in R∗ g by wavefront of g at (ϕ0 +π, −p0 , −ω0 (αdϕ+dp)). However, because the functions g we consider are symmetric as in (A.9) and R∗ is also symmetric, we can view R∗ as a FIO with domain E ′ (([0, 2π] × R)/R) where R is the relation (ϕ, p) ≡ (ϕ + π, −p). With this identification, then πL is an injective immersion and so R∗ satisfies the microlocal Bolker assumption (e.g., [36]) and one can use arguments in the proof in [33], for example, to infer (A.10). Then, a calculation using the expression (A.7) proves (A.11) and (A.12). A.2. Proof of Proposition 2.4. Recall that (2.9) is WFL(ϕ0 ,p0) (f ) 6= ∅ if and only if WF(ϕ0 ,p0 ) (Rf ) 6= ∅.
To prove this, we note that covectors in the left side of (A.5) are in WFL(ϕ0 ,p0 ) (f ) and covectors in the right side of (A.5) are in WF(ϕ0 ,p0 ) (Rf ). Therefore, the equivalence between covectors in (A.5) implies the equivalence (2.9), and this finishes the proof. A.3. Proof of Theorem 3.1. Let A ∈ [0, 2π] × R be closed and satisfy symmetry condition (3.1) and f ∈ L2 (D). Let (x0 , ξ0 dx) ∈ T ∗ (R2 ) \ 0 and let (ϕ0 , p0 ) ∈ [0, 2π] × R and ω0 6= 0 such that p0 = x0 · θ(ϕ0 ) and ξ0 = ω0 θ(ϕ0 ). Finally, let α = x0 · θ ⊥ (ϕ). To prove part (i), assume (x0 , ξ0 dx) ∈ WF(LA f ). Using the definition of composition and the expression (A.7) one sees that there are exactly two solutions λ ∈ T ∗ ([0, 2π] × R) to C t ◦ {λ} = {(x0 , ξ0 dx)} and they are the covectors (A.15)
λ1 = (ϕ0 , p0 , ω0 (−αdϕ + dp)) R∗ g
λ2 = (ϕ0 + π, −p0 , −ω0 (αdϕ + dp))
Let g = Λ1A Rf . Since LA f = we can use Theorem A.3 (A.10) and then (A.13) to conclude that both λ1 and λ2 are in WF(g). Because Λ is a pseudodifferential operator and g = Λ(1A Rf ), λ1 and λ2 are both in WF(1A Rf ). Since sing supp(1A Rf ) ⊂ sing supp(1A ) ∪ sing supp(Rf ), either and either
WF(ϕ0 ,p0 ) (1A ) 6= ∅ or WF(ϕ0 ,p0) (Rf ) 6= ∅,
WF(ϕ0 +π,−p0 ) (1A ) 6= ∅ or WF(ϕ0 +π,−p0 ) (Rf ) 6= ∅ as well. This proves part (i). ∗ Part (ii) follows from (A.10) and the fact that if λ ∈ T(ϕ ([0, 2π] × R) is a singularity above 0 ,p0 ) t ∗ (ϕ0 , p0 ), then C ◦ {λ} ∈ N (L(ϕ0 , p0 )) \ 0. Remark A.4. Here we provide an example in which there might not be artifacts along L(ϕ0 , p0 ) even though f has singularities conormal to the line. The basic idea is that those singularities come from points outside the scanned region (i.e., union of lines in the data set). Part (i) of Theorem 3.5 asserts that, if f has a singularity conormal to L(ϕ0 , p0 ), then there could be artifacts along 27
L(ϕ0 , p0 ) (WFL(ϕ0 ,p0 ) (LA f ) ⊂ N ∗ (L(ϕ0 , p0 )) \ 0). However, we cannot in general assert more since artifacts might not be added. To show this, we construct a function f with an artifact conormal to a line on the boundary of the data set but for which the data, RA f , is identically zero. Therefore, no artifacts are created. Let A = {(ϕ, p) : ϕ ∈ [0, π/2], p ≤ 0} Then the set of lines in A do not meet the open first quadrant of the coordinate plane. If f is the characteristic function of a disk of radius 1 at the point (1, 1) then ((0, 1), dx1 )) ∈ WFL(0,0) (f ) but LA f has no artifacts since LA f = 0 since the data RA f is identically zero as f is supported in the first quadrant. To summarize, the singularity of f that could cause an added artifact is not visible in the data set since f is zero on all lines in the scanned region (i.e., all (ϕ, p) in the data set). A.4. Proof of Theorem 3.3. Assume (ϕ0 , p0 ) ∈ / bd(A). Then, there are two cases: either (ϕ0 , p0 ) ∈ int(A), the interior of A or (ϕ0 , p0 ) ∈ ext(A), the exterior of A. If (ϕ0 , p0 ) ∈ int(A) then, by definition of interior, there is a neighborhood U of (ϕ0 , p0 ) that is contained in A. By Lemma A.1 part (i), since Rf and 1A Rf agree in U , WF(ϕ0 ,p0 ) (Rf ) = WF(ϕ0 ,p0 ) (1A Rf ) and since Λ is a pseudodifferential operator, WF(ϕ0 ,p0 ) (ΛRf ) = WF(ϕ0 ,p0 ) (Λ1A Rf ). Using (in order from left to right) the filtered backprojection inversion formula for R, f = and (A.11), the equality just above, and (A.11) again,
1 ∗ 4π R ΛRf ,
WFL(ϕ0 ,p0 ) (f ) = C t ◦ WF(ϕ0 ,p0 ) (ΛRf ) = C t ◦ WF(ϕ0 ,p0) (Λ1A Rf ) = WFL(ϕ0 ,p0 ) (LA f ).
This proves part (i). Now, assume (ϕ0 , p0 ) ∈ ext(A), then 1A Rf is zero in a neighborhood U of (ϕ0 , p0 ). Because Λ is a pseudodifferential operator, Λ1A Rf is smooth in U so WF(ϕ0 ,p0 ) (Λ1A Rf ) = ∅. Now, we use statement (A.11) to show WFL(ϕ0 ,p0) (LA f ) = ∅. This proves part (ii). 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). Remark A.5. In order to multiply distributions u and v as distributions one needs the noncancellation condition (A.16)
∀ (ϕ, p, ξ) ∈ WF(u), (ϕ, p, −ξ) ∈ / WF(v)
to hold, then uv is a distribution and an upper bound for WF(uv) can be calculated in terms of WF(u) and WF(v) ( [18, Theorem 8.2.10]). However, this non-cancellation condition does not hold, in general for 1A and Rf when 1A either is smooth with finite slope or is not smooth at (ϕ0 , p0 ). That is why we consider functions f ∈ L2 (D) since 1A Rf will be a function in L2 ([0, 2π] × R) even if [18, Theorem 8.2.10] does not apply. Furthermore, multiplying 1A and Rf will not create singularities at places besides where these two functions have singularities because sing supp(1A Rf ) ⊂ sing supp(1A ) ∪ sing supp(Rf ). A.5. Proof Theorem A.6, The Complete Artifact Characterization. We state our artifact characterization as one theorem to emphasize that it addresses all cases. We divided it up in section 3 so we could focus on each individual type of artifact. Theorem A.6 (Added Artifacts in the Reconstruction). Let f ∈ L2 (D) and let A be a closed subset of [0, 2π] × R that satisfies the symmetry condition (3.1). Let (ϕ0 , p0 ) ∈ bd(A). Streak artifacts or point artifacts can occur in the reconstruction LA f conormal the line L(ϕ0 , p0 ) under the following conditions. I. Assume bd(A) is smooth at (ϕ0 , p0 ). 28
A. Assume bd(A) has finite slope at (ϕ0 , p0 ) and let the smooth function p = p(ϕ) parameterize this part of bd(A) for ϕ in some open neighborhood I of ϕ0 . For ϕ ∈ I define (A.17)
xb = xb (ϕ) = p(ϕ)θ(ϕ) + p′ (ϕ)θ ⊥ (ϕ). (i) If f is not smooth conormal to L(ϕ0 , p0 ), then WFL(ϕ0 ,p0 ) (LA f ) ⊂ N ∗ (L(ϕ0 , p0 )) \ 0,
(A.18)
so there can be a streak artifact along L(ϕ0 , p0 ). (ii) If f is smooth conormal to L(ϕ0 , p0 ), then LA f is smooth conormal to L(ϕ0 , p0 ) except perhaps above xb : WFL(ϕ0 ,p0 ) (LA f ) ⊂ (xb (ϕ0 ), ωθ(ϕ0 )) : ω 6= 0 .
If, in addition, Rf (ϕ0 , p0 ) 6= 0, then equality holds in (A.18). B. Now, assume bd(A) has a vertical tangent at (ϕ0 , p0 ) (infinite slope). (i) If f is not smooth conormal to L(ϕ0 , p0 ), then WFL(ϕ0 ,p0 ) (LA f ) ⊂ N ∗ (L(ϕ0 , p0 )) \ 0.
This means that there can be a streak artifact all along L(ϕ0 , p0 ). (ii) If f is smooth conormal to L(ϕ0 , p0 ) then LA f is smooth conormal to L(ϕ0 , p0 ): WFL(ϕ0 ,p0 ) (f ) = ∅ implies that WFL(ϕ0 ,p0) (LA f ) = ∅. There are no streak artifacts or point artifacts conormal to L(ϕ0 , p0 ). II. Assume that bd(A) is not smooth at (ϕ0 , p0 ). Then, (A.19)
WFL(ϕ0 ,p0 ) (LA f ) ⊂ N ∗ (L(ϕ0 , p0 )) \ 0. Thus, LA f can have a streak artifact on L(ϕ0 , p0 ). If f is smooth conormal to L(ϕ0 , p0 ), Rf (ϕ0 , p0 ) 6= 0, and bd(A) has a corner at (ϕ0 , p0 ), then LA f does have a streak artifact all along L(ϕ0 , p0 ) – equality holds in (A.19).
Our next remark will be used in ellipticity proofs that follow. Remark A.7. 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 ∈ L2 (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 ∈ L2 (D), WF(ΛRf ) = WF(Rf ). By (A.8), even if (ϕ, ωdϕ) ∈ WF(1A Rf ), that covector will not affect the calculation of C t ◦ WF(Λ1A 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. A.6. Proof of Theorem A.6. Let f ∈ L2 (D), let A be a symmetric closed subset of [0, 2π] × R, and let g = Λ1A Rf . First note that for (ϕ0 , p0 ) ∈ [0, 2π] × R (A.20)
WFL(ϕ0 ,p0 ) (LA f ) =C t ◦ WF(ϕ0 ,p0 ) (g) ∗ ⊂ C t ◦ T(ϕ ([0, 2π] × R) \ 0 = N ∗ (L(ϕ0 , p0 )) \ 0 ,p ) 0 0
by (A.11), (A.7), and (A.8). To prove part I., we assume bd(A) is smooth at (ϕ0 , p0 ). This part of the proof uses ideas from the paradigm in [11, Sect. 3] except that, in this proof, the non-cancellation condition (A.16) does not always hold, so the proof is more intricate. We prove part I.A. and assume bd(A) has finite slope α at (ϕ0 , p0 ). Recall that the “bad point” (A.17) is given by xb = p0 θ(ϕ0 ) + αθ ⊥ (ϕ0 ). 29
Since (1, α) is tangent to bd(A) at (ϕ0 , p0 ), (−α, 1) is normal to this curve so (A.21)
WF(ϕ0 ,p0) (1A ) = {ϕ0 , p0 , ω(−αdϕ + dp)) : ω 6= 0} .
Here we use that WF(1A ) is conormal to bd(A) at points where bd(A) is smooth. Consider case I.A.(i) and assume f is not smooth conormal to L(ϕ0 , p0 ). Then, by (A.20), WFL(ϕ0 ,p0 ) (LA f ) ⊂ N ∗ (L(ϕ0 , p0 )) \ 0. Now, consider part I.A.(ii). Since we are assuming f is smooth conormal to L(ϕ0 , p0 ), WF(ϕ0 ,p0 ) (Rf ) = ∅ by Proposition 2.4. This means that Rf is a smooth function in some neighborhood U of (ϕ0 , p0 ) by Lemma A.1 part (iii). Therefore, WF(ϕ0 ,p0 ) (1A Rf ) ⊂ WF(1A ) by this same lemma, parts (i) and (ii). Now using Theorem A.3 and that Λ is elliptic in the sense of Remark A.7 WFL(ϕ0 ,p0 ) (LA f ) = C t ◦ WF(ϕ0 ,p0 ) (Λ1A Rf ) = C t ◦ WF(ϕ0 ,p0 ) (1A Rf )
⊂ C t ◦ WF(ϕ0 ,p0 ) (1A )
(A.22)
= C t ◦ {ϕ0 , p0 , ω(−αdϕ + dp)) : ω 6= 0} = xb (ϕ0 ), ωθ(ϕ0 ) : ω 6= 0 .
Next, assume Rf (ϕ0 , p0 ) 6= 0, then, because 1A Rf is the product of 1A with a smooth nonzero function (at least on a neighborhood of (ϕ0 , p0 )), WF(ϕ0 ,p0) (1A Rf ) = WF(ϕ0 ,p0) (1A ) by Lemma A.1 part (ii). The one containment in the calculation (A.22) just above becomes an equality and this proves WFL(ϕ0 ,p0 ) (LA f ) = xb (ϕ0 ), ωθ(ϕ0 ) : ω 6= 0 . This finishes the proof of part I.A.(ii). The proof of part I.B. is similar, but the difference is that the normal to bd(A) at (ϕ0 , p0 ) is horizontal because of the vertical tangent and so WF(ϕ0 ,p0 ) (1A ) = {(ϕ0 , p0 , ωdϕ) : ω 6= 0}. Then, the proof proceeds using the paradigm in [11] since the non-cancellation condition (A.16) below holds between 1A and Rf and so 1A Rf is defined as a distribution and an upper bound for WF(1A Rf ) can be calculated using [18, Theorem 8.2.10]. Basically, the proof follows in the same way as the proof of part I.A. but without the point xb . This finishes the proof of part I. To prove part II., assume bd(A) is not smooth at (ϕ0 , p0 ). Using (A.20), we see WFL(ϕ0 ,p0 ) (LA f ) ⊂ N ∗ (L(ϕ0 , p0 )) \ 0. Now, assume f is smooth conormal to L(ϕ0 , p0 ) and Rf (ϕ0 , p0 ) 6= 0. Then, by Proposition 2.4, WF(ϕ0 ,p0 ) (Rf ) = ∅ so Rf is smooth near (ϕ0 , p0 ) by Lemma A.1 part (iii). Since Rf is nonzero at (ϕ0 , p0 ), WF(ϕ0 ,p0 ) (1A ) = WF(ϕ0 ,p0) (1A RF ) by part (i) (to replace Rf by a smooth function of compact support equal to Rf on a neighborhood of (ϕ0 , p0 )) and then part (ii) of that lemma. Because bd(A) has a corner at (ϕ0 , p0 ), then ∗ WF(ϕ0 ,p0 ) (1A Rf ) = WF(ϕ0 ,p0) (1A ) = T(ϕ ([0, 2π] × R) \ 0. 0 ,p0 ) ∗ Because Λ is elliptic in the sense of Remark A.7, WF(ϕ0 ,p0 ) (g) = T(ϕ ([0, 2π] × R) \ 0 (up 0 ,p0 ) to covectors of the form (ϕ0 , p0 , αdϕ) but such vectors will not affect the composition with C t ). Therefore, by (A.20), we see ∗ WFL(ϕ0 ,p0 ) (LA f ) = WFL(ϕ0 ,p0 ) (R∗ (g)) = C t ◦ T(ϕ ([0, 2π] × R) \ 0 = N ∗ (L(ϕ0 , p0 )) \ 0. 0 ,p0 )
This finishes the proof of part II. 30
A.7. Proof of Theorem 5.2. We first prove a lemma giving the correspondence between Sobolev wavefront set and R∗ . Theorem A.8 (Sobolev Wavefront Correspondence for R∗ ). Let (ϕ0 , p0 ) ∈ [0, 2π] × R and let ω0 6= 0 and α ∈ R. Let η0 = ω0 (−αdθ + dp),
x0 = p0 θ(ϕ0 ) + αθ ⊥ (ϕ0 ), and ξ0 = ω0 θ(ϕ0 )dx.
If g is in Hs at (ϕ0 , p0 , η0 ) then R∗ g is in Hs+1/2 at (x0 , ξ0 ). On the other hand, if (ϕ0 , p0 , η0 ) ∈ WFs (g) then (x0 , ξ0 ) ∈ W Fs+1/2 (R∗ g). Proof. The easier part of the proof follows the arguments for the analogous theorem for R, [34, Theorem 3.1]. Assume g is in Hs at (ϕ0 , p0 , η0 ). Then, by [32, Theorem 6.1, p. 259], g can be written as the sum of two distributions, g = g1 + g2 where g1 ∈ Hs ([0, 2π] × R) and g2 ∈ D ′ ([0, 2π] × R) with (ϕ0 , p0 , η0 ) ∈ / WF(g2 ). Because R∗ is a FIO of order −1/2 associated to a local canonical graph (and ∗ R is defined on all distributions), R∗ is continuous from Hs,loc([0, 2π] × R) to Hs+1/2,loc (R2 ) by [17, Theorem 4.3.1] or [7, Corollary 4.4.5]. Therefore R∗ g1 is in H(s+1/2),loc (R2 ). Since (ϕ0 , p0 , η0 ) ∈ / ∗ ∗ ∗ WF(g2 ), then (x0 , ξ0 ) ∈ / WF(R g2 ) by Proposition A.3 equation (A.12). Therefore, R g1 + R g2 = R∗ g is in Hs+1/2 at (x0 , ξ0 ). To prove the second statement in the lemma, assume R∗ g is in Hs+1/2 at (x0 , ξ0 ). We will prove g is in Hs at (ϕ0 , p0 , η0 ). We start by writing R∗ g = h1 +h2 where h1 ∈ Hs+1/2 and (x0 , ξ0 ) ∈ / WF(h2 ). Let E be an open disk centered at the origin containing x0 , and let ψ ∈ D(R2 ) be equal to 1 on E. Then, ψh1 and ψh2 are distributions of compact support so R(ψhj ) is a distribution of compact support and gej = ΛR(ψhj ), for j = 1, 2 is a distribution. By the Radon inversion formula,
R∗ (e g1 + e g2 ) = ψ(h1 + h2 ) = ψR∗ (g).
Let E ′ be a disk centered at the origin containing supp(ψ). Because ψh1 is in Hs+1/2 at (x0 , ξ0 ), Rψh1 is in Hs+1 at (ϕ0 , p0 , η0 ) by the Sobolev wavefront correspondence for R, [34, Theorem 3.1]. Now, because Λ is a pseudodifferential operator of order one, ge1 is in Hs at (ϕ0 , p0 , η0 ).
Note that (x0 , ξ0 ) ∈ / WF(ψh2 ) as (x0 , ξ0 ) ∈ / WF(h2 ) (see Lemma A.1, part (ii)). By (A.12) in Proposition A.3 and the fact that Λ is a pseudodifferential operator, (ϕ0 , p0 , η0 ) ∈ / WF(e g2 ). Let e g = ge1 + ge2 . The following argument shows that (ϕ0 , p0 , η0 ) ∈ / WF(g − ge):
R∗ (g − ge) = R∗ g − R∗ (e g1 ) − R∗ (e g2 ) = R∗ g − ψ(h1 + h2 ) = (1 − ψ)(R∗ g) = 0 for x ∈ E
as ψ = 1 on E. Therefore, (x0 , ξ0 ) ∈ / WF(R∗ (g − (e g1 + ge2 ))). Again, by Proposition A.3,
(A.23)
(ϕ0 , p0 , η0 ) ∈ / WF(g − (e g1 + e g2 )).
Because ge1 is in Hs at (ϕ0 , p0 , η0 ), (ϕ0 , p0 , η0 ) ∈ / WF(e g2 ), and by (A.23), g is in Hs at (ϕ0 , p0 , η0 ). This finishes the proof of the theorem.
For g ∈ D ′ ([0, 2π] × R), we let (WFs )(ϕ0 ,p0) (g) denote the Sobolev wavefront set of g above (ϕ0 , p0 ) (in analogy to (2.8)). Now, we prove Theorem 5.2. Let f ∈ L2 (D) and let A be a closed symmetric subset of [0, 2π] × R. Let (ϕ0 , p0 ) ∈ bd(A) and assume Rf (ϕ0 , p0 ) 6= 0 and f is smooth conormal to L(ϕ0 , p0 ), 31
i.e., WFL(ϕ0 ,p0 ) (f ) = ∅. Because f is smooth conormal to L(ϕ0 , p0 ), WF(ϕ0 ,p0 ) (Rf ) = ∅. Since Rf (ϕ0 , p0 ) 6= 0, for each s, (A.24)
(WFs−1 )(ϕ0 ,p0 ) (Λ1A Rf ) = (WFs )(ϕ0 ,p0) (1A Rf ) = (WFs )(ϕ0 ,p0 ) (1A ) ;
the left-hand equality is true because Λ is an elliptic pseudodifferential operator of order one in these directions and, the right-hand equality is true by arguments similar to the proof of part (ii) of Lemma A.1 (or realizing that multiplication by Rf is, locally near (ϕ0 , p0 ) an elliptic order zero pseudodifferential operator). To prove part (i) of the theorem, assume bd(A) is smooth at (ϕ0 , p0 ). First, assume the boundary has finite slope at (ϕ0 , p0 ). Let α be the slope of the curve bd(A) at (ϕ0 , p0 ), and let η0 = −αdϕ + dp. We claim that (ϕ0 , p0 , ±η0 ) ∈ WF1/2 (1A ) and, for s < 1/2, 1A is in Hs at (ϕ0 , p0 , ωη0 ). Furthermore 1A is smooth in every other direction above (ϕ0 , p0 ). The proofs of these two statements are now outlined. Use a diffeomorphism to transform bd(A) locally to a horizontal line (α = 0). Use a product cutoff function ψ = ψ1 (ϕ)ψ2 (p) to calculate F(ψ1A ). Let (ν, τ ) be the Fourier dual variables to (ϕ, p). Using integrations by parts, one shows that this localized Fourier transform is of the form S(ν)T (τ ) where S is a smooth, rapidly decreasing function and T is O(1/ |τ |) is in the vertical direction. Therefore S(ν)T (τ ) is rapidly decaying in all directions but the vertical. Changing coordinates back, this implies that 1A is in Hs for s < 1/2 at (ϕ0 , p0 , ±η0 ) and (ϕ0 , p0 , ±η0 ) ∈ WF1/2 (1A ) (here we use that the wavefront set transforms as the cotangent space does under diffeomorphisms). This also shows that this localized Fourier transform is rapidly decaying in all directions except ±(−αdϕ + dp). Now, using (A.24) one sees that (ϕ0 , p0 , ±η0 ) ∈ WF−1/2 (Λ1A Rf ); Λ1A Rf is in Hs for s < −1/2 at (ϕ0 , p0 , ±η0 ); and (ϕ0 , p0 , η) ∈ / WF(Λ1A Rf ) for any η not parallel η0 . Now, by Theorem A.8, LA f = R∗ Λ1A Rf is in Hs at (xb , ±θ(ϕ0 )) for s < 0 and (xb , ±θ(ϕ0 )) ∈ W F0 (LA f ) where xb is given by (A.17). Using this theorem again, one sees that for any x ∈ L(ϕ0 , p0 ), if x 6= xb , (x, ±θ(ϕ0 )dx) ∈ / WF(LA f ).
Therefore, the only covectors in [N ∗ (L(ϕ0 , p0 )) \ 0] ∩ WF(LA f ) are (xb , αθ(ϕ0 )dx) for α 6= 0. If bd(A) is vertical at (ϕ0 , p0 ) then there is no point xb and the proof just given shows that LA f is smooth conormal to L(ϕ0 , p0 ). This finishes the proof of part (i). To prove part (ii), assume bd(A) has a corner at (ϕ0 , p0 ). Let α1 and α2 be the slopes at (ϕ0 , p0 ) of the two parts of bd(A). Let (A.25)
ηj = −αj dϕ + dp,
xbj = p0 θ(ϕ0 ) + αj θ ⊥ (ϕ0 ),
j = 1, 2.
an argument similar to the diffeomorphism/integration by parts argument in the last part of the proof is used. First a diffeomorphism is used to transform the corner so, locally A becomes A = {(ϕ, p) : ϕ ≥ 0, p ≥ 0}. Then one uses a product cutoff ψ = ψ1 (ϕ)ψ2 (p) at (0, 0). Then, the Fourier transform can be written F (ψ1A ) = S(ν)T (τ ) where S(ν) = O(1/ |ν|) and T (τ ) = O(1/ |τ |). So, the localized Fourier transform is decreasing of order −1 in the vertical and horizontal directions and −2 in all other directions. Transforming coordinates back, one gets that (ϕ0 , p0 , ±ηj ) ∈ WF−1/2 (Λ1A Rf ) and, for s < −1/2, Λ1A Rf is in Hs at (ϕ0 , p0 , ηj ). Other covectors are in WF1/2 (Λ1A Rf ). Part (ii) now follows from the Sobolev correspondence under R∗ in Theorem A.8. Part (iii) is valid by the relation between Sobolev wavefront set of f and Rf , [34, Theorem 3.1], since (ϕ0 , p0 ) ∈ int(A), 1A Rf = Rf in a neighborhood of (ϕ0 , p0 ), Λ is elliptic as in Remark A.7, and by the relation between Sobolev wavefront set of g and R∗ g, Theorem A.8. 32
References [1] L. L. Barannyk, J. Frikel, and L. V. Nguyen. On artifacts in limited data spherical radon transform: curved observation surface. Inverse Problems, 32:015012, 2016. [2] R. Bates and R. Lewitt. Image reconstruction from projections: I: General theoretical considerations, III: Projection completion methods (theory), IV: Projection completion methods (computational examples). Optik, 50:I: 19–33, II: 189–204, III: 269–278, 1978. [3] F. E. Boas and D. Fleischmann. CT artifacts: Causes and reduction techniques. Imaging Med., 4:229–240, 2012. [4] L. Borg, J. S. Jørgensen, J. Frikel, and J. Sporring. Reduction of variable-truncation artifacts from beam occlusion during in situ X-ray tomography. Meas. Sci. Tech., 28:19pp, 2017. [5] L. Borg, J. S. Jørgensen, and J. Sporring. Towards characterizing and reducing artifacts caused by varying projection truncation. Technical report, Department of Computer Science, University of Copenhagen, 2017/1. 42pp, http://static-curis.ku.dk/portal/files/174175124/teknisk rapport leise borg1.pdf. [6] J. K. Choi, H. S. Park, S. Wang, Y. Wang, and J. K. Seo. Inverse problem in quantitative susceptibility mapping. SIAM J. Imaging Sci., 7:16691689, 2014. [7] J. J. Duistermaat. Fourier integral operators, volume 130 of Progress in Mathematics. Birkh¨ auser, Inc., Boston, MA, 1996. [8] A. Faridani, E. L. Ritman, and K. T. Smith. Local tomography. SIAM J. Appl. Math., 52:459–484, 1992. [9] D. V. Finch, I.-R. Lan, and G. Uhlmann. Microlocal Analysis of the Restricted X-ray Transform with Sources on a Curve. In G. Uhlmann, editor, Inside Out, Inverse Problems and Applications, volume 47 of MSRI Publications, pages 193–218. Cambridge University Press, 2003. [10] J. Frikel and E. T. Quinto. Characterization and reduction of artifacts in limited angle tomography. Inverse Problems, 29:125007, 2013. [11] J. Frikel and E. T. Quinto. Artifacts in incomplete data tomography with applications to photoacoustic tomography and sonar. SIAM J. Appl. Math., 75:703–725, 2015. [12] J. Frikel and E. T. Quinto. Limited data problems for the generalized Radon transform in Rn . SIAM J. Math. Anal., 48:2301–2318, 2016. [13] A. Greenleaf and G. Uhlmann. Non-local inversion formulas for the X-ray transform. Duke Math. J., 58:205–240, 1989. [14] V. Guillemin. On some results of Gelfand in integral geometry. Proc. Symp. Pure Math., 43:149–155, 1985. [15] V. Guillemin and S. Sternberg. Geometric Asymptotics. American Mathematical Society, Providence, RI, 1977. [16] M. G. Hahn and E. T. Quinto. Distances between measures from 1-dimensional projections as implied by continuity of the inverse Radon transform. J. Wahrscheinlichkeitstheorie verw Gebiete, 70:361–380, 1985. [17] L. H¨ ormander. Fourier Integral Operators, I. Acta Math., 127:79–183, 1971. [18] L. H¨ ormander. The analysis of linear partial differential operators. I. Classics in Mathematics. Springer-Verlag, Berlin, 2003. Distribution theory and Fourier analysis, Reprint of the second (1990) edition [Springer, Berlin]. [19] A. I. Katsevich. Local tomography for the limited-angle problem. J. Math. Anal. Appl., 213:160–182, 1997. [20] V. P. Krishnan and E. T. Quinto. Microlocal Analysis in Tomography. In O. Scherzer, editor, Handbook of Mathematical Methods in Imaging, chapter 18, pages 847–902. Springer Verlag, New York, Berlin, 2015. [21] A. K. Louis. Ghosts in tomography, the null space of the Radon transform. Math. Methods Appl. Sci., 3:1–10, 1981. [22] A. K. Louis. Incomplete data problems in X-ray computerized tomography I. Singular value decomposition of the limited angle transform. Numer. Math., 48:251–262, 1986. [23] F. Natterer. Efficient implementation of ‘optimal’ algorithms in computerized tomography. Math. Methods Appl. Sci., 2:545–555, 1980. [24] F. Natterer. The mathematics of computerized tomography. Classics in Mathematics. Society for Industrial and Applied Mathematics, New York, 2001. [25] F. Natterer and F. W¨ ubbeling. Mathematical methods in image reconstruction. SIAM Monographs on Mathematical Modeling and Computation. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2001. [26] L. V. Nguyen. How strong are streak artifacts in limited angle computed tomography? Inverse Problems, 31:055003, 2015. [27] L. V. Nguyen. On artifacts in limited data spherical radon transform: flat observation surfaces. SIAM J. Math. Anal., 47:2984–3004, 2015. [28] L. V. Nguyen. On The Strength Of Streak Artifacts In Filtered Back-projection Reconstructions For Limited Angle Weighted X-Ray Transform. J. Fourier Anal. Appl., 23:712–728, 2017. [29] B. Palacios, G. Uhlmann, and Y. Wang. Reducing Streaking Artifacts in Quantitative Susceptibility Mapping. SIAM J. Imaging Sci., 10:1921–1934, 2017. 33
[30] X. Pan, E. Y. Sidky, and M. Vannier. Why do commercial CT scanners still employ traditional, filtered backprojection for image reconstruction? Inverse Problems, 25:123009, 2009. [31] H. S. Park, J. K. Choi, and J. K. Seo. Characterization of metal artifacts in x-ray computed tomography. Comm. Pure Appl. Math., 70:2191–2217, 2017. [32] B. Petersen. Introduction to the Fourier Transform and Pseudo-Differential Operators. Pittman, Boston, 1983. [33] E. T. Quinto. The dependence of the generalized Radon transform on defining measures. Trans. Amer. Math. Soc., 257:331–346, 1980. [34] E. T. Quinto. Singularities of the X-ray transform and limited data tomography in R2 and R3 . SIAM J. Math. Anal., 24:1215–1225, 1993. [35] E. T. Quinto. Exterior and Limited Angle Tomography in Non-destructive Evaluation. Inverse Problems, 14:339– 353, 1998. ´ [36] E. T. Quinto. An Introduction to X-ray tomography and Radon Transforms. In G. Olafsson and E. T. Quinto, editors, The Radon Transform and Applications to Inverse Problems, volume 63 of AMS Proceedings of Symposia in Applied Mathematics, pages 1–23, Providence, RI, USA, 2006. American Mathematical Society. [37] G. Rigaud. On Analytical Solutions to Beam-Hardening. Sens. Imaging, 18:17pp, 2017. [38] W. Rudin. Functional analysis. McGraw-Hill Book Co., New York, 1973. McGraw-Hill Series in Higher Mathematics. [39] L. Shen, E. T. Quinto, S. Wang, and M. Jiang. Simultaneous reconstruction and segmentation with the MumfordShah functional for electron tomography. In 38th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC), pages 5909 – 5912, 2016. [40] F. Tr`eves. Introduction to Pseudodifferential and Fourier Integral Operators. Volume 2: Fourier Integral Operators. Plenum Press, New York and London, 1980. [41] J. Vogelgesang and C. Schorr. Iterative Region-of-Interest Reconstruction from Limtied Data using Prior Information. Sens. Imaging, 18:21pp, 2017. [42] F. W. Warner. Foundations of differentiable manifolds and Lie groups, volume 94 of Graduate Texts in Mathematics. Springer-Verlag, New York-Berlin, 1983. Corrected reprint of the 1971 edition. [43] Y. Yang, S. S. Hakim, S. Bruns, K. N. Dalby, K. Uesugi, S. L. S. Stipp, and H. O. Sørensen. Wormholes grow along paths with minimal cumulative surface. Technical report, University of Copenhagen and SPring-8, 2016. 22pp, https://arxiv.org/abs/1704.01062.
34
