Artifacts in the inversion of the broken ray transform in the plane

We study the integral transform over a general family of broken rays in $\mathbb{R}^2$. It is natural for broken rays to have conjugate points, for example, when they are reflected from a curved boundary. If there are conjugate points, we show that the singularities cannot be recovered from local data and therefore artifacts arise in the reconstruction. We apply these conclusions to two examples, the V-line Radon transform and the parallel ray transform. In each example, a detailed discussion of the local and global recovery of singularities is given and we perform numerical experiments to illustrate the results.


INTRODUCTION
The purpose of this work is to study the integral transform over a general family of broken rays in the plane. A broken ray in the Euclidean space is usually defined as the linear path following the law of reflection on the boundary, which will be an important example of our setting in the following. In fact, a major motivation of this work is the reconstruction of an unknown function from integrals over such broken rays in medical imaging. One promising application is the Single Photon Emission Computed Tomography (SPECT) with Compton cameras.
We define a more general family of broken rays. Suppose f is a distribution with compact support. A broken ray is defined as the union of two rays l 1 and l 2 that do not intersect in the support of f and are related by a diffeomorphism, as in Figure 1. The broken ray transform Bf is the integral of f along ν Bf (ν) = f (ν(t))dt.
One way to think about this is to imagine that there is a curve smoothly connecting these two half-lines. Then it becomes an X-ray type transform over smooth curves. The connecting curve plays no rule in the analysis, if we always assume that the distribution f is compactly supported away from it. The goal is to understand what part of the wave front set WF(f ) can be recovered from the transform Bf . The wave front set allows us to recover, in particular, the location and the size of jumps of f . Conjugate points naturally exist for broken rays, for instance, when rays are reflected from a curved boundary. One would expect and we confirm that recovery of singularities are affected by the existence of conjugate points on ν. Much work has been done for the class of X-ray type transform with conjugate points [32,31,22,11]. In the case of transform for a generic family of smooth curves [8], if there are no conjugate points, the Partly supported by NSF Grant DMS-1600327. localized normal operator is an elliptic pseudodifferential operator of order −1. Injectivity and the stability estimates are established, which in particular implies that we can recover the singularities uniquely. When conjugate points exist, however, artifacts may arise, and in some situations they cannot be resolved. A similar situations occurs in synthetic aperture radar imaging [32]; it is impossible to recover WF(f ) if the singularities hit the trajectory only once, because of the existence of mirror points. On the other hand, if the trajectory is the boundary of a strictly convex domain and we have the priori that f has singularities in a compact set, then we can recover WF(f ) from the global data. However, this is a global procedure and we cannot expect a local reconstruction. In the case of X-ray transforms over geodesic-like families of curves with conjugate points of fold type, a detailed description of the normal operator is given in [31]. The analysis of the normal operator for general conjugate points is in [12]. Further, [22] shows that regardless of the type of the conjugate points, the geodesic ray transform on Riemannian surfaces is always unstable and we have loss of all derivatives, which leads to the artifacts in the reconstruction. It is also proved that the attenuated geodesic ray transform is well posed under certain conditions. Most recently, [11] provides a thorough analysis of the stability of attenuated geodesic ray transform and shows what artifacts we can expect when using the Landweber iterative reconstruction for unstable problems.
As mentioned above, an important example of this setting is the broken ray transform in the usual sense, which is also called V-line Radon transform. As is shown in Figure 1, the diffeomorphism is given by the law of reflection. The inversion of the V-line Radon transform arises in Single Photon Emission Computed Tomography (SPECT) with Compton cameras in two dimensions. SPECT based on Anger camera is a widely used technique for functional imaging in medical diagnosis and biological research. The using of Compton camera in SPECT is proposed to greatly improve the sensitivity and resolution [36,6,29]. The gamma photons are emitted proportionally to markers density and then are scattered by two detectors. Photons can be traced back to broken lines. The mathematical model is the cone transform (or conical Radon transform) of an unknown density. Various inversion approaches for certain cases are proposed in [2,27,35,30,21,9,33,20,18,28,23,34,24,10]. The V-line Radon transform is a special case in two dimensions where each vertex is restricted on a curve and is associated with a single axis. There are also some injectivity and stability results when we allow the rays to reflect on the boundary more than once [16,14,17,15]. These reconstructions are from full data and most of them assume specific boundaries at least for the reflection part, for example a flat one or a circle. It also should be mentioned that the broken ray transform or the V-line Radon transform sometimes refers to a different transform from the one we consider in this work, see [25,7,1,19]. In their settings, the V-line vertices are inside the object with a fixed axis direction. The integral near the vertices in the support of f makes it possible to recover singularities there. In this work, however, the vertices are always away from support of f and the axis direction changes, which make the recovery more difficult.
Another motivation is the application of parallel ray transform in X-ray luminescence computed tomography (XLCT). A multiple pinhole collimator based on XLCT is proposed in [37] to promote photon utilization efficiency in a single pinhole collimator. In this method, multiple X-ray beams are generated to scan a sample at multiple positions simultaneously, which we mathematically model by the parallel ray transform, see Section 5.
We are inspired by the spirit of [32,22,11]. In Section 2 and 3, we introduce the definition of conjugate points and conjugate covectors along broken rays and give a characterization of them. Then we consider the local problem, that is, the data Bf is known in a small neighborhood of a fixed broken ray. We show that B is an FIO and the image of two conjugate covectors under its canonical relation are identical. Singularities can be canceled by these conjugate covectors. This implies that we can only reconstruct f up to an error in the microlocal kernel. We also provide a similar analysis for the numerical result as in [11], if the Landweber iteration is used to reconstruct f . The main results are listed in the following (1) The local problem is ill-posed if there are conjugate points. Singularities cannot be recovered uniquely. (2) The global problem might be well-posed in some cases for most singularities, depending on a dynamical system inside the domain.
In Section 4 and Section 5, we apply these conclusions to two special examples, the V-line Radon transform and the parallel ray transform, as mentioned above. The conjugate points appearing in the V-line Radon transform coincide with the caustics in geometrical optics, see [3]. Additionally, when the boundary is a circle, we show that there exists conjugate points of fold type as well as cusps. Geometrically, the caustic inside a circle is an interesting problem itself, which can be traced back to the middle of 19th century [4]. We discuss the local and global recovery of singularities and we perform numerical experiments to illustrate the results. In particular, for the reflection case in a circle, we connect our analysis with the inversion formula derived in [24].
Some assumptions and notation are introduced in the following. Throughout this work, we assume f is a distribution supported in compact set in R 2 . As mentioned before, a broken ray ν is the union of two directed lines l 1 and l 2 related by a diffeomorphism χ. The two lines do not intersect with each other in the support of f . We always suppose there is a smooth family of them. Let v(α) = (cos α, sin α) and w(α) = (− sin α, cos α). We use (s j , α j ) to represent a directed line {x ∈ R 2 |x · w(α j ) = s j } with the direction v(α j ) and the unit normal w(α j ), for j = 1, 2. Suppose l 1 is the incoming part represented by (s 1 , α 1 ) and l 2 is the reflected part represented by (s 2 , α 2 ). They are related by χ : (s 1 , α 1 ) → (s 2 , α 2 ).
Acknowledgements. The author would like to express the acknowledgments to Prof. Plamen Stefanov for suggesting this problem and for the patient guidance and helpful suggestions he has provided throughout this project.

CONJUGATE POINTS
On a Riemannian manifold, the conjugate vectors of a fixed point p are vectors such that the differential of the exponential map exp p (v) with respect to v is not an isomorphism. The conjugate points are the image of these vectors under the exponential map. In [22], it is shown that singularities can be canceled by conjugate points in the geodesic ray transform case. Conjugate points also exist in the case of broken ray transform, for exmaple, the caustics in geometrical optics, see [4,3]. The light rays reflected or refracted by a curved surface form an envelope, which are conjugate points of the light source.
Suppose the incoming ray l 1 starts from p with the angle α 1 . Consider the conjugate points of p belonging to l 2 . We show below the conjugate points on l 2 do not depend on what kind of connecting curve we choose, if l 1 and l 2 are given. Notice the X-ray parametrization (p, α 1 ) gives us a Radon parameterization (s 1 , α 1 ) by s 1 = p, w(α 1 ) . Fix p, for each α 1 , the diffeomorphism χ gives another ray (s 2 , α 2 ). Suppose the reflected ray starts at q 0 at time t 2 . Here q 0 = q 0 (s 1 , α 1 ) is chosen in a smooth way and should always satisfy q 0 , w(α 2 ) = s 2 . Then the exponential map is given in the following To calculate the differential of the exponential map, we choose α 1 and t as the parameters, that is, we change the starting angle and the length. We have Then the Jacobi matrix ∂l 2 (t) ∂(t,α 1 ) has the following determinant det(dexp p (tv(α 1 ))) = det ∂l 2 ∂t ∂l 2 ∂α 1 The last equality comes from the observation that det[v(α 2 ) w(α 2 )] = 1. Thus, the determinant vanishes if and only if (t − t 2 ) dα 2 dα 1 = − dq 0 dα 1 , w(α 2 ) . Notice if we assume q 0 is the starting point of l 2 , this equation makes sense when t − t 2 > 0, that is, when dα 2 Then differentiating q 0 , w(α 2 ) = s 2 with respect to α 1 shows If q is the point on l 2 at t such that dexp p (tv 1 ) is not an isomorphism, then we have t − t 2 = q − q 0 , v(α 2 ) . By (1)(2), q should satisfy On the contrary, if there exists q on l 2 such that the equation (3) is true, then the determinant of the differential will be zero. Notice dα 2 dα 1 and ds 2 dα 1 cannot vanish at the same time with the assumption that χ is a diffeomorphism. This proves the following. Proposition 1. Suppose the incoming ray l 1 starting from p represented by (s 1 , α 1 ) and the reflected ray l 2 starting from q 0 represented by (s 2 , α 2 ) are related by a diffeomorphism χ. The point q 0 is chosen smoothly depending on (s 1 , α 1 ). Then (a) p has a conjugate point q belonging to l 2 if and only if Remark. If we consider the whole straight line which l 2 belongs to instead of the ray, then we can always find one and the only one conjugate point q satisfying (b), unless dα 2 dα 1 = 0. The condition (a) is to check whether this q belong to the reflected ray that we define. Additionally, if we perturb q 0 a little bit, that is, let q 0 = q 0 + (α 1 )v(α 2 ). Then This shows a small enough perturbation of q 0 doesn't change the sign of dq 0 dα 1 , w(α 2 ) . Therefore the existence of conjugated points is not affected by the choice of q 0 in a small neighborhood.
The proof above implies p and q are conjugate points to each other in some sense by the following reasons. Suppose we know that p and q belong to ν. The point q is the conjugate point of p if and only if .
Solving p, v(α 1 ) out, we have . Now let ν be a broken ray coinciding with ν but in the opposite direction. Actually ν is one of the family of broken rays that are associated with χ −1 . We list the Jacobian matrix in the following dχ = ∂s 2 ∂s 1 Notice equation (4) exactly means p is the conjugate point of q along ν .

MICROLOCAL ANALYSIS FOR LOCAL PROBLEMS
Recall the definition of a broken ray in Section 1. We define the broken ray transform Bf as the integral of f along ν where a(y, α) is a smooth nonzero weight. Notice χ could be the reflection operator, in which case we have a broken ray transform in the classical sense. Suppose f has support in a compact subset away from the connecting part. Then the support of f implies the transform can be interpreted as the sum of Radon transforms over two lines. We can only expect to recover the singularities in their conormal bundle. Thus, for a fixed broken ray ν 0 , we consider (x 1 , ξ 1 ) and (x 2 , ξ 2 ) in its incoming and outgoing part respectively, with ξ 1 and ξ 2 conormal to them. Let Γ(ν 0 ) be a small neighborhood of ν 0 , and U i be disjoint small open neighborhoods of x i , i = 1, 2. We choose these neighborhoods small enough, such that U 1 is disjoint from all outgoing part and U 2 is disjoint from all incoming part in Γ(ν 0 ). We extend U i to some small conic neighborhood V i of (x i , ξ i ) in the conormal bundle, for i = 1, 2.
Here we use the notation R i because in the small neighborhood Γ(ν 0 ), Bf 1 is the Radon transform of f 1 and Bf 2 can be regarded as the Radon Transform performing along the line (s 2 , α 2 ). More precisely, the restricted operators R 1 and R 2 have the following form: where φ(s, α) is a smooth cutoff function with supp φ ⊂ Γ(ν 0 ); ϕ i are cutoff pseudodifferential operators with essential support in V i , i = 1, 2; the pull back χ * g(s, α) = g(χ(s, α)) is induced by the diffeomorphism χ. We should note that outside Γ(ν 0 ), there might be another broken ray which carries the singularities (x 1 , ξ 1 ) but with it in the outgoing part. Thus, we actually multiply φ to B itself as well to make equation (6) valid.
To analyze the canonical relation of R 1 and R 2 , we need some conclusions of Radon transform. The weighted Radon transform R is defined in the following: where w(α) = (− sin α, cos α) as before, and a(y, α) is a smooth weight function.

Proposition 2. The Radon transform R is an FIO with the canonical relation
where v(α) = (cos α, sin α) and w(α) are defined as before. Specifically, C R has two components, corresponding to the choice of the sign of λ. Each component is a local diffeomorphism. The inverse is also a local diffeomorphism.
Proof. We write the Radon transform as The characteristic manifold is Z = {(s, α, y)|Φ(s, α, y) = λ(s − y · w) = 0}. Then the Lagrangian Λ is given by Therefore, the Radon transform is an FIO associated with Λ and the canonical relation C R is obtained by twisting the Lagrangian. The sign of λ is chosen corresponding to the orientation of η with respect to w. It is elliptic at (y, η) if and only if a(y, α) = 0 for α such that w(α) is colinear with η.
The composition operator χ * R is also an FIO of which the canonical relation C χ * R = C χ * • C R is a local diffeomorphism. Additionally, since the multiplication of cutoff functions does not influence the Lagrangian, therefore, the restricted operator R i has the same canonical relation as above, call it C i , in the small neighborhood of ν 0 , i = 1, 2.
Suppose (s 1 , α 1 , s 1 , α 1 ) and (s 2 , α 2 , s 2 , α 2 ) are images of (x 1 , ξ 1 ) and (x 2 , ξ 2 ) under C R . That is, with s i and α i given by (8), we have Then from the analysis above, The first equality says there is a broken ray ν of which (s 1 , α 1 ) and (s 2 , α 2 ) are the incoming and outgoing part. The second condition is equivalent to Notice ∂α 2 ∂α 1 − x 1 , v(α 1 ) ∂α 2 ∂s 1 and ∂s 2 ∂α 1 − x 1 , v(α 1 ) ∂s 2 ∂s 1 are exactly dα 2 dα 1 and ds 2 dα 1 if we fixed x 1 and consider s 2 , α 2 as functions of one variable α 1 . Therefore (11) is true if and only if if and only if there is a broken ray ν joining x and y such that (a) x and y are conjugate points along ν.
where α 1 is the angle of the incoming part and α 2 is the angle of the reflected part of ν.
For (y, η) satisfying this theorem, we call it the conjugate covectors of (x, ξ). Since C 1 is a local diffeomorphism, it maps a small neighborhood of (x, ξ) to a small neighborhood of (s 1 , α 1 , s 1 , α 1 ). The similar is true with C 2 . Then by shrinking V 1 and V 2 a bit, we can assume We can also define F 21 and C 21 in a similar way. This proves the following. Then if and only if Thus, given a distribution f 1 singular in V 1 , there exists a distribution f 2 singular in V 2 such that B(f 1 + f 2 ) is smooth. One possible choice is f 2 = −F 12 f 1 . It is also the only choice if we consider it up to smooth functions. If we introduce the definition of the microlocal kernel as in [11], then for any h with the wave front set in V 1 , h − F 12 h is in the microlocal kernel of B. This implies the reconstruction of f = f 1 always has some error in form of h − F 12 h for some h. In other words, the singularities of f cannot be resolved from the singularities of B(f ) and it can only be recovered up to an error in the microlocal kernel. A more detailed description of the kernel is in [11].
With the notation above, we are going to find out the artifacts arising when we use the backprojection B * B to reconstruct f . Without loss of generality, we assume the weight a(y, α) = 1 in the following. Suppose ν is the broken ray in Theorem 1. In the small neighborhood of ν, we have Recall R 1 and R 2 are defined microlocally. On the one hand, the assumption on f i play the same role as restricting the operator on U i , i = 1, 2. For simplification, we just ignore them. On the other hand, if we concentrate on the small neighborhood of ν, then we exclude the broken ray that carries (x, ξ) on its outgoing part. Microlocally R 1 is equivalent to the Radon transform operator near (x, ξ), which indicates R * 1 R 1 is an elliptic pseudodifferential operator of order −1 for distributions singular near (x, ξ). Especially, it has the principal symbol 4π |ξ| . The similar is true for R * 2 R 2 .
One can follow the same argument in [22,11] to show the properties of the normal operators. Additionally, similar to Radon transform, we can apply a filter to the backprojection to get a better reconstruction. Thus, we have i is the inverse of that of R i , and therefore by Egorov's theorem [13], R * i ΛR i is a pseudodifferential operator of order 0 with principal symbol σ(R where σ(P ) refers to the principal symbol of P and g i is the canonical transformation corresponding to C R i for i = 1, 2. Recall Proposition 2, we have σ(Λ) • g i = 1 4π |ξ|, which implies This also coincides with the inversion formula for Radon transform. Then with the observation R * 1 ΛR 2 F 21 = R * 1 ΛR 1 and F 21 R * 1 ΛR 2 ≡ I, we have R * 1 ΛR 2 ≡ F 12 up to a lower order. The same is true with R * 2 ΛR 1 . Notice he following calculations are all microlocal and up to order −1. Hence, we have where we follow the convention in [11] to think f = f 1 + f 2 as vector functions. This implies when performing the filtered backprojection, the reconstruction has two parts of artifacts, F 12 f 2 in V 1 and F 21 f 1 in V 2 . In [11], it is also shown that F 12 and F 21 are principally unitary in H − 1 2 , and the artifacts have the same strength.
Next, consider the numerical reconstruction by using the Landweber iteration. We still focus on the local problem, that is, we consider Bf in the small neighborhood Γ(ν 0 ) of fixed ν 0 . With the notation above, we use a slightly different Landweber iteration to solve the equation Bf = g, with g being the local data and in the range of B. Here, we set L = Λ 1 2 B to have (14) (Id − (Id − γL * L))f = γL * Λ 1 2 g.
Then with a small enough and suitable γ > 0, it can be solved by the Neumann series Suppose the original function is f = f 1 + f 2 . We track the terms of highest order, that is, order zero, to have the approximation sequence With the observation M k = 2 k−1 M for k ≥ 1, a straightforward calculation shows The numerical solution is Therefore, the error equals to 1 2 (f 1 − F 12 f 2 ) + 1 2 (f 2 − F 21 f 1 ), which belongs to the microlocal kernel.

THE V-LINE RADON TRANSFORM
In this section we are going to apply the conclusions to the V-line Radon transform, that is, the case when χ is a reflection and obeys the law in geometric optics. We consider the situation when the weight a(y, α) = 1. First we verify the reflection operator is a diffeomorphism. Then it is followed by the potential cancellation of singularities due to the existence of conjugates points. We derive an explicit formula to illustrate when conjugate points exist in this case.
FIGURE 3. A sketch of a broken ray reflected on a smooth boundary and the notation.

Conjugate Points.
The incoming ray l 1 (t) and reflected ray l 2 (t) are given in the following where q 0 = γ(τ 0 ) is the intersection point on the boundary. Compared with (2), now q 0 connects l 1 and l 2 and t 1 = t 2 . We use t 1 instead of t 2 in the following. By equation (15)(17), a straightforward calculation shows where t 1 = q 0 − p, v(α 1 ) is the time or length from p to q 0 . Plugging these back into (1), we get the determinant is (t − t 1 ) dα 2 dα 1 − t 1 . Especially, the matrix is in the following, (b) If this happens, q is uniquely determined by ∆t 2 = dα 2 dα 1 −1 ∆t 1 , where∆t 1 = t 1 is the time or length from p to γ(τ 0 ) and ∆t 2 = q − γ(τ 0 ), v(α 2 ) is the time or length from q to γ(τ 0 ).
The statement (a) comes from the observation that the other factor dq 0 dα 1 , w(α 2 ) in Proposition 1(a) is always negative in the reflection case, as shown in Figure 3. This statement has a straightforward geometrical explanation, see Figure 4. There are conjugate points if and only if α 2 increases as α 1 increases.
For negatively oriented smooth curve which is the boundary of a convex set, the curvature κ < 0 and the inner product w(α 1 ),γ < 0. The inequality actually says |κ(τ 0 )| > | w(α 2 ),γ | 2t 1 . Additionally, since w(α 1 ),γ = − cos β where β is the incident and the reflected angle, each component involved in the criterion is geometrical and therefore is invariant regardless of what kind of parameterization we choose for the boundary. We should mention the equation in (b) coincides with the Generalized Mirror Equation in [3] but is in different form and is derived from the perspective of the exponential map.  . Then the incoming ray has the direction along − → pq 0 , and w(α 1 ),γ(x 0 ), κ(x 0 ), t 1 can be calculated directly by definition. After simplification, the criterion is equivalent to We have the following three cases. case 1: If d > a, p has conjugate points if and only if the incoming ray hits the boundary at the region x 2 < 4 3 ad, as is shown in Figure 5(a). case 2: If d < a, p has conjugate points if and only if the incoming ray hits the boundary at the region x 2 > 4 3 ad, as is shown in Figure 5(b). case 3: If d = a, p has no conjugate points for all directions, which coincides with the fact that all rays of light emitting from the focus reflect and travel parallel to the y axis, as is shown in Figure 5  Example 2. The second example is to illustrate that we have different type of conjugate points, specifically fold and cusps, if we have a circular mirror with the light source inside. For simplification, we use α to replace α 1 in this example. We follow the notations in the paper [22]. The tangent conjugate locus S(p) is the set of all vectors V such that the differential of the exponential mapat p is not a isomorphism. The kernel of dexp p (V ) is denoted by N p (V ). By the previous results, We fix V 0 ∈ S(p). By (18), the differential dexp p (V 0 ) has the matrix form [v(α 2 ), 0], which shows N p (V 0 ) is spanned by ∂ ∂α . The conjugate vector V 0 is called of fold type, if ∂F ∂α = 0 for all (α, t) that satisfies F (α, t) = 0. Otherwise, we may have cusps. We will show in the following that the cusps exist in some cases.
We assume the origin O is at the center of the circle and the source is not there. Suppose the mirror has radius 1, then κ = −1. The tangent conjugate locus is the zero level set of F (α, t) given by Recall t 1 and cos β are smooth functions of α. A straightforward calculation shows Suppose we have conjugate points, then 2t 1 − cos β > 0. The incidence angle β ∈ (− π 2 , π 2 ) so we at most have two zeros for ∂F ∂α , • β = 0, which means the incoming ray and reflected ray coincide. This is a simple zero, because d dα sin β = t 1 − cos β = t 1 − 1 = 0.
• cos β − t 1 = 0 is true for some α 0 . This happens when pO is perpendicular to the incoming ray. We check d dα (cos β − t 1 ) = sin β = 0. This is also a simple zero. 4.3. Numerical Examples. This subsection aims to illustrate the artifacts arising in the reconstruction of V-line Radon transform by numerical experiments. We use Γ to denote the smooth family of all broken rays chosen for tomography. We say (x, ξ) is visible if there is a broken ray γ in the family of tomography such that (x, ξ) is in the conormal bundle of γ excluding the connecting part. The fact that (x, ξ) is visible does not necessarily imply that (x, ξ) is recoverable.
Example 3. In this example, we use filtered backprojection to recover f , which usually serves as the first attempt of reconstruction. We choose the domain as a disk with radius R and suppose the boundary is negatively oriented. The set of tomography Γ contains all broken rays whose incoming part has positive projection onto the tangent of the boundary. We choose f 1 to be a Gaussian concentrated near a single point, as an approximation of a delta function and f 2 to be zero. The support of f = f 1 + f 2 is in this disk. In the code, Bf is parameterized in the coordinate (x p , α) ∈ [−R, R] × [0, 2π). Here (x p , α) refers to the incoming part of a broken ray and we use it to represent the broken ray. This parameterization follows the convention in Radon transform in MATLAB. The radial coordinate x p is the value along the x -axis, which is oriented at α degree counterclockwise from the x-axis. We use the function radon to numerically construct our operator B by the following formula where (x p , α ) is given by the reflection. Since numerically Rf is known on discrete values of (x p , α), we use interpolation methods to approximate Rf (x p , α ). Similarly, B * is numerically constructed by the function iradon and interpolation methods. To better recover f , we apply the filter Λ to the data before applying the adjoint operator. The plots are shown in the Figure 6. We can clearly see the artifacts appear exactly in the location of conjugate points, compared with the caustics caused by a light source. Furthermore, they are expained by equation (12). Example 4. This example is to illustrate the reconstruction from local data by Landweber iteration. Assume for each (x, ξ) in WF(f ), it is visible and is perceived by only one broken ray. Then it has at most one conjugate point. To make it true, we use part of the circle as the reflection boundary. The tomography family Γ is the set of all broken rays which comes from the left side with vertices on the boundary.
By [11] , we choose f to be a modified Gaussian with singularities located both in certain space and in direction, that is, a coherent state, as is shown in Figure 7(a). We use the Landweber iteration to reconstruct f . The artifacts are still there after 100 iterations and the error becomes stable. Then we rotate f or move it to see what happens to the artifacts. Specifically, in Figure 7(c) and (d), f remains in the location but is rotated by some angles. In (e) and (f), we move f closer to the center and rotate it a bit. As the wave front set of f changes, the artifacts changes and always appear in the location of their conjugate vectors.  4.4. Global Problems. In this section, we consider the artifacts when we use full data for reconstruction. Suppose (x 0 , ξ 0 ) in WF(f ). There are two broken rays in Γ that carry this singularity. One broken ray ν 0 represented by (s 0 , α 0 ) has it in the incoming part, and the other one ν −1 represented by (s −1 , α −1 ) has it in the reflected part. Suppose (x 1 , ξ 1 ) and (x −1 , ξ −1 ) are its conjugate covectors along ν 0 and ν −1 , if they exist. We have the following cases.
If both (x 1 , ξ 1 ) and (x −1 , ξ −1 ) exist, then the singularities might be canceled by them. We continue to consider ν 1 , ν −2 and so on. Then we get a sequence of broken rays and conjugate covectors. We define the set of all conjugate covectors related to (x 0 , ξ 0 ) in the following If M(x 0 , ξ 0 ) contains finitely many (x i , ξ i ) whose index i is positive(or negative), we say it is incomplete in positive(or negative) direction. Otherwise, we say it is complete.

Example 5.
As is shown in the Figure 8, we use the same domain and family of tomography as in Example 2. Especially, we suppose the disk is centered at the origin for simplification. Considering a point (p 0 , ξ 0 ), . Inside a circular mirror, a sequence of broken rays and conjugate points on them.
we have a sequence of broken rays as well as the set M(p 0 , ξ 0 ).
has a solution in (0, 2 cos β). To simplify, we change the variable that d i = cos β(a i + 1). Thus, The requirement that p i is inside the domain means we are finding solutions for a i ∈ (−1, 1). case 1. a 0 = 0, which is followed by a i = 0 for any integer i. This is the case when we have p 0 at the midpoint of some chord and ξ 0 is the conormal of the chord. The same is true with all (p i , ξ i ). We have a complete M(p 0 , ξ 0 ). case 2. a i = 0. Then (19) can be reduced to the following iteration formula 1 a i+1 = 1 a i + 2 Suppose we start from some a 0 . Each time, the next 1 a i increases or decreases by 2. With 1 a 0 ∈ (−∞, −1) (1, ∞), finally we must have some 1 a i belonging to the interval (−1, 1), which mean p i goes out of the domain. In this case, M(p 0 , ξ 0 ) is always incomplete.
Next, let V i be a small conic neighborhoods of a fixed (x i , ξ i ) ∈ M(x 0 , ξ 0 ) and U i = π(V i ). Let f i be the restriction of f on U i . By shrinking V i carefully, we have C( Then the cancellation of singularities shows, For M(x 0 , ξ 0 ) that is finite in positive direction, finally we have By applying the diffeomorphism χ * and forward substitution, we can show all f i must be smooth. It is similar if M(x 0 , ξ 0 ) is incomplete in negative direction. This proves when M(x 0 , ξ 0 ) is incomplete, (x 0 , ξ 0 ) is a recoverable singularity.

Corollary 2. Suppose everything as in Example 5. Then
Example 6. With the same set up as above, we first choose f 1 to be a modified Gaussian of coherent state whose singularities are not radial. To compare, then we choose f 2 to be with radial singularities. As is shown in Figure 9, after performing Landweber iteration of 100 steps, all artifacts fade out and the reconstruction has a small error if f has non-radial singularities. On the contrary, if f has radial singularities, the error still decreases as the iteration but in a much slower speed. In these two cases, since f is only supported in a small set, the artifacts arising in the reconstruction may seem not so obvious. However, when f is more complicated, the artifacts might be unignorable. In the following we choose f 3 to be a Modified Shepp-Logan phantom. The error plots of these three cases are in Figure 11 to better illustrate the difference between radial and non-radial singularities. They also show where the artifacts appear( for more details, see 4.5). It is clear to see the error of reconstruction is much smaller when we have non-radial singularities than radial ones. Example 7. In this example we consider the reconstruction of the V-line Radon transform in an elliptical domain Q from global data. By [26], the billiard trajectory in an elliptical table has the following cases. If the trajectory crosses one of the focal points, then it converges to the major axis of Q. If the trajectory crosses the line segment between the two focal points, then it is tangent to a unique hyperbola, which is determined by the trajectory and shares the same focal points with Q. If it does not cross the line segment between the two focal points, then it is tangent to a unique ellipse, which shares the same focal points with Q.
In the following numerical experiments, we choose f as a coherent state. It is located and rotated such that the trajectory carrying its singularities falls into the last two cases above. We use Landweber iterations to reconstruct f by iterating 100 steps. As in Figure 12, in the reconstruction of the first coherent state, the artifacts disappear as we iterate, since some conjugate points are outside the domain. On the contrary, with conjugate points staying in the domain at least for the first reflection, there is a relative larger error in the reconstruction of the second one. A more complete analysis of the ellipse case is behind the scope of this work.

4.5.
Comparison with previous results when the boundary is a circle. This subsection is to connect our analysis to the results in [24]. By expanding f and the data Bf as Fourier series with respect to the angular variable, [24] gives an inversion formula (2.8) for V-line Radon transform with vertices on the circle. The denominator inside the integral has zeros for certain radius r and with noises it could be very unstable. This indicates we can expect certain patterns of the artifacts in the reconstruction. We are going to show these artifacts predicted by (2.8) coincides with the conjugate covectors of radial singularities in the following.
When (x 0 , ξ 0 ) is radial, M(x 0 , ξ 0 ) is complete and we have two cases. One is the case that M(x 0 , ξ 0 ) is a periodic set with period p. That is, the broken rays that carry (x 0 , ξ 0 ) after several reflections form a regular polygon of p edges, a convex or star one. The set P of all possible regular polygons can be described by the Schläfli symbol [5], P = {(p/q), p, q ∈ N, 2 ≤ 2q < p, gcd(p, q) = 1}.
Here (p/q) refers to a regular polygon with p sides which winds q times around its center. When q = 1, it is a convex regular one; otherwise it is a star one. For the polygon (p/q), the internal angle equals to π(p−2q) p . This implies |x| = cos qπ p , where x is the midpoint of one edge. Suppose Bf is smooth. We have By forward substitution, we get When p is odd, f 0 must be smooth, which implies f is smooth and therefore (x 0 , ξ 0 ) is recoverable. When p is even, it possibly causes artifacts. These artifacts are located at radius |x| = cos (2k+1)π 2n , where p = 2n and q = 2k + 1 with 0 ≤ 2q < p. These radius are exactly the positive solution of s such that cos(n(arcsin(s) − π/2)) = 0 in Formula (2.8) in [24].
In the following example, we use the same function as in Figure 9 but move them closer to the origin. The plot of error shows the artifacts are centered at the midpoint of each edges of regular stars.
We should mention that in the numerical reconstruction in [24], the regularization (2.12) is used to remove the instabilities caused by these zeros. Therefore the artifacts are removed but on the other hand some true singularities are removed as well. In [28], the regularization is also used in the numerical reconstruction of a smiley phantom but we can still see some artifacts caused by the radial singularities (see Figure 2 in [28]). , error for f with radial singularities after 100 iterations. The relative error e is defined as before.

THE PARALLEL RAYS TRANSFORM
We define the parallel X-ray transform as an integral transform over two or more equidistant parallel rays. The simplest case is the one over two parallel rays and is defined in the following It can be regarded as one example of the broken ray transform that we defined in Section 1, if we suppose the two rays are connected by a smooth curve outside the support of f or simply at the infinity. Additionally, the diffeomorphism χ is the translation which maps (s, α) to (s + d, α). Following the previous notations and calculation, we have dα 2 dα 1 = 1, ds 2 ds 1 = − p, v(α 1 ) .
By Proposition 1, if p in the incoming ray (s 1 , α 1 ) has a conjugate point q in the outgoing ray (s 2 , α 2 ), then q is determined by q, v(α 2 ) = p, v(α 1 ) , which implies q = p + w(α 1 )d. Then Theorem 1 shows a singularities (x, ξ) can be canceled by (y, η) if and only if x and y are conjugate points and ξ = η. It is shown in Figure 3 that the artifacts arising when we use the backprojection as the first attempt to recover f . Now we consider the reconstruction by iteration process. Suppose (x 0 , ξ 0 ) ∈ WF(f ) belongs to the ray (s, α). It has two conjugate vectors, that is, (x 0 ±w(α)d, ξ 0 ). We follow the same analysis as in the previous section to have M(x 0 , ξ 0 ) = {(x 0 + iw(α)d, ξ 0 ), i = 0, ±1, ±2, . . .}. The typology of M(x 0 , ξ 0 ) is quite clear. It is a discrete set of points which has equal distance. Assume Pf is smooth. Then (x 0 , ξ 0 ) ∈ WF(f ) implies M(x 0 , ξ 0 ) ⊂ WF(f ), by the same argument as before. Thus, we have the following result, see also [32]. In particular, with the prior knowledge that WF(f ) is in a compact set, the singularities are recoverable.
Corollary 3. Suppose f ∈ E (R 2 ) and assume Pf is smooth. Then f is smooth.
In the numerical experiment, we use the Landweber iteration to reconstruct f . With the assumption that f ∈ E (R 2 ), a cutoff operator is performed at every step. After 100 iterations, we get a quite good reconstruction (with f (100) − f ∞ = 0.003).
It should be mentioned that Corollary 3 shows f with singularities in a compact set could be recovered from the global data. This implies when performing the transform, we move the parallel rays around until all of them leave the compact set. In fact, from our analysis above, the condition that the rays leaving at least one side the compact set is enough. On the other hand, the local problem (illumination of a region of interest only) could create artifacts.