We generalize Y. Nievergelt's inversion method for the Radon transform on
lines in the 2-plane to the $k$-plane Radon transform of continuous and $L^p$
functions on $R^n$ for all $1\leq k<n$.Comment: 9 page

Ournycheva, Elena

Rubin, Boris

Fractional Calculus and Applied Analysis

English

arXiv

Mathematics Subject Classification 2010: 42C40, 44A12.In 1986 Y. Nievergelt suggested a simple formula which allows to reconstruct a continuous compactly supported function on the 2-plane from its Radon transform. This formula falls into the scope of the classical convolution-backprojection method. We show that elementary tools of fractional calculus can be used to obtain more general inversion formulas for the k-plane Radon transform of continuous and L^p functions on R^n for all 1 ≤ k < n. Further generalizations and open problems are discussed

Ournycheva, E.

Rubin, B.

Bulgarian Digital Mathematics Library at IMI-BAS

ON Y. NIEVERGELT’S INVERSION FORMULAFOR THE RADON TRANSFORME. Ournycheva ∗, B. Rubin ∗∗AbstractIn 1986 Y. Nievergelt suggested a simple formula which allows to re-construct a continuous compactly supported function on the 2-plane fromits Radon transform. This formula falls into the scope of the classicalconvolution-backprojection method. We show that elementary tools of frac-tional calculus can be used to obtain more general inversion formulas forthe k-plane Radon transform of continuous and Lp functions on Rn for all1 ≤ k < n. Further generalizations and open problems are discussed.Mathematics Subject Classification 2010 : 42C40, 44A12Key Words and Phrases: The k-plane Radon transform, Nievergelt’sinversion formula, the convolution-backprojection method.1. IntroductionInversion formulas for Radon transforms of different kinds are of greatimportance in mathematics and its applications; see, e.g., [2, 4, 7, 10, 11,18, 20, 21, 24, 29], and references therein. Since many of them are prettyinvolved or applicable under essential restrictions, every more or less ele-mentary inversion algorithm deserves special consideration. In 1986, YvesNievergelt came up with the paper [12], entitled “Elementary inversion ofRadon’s transform”. His result can be stated as follows.c© 2010, FCAA. All rights reserved.44 E. Ournycheva, B. RubinTheorem 1. LetGa(t) =1/(pia2) if |t| ≤ a,1pia2(1− 1√1− a2/t2)if |t| > a;(1)a > 0. Any continuous compactly supported function on the 2-plane canbe reconstructed from the Radon transform over lines in this plane by theformulaf(x, y) = lima→01pi∫ pi0∫ ∞−∞(Rf)(t− xcosα− y sinα, α)Ga(t) dt dα, (2)where the double integral on the right-hand side equals the average of fover the disc of radius a centered at (x, y).Formulas (1) and (2) indeed look elementary. The following questionsarise:1. What is the basic idea of the Nievergelt’s method from the point ofview of modern developments?2. Whether this method is applicable to more general k-plane Radontransforms on Rn for all 1 ≤ k < n and arbitrary continuous or Lp functions,satisfying minimal assumptions at infinity?In the present article we apply elementary facts from fractional calculus[25] to answer these questions and indicate possible generalizations.It is worth noting that Y. Nievergelt’s inversion formula is based on therelationWa ∗ f = R∗(wa ∗Rf), (3)where ∗ is the convolution operator, R∗ is the backprojection, wa is a delta-sequence, andWa = R∗(wa). As it was explained in [11], the relation (3) is abasis of various inversion algorithms, cf. [8, 9]. Our interest in generalizationof (2) is purely theoretical and no claim about practical importance of theresults is made. However, it is felt that the elementary use of operators offractional integration might appeal to the applied mathematician as well.Notation and main results. Let Gn,k and Gn,k be the affine Grass-mann manifold of all non-oriented k-planes τ in Rn, and the ordinary Grass-mann manifold of k-dimensional subspaces ζ of Rn, respectively. Givenζ ∈ Gn,k, each vector x ∈ Rn can be written as x = (x′, x′′) = x′+x′′ whereON Y. NIEVERGELT’S INVERSION FORMULA . . . 45x′ ∈ ζ and x′′ ∈ ζ⊥, ζ⊥ being the orthogonal complement to ζ in Rn. Eachk-plane τ is parameterized by the pair (ζ, x′′) where ζ ∈ Gn,k and x′′ ∈ ζ⊥.The manifold Gn,k will be endowed with the product measure dτ = dζdx′′,where dζ is the SO(n)-invariant measure onGn,k of total mass 1, and dx′′ de-notes the usual volume element on ζ⊥. We write C0 = C0(Rn) for the spaceof continuous functions on Rn vanishing at infinity; σn−1 = 2pin/2/Γ(n/2)denotes the area of the unit sphere Sn−1 in Rn.The k-plane transform of a function f on Rn is a function fˆ on Gn,kdefined byfˆ(τ) =∫ζf(x′ + x′′) dx′, τ = (ζ, x′′) ∈ Gn,k. (4)This expression is finite for all τ ∈ Gn,k if f is continuous and decays likeO(|x|−λ) with λ > k. Moreover [24, 27, 29], if f ∈ Lp(Rn), 1 ≤ p < n/k,then fˆ(τ) is finite for almost all planes τ ∈ Gn,k. The above-mentionedbounds for λ and p are best possible.Following [23], we define the wavelet-like transform(W ∗aϕ)(x) =1an∫Gn,kϕ(τ)w( |x− τ |a)dτ, a > 0, (5)where |x− τ | denotes the Euclidean distance between the point x ∈ Rn andthe k-plane τ .Theorem 2. [23, Th. 3.1] Let f ∈ Lp, 1 ≤ p < n/k, and let ψ(| · |) bea radial function on Rn, which has an integrable decreasing radial majorant.If w is a solution of the Abel type integral equationcr2−n∫ r0sn−k−1w(s)(r2−s2)k/2−1 ds=ψ(r), c= σk−1 σn−k−1σn−1, (6)then(W ∗a fˆ)(x) =∫Rnf(x− ay)ψ(|y|) dy, (7)and therefore,lima→0(W ∗a fˆ)(x) = λf(x), λ =∫Rnψ(|x|) dx. (8)The limit in (8) is understood in the Lp-norm and in the almost everywheresense. If f ∈ C0 ∩ Lp for some 1 ≤ p < n/k, then (8) holds uniformly onRn.46 E. Ournycheva, B. RubinThis theorem is a core of the convolution-backprojection method for thek-plane Radon transform, and the most difficult task is to choose relativelysimple functions w and ψ satisfying (6). The crux is that the left-handside of (6) has, in general, a bad behavior when r→∞. Hence, to achieveintegrability of ψ, the solution w must be sign-changing.Our first observation is that the essence of Y. Nievergelt’s Theorem 1can be presented in the language of Theorem 2 as follows.Theorem 3. Let k = 1, n = 2,ψ(r) ={1 if 0 ≤ r ≤ 1,0 if r > 1.(9)Then (6) has a solutionw(r) ={1 if 0 ≤ r ≤ 1,1− r√r2 − 1 if r > 1,(10)such that for every compactly supported continuous function f on R2,(W ∗a fˆ)(x) =∫|y|<1f(x− ay) dy, (11)where W ∗a is the wavelet-like transform (5) generated by w.For the convenience of presentation, we will keep to the following con-vention.Definition 4. The convolution-backprojection algorithm in Theorem2 will be called the Nievergelt’s method if ψ is chosen according to (9).Of course, Theorem 2 deals with essentially more general classes of func-tions than Theorem 1, however, the main focus of our article is different:we want to find auxiliary functions ψ and w, having possibly simple analyticexpression.Theorem A.(i) The Nievergelt’s method is applicable to the X-ray transform (thecase k = 1) in any dimension. Namely, if ψ is chosen according to (9), then(6) has a solutionw(r) =1 if 0 ≤ r ≤ 1,−Γ(n−12)r3−n2√pi Γ(n2) ∫ 10vn/2−1 (r2−v)−3/2 dv if r > 1,(12)ON Y. NIEVERGELT’S INVERSION FORMULA . . . 47and inversion formula (8) holds.(ii) If k > 1, then Nievergelt’s method is inapplicable.An integral in (12) can be expressed through the hypergeometric func-tion and explicitly evaluated in some particular cases; see Section 2.2. Forinstance, if n = 2, then (12) is the Nievergelt’s function (10).To include all 1 ≤ k < n, we modify the Nievergelt’s method by choosingψ(r) in a different way as follows.Theorem B. Letψ(r) =0 if 0 ≤ r ≤ 1,(r2 − 1)`rn+2`+1if r > 1;(13)` ≥ 0. Then the corresponding function ψ(| · |) has a decreasing radialmajorant in L1(Rn) and (6) has the following solution:(i) In the case k = 2`; ` = 1, 2, . . . :w(r) =0 if 0 ≤ r ≤ 1,c` r2+2`−n(12rddr)` [(r2 − 1)`r2`+3]if r > 1,(14)c` =Γ(n/2− `)Γ(n/2).(ii) In the case k = 2`+ 1; ` = 0, 1, 2, . . . :w(r)=0 if 0≤r≤1,c˜` r3+2`−n(12rddr)`+1 [(r2−1)`+1/2r2`+2]if r > 1,(15)c˜` =Γ((n− 1)/2− `) `!Γ(n/2) Γ(`+ 3/2).In both cases the inversion result in Theorem 2 is valid.Theorems A and B are proved in Sections 2 and 3, respectively.48 E. Ournycheva, B. RubinPossible generalizations.1◦. The convolution-backprojection method is well-developed in thegeneral context of totally geodesic Radon transforms on spaces of constantcurvature. Apart of Rn, the latter include the n-dimensional unit sphereSn and the hyperbolic space Hn; see [1, 22, 23]. As above, the key rolein this theory belongs to a certain Abel type integral equation and therelevant sign-changing solution w. Moreover, passage to the limit in (8)as a → 0, can be replaced by integration in a from 0 to ∞ against thedilation-invariant measure da/a. This leads to inversion formulas, whichresemble the classical Caldero´n’s identity for continuous wavelet transforms[3]. The corresponding wavelet function is determined as a solution of asimilar Abel type integral equation; see [1, 22, 23] for details. In all thesecases analogues of Theorems A and B can be obtained. We leave thisexercise to the interested reader.2◦. Unlike the classical k-plane transforms on Rn, the correspondingtransforms on matrix spaces [5, 15, 16, 17] are much less investigated. Tothe best of our knowledge, no pointwise inversion formulas (i.e., those, thatdo not contain operations in the sense of distributions) are available for thesetransforms if the latter are applied to arbitrary continuous or Lp functions(p 6= 2). One of the reasons of our interest in Nievergelt’s idea is that itmight be applicable to the matrix case. Moreover, as in 1◦, it may pave theway to implementation of wavelet-like transforms in the corresponding re-construction formulas. We plan to study these questions in our forthcomingpublication.2. The case k = 12.1. Proof of Theorem AWe will be dealing with Riemann-Liouville fractional integrals(Iαa+g)(u) =1Γ(α)∫ ua(u− v)α−1g(v) dv, α > 0. (16)Changing variables, we transform the basic integral equationcr2−n∫ r0sn−k−1w(s)(r2−s2)k/2−1 ds=ψ(r), c= σk−1 σn−k−1σn−1, (17)ON Y. NIEVERGELT’S INVERSION FORMULA . . . 49(cf. (6)) to the formIk/20+ w˜ = ψ˜, (18)wherew˜(u) = u(n−k)/2−1w(√u), ψ˜(u) =2un/2−1cΓ(k/2)ψ(√u). (19)Suppose that ψ is defined by (9). If 0 < u ≤ 1, then, by homogeneity,w˜(u) = c′u(n−k)/2−1, c′ =2Γ(n/2)cΓ(k/2) Γ((n−k)/2) = 1, (20)c being the constant from (17). Hence, for u > 1, we necessarily have(Ik/21+ w˜)(u) = −1Γ(k/2)∫ 10(u− v)k/2−1v(n−k)/2−1 dv. (21)If k ≥ 2, this equation has no solution w˜ ∈ L1loc(1,∞), because, otherwise,we getlimx→1+(l.h.s) = 0, limx→1+(r.h.s) = const 6= 0.This proves the second statement in Theorem A.Consider the case k = 1. If 0 < u ≤ 1, then w˜(u) = u(n−1)/2−1. If u > 1,then, setting w0(u) = u(n−1)/2−1, from (21) we get(I1/21+ [w˜ − w0])(u) = −1√pi∫ u0(u− v)−1/2v(n−1)/2−1 dv= −cn un/2−1, cn = Γ((n− 1)/2)Γ(n/2) .This givesw˜(u) = w0(u)− cn√piddu∫ u1(u− v)−1/2vn/2−1 dvor (set∫ u1 =∫ u0 −∫ 10 )w˜(u) = − cn2√pi∫ 10(u− v)−3/2vn/2−1 dv.50 E. Ournycheva, B. RubinThus, if I1/20+ w˜ = ψ˜, then, necessarily,w˜(u) =u(n−1)/2−1 if 0 < u ≤ 1,− cn2√pi∫ 10(u− v)−3/2vn/2−1 dv if u > 1.(22)One can readily see that function (22) is locally integrable on (0,∞). Letus prove that it satisfies (I1/20+ w˜)(u) = ψ˜(u) for all u > 0. It suffices to showthat (I1/20+ w˜)(u) ≡ 0 when u > 1. We have I1/20+ w˜ = I1 − I2, whereI1 =1√pi∫ 10(u− v)−1/2v(n−1)/2−1 dv,I2 =cn2pi∫ u1(u− v)−1/2 dv∫ 10(v − s)−3/2sn/2−1 ds.Both integrals can be expressed in terms of hypergeometric functions. ForI1, owing to 3.197 (3) and 9.131 (1) from [6], we obtainI1 =u−1/2√piB(n− 12, 1)F(12,n− 12;n+ 12;1u)=Γ((n− 1)/2)√pi Γ((n+ 1)/2)(u− 1)−1/2 F(1,12;n+ 12;11− u). (23)For I2, changing the order of integration and using [6, 3.238 (3)], we haveI2 =cn (u− 1)1/2pi∫ 10sn/2−1 (1− s)−1/2u− s ds.By [6, 3.228 (3)] this expression coincides with (23).To complete the proof, we recall that w(r) = r3−nw˜(r2), which givesw(r) =1 if 0 ≤ r ≤ 1,−cn r3−n2√pi∫ 10(r2 − v)−3/2vn/2−1 dv if r > 1.This coincides with (12).ON Y. NIEVERGELT’S INVERSION FORMULA . . . 512.2. ExamplesLet us give some examples of functions w(r) defined by (12) in the caser > 1. By [6, 3.197(3)],w(r) = −Γ((n− 1)/2) r−n2√pi Γ(n/2 + 1)F(32,n2;n2+ 1;1r2).Keeping in mind that F (a, b; c; z) = F (b, a; c; z), and using formulas 156,203, and 211 from [19, 7.3.2], we obtain:For n = 2:w(r) = − 12r2F(32, 1; 2;1r2)= 1− r√r2 − 1 . (24)For n = 3:w(r)=− 12r3√pi Γ(5/2)F(32,32;52;1r2)=2pi(arcsin1r− 1√r2−1). (25)For n = 4:w(r) = − 18r4F(32, 2; 3;1r2)= 1− 2r2 − 12r√r2 − 1 . (26)3. The general caseAs in Section 2.1, our main concern is integral equation (17), which isequivalent toIk/20+ w˜ = ψ˜, (27)w˜(u) = u(n−k)/2−1w(√u), ψ˜(u) =2un/2−1cΓ(k/2)ψ(√u). (28)We want to find relatively simple functions w and ψ, which are admissiblein the basic Theorem 2 and such that the corresponding functions w˜ andψ˜ obey (27). It is convenient to consider the cases of k even and k oddseparately.3.1. The case of k evenLet k = 2`; ` = 1, 2, . . . . We chooseψ(r) =0 if 0 ≤ r ≤ 1,(r2 − 1)`rn+2`+1if r > 1.(29)52 E. Ournycheva, B. RubinThe corresponding function x → ψ(|x|) obviously has a decreasing radialmajorant in L1(Rn). By (28), equation (27) is equivalent to(I`1+w˜)(u) =c` (u− 1)`u`+3/2, c` =Γ(n/2− `)Γ(n/2), u > 1. (30)The `th derivativew˜(u) =[c` (u− 1)`u`+3/2](`)= c`∑`j=0(`j)[(u− 1)`](j) [u−`−3/2](`−j) (31)is integrable on (1,∞). The fact that (31) satisfies (30) can be easily checkedusing integration by parts. Thus, the pair of functions w(r) and ψ(r),defined byw(r) =0 if 0 ≤ r ≤ 1,c` r2+2`−n(12rddr)` [(r2 − 1)`r2`+3]if r > 1,(32)and (29), falls into the scope of Theorem 2 and the “even part” of TheoremB is proved.3.2. The case of k oddLet k = 2`+1; ` = 0, 1, 2, . . . . We define ψ(r) by (29), as above. Then,instead of (30), we have(I`+1/21+ w˜)(u) =c′` (u− 1)`u`+3/2, c′` =Γ((n− 1)/2− `)Γ(n/2), u > 1. (33)This gives w˜(u) = g(`+1)(u), whereg(u) =c′`√pi∫ u1(u− v)−1/2(v − 1)`v`+3/2dv = c˜`(u− 1)`+1/2u`+1,c˜` =Γ((n− 1)/2− `) `!Γ(n/2) Γ(`+ 3/2),(use [6, 3.238(3)]). Let us show that w˜(u) ≡ g(`+1)(u) satisfies (33). Inte-grating by parts, we have I`+1/21+ w˜ = I1/21+ g′, whereg′(u) =dg(u)du= c˜`[(`+ 1/2) (u− 1)`−1/2u`+1− (`+ 1) (u− 1)`+1/2u`+2].ON Y. NIEVERGELT’S INVERSION FORMULA . . . 53This gives I`+1/21+ w˜ = c˜` (I1 − I2) where (use [6, 3.238(3)] again)I1 =`+ 1/2√pi∫ u1(u− v)−1/2(v − 1)`−1/2v`+1dv =Γ(`+ 3/2)`!(u− 1)`u`+1/2,I2 =`+ 1√pi∫ u1(u− v)−1/2(v − 1)`+1/2v`+2dv =Γ(`+ 3/2)`!(u− 1)`+1u`+3/2.Hence,I1 − I2 = Γ(`+ 3/2)`!(u− 1)`u`+3/2and therefore,(I`+1/21+ w˜)(u) =c˜` Γ(`+ 3/2)`!(u− 1)`u`+3/2= c′`(u− 1)`u`+3/2,as desired. Thus, functions w(r) and ψ(r), defined byw(r) =0 if 0 ≤ r ≤ 1,c˜` r3+2`−n(12rddr)`+1 [(r2 − 1)`+1/2r2`+2]if r > 1,(34)and (29), obey Theorem 2. This completes the proof of Theorem B.Acknowledgements. We are grateful to Professor Yves Nievergelt,who read our article and paid attention to close papers [13, 14, 26]. Ourapproach differs from that in these works.The second author (B. Rubin) was partially supported by the NSFgrants PFUND-137 (Louisiana Board of Regents) and DMS-0556157.References[1] C.A. Berenstein and B. Rubin, Totally geodesic Radon transform ofLp-functions on real hyperbolic space. Fourier Analysis and Convexity,Series: Applied and Numerical Harmonic Analysis, Ed. L. Brandolini,L. Colzani, A. Iosevich and G. Travaglini, No VIII (2004), 280.[2] L. Ehrenpreis, The Universality of the Radon Transform. Oxford Uni-versity Press, Oxford (2003).54 E. Ournycheva, B. Rubin[3] M. Frazier, B. Jawerth and G. Weiss, Littlewood-Paley theory andthe study of function spaces. CBMS Regional Conference Series inMathematics, vol. 79, AMS, Providence, R.I. (1991).[4] I.M. Gelfand, S.G. Gindikin and M.I. Graev, Selected Topics in In-tegral Geometry. Translations of Mathematical Monographs, v. 220,AMS, Providence, RI (2003).[5] F.B. Gonzalez and T. Kakehi, Invariant differential operators and therange of the matrix Radon transform. J. Funct. Anal. 241 (2006),232-267.[6] I.S. Gradshteyn and I.M. Ryzhik, Table of Integrals, Series and Prod-ucts. Academic Press, NY (1980).[7] S. Helgason, The Radon Transform. Birkha¨user, Boston, Second Ed.(1999).[8] A.V. Lakshminarayanan and G.N. Ramachandran, Three-dimensionalreconstruction from radiographs and electron micrographs: Applica-tion of convolutions instead of Fourier transforms. Proc. Nat. Acad.Sci. U.S.A. 68 (1971) 2236-2240.[9] B.F. Logan and L. Shepp, The Fourier Reconstruction of a Head Sec-tion. IEEE Trans. on Nuclear Science 21, (1974), 21-43.[10] A. Markoe, Analytic Tomography. Encyclopedia of Mathematics andits Applications, vol. 106, Cambridge University Press, Cambridge(2006).[11] F. Natterer, The Mathematics of Computerized Tomography. SIAM,Philadelphia (2001).[12] Y. Nievergelt, Elementary inversion of Radon’s transform. SIAM Rev.28, No 1 (1986), 79-84.[13] Y. Nievergelt, Exact reconstruction filters to invert Radon transformswith finite elements. J. Math. Anal. Appl. 120, No 1 (1986), 288-316.[14] Y. Nievergelt, Elementary inversion of the exponential X-ray trans-form. IEEE Trans. Nuclear Sci. 38, No 2, Part II (1991), 873-876.[15] G. O´lafsson, E. Ournycheva and B. Rubin, Higher rank wavelet trans-forms, ridgelet transforms, and Radon transforms on the space of ma-trices. Appl. Comput. Harmon. Anal. 21 (2006), 182-203.ON Y. NIEVERGELT’S INVERSION FORMULA . . . 55[16] E. Ournycheva and B. Rubin, Method of mean value operators forRadon transforms in the space of matrices. Intern. J. Math. 19 (2008),245-283.[17] E. Ournycheva and B. Rubin, Semyanistyi’s integrals and Radon trans-forms on matrix spaces. The J. of Fourier Analysis and Appl. 14, No1 (2008), 60-88.[18] V. Palamodov, Reconstructive Integral Geometry. Monographs inMathematics, vol. 98, Birkha¨user Verlag, Basel (2004).[19] A.P. Prudnikov, Y.A. Brychkov and O.I. Marichev, Integrals and Se-ries, Vol. 3: More Special Functions. Gordon and Breach, Amsterdam(1990).[20] E.T. Quinto, M. Cheney and P. Kuchment (Ed-s), Tomography,Impedance Imaging, and Integral Geometry. Lectures in Appl. Math.,30, AMS, Providence, RI (1994).[21] A.G. Ramm and A.I. Katsevich, The Radon Transform and Local To-mography. CRC Press, Boca Raton (1996).[22] B. Rubin, Inversion formulas for the spherical Radon transform andthe generalized cosine transform. Advances in Appl. Math. 29 (2002),471-497.[23] B. Rubin, Convolution-backprojection method for the k-plane trans-form, and Caldero´n’s identity for ridgelet transforms. Appl. Comput.Harmon. Anal. 16 (2004), 231-242.[24] B. Rubin, Reconstruction of functions from their integrals over k-planes. Israel J. of Math. 141 (2004), 93117.[25] S.G. Samko, A.A. Kilbas and O.I. Marichev, Fractional Integrals andDerivatives. Theory and Applications. Gordon and Breach Sc. Publ.,New York (1993).[26] D. Schuster, Der Faltungsalgorithmus fu¨r dreidimensionale Radon-Probleme. Numer. Funct. Anal. Optim., 10 (1989), 1203-1212.[27] D.C. Solmon, A note on k-plane integral transforms. Journal of Math.Anal. and Appl. 71 (1979), 351-358.[28] E.M. Stein, Singular Integrals and Differentiability Properties of Func-tions. Princeton Univ. Press, Princeton, NJ (1970).56 E. Ournycheva, B. Rubin[29] R.S. Strichartz, Lp-estimates for Radon transforms in Euclidean andnon-euclidean spaces. Duke Math. J. 48 (1981), 699-727.∗ University of Pittsburgh at Bradford300, Campus Drive16701, Bradford, PA, USAe-mail: elo10@pitt.edu∗∗ Department of MathematicsLouisiana State UniversityBaton Rouge, LA, 70803, USA Received: December 9, 2009e-mail: borisr@math.lsu.edu

On Y. Nievergelt's inversion formula for the Radon transform

Abstract

Similar works

Full text

Available Versions

Bulgarian Digital Mathematics Library at IMI-BAS