Inversion formula for an integral geometry problem over surfaces of revolution

An integral geometry problem is considered on a family of n$n$ ‐dimensional surfaces of revolution whose vertices lie on a hyperplane and directions of symmetry axes are fixed and orthogonal to this plane, in Rn+1$\mathbb {R} ^{n+1}$ . More precisely, the reconstruction of a function f(x,y)$f(\mathbf {x,}y)$ , x∈Rn$\mathbf {x}\in \mathbb {R} ^{n}$ , y∈R$y\in \mathbb {R}$ , from the integrals of the form f(x,y)dx$f(\mathbf {x,}y) d\mathbf {x}$ extended over a chosen side of all surfaces of revolution of a given family is investigated. Unlike the usual Radon transform, the integrals considered here are not taken with respect to the surface area element. A Fourier slice identity and a backprojection‐type inversion formula are obtained with a method based on the Fourier and Hankel transforms. The reconstruction procedure and some analytical and numerical implementations of the obtained inversion formulas in the cases of n=1$n=1$ and n=2$n=2$ are provided.


INTRODUCTION
The Radon transform associates a function defined in -dimensional space with its integrals over all hyperplanes in the same space. 1The problem of inverting this transform can be considered as the basis and a special case of the problems of integral geometry, where the integrals are considered over the manifolds of a given family of manifolds. 2-4The Radon transform and its generalizations provide mathematical background for the reconstruction problems that arise in several important applications, especially imaging applications such as medical, seismic, or optical imaging (see, e.g., Refs.[5][6][7][8]).The importance of these transformations and their inversions lies mainly in the fact that they allow the internal structure of an object to be determined from the Radon projections without damaging this object.As a well-known imaging application, in transmission computed tomography the object is probed with some radiation using the scanning geometries, such as parallel beam, fan beam, or cone beam geometry, and the measured projections are represented by the Radon-type transforms.In practice, there is no common scanning geometry and inversion procedure, and developing efficient and feasible reconstruction algorithms for a particular industrial application depends on the complexity of the model and physical limitations.Specifically, backprojection-type inversion algorithms are based on the knowledge of the exact analytical inversion formulas.In this study, we introduce an integral geometry problem on a family of -dimensional surfaces of revolution and provide an analytical formula for the solution of this problem with some examples on the use of this formula.
Let (,) be a smooth real-valued function on ℝ +1 with compact support in ℝ  × (, ∞).We investigate the problem of determining the function  from its integrals for  ∈ ℝ  ,  > 0. Here, the integral in the first equality is the surface integral of the form (,) extended over the upper side, the side where the angle made by the normal of the surface with the positive direction of the -axis is acute, of the surface   , , and  represents the volume of the projection of the surface element on the plane  = 0.This integral form is called a surface integral of the second type or an integral over a projection onto a coordinate plane, and it can be reduced to an ordinary multiple integral (see, e.g., Refs.[9, section 8.5.2] and [10, section 22.4]).Considering that the projection of a surface   , onto the plane  = 0 is the whole ℝ  and denoting the standard Euclidean volume element on ℝ  by , the second integral in (2) is the ordinary multiple integral reduced from the first one. maps a function to the integrals (2) over all surfaces of a given family {  , } in the sense described above, and it can be considered as a generalized Radon transform.It is worth noting that, the integrals in (2) are not taken with respect to the surface area element of   , , that is, they are not the surface integrals of first type, as in many Radon-type transforms.Inversion of generalized Radon transforms over some families of plane curves with the vertical shifting and scaling parameterization is studied in Refs.[11,12].Similar to the surface integral of the second type, the line integrals in these seismic-type Radon transforms are taken with respect to the standard element  on ℝ, not with respect to the arc length elements of the curves.In Ref. [13], an inversion formula is obtained for the problem of reconstructing a function from the second-type surface integrals of it over -dimensional paraboloids  −  = | − | 2 , by a method based on reducing the problem to an equivalent boundary value problem.Inversion of generalized Radon transforms over -dimensional paraboloids with fixed focal length and fixed direction of symmetry axis, without any restrictions on their vertices are also considered in Ref. [14] with a method based on the Fourier transform and the solution of Volterra integral equations of the first kind.The parameters of the surfaces considered here differ from those just mentioned above, and the inversion result in this study is mainly based on employing the Fourier and Hankel transforms, in a similar way that is used in, for example, Refs.[15][16][17][18].In Ref. [16], the problem of reconstructing a function from its -dimensional spherical averages is considered and an inversion formula is obtained, and under more general conditions on the reconstructed function and its spherical averages this inversion formula is provided in Ref. [15].In Refs.[17] and [18], inversions of -dimensional conical Radon transform and -dimensional elliptical Radon transform are derived, respectively.If the family { , } consists of conical surfaces, then  can be considered as the conical Radon transform of , and our inversion formula reduces to the one given in Ref. [17].
From the applied point of view, the conical Radon transform is related with emission imaging based on Compton scattered radiation. 19In emission computed tomography, after injection of some radioactive tracers to an object of interest and using the external projection measurements of emitting photons, one determines the distribution of radiotracer inside the object.These projections can be represented by the Radon transform in the case of nonscattered radiation, while scattered radiation gives rise to the generalized Radon transforms; conical projections lead to the conical Radon transforms.It is worth mentioning that, to investigate the reconstruction problems in practice, such as tomography, one must consider the transforms not only for smooth functions but also for piecewise smooth functions and analyze the inversion formula in a neighborhood of a point of nonsmoothness, as indicated in Ref. [2, p. 4].
The rest of this paper is organized as follows.In Section 2, an analog of the Fourier slice identity is derived.Using this identity, an explicit inversion formula for , which is the main result of this paper, is provided in Section 3. Finally, in Section 4, some analytical and numerical implementations of the obtained inversion formula are presented for the cases of  = 1 and  = 2.

FOURIER SLICE IDENTITY
Let us introduce the notation that will be used throughout the paper.The Fourier transform and its inverse with respect to the first argument are denoted by Denoting by  (−) the Riesz fractional derivative operator, we have the relations and  (−) yields the (left) inverse to the Riesz potential operator  where   is the th-order Bessel function of the first kind. 22Some sufficient conditions for the validity of the formulas ( 4) and ( 8) can be seen in, for example, Ref.
can be written as or, with the change of variables  = ∕, as where  is surface measure on  −1 .
Then for (, ) ∈ ℝ  × (0, ∞), the following identity holds where Proof.By taking the Fourier transform of ( 10) with respect to the first component,  , followed by interchanging the order of integration, and using the shifting property of the Fourier transform, we obtain The inner integral above can be evaluated as where  ∕2−1 is the Bessel function of the first kind of order (∕2 − 1).The proof of the identity (12) can be found in Refs.[25, section 4.3] or [14], and a general form of this formula, based on the Funk-Hecke theorem on spherical harmonics, is given in Ref. [6, p. 198].So, from (12) it turns out that Introducing the notation and using the definition of the Hankel transform (7), we obtain the Fourier slice identity

INVERSION
Using the above definitions, notations, and the Fourier slice identity (11), we now present the main result of this paper: Theorem 2 provides the inversion formula for  defined in (2) over the surfaces of revolution of the family {  , } determined by a function .
Proof.Since  ∈  ∞ (ℝ +1 ) have compact support in ℝ  × (, ∞), considering the definitions of the inverse Hankel transform (8), with respect to the second argument, and the inverse Fourier transform (4), with respect to the first argument, and denoting (, ) = ( −1 (, ))(, ), we have By using the identity (11), followed by the identity (12), (, ) becomes Considering relation (6) and the definition of the inverse Fourier transform, or directly from (5), we have Equation ( 13) yields thus we obtain Finally, the substitution  = 1  −1 () followed by the one  =  +  yields the desired inversion formula (14), that is, .□ For  = 1, the family { , } consists of the plane curves  = (| − |∕) and we have the relation  (−1) =     , where   is the Hilbert transform operator with respect to the first argument and   is the differential operator with respect .So in the special case  = 1, we have the following consequence of Theorem 2.
For  = 2, the family {  , } consists of 2-dimensional surfaces of revolution in ℝ 3 and considering the relation  (−2) = −Δ  , where is the Laplace operator with respect to  = ( 1 ,  2 ), we have another consequence of Theorem 2.

IMPLEMENTATIONS
In this section, we present the validation of the exact inversion formulas (15) and (16) with some examples, it is not aimed to provide a comprehensive numerical inversion procedure and its convergence behavior.Here, we provide some analytical and numerical implementations using some convenient tools of the numerical analysis.
As is well known, the performance and feasibility of a numerical approach can be affected by several factors such as available projection data, noise, the evaluation of the  (−) , specifically      for  = 1 and Δ   for  = 2, and the evaluation of the arising integrals.
To implement the above inversion formulas, in each example given below, we define a test function  and a family of surfaces of revolution generated by rotating the curves  = ( 1 ∕) on ( 1 , )-plane about the line  = .In Examples 1 and 2 below, we consider some smooth functions, which are not compactly supported but satisfy some decaying conditions at infinity, and we consider some smooth compactly supported functions in Examples 3 and 4.
In Examples 1 and 2, to validate the inversion formulas, we first evaluate the projection function  and then the integrands in ( 15) and ( 16) analytically.The arising integrals are computed by employing the built-in functions "integral.m"and "integral2.m",based on the global adaptive quadrature method, of the MATLAB software using the default absolute error tolerance of 10 −10 and relative error tolerance of 10 −6 . 26,27In Examples 3 and 4, the above steps are performed fully numerically.
In the following example, we consider the case of  = 1, that is, the reconstruction of a function from the integrals (2) over the plane curves that belong to a family { , } in ℝ 2 .
on the restricted domain  = ℝ × (0, ∞), then ()(, ) can be evaluated as Now we reconstruct (, ) by using the inversion formula (15).Since and see, for example, [ 28 , p. 465], where erf i is the imaginary error function, we have with  −1 () = √ , from the inversion formula (15), we have and the evaluation of the integral in (18) yields the desired values of (, ).With the absolute error tolerance of 10 −10 and relative error tolerance of 10 −6 to limit the estimates of the absolute error or the relative error of the global adaptive quadrature method used in the standard function "integral.m" of MATLAB, 26 the integrals are numerically evaluated for 64 × 64 values of (, ) on [−2, 2] × (0, 4].The computed values are exact up to 10 decimal digits, with the default error tolerances, and Figure 1 shows the reconstructed . The following example provides an inversion for the case of  = 2, that is, the reconstruction of a function from the integrals (2) over 2-dimensional surfaces of revolution that belong to a family {  , } in ℝ 3 .
Remark 2. Examples 1 and 2 show that for  ∈  ∞ (ℝ  ), having compact support is a sufficient condition but not a necessary condition on the proposed inversion formulas to hold.Moreover, these examples provide the implementation of the inversion formulas analytically and are helpful in checking and improving the steps of the numerical inversion procedure, for example, in Examples 3 and 4.Moreover, in Step 1 below, we generate the discrete data by using the definition of  and then the reconstruction process follows with these data and the inversion formula (15) or (16).In practice, various noise models exist and it is often not possible to obtain or measure all the exact data.So, in Step 1, one needs to consider the noise model and limited data of a possible particular application for a comprehensive numerical procedure.Hence, the below given process is just the beginning step toward this goal.
The procedure of reconstructing a test function  from the integrals  by using the formulas ( 15) and ( 16), used in Examples 3 and 4 below, consists of the following steps.
Step 1 Generate the discrete test data of ()( , ).Choose a test function  and construct a rectangular grid in a set Ω.Then, numerically evaluate  at the nodes ( , ) of the grid in Ω for the given family {  , }.
• Case  = 1: Using the discrete data obtained in Step 1, compute the approximate values of the derivative (  )(, ), for example, by employing a finite difference formula.Then, using these approximate values, compute the Hilbert transform of   , that is, the values of (    )(, ), for example, by employing an appropriate numerical integration method or the standard function "hilbert.m" of the MATLAB program.
• Case  = 2: Using the discrete data obtained in Step 1, compute the approximate values of (Δ  )( , ), for example, by employing a finite difference formula.
Step 3 Reconstruct  by using the inversion formula (15) or (16).Using the values of ( (−) )( , ) obtained in Step 2, reconstruct  by evaluating the integral in (15) or ( 16) at the discrete points (, ), for ( , Here, the values of ( (−) )( , ) can be determined approximately via an appropriate interpolation method from the computed values of  (−)  at the nodes ( , ) of the grid in Ω.
Example 3. Let us consider the transform (2) defined on the family { , } given in Example 1.We consider the smooth compactly supported function  defined by as the test function, and we present the results for the sample image, given in Figure 3, defined by .
In Step 1 of the above given reconstruction procedure, an  ×  rectangular grid in Ω = [  ,   ] × [  ,   ] is constructed with the linear scaling for the values of  on [  ,   ] and logarithmic scaling for the values of  on [  ,   ], and the integrals are computed using the MATLAB function "integral.m,"with the default error tolerances.Here, the approximations by using three different data sets, generated for Ω and  given in Table 1, are considered.In Step 2, the approximate values of (  )(, ) are computed using the centered difference formula for ℎ = 0.2, and using these values, the Hilbert transform is computed with the function "hilbert.m" of MATLAB.
By scaling the colormap to the range of , 2D images representing the exact image and the approximation results on reconstructing  using the formula (15) are given in Figure 4.
In ) are determined from the values (Δ  )( , ) via the nearest neighbor interpolation.
In Figures 4 and 5, it can be seen that images similar to the exact image can be obtained using less data and smaller Ω.Moreover, it is observed that the exact values can be better approached with a larger set Ω and a sufficient number of grid points.The reason for taking relatively larger sets in Example 3 is mainly based on the computation of the Hilbert transform, which is not a local operation, of the discrete data.

CONCLUSION
The problems of integral geometry are important from both theoretical and practical point of view.
Obtaining an analytical inversion formula, that is, expressing the unknown function in terms of its given integrals over the curves or the surfaces of a given family, is one of the fundamental issues in the study of these problems.
There have been considerable results on the integral geometry problems in which families of surfaces consist of spheres, cylinders, ellipsoids, paraboloids, or cones.In this study, some families of -dimensional surfaces of revolution in ℝ +1 have been taken as the families of surfaces, and we have derived a Fourier slice identity and a backprojection-type inversion formula.The families of circular cones and circular paraboloids are two special cases of the families we have considered, and our results on inversion are valid for the integral geometry problems over these families.
The special forms of the inversion formula ( 14) have been given in the cases of  = 1 and  = 2, and some analytical and numerical reconstruction examples have been presented to implement these inversion results.Based on the obtained analytical inversion formulas and provided implementations in this study, it is possible to improve the numerical approaches for the practical applications considering the real data and noise.For a comprehensive numerical approach, feasible in practice, the incomplete data and appropriate noise modeling need to be taken into account.Moreover, the inversion needs to be analyzed in a neighborhood of a point of nonsmoothness of piecewise smooth functions.The intrinsic limitations in a reconstruction procedure for an underlying practical problem can be understood by investigating the visibility and location of singularities.The microlocal analysis on reconstructing singularities is provided in Ref. [30] (A)   2, are used, respectively.and some elementary introduction to microlocal methods in integral geometry and tomographic applications can be seen, for example, in Refs.[31][32][33].

A C K N O W L E D G M E N T
The authors thank the anonymous reviewers for their many insightful comments and suggestions.

C O N F L I C T O F I N T E R E S T S TAT E M E N T
The author declares no potential conflicts of interest.

2 F I G U R E 2
The cross section of reconstructed  by the inversion formula(16) at  = 1.

Example 4 .
Let us consider the transform (2) defined on the family of circular paraboloids {  , } given in Example 2. Here, we present the results for the sample image defined by the smooth compactly supported function