INVERSION OF THE SPHERICAL RADON TRANSFORM ON SPHERES THROUGH THE ORIGIN USING THE REGULAR RADON TRANSFORM

. A spherical Radon transform whose integral domain is a sphere has many applications in partial diﬀerential equations as well as tomography. This paper is devoted to the spherical Radon transform which assigns to a given function its integrals over the set of spheres passing through the origin. We present a relation between this spherical Radon transform and the regular Radon transform, and we provide a new inversion formula for the spherical Radon transform using this relation. Numerical simulations were performed to demonstrate the suggested algorithm in dimension 2.

1. Introduction. The Radon transform, an integral transform which maps a given function onto its integrals over the lines in two-dimensional space, was introduced in 1917 by Radon. The Radon transform has many applications to partial differential equations [13,17], as well as to science, engineering, and medicine [2,3,30]. Once these application were perceived, many mathematicians started to study generalizations of the Radon transform to integrations over various domains: spheres, ellipsoids, broken rays, and so on.
Among these, a Radon transform whose a integral domain is a sphere is called a "spherical Radon transform." More precisely, the spherical Radon transform assigns to a given function its integrals over a set of spheres. This transform is important in mathematics [11,13,27,33], as well as in applications to tomography, including SONAR [5,12,23], seismic testing, and RADAR. There are many types of spherical Radon transforms: the centers of the spheres of integration are centered at any point in the whole space with a fixed radius [33], the centers are on a hyperplane and the radius is variable [1,5,18,20,24,28,32], or the centers are on a sphere and the radius is variable [6,7,8,13,14,15,16]. Here, we study another type of spherical Radon transform: the integrals over the set of all spheres passing through the origin in R n . This transform provides the integral of a given function over each sphere. Cormack and Quinto discovered an inversion formula for this transform by using spherical harmonics, and this inversion formula implies a support restriction in [2,4]. Rhee found a different inversion formula in odd dimensions [29]. Also, Quinto discussed the null space and range for this spherical Radon transform in [25] and he presented how to recover a given function from the spherical Radon transform using a singular value decomposition in [26]. Yagle showed the spherical Radon transform is related to the back projection of the regular Radon transform, and using this fact, found a new inversion formula [31].
We introduce a map which transforms straight lines into circles through the origin. Using this map, we show that the spherical Radon transform can be reduced to the regular Radon transform. Using the relation, we provide a new inversion formula. Additionally, this spherical Radon transform becomes a circular-arc Radon transform in the 2-dimensional case when the given function has compact support in the upper half plane. In this case, it is related to Compton scattering tomography [21]. Some works [21,22] have studied this circular-arc Radon transform.
The rest of this paper is organized as follows: in section 2, we formulate the spherical Radon transform and show how to reduce this spherical Radon transform to the regular Radon transform. We present numerical simulations to demonstrate the suggested 2-dimensional algorithm in section 3.
• The map m n : x n x , |x| 2 .
If an n−1-dimensional sphere intersects the x n = 0 plane only at the origin, then m −1 n transforms it into a hyperplane. On the other hand, if an n − 1-dimensional sphere intersects the x n = 0 plane in an n − 2-dimensional sphere, m −1 n transforms it into a hyperplane missing an n − 2-dimensional plane. Changing variables using this map m n plays a critical role in reducing the spherical Radon transform to the regular Radon transform.
Proof. We can easily check that Since |x − u/2| = |u/2|, |x| 2 − x · u = 0 and thus (1) becomes Let us define a function k(y) by where y = (y , y n ) ∈ R n . Then we have for x n = 0, The regular n-dimensional Radon transform Rk(e θ , t) is defined by The next theorem shows that the spherical Radon transform R S f can be reduced to the regular Radon transform Rk.
Theorem 2. Let f ∈ C(R n ) have compact support in R n . Then we have for u n = 0, Proof. Using the Dirac delta function, R S f can be written as Changing the variables u/2 + |u|x/2 → x gives us We change the variable m −1 n (x) = y, i.e.,

SUNGHWAN MOON
Then we have where the Jacobian J(y) is By the definition of k, we have |u| n−2 R S f (u) = 2 n−1 R n f (m n (y))δ (y n − y · u − u n ) |y n | n−2 (1 + |y | 2 ) n−1 dy We recognize the right hand side as the integral along the hyperplane perpendicular to (−u , 1)/ 1 + |u | 2 with (signed) distance u n / 1 + |u | 2 from the origin. In this case, the measure for the hyperplane becomes 1 + |u | 2 dy . Setting e θ = (−u , 1)/ 1 + |u | 2 and t = u n / 1 + |u | 2 , we have our assertion. Now we can prove our main theorem.
Theorem 3. Let f ∈ C ∞ (R n ) have compact support in R n . If g(u) = |u| n−2 R S f (u), then we have for x n = 0, where P.V. means the Cauchy principal value.
(Note that e θ = (e θ sin θ n−1 , cos θ n−1 ) ∈ S n−1 .) By the projection slice theorem for the regular Radon transform, we have =ĝ(e θ tan θ n−1 , σ cos θ n−1 ), wherek is the n-dimensional Fourier transform of k andĝ is the 1-dimensional Fourier transform of g with respect to u n . Thus we have We have for x n = 0, 3. Numerical implementation. Here we discuss the results of 2-dimensional numerical implementations. Although many inversion formulas were derived in [2,4,25,26,31], neither these works nor any others in the literature present numerical results for these formulas.
In the experiments presented here we use the phantom shown in Figure 1  This implies that the function k has support in the unit ball and this makes it sufficient to consider the range [−1, 1] in t.) The 128 × 128 images are used in Figure 1. To reconstruct the image in Figure 1 (b), we use 256 × 256 projections for θ and t in 2Rk(e θ , t) = | sec θ|R S f (− tan θ, t sec θ) obtained by Theorem 2. After obtaining the function k using the inversion code for the regular Radon transform, we reconstruct the function f using (2). All computations have been performed in Matlab. While Figure 1 (b) demonstrates the image reconstructed from the exact data, Figure 1 (c) and (d) show the results of reconstruction from noisy data. The noise is modelled by normally distributed random numbers. In Figure 1 (c), the noisy data is modelled by adding the noise values scaled to 5% of the norm of the exact data to the exact data. In Figure 1 (d), the noisy data is modelled by adding to the exact data the noise value scaled to 10%. These reconstructions from the noisy data show the stability of our algorithm in the presence of noise. In Figure 2 the surface plots of all images in Figure 1 are provided.
Herem n (x) = x/|x| 2 is reflection with respect to the unit circle. Hence, if f has a nonzero value near the origin, the support of k is very large. This fact is a disadvantage in implementation because a large range of Rk in t is required. In Figure 3, we reconstruct f using the relation (4). All images in Figure 3 are 128×128 as in Figure 1. The phantom, supported within the rectangle [−0.5, 0.5] × [0, 1], is the sum of multiples of characteristic functions of disks centered at (0.4, 0.4), (0.1, 0.6), and (0, 0.61) with radii 0.03, 0.05, and 0.1, whose values are 2, 1, and 0.5, respectively. Here although we use disks with the same size as in Figure 1, the center of one disk is changed from (0, 0.3) to (0.4, 0.4) because of the support ofk. If we had chosen the disk centered at (0, 0.3), the support ofk would have been inside the disk centered at the origin with radius 100/33, so the required range of Rk in t is [−100/33, 100/33]. This is large, so to reduce the range we moved the center of one disk. Now it is enough to consider the range [−2, 2] in t. To reconstruct the image in Figure 3 (b), we use 256 × 256 projections for θ and t in Rk(e θ , t) = R S f (e θ /t)/|t| (n = 2), obtained from (4), again. After obtaining the functionk, we reconstruct the function f using f (x) =k(x/|x| 2 )/(2|x|) (n = 2).
In fact, although we use 256 × 256 projections for θ and t as in Figure 1, the distance between grid points in t is twice as large as in Figure 1 because t ranges from -2 to 2 rather than from -1 to 1 as before. Hence, to get the same resolution, we have to use 256×512 projections for θ and t and this increases the computational cost. We present the image reconstructed using 256 × 512 projections in Figure 3 (c). Figure 3 (d) and (e) show the reconstructions from noisy data with the noise value scaled 10%. Although the numerical results look similar to ours, the closer the phantom is to zero, the bigger the range of t that is required, and thus the more computational cost is required to get the same resolution in this approach using the relation (4). In Figure 4 the surface plots of all images in Figure 3 are provided.

4.
Conclusion. This paper is devoted to studying the spherical Radon transform, which is one of the classic problems in computed tomography: the spherical Radon transform mapping a given function to its integrals over the set of spheres passing through the origin. We suggest the reduction of this spherical Radon transform to the regular Radon transform and provide a new inversion formula for this transform. Also, we show numerical simulations to demonstrate our algorithm.