Fast Backprojection Techniques for High Resolution Tomography
FFast Backprojection Techniques for HighResolution Tomography
Nikolay Koshev, Elias S. Helou, Eduardo X. MiquelesOctober 15, 2018
Abstract
Fast image reconstruction techniques are becoming important with theincreasing number of scientific cases in high resolution micro and nanotomography. The processing of the large scale three-dimensional data de-mands new mathematical tools for the tomographic reconstruction taskbecause of the big computational complexity of most current algorithms asthe sizes of tomographic data grow with the development of more power-ful acquisition hardware and more refined scientific needs. In the presentpaper we propose a new fast back-projection operator for the process-ing of tomographic data and compare it against other fast reconstructiontechniques.
Tomographic imaging is a very powerful instrument of non-destructive researchand control of the internal structure of non-opaque objects. An importantbranch of tomographic techniques is transmission tomography , which can beused at nano, micro and macro resolution levels. For further consideration wedescribe in general the basic principles of transmission tomography from parallelrays, and define some notations.Physically, all types of transmission tomography are based on registeringthe energy loss or/and intensity loss of the incoming electromagnetic wave ( x -rays for instance), after passing through the object under investigation alsoreferred to here as sample ). In our case, we consider that x -rays generatedfrom a synchrotron light source hit the object under investigation determininga projection image (also referred as frame ) at a ccd (charge coupled device)camera. A typical dataset is shown in Figure 1.A, where a high-resolution frame P gathered using the x -rays source is shown, with dimensions 2048 × sinogram , and that will be used as input to an appropriate inversion algorithmin order to reconstruct the slice of the sample.1 a r X i v : . [ m a t h . NA ] A ug igure 1: (A) Projection (or frame ) for a cylindrical sample obtained with a ccd camera (B) Three-dimensional representation of the measured data: P isthe measured frame and G is the sinogram image at a given row of the areadetector.We introduce the cartesian coordinate system in the plane of a given slice ofthe object. Let the function f ( x ) ∈ U be the feature function , i.e., a functionwhich depends on the internal structure of the object in the plane of the sliceand which defines the linear absorption coefficient of the sample. Set U , referredto here as the feature space , is a Schwartz space S ( R ).A given frame (see Figure 1.A) represents the integral of f ( x ) over straightlines passing through the sample and perpendicular to the detector’s plane. Onerow of each of these frame images contains the integrals relevant to a slice ofthe object, and orthogonal to the rotation axis. Let us introduce an axis t overthe detector’s row. It is clear that for each angle θ (see Figure 2), such a row ismathematically determined by g ( θ, t ) ≡ g θ ( t ) = (cid:90) L ( θ,t ) f ( x )d s = (cid:90) R f ( x ) δ ( x · ξ θ − t )d x , (1.1)where L ( θ, t ) is a straight line defining the x-ray path, L ( θ, t ) = (cid:8) x ∈ R : x · ξ θ = t (cid:9) , ξ θ = (cos θ, sin θ ) T . (1.2)From (1.1) we have a linear operator acting on the feature function f , i.e., R : f ∈ U (cid:55)→ g ∈ V , which is called the Radon transform . Space V is theSchwartz space S ( R + × [0 , π ]). The operator B : V → U defined as b ( x ) = B g ( x ) = (cid:90) [0 ,π ] g ( x · ξ θ , θ )d θ, (1.3)2igure 2: Geometry of incoming x-rays for parallel tomography.is defined as the backprojection operator, and is the adjoint of R in the followingsense (cid:90) R + × [0 ,π ] R f ( t, θ ) g ( t, θ )d t d θ = (cid:90) R f ( x ) B g ( x )d x , (1.4)More about the theory of the integral operators { R , B } can be found on [1, 2,3, 4].At this point, it is convenient to introduce some notations. We first introducethe notations for the representation of feature function f : R → R in differentcoordinate systems, and their respective jacobians: (a) Pr¨ufer coordinates (see [5]) : x = p ( µ ) ξ θ , d x = | p (cid:48) ( µ ) p ( µ ) | d µ d θ . Therepresentation is denoted by [ f ] Pr ( µ, θ ). Function p will always be welldefined within the context by special notation as follows; (b) Log-polar coordinates : particular case of Pr¨ufer coordinates when p ( µ ) = e µ . Here, d x = e µ d µ d θ . The representation is denoted by [ f ] L ( µ, θ ); (c) Semi-polar coordinates : particular case of Pr¨ufer coordinates when p ( µ ) = µ . Here, d x = µ d µ d θ . The representation is denoted by [ f ] P ( µ, θ ). (d) Sinogram coordinates are similar to the semi-polar coordinates and, in fact,can be obtained by flipping the angles θ ∈ [ π, π ) to the negative part ofthe t -axis , so that t ∈ [ − , g in (1.1) can be written as [ g ] P ( t, θ ) in orderto indicate semi-polar coordinates. An example, using the well-known Shepp-Logan phantom [6] is presented on Figure 3. The sinogram of the Shepp-loganfeature function f is presented in the sinogram coordinate system mentionedabove. 3a) (b) | x | ≤ | x | ≤ θ ∈ [0 , π ] t ∈ [ − √ , √ ] (c) (d) θ ∈ [0 , π ] t ∈ [ , √ ] θ ∈ [0 , π ] ρ ∈ ( − ∞ , ] Figure 3: (a) Shepp-Logan feature function f ( x ) and his associated sinograms, indifferent coordinate systems: (b) Semi-polar coordinates [ g ] s , (c) Polar coordinates[ g ] P and (d) Log-polar coordinates [ g ] L . b ∈ U for a given sinogram g ∈ V . The computationof the tomographic image from sinogram data depends on the backprojectionoperator B , which bears the major computational cost of reconstruction meth-ods: O ( N ) for images with N pixels and N projections.For a high-resolution tomographic synchrotron experiment, the amount ofdata at a micro-tomography setup is considerably large for today’s computa-tional standards, mainly because of this asymptotic floating point operations(flops) count. Indeed, at the Brazilian National Synchrotron Light Source ( lnls )one wishes to obtain 2048 reconstructions images with 2048 × × B represent the main bottleneck of the reconstruction process. If certainuseful mathematical properties of B are exploited, the computational effort canbe significantly reduced to O ( N log N ) flops [7, 8].Several techniques were developed aiming at a reduction to O ( N log N ) flopsfor computing B . One approach was established in [9], where the computationof B g is performed after a change from cartesian to log-polar coordinates in thedata. This approach leads to a convolution, which is computable through FastFourier Transform ( fft ) algorithms. Although elegant, the methods suffer fromthe ill-conditioning of the Log-polar transform at the “fovea”. Nevertheless, itis possible to translate the fovea to different regions of the cartesian plane,in order to enclose the reconstruction region. This leads to the concept of partial-backprojection which can be easily implemented in a parallel form. Othermethods for fast computation of B were presented in [10, 11, 12, 13], using adivide and conquer strategy based on hierarchical decompositions of the fullbackprojection which are simpler than the full backprojection. Hierarchicaldecompositions can be created both in the image [12] or in the data [10, 11].Yet another approach is based on Non-uniform Fast Fourier Transform ( nfft )algorithms (see [12, 14, 15]) and the so-called Fourier Slice Theorem ([2, 16].In this paper, we propose another fast method for the computation of B g , alsobased on Fourier transforms. We claim that the backprojection of g ∈ V canbe easily done by filtering the lines of the ˜ g one by one, where ˜ g is the polarrepresentation of g in S + = R × [0 , π ].This manuscript is organized as follows: Section 2 presents a discussionof low-complexity algorithms for the computation of B . Our low-complexityformula is presented in Section 3 and a discussion of the implementation ispresented in Section 5. Further comparison of all algorithms is presented inSection 7 and a discussion of the results is shown in Section 8. Remark:
In this manuscript we use the the integral operator, sometimes withd x placed before the integrand, as it is more convenient to make explicit thevariables being considered. Whenever the integrand is short, we adopt theclassic notation (cid:82) f ( x )d x . 5 Class of Algorithms for the Backprojection
Let g ∈ V be a given sinogram. A na¨ıve implementation of the typical backpro-jection formula (1.4) has to be done using nested loops. Indeed, for each pixel x lying on a predefined meshgrid within the square (cid:107) x (cid:107) ∞ ≤
1, the approximationof b ( x ) = B g ( x ) is given by b ( x ) ≈ ∆ θ N θ (cid:88) k =1 g ( x · ξ θ k , θ k ) . (2.1)It is easy to realize that the above approximation has a computational costof O ( N θ ) for each pixel x , where N θ is the total number of sampled angles.For a high resolution frame (see Figure 1.A), a linear interpolation for x · ξ θ on the grid of − ≤ t ≤ b isrepresented by a square image of order N × N , the total cost for computing thefinal backprojected image b is O ( N N θ ). In practice, N θ has almost the samemagnitude of N , and thus we can state that the asymptotic cost to obtain b is O ( N ). Such an algorithm is impractical for high-resolution images.There are at least three other types of backprojection algorithms which candramatically reduce the computing time of the backprojected image b , for largedatasets: (i) A fast slant-stack based approach [17] was proposed by Averbuch et al .Although this is an elegant and fast approach, it will not be covered inthis manuscript; (ii)
Hierarchical decomposition [10, 12, 11]: Two different approaches that ap-ply the divide-and conquer paradigm to the backprojection computation,splitting it recursively into smaller and simpler subproblems; (iii)
NFFT [18]: The Fourier Slice Theorem sets the Fourier Transform asbridge between the Radon Transform R f ( θ, t ) and the original image f .However, tomographic data does not provide an evenly distributed sam-pling of the Fourier space, as required by traditional fft techniques (see[19]). Use of this Fourier approach was enabled by research on nfft al-gorithms (see [12, 14, 15]); (iv) Anderson’s formula [9]: Such a formula is obtained with an appropriatechange of variables on the classical equation of the backprojection formula(1.4). The main idea is to convolve the sinogram in log-polar coordinateswith an ideal kernel using fft algorithms.In this paper, we focus mainly on the description of our algorithm and algorithm(iii).
A fast method for obtaining the backprojection was derived by Andersson [9].His approach is based on a representation of the Backprojection/Radon trans-6orm as a convolution, by casting the computation in a log-polar coordinatesystem. In this section we propose a different proof for his formula.Let g = R f be a given sinogram, i.e., the Radon transform of a compactlysupported function f . Using the coordinate system notation of the previoussection, where [ · ] L denotes the log-polar representation of some function, themain formula of the log-polar backprojection is written as[ B g ] L ( ρ, θ ) = [ g ] L ∗ [ K ] L ( ρ, θ ) , (2.2)where ∗ stands for the two-dimensional convolution, and K is the convolutionKernel [ K ] L ( ρ, θ ) = δ (1 − e ρ cos θ ) . (2.3)Using above formula and the convolution theorem, we obtain[ B g ] L = F − (cid:0) F [ g ] L · F [ K ] L (cid:1) . (2.4)Let us give a simple proof of the above equation, assuming that f lies in aSchwartz space S ( R ), and g ∈ S ( R + × [ − π, π ]). Proof:
We start with the integral representation of the backprojection operator,given in (A.7) (See Appendix A). Now, formula (2.2) is derived in four steps: (a)
Changing the integral (A.7) from cartesian coordinates y ∈ R to Pr¨ufercoordinates, i.e., y ≡ y µ,θ = p ( µ ) ξ φ we get d y = | p (cid:48) ( µ ) p ( µ ) | d µ d θ and B g ( x ρ,θ ) = (cid:90) S + g ( y µ,φ ) δ (cid:0) κ x ρ,θ ( y µ,φ ) (cid:1) | p (cid:48) ( µ ) p ( µ ) | d µ d φ (2.5) (b) The support of the Delta distribution in (2.5) is κ x ρ,θ ( y µ,φ ) = p ( µ ) (cid:20) − p ( ρ ) p ( µ ) ξ φ · ξ θ (cid:21) = p ( µ ) (cid:20) − p ( ρ ) p ( µ ) cos( φ − θ ) (cid:21) (2.6) (c) Let [ · ] Pr be the representation in Pr¨ufer coordinates. From (2.6) and (2.5)we arrive at[ B g ] Pr ( ρ, θ )= (cid:90) S + [ g ] G ( µ, φ ) δ (cid:18) p ( µ ) (cid:20) − p ( ρ ) p ( µ ) cos( φ − θ ) (cid:21)(cid:19) | p (cid:48) ( µ ) p ( µ ) | d µ d φ = (cid:90) S + [ g ] G ( µ, φ ) δ (cid:18) − p ( ρ ) p ( µ ) cos( φ − θ ) (cid:19) | p (cid:48) ( µ ) p ( µ ) | p ( µ ) d µ d φ (2.7)where S + = R + × [ − π, π ] 7 d) A convolution is obtained in (2.7) only if p is such that p ( ρ ) = p ( µ ) p ( ρ − µ ),which in turn implies that p is an exponential function. Hence, Pr¨ufer co-ordinates reduce to log-polar coordinates, which we denote by [ · ] L . Finally,we obtain[ B g ] L ( ρ, θ ) = (cid:90) S + [ g ] L ( µ, φ ) δ (cid:0) − e ρ − µ cos( φ − θ ) (cid:1) d µ d φ (2.8)which is the final convolution formula. Although Anderson’s approach is asymptotically fast, it has a few drawbacks.Firstly, the gain of speed using Fourier transforms to compute the convolu-tion is reduced with forward/backward log-polar transformations. Also, theseinterpolations can produce errors, especially near the origin, due to a strongnon-uniformity of the log-polar mesh in that region. To avoid these factors, an-other approach for the calculation of the backprojection operator can be used.This approach is based on the following theorem:
Theorem ( Backprojection Slice Theorem (BST)).
Let g = g ( t, θ ) ∈ V agiven sinogram and (cid:98) · denotes the Fourier transform operation. It follows thatthe backprojection B satisfies (cid:100) B g ( σ ξ θ ) = ˆ g ( σ, θ ) σ (3.1) with σ > ∈ R and θ ∈ [0 , π ] .Proof. Using the sifting property of the δ -distribution, the backprojection (1.3)can be presented in the following form B g ( x ) = (cid:90) π g ( x · ξ θ , θ )d θ = (cid:90) π (cid:90) R g ( t, θ ) δ ( t − x · ξ θ )d t d θ Considering the two-dimensional Fourier transform of B g , i.e., F : B g (cid:55)→ (cid:100) B g and using representation (A.7) (see Appendix), (cid:100) B g ( ω ) = (cid:90) R B g ( x ) e − i ω · x d x = (cid:90) R d x (cid:90) R d y [ g ] c ( y ) δ ( y · ( y − x )) e − i ω · x = (cid:90) R d y [ g ] c ( y ) (cid:90) R d x δ ( y · ( y − x )) e − i ω · x ≡ (cid:90) R d y [ g ] c ( y ) T ( y , ω )where y , ω ∈ R and T ( y , ω ) = (cid:90) R d x δ ( h y ( x )) e − i ω · x , h y ( x ) = y · ( y − x ) (3.2)8ince the distribution (3.2) is supported in the set h − y (0) = { x ∈ R : h y ( x ) =0 } , it follows from (A.1) (See Appendix A) and ∇ h y = − y that T ( y , ω ) = 1 (cid:107) y (cid:107) (cid:90) h − y (0) e − i ω · x d s ( x ) = (cid:90) h − y (0) e − i ω · x ( q ) d q (3.3)The set h − y (0) determines a straight line passing through y and with normalvector y . Thus, h − y (0) = y + span { S y } , being S y ⊥ y and S a π -rotationmatrix. Therefore, x ( q ) ∈ h − x (0) is on the form x ( q ) = y + qS y and theintegral in (3.3) can be written as: T ( y , ω ) = (cid:90) R e − i ω · [ y + qS y ] d q = e − i ω · y (cid:90) R e − iq ω · ( S y ) d q = e − i ω · y δ ( ω · S y )(3.4)Hence, the Fourier transform of B g becomes (cid:100) B g ( ω ) = (cid:90) R [ g ] c ( y ) δ ( ω · S y ) e − i ω · y d y (3.5)For ω fixed, { y ∈ R : ω · ( S y ) = 0 } = span { ω } , with S a π -rotation matrix.Indeed, since S y ⊥ w and S y ⊥ y , it follows ω (cid:107) y . Once again, using therepresentation (A.1) (see Appendix A) for (3.5) we arrive at (cid:100) B g ( ω ) = (cid:90) R [ g ] c ( q ω ) (cid:107) S ω (cid:107) e − i ω · ( q ω ) d s ( ω ) (3.6)Since (cid:107) S ω (cid:107) = (cid:107) ω (cid:107) and d s ( ω ) = (cid:107) ω (cid:107) dq , we finally obtain (cid:100) B g ( ω ) = (cid:90) R [ g ] c ( q ω ) e − iq (cid:107) ω (cid:107) d q (3.7)From the above equation, the backprojection is a polar convolution. Indeed,switching the frequency domain to polar coordinates, i.e., ω = σ ξ θ (with σ ∈ R + and θ ∈ [0 , π ]) we get (cid:100) B g ( σ ξ θ ) = (cid:90) R [ g ] c ( qσ ξ θ ) e − iq (cid:107) σ ξ θ (cid:107) d q = (cid:90) R [ g ] c ( u ξ θ ) σ e − iuσ d u. (3.8)Now, letting [ · ] s be the representation in semi-polar coordinates, it is true that[ g ] c ( u ξ θ ) = g ( u, θ ) is the input sinogram g ( u, θ ). From (3.7) and (3.8), usingpolar coordinates[ (cid:100) B g ] p ( σ, θ ) = (cid:100) B g ( σ ξ θ ) = 1 σ (cid:90) R g ( u, θ ) e − iuσ d u (3.9)Identity (3.9) is our backprojection-slice Theorem (3.1) for computing the op-erator B . 9ndeed, at each radial line θ in the frequency domain, the two-dimensionalFourier transform of B equals the one-dimensional radial Fourier transform ofthe projection g ( t, θ ) multiplied by the kernel 1 /σ for σ > Remark 1:
The mathematical proof outlined above provides a direct formulafor the computation of a backprojected image, i.e., given a sinogram g , theexplicit steps to compute the backprojection in the frequency polar coordinatesresults in formula (3.1). In practice, there are several iterative methods thatdepends explicitly on the computation of the backprojection of any sinogram.In the other hand, analytical formulas usually handle with the backprojection ofa filtered sinogram, from where standard formulas like the filtered backprojection or the filter of the backprojection are established. To validate our backprojectionresult we remark the following items: (i) It is a well known fact [1, 4, 2] that, for a given feature function f ∈ U , thefollowing property holds BR f ( x ) = ( f ∗ h )( x ) , h ( x ) = 1 (cid:107) x (cid:107) (3.10)which, in the frequency domain, is written as (cartesian and polar representa-tion, respectively) (cid:92) BR f ( w ) = ˆ f ( w ) 1 (cid:107) w (cid:107) ⇔ (cid:92) BR f ( σξ θ ) = ˆ f ( σξ θ ) 1 σ (3.11)due to the fact that F : (cid:107) x (cid:107) (cid:55)→ (cid:107) w (cid:107) . Now, replacing the backprojection slicetheorem (3.1) into (3.11), we obtain (cid:100) R f ( σ ξ θ ) 1 σ = ˆ f ( σξ θ ) 1 σ ⇒ (cid:100) R f ( σ ξ θ ) = ˆ f ( σξ θ ) (3.12)which is the celebrated Fourier Slice-Theorem [3]. (ii) From the classical inversion of the Radon transform, i.e., the filtered-backprojection algorithm, it is true that B F g ( x ) = f ( x ) , g = R f (3.13)where F is a low-pass filtering operator, that is (cid:99) F g ( ν, θ ) = ˆ g ( ν, θ ) | ν | , for ν ∈ R .In the polar frequency domain, (3.13) reads (cid:91) B F g ( σ ξ θ ) = ˆ f ( σ ξ θ ). From thebackprojection slice Theorem (3.1), such equation becomes1 σ (cid:99) F g ( σ ξ θ ) = ˆ f ( σ ξ θ ) ⇒ σ ˆ g ( σ ξ θ ) σ = ˆ f ( σ ξ θ ) , σ ∈ R + (3.14)Once again, the above equation yields the Fourier Slice-Theorem.10 emark 2: The dc -component of the Backprojection of some function g lyingin the sinogram space is defined by (cid:100) B g ( ) = (cid:90) R B g ( x )d x (3.15)= (cid:90) R d x (cid:90) π d θ g ( x · ξ θ , θ ) (3.16)= (cid:90) π d θ (cid:90) R d t (cid:90) R d s g ( t, θ ) ≡ M (3.17)where we have used d x = d t d s to make explicit the change of variables from x to( t, s ), being s the variable along the direction ξ ⊥ θ . The dc of an arbitrary g ∈ V provide M = ∞ . In this sense, (cid:100) B g behaves like a tempered distribution since B g lies in a Schwartz space, where the Fourier transform is an automorphism.Also, it is easy to note thatˆ g ( σ, θ ) σ = i ˆ h ( σ, θ ) , h ( t ) = (cid:90) t ∞ g ( t, θ )d t (3.18)i.e., h is a primitive of g . Hence, using (3.18) as σ →
0, the limit of the ratioˆ g ( σ, θ ) /σ diverge in σ = 0. Finally, bst formula can be easily applied for some g ∈ V with a nonzero dc -component. In fact, setting p ( t, θ ) = g ( t, θ ) − ˆ g (0 , θ ),it is true that ˆ p (0 , θ ) = 0 and the backprojection of g follows with B g ( x ) = B p ( x ) + ˆ g (0 , θ ). In this section, we provide two examples where the analytical computation ofthe backprojection operator is possible. This is important to validate furthernumerical simulations. The first example given is for a point source function,both in log-polar and cartesian coordinates. The second one, for a symmetricalcircular function. In what follows, we consider that g = g ( t, θ ) is the Radontransform of f = f ( x ), while b = b ( x ) is the final backprojected image. Example I:
Following Andersson’s formula , the backprojection of any sino-gram g is written as a convolution in Log-polar coordinates b ( e ρ ξ θ ) = B g ( e ρ ξ θ ) = (cid:90) d u (cid:90) d β g ( e u , β ) δ (1 − e ρ − u cos( θ − β )) (4.1)It is a well known fact that the Radon transform of a single point source, locatedat x = a is f ( x ) = δ ( x − a ) ⇒ g ( t, θ ) = δ ( t − a · ξ θ ) (4.2) Even if g is the sinogram of a compactly supported function feature function on the unitdisk (cid:107) x (cid:107) ≤
1, we have ˆ g (0 , θ ) = constant, although with M = ∞ . a = e A ξ φ as the log-polar representation of the source point a , we use(4.2) and (4.1) to obtain b as b ( e ρ ξ θ ) = (cid:90) d β (cid:90) d u δ ( a · ξ β − e u ) δ (1 − e ρ − u cos( θ − β )) (4.3)= 1 (cid:112) (cos( θ − φ ) e ρ ) + ( e A − sin( θ − φ ) e ρ ) (4.4)where A = ln (cid:107) a (cid:107) . The details of the log-polar representation (4.4) are presentedin the Appendix B. To obtain a cartesian representation we use x = e ρ ξ θ andcos( θ − φ ) = ξ θ · ξ φ , sin( θ − φ ) = ξ θ · ξ ⊥ φ (4.5)Now, (4.4) becomes b ( x ) = 1 (cid:113) [( ξ θ · ξ φ ) e ρ ] + [ e A − ( ξ θ · ξ ⊥ φ ) e ρ ] (4.6)= 1 (cid:113) [( e ρ ξ θ ) · ( e A ξ φ ) e − A ] + [ e A − ( e ρ ξ θ ) · ( e A ξ ⊥ φ ) e − A ] (4.7)= 1 (cid:112) [( x · a ) e − A ] + [ e A − ( x · ˆ a ) e − A ] (4.8)where ˆ a = e A ξ ⊥ φ is a counterclockwise rotation by π of the point source a .Finally, since e A = (cid:107) a (cid:107) we obtain b ( x ) = (cid:107) a (cid:107) (cid:112) [ x · a ] + [ (cid:107) a (cid:107) − x · ˆ a ] (4.9)which is the cartesian representation of the backprojection of the sinogram g given in (4.2). Example II:
Considering f ( x ) = circ( (cid:107) x (cid:107) ), i.e., f ( x ) = , (cid:107) x (cid:107) ≤ , (cid:107) x (cid:107) > , (cid:107) x (cid:107) = 1 (4.10)it is known [20] that ˆ f ( ω ) = J ( (cid:107) ω (cid:107) ) (cid:107) ω (cid:107) (4.11)with J the order 1 Bessel function of the first kind. From the Fourier-Slice-Theorem ˆ f ( σ ξ θ ) = ˆ g ( σ, θ ) and (4.11) it is easy to obtain g satisfyingˆ g ( σ, θ ) = J ( (cid:107) σ ξ θ (cid:107) ) (cid:107) σ ξ θ (cid:107) = J ( σ ) σ bst formula (3.1), the Fourier representation for b becomesˆ b ( σ ξ θ ) = ˆ g ( σ, θ ) σ = J ( σ ) σ ⇐⇒ ˆ b ( ω ) = J ( ω ) (cid:107) ω (cid:107) (4.12)The above equation provides a testing algorithm (in the Fourier domain) forour numerical strategies. In fact, either bst or Andersson’s algorithm presentsdifficulties as ω → , as discussed in the next section. All the notation needed for the implementation of the bst formula is presentedin Table 1. In this section we assume to have the sinogram presented at thenodes of the uniform sinogram grid G s . Quantification of the sinogram function g ( θ, t ) over G s will be denoted by g ij . Meshsequence Descriptionof coordinates MeshSize Stepsize { t k } = {− . , . . . , . } sinogram | t k | = N t ∆ t = 2 /N t { s k } = { , . . . , . } polar | s k | = N s = N t / s = 2 /N t { ρ k } = { ρ , . . . , ln 1 } log-polar | ρ k | = N ρ ∆ ρ = ( − ρ + ln 1) /N ρ { θ k } = { , . . . , π } angles within [0 , π ] | θ k | = N θ ∆ θ = π/N θ { φ k } = { , . . . , π } angles within [0 , π ] | φ k | = 2 N θ ∆ φ = π/N θ G s = { θ k } × { t i }
2D sinogram mesh | G s | = ( N t , N θ ) - G P = { φ k } × { s i }
2D polar mesh | G P | = ( N t / , N θ ) - G L = { φ k } × { ρ i }
2D log-polar mesh | G L | = ( N ρ , N θ ) - Table 1: Glossary of symbols used for implementation details.
In the implementation of the simplest algorithm, based on the Andersson’s ap-proach, we don’t take into account the irregularity in the origin. The algorithmcan be summarized as follows:
Step 1.
Interpolate the sinogram to log-polar coordinates, [ g ] s → [ g ] L : Thedomain of the experimentally obtained sinogram is S + = [0 , π ] × [ − , polar grid G P (see Table 1). The change from semi-polarto log-polar coordinates can be easily done using the linear interpolationalong the ray for every θ = Const . Let us denote the number of pointsalong every ray as N ρ . Since ln( (cid:15) ) → −∞ as (cid:15) → ρ < G s to the log-polar grid G L (see Table 1). In fact,13or all change of coordinates computed in this work, we use simple linearinterpolation. The problem of selection of the first node ρ of log-polarmesh considered below, in Sec. 5.2.One of the specific features of Log-Polar mesh is its non-uniformity. Toobtain clear interpolation without losing information, we need to makethe biggest step of log-polar mesh equal to radial step ∆ s of the originalsinogram. It can be easily done by finding the mesh size N ρ from theequation: exp(0) − exp( − ∆ ρ ) = ∆ s (5.1)The above equation comes from the definition of the mesh G L since { ρ k , ρ k − } is in fact { , − ∆ ρ } with k = N ρ −
1. Since ∆ ρ = ρ /N ρ we finally obtainthe number of points in the log-polar system N ρ = − ρ ln(1 − ∆ s ) . (5.2)In the next Section we consider a different choice for the parameter ρ . Step 2.
Calculate the kernel K by formula (2.3) : The kernel represented withthe formula (2.3): [ K ] L ( ρ, θ ) = δ (1 − e ρ cos θ ) and can be approximatedon the mesh G L using the condition: k ji = (cid:40) / ∆ ρ, if | e ρ i cos θ j − | ≤ ∆ ρ, , otherwise . (5.3)Here ∆ ρ = − ρ N ρ is the step on the mesh for variable ρ . The kernelin log-polar coordinates and the absolute value of its Fourier transformare shown on Figure 4. The approximation (5.3) could be numericallyimproved using appropriate strategies for the evaluation of a Delta distri-bution concentrated on the zero level set of a function, see [21]. Step 3.
Calculate the convolution [ K ] L (cid:63) [ g ] L : Using the uniform grids, theconvolution is calculated using the Fast Fourier Transform pair through fftw3 software library. Step 4.
Interpolate the result from previous step back to Cartesian coordinatesystem . We are using bilinear interpolation, which already has been de-scribed in Step 1.
Due to the log-polar representation, s → ρ → −∞ . Of course, inreal calculations it is not possible to get a proper interpolation to this grid.This practical problem can be solved in two ways. The first approach is clearlymathematical, and was proposed by Andersson in his work [9], called partialbackprojection method . This method is based on moving the origin outside theregion of interest, which allows us to make a clear interpolation of all points of14a) (b) θ ∈ [ − π, π ] ρ ∈ [ ρ , ] K ( ρ, θ )Figure 4: Kernel in log-polar representation (a) and its Fourier image (b) the sinogram with non-zero values. The second approach is to select a proper ρ , which adapts properly to the resulting cartesian grid. As we will show, thisway also gives good results. Adaptive selection of ρ . For clear interpolation we have to use ρ as a verybig negative number, from which we start the approximation to the log-polarmesh. But, in fact, this number is connected to the mesh which is chosen forcartesian representation of the result. Assume the cartesian mesh with the direc-tion steps (∆ x, ∆ y ). In this case, to avoid the loss of information near the origin,we have to set ρ < ln(min(∆ x, ∆ y )). The result of using the Anderson’s withthe rightly chosen ρ is presented Fig.6.(a) in the “Numerical results” section. Infact, it is important to note that the performance of Log-polar backprojectiondepends on the desired resolution of the resulting image. The oversampling oflog-polar mesh grows up fast as the number of pixel increases in the cartesiangrid, as it shown below in this section.Now let g ( t, θ ) be some sinogram and [ g ] L ( ρ, θ ) his log-polar representationfor ρ ∈ ( −∞ , g ] L with the following com-pactly supported function: g − L ( ρ, θ ) = (cid:40) [ g ] L ( ρ, θ ) , ρ ∈ [ ρ , , , ρ < ρ (5.4)where ρ (cid:28) B g and B g − , i.e., (cid:107) [ B g ] L − [ B g − ] L (cid:107) L = π (cid:82) (cid:82) −∞ (cid:0) [ B g ] L ( ρ, θ ) − [ B g − ] L ( ρ, θ ) (cid:1) e ρ d ρ d θ = π (cid:90) (cid:90) ρ + π (cid:90) ρ (cid:90) −∞ (cid:0) [ B g ] L ( ρ, θ ) − [ B g − ] L ( ρ, θ ) (cid:1) e ρ d ρ d θ (5.5)From (5.4) it follows that (cid:107) [ B g ] L − [ B g − ] L (cid:107) L = π (cid:90) ρ (cid:90) −∞ [ B g ] L ( ρ, θ ) e ρ d ρ d θ ≤ π c e ρ , (5.6)where c = max ρ ≤ ρ [ B g ] L is a constant, which, in practice, refers to the value ofBackprojection in the origin. This value can be easily estimated. Equation (5.6)give us a bound for the error, when we remove the origin in the computation ofthe backprojected image.To obtain a good reconstruction, it is easy to obtain the number N ρ fromusing the following formula N ρ ≈ ln(min(1 /N x , /N y )ln(1 − ∆ s ) . (5.7)For example, assuming that N s = 1024 and N x = N y = 1024, then N ρ = 3546.Therefore, an oversampling of the input data is usually needed (about ≈ Partial Backprojections.
The Partial Backprojections method is based onthe shifting property of the Radon transform, defined as u ( x ) = f ( x − ∆ x ) ⇒ R u ( t, θ ) = g ( t − ξ θ · ∆ x , θ ) , g = R f (5.8)Using this formula we can transform the original sinogram, presented in semi-polar coordinates and take the part, which is located far away from zero. Wechoose some angle β and consider the sector of the original sinogram θ ∈ [ θ , θ ],where θ = θ + 2 β . Rescaling the original sinogram to the size a r , as it shownon Fig.5.(a) and rotating it in the way that sector under investigation will belocated in θ ∈ [ − β, β ]. Now, it is easy to obtain the values of the distancesbetween old and new origins (1 − a r ) and between the sector under investigationand new origin 1 − a r , which defines the minimal s and ρ for the interpolationfrom semi-polar to log-polar coordinates. This method is described in details in[9]. 16a) (b) x x θ ∈ [ − π, π ] ρ ∈ [ ρ , ] Figure 5: (a) The scheme of the reconstruction with sectoral method; (b) transformedsinogram and sector outside of zero (inside red lines)
Using partial backprojection is convenient because it is possible to highlightany sector of interest of the sinogram without any information loss. Also, itis not necessary to process all the sinogram at once - we can process only theparts that we are interested in. The first disadvantage of this method is thatthe Fourier image of the kernel is singular at the origin (see Fig.4.(b)), whichcauses artifacts on the result of sector backprojection. We note that, mathe-matically, this problem also exists in the previous method (adaptive choosingthe ρ ), but, since we exclude the origin from the calculations, we just ”skip”this singularity. The second disadvantage of the partial method is that furthermathematical operations are needed, e.g., translation of the origin in differentcoordinate systems (log-polar and Cartesian) and an extension of angles of thesectors under reconstructions to avoid lost of information at the boundaries ofeach sector. Also, application of partial backprojections for whole object canincrease the calculation time. However, partial backprojection algorithms doesnot need a large oversampling to obtain good resolution at a small region ofinterest and far from the origin. The second algorithm we consider in this paper is based on the BackprojectionSlice Theorem. Due to the fact that the log-polar coordinates interpolation isnot needed for this reconstruction, the algorithm is simpler than the above one.Also we note that straight usage of the Fourier transform may produce ratherbig artifacts near the origin (boundary effect), caused by the fact that the valueson sinogram on the line s = 0 are not equal to 0. This problem can be solved17ith usage of short-time Fourier transform:ˆ f ( t, σ ) = + ∞ (cid:90) −∞ f ( x ) w ( x − t ) e iσx dx (5.9)with t = 0. In this work as a function w we are using the Kaiser-Bessel window[22], which can be defined on the mesh G s - see (1) - using the following formula: w ( β ) = (cid:12)(cid:12)(cid:12) I (cid:16) β (cid:113) − (cid:0) i − N s +1 N s − (cid:1) (cid:17)(cid:12)(cid:12)(cid:12) | I ( β ) | , (5.10)where I ( · ) is a modified Bessel function of the zeroth order.In this case, the sequence of steps to retrieve the backprojected image follows: Step 1.
Transform the source image from sinogram coordinates to semi-polar:[ g ] s → [ g ] P . Step 2.
For each constant θ :2.1. Multiply the s − axis (polar domain, s ≥
0) of a sinogram withthe window function w ( s ), defined by (5.10): [ (cid:101) g ] P ( s, θ ) = [ g ] P ( s, θ ) · w ( s )2.2. Derive for the obtained s − axis;2.3. Multiply the obtained sinogram with the kernel K σ = σ in thefrequency domain or its approximation (to avoid division by zero). Step 3.
Interpolate the resulting image to cartesian coordinates, in the fre-quency domain.
Step 4.
Apply two-dimensional inverse Fourier transform to obtain the finalbackprojected image.Considering the bst formula, we can notice that the method also has anirregularity at the origin. This irregularity, caused by the division on zero inthe frequency domain, can be a problem in calculations. The simplest way toavoid this problem is to exclude the origin from the calculations. The other wayis to approximate this division, changing σ = 0 with some α , where α > K σ = (cid:40) σ , if σ (cid:54) = 0 , σ , otherwise , (5.11)where ∆ σ is the step of the mesh on σ .18 Regularized FBP: An Application
The bst formula (3.1) can be used to obtain an analytical solution of the stan-dard Tikhonov regularization problem in the feature space U minimize f ∈ U (cid:107) R f − g (cid:107) L + λ (cid:107) f (cid:107) L (6.1)In fact, the Euler-Lagrange equations provide the optimality condition for theabove optimization problem, i.e., f minimizes (6.1) if and only if [23]( R ∗ R + λ I ) f ( x ) = R ∗ g ( x ) (6.2)with R ∗ standing for the adjoint operator of the Radon transform and I theidentity operator in U . In fact, (6.2) are the so-called normal equations in theHilbert spaces U and V . Since R ∗ = B in the usual inner-product for L , theabove equation becomes ( BR + λ I ) f ( x ) = B g ( x ) (6.3)Applying the Fourier transformation on (6.3) and using property (3.10), weobtain the following standard resultˆ f ( ω ) 1 (cid:107) ω (cid:107) + λ ˆ f ( ω ) = (cid:100) B g ( ω ) ⇐⇒ ˆ f ( ω ) (cid:18) λ (cid:107) ω (cid:107) (cid:107) ω (cid:107) (cid:19) = (cid:100) B g ( ω ) (6.4)From (6.4) is easy to obtain f as a convolution of B g with an specific two-dimensional filter. If λ = 0 the analytical formula obtained is exactly the‘rho-filter layergram’ proposed in [24] consisting in a post-processing of thebackprojection (also mentioned earlier in this manuscript as filter of the back-projection ).The novelty here is that, if we change (6.4) to polar coordinates, we canimmediately apply the bst formula (3.1). Indeed, since ω = σ ξ θ , the pointwiseproduct becomes ˆ f ( σ ξ θ ) (cid:18) λσσ (cid:19) = (cid:100) B g ( σ ξ θ ) (6.5)which is essentially the same asˆ f ( σ ξ θ ) = (cid:18)
11 + λσ (cid:19) ˆ g ( σ, θ ) (6.6)The above equation is a regularized version of the Fourier-Slice-Theorem andcan be used to obtain f explicitly through any gridding strategy [8].Applying (6.6) in the change of variables of the Fourier representation of f ( x ) we finally obtain a new representation for the reconstructed image f , f ( x ) = (cid:90) R ˆ f ( ω ) e i ω · x d ω (6.7)= (cid:90) R d σ (cid:90) π d θf ( σ ξ θ ) | σ | e iσ x · ξ θ (6.8)= (cid:90) R d σ (cid:90) π d θ (cid:18)
11 + λ | σ | (cid:19) ˆ g ( σ, θ ) | σ | e iσ x · ξ θ (6.9)19quation (6.9) provides exactly the same reconstruction pattern as a typicalfiltered backprojection reconstruction algorithm, but with a different filter. Infact, we can generalize our regularized strategy in the following representation f λ ( x ) = B F λ g ( x ) (6.10)Now, { f λ } is a family of solutions of the optimization problem (6.1), dependingon the regularization parameter λ . The filter function F λ , in the frequencydomain reads (cid:99) F λ ( σ ) = | σ | λ | σ | (6.11)Our regularized solution (6.10) depends explicitly on the computation of theBackprojection operator B , and either the bst or Andersson’s formula can beused. All the algorithms were implemented using the fast Fourier framework fftw3 [7]. We validate our approach using five datasets: two real sinograms and threesimulated. The experimental sinograms (a slice from a wood-fiber and a porousrock) were obtained at the imaging beamline of the Brazilian Synchrotron lightsource and are high-resolution images with 2048 × × angles). There-fore, the feature images (either backprojected or filtered-backprojected) wererestored with 2048 × f ( x ) = (cid:88) j =1 δ ( x − a j ) (7.1)where { a j } are points randomly spanned over the domain [ − . , . × [ − . , . ρ ; on Fig.6 - logpolar reconstruction with usage of Partial Backprojections. Also we note thaton practice the first (adaptive) algoritm works 1.5-2 times faster.In futher tests we compare the results of bst with Log-Polar reconstructionwith an adaptive ρ selection, and the nfft approach [18] for the Fourier-Slice-Theorem.The regularized filtered-backprojection algorithm described in Section 6 wasapplied to the noisy data gathered for the wood-fiber and the rock sample. Theresults are shown in Figure 9 for the wood-fiber and the rock sample using onlythree values for the regularization parameter λ . In fact, an algorithm for the20a) (b) (a)(b) (c)(d) Figure 6: The results of the log-polar reconstuction using two approaches tocope with origin irregularity: (a),(b) - adaptive ρ selection; (c),(d) - partialBackprojections. 21a) (b) (c)Figure 7: Comparison between backprojected images. Column (a) shows the resultsobtained with bst , (b) with Andersson’s algorithm and (c) using nfft . See text fordetails. bst , (b) with Andersson’s algorithm and (c) using nfft . 23a) (b) (c)Figure 9: Reconstruction with real data and the regularized filtered-backprojection described in Section 6 for different λ (cid:48) s and using bst . (a) λ = 0 . λ = 0 .
02 and (c) λ = 0 . λ ∈ { . , . , . } .As it is known from Tikhonov regularization schemes, the bigger is λ , thesmoother the resulting image will be. This is clearly visible in Figure 9, whatindicates that such an approach could be used to increase the constrast in thefinal reconstructed image.All algorithms are fast due to usage of convolutions and Fast Fourier tech-niques. Computational complexity - that we denote by Ω - of Andersson’sapproach is similar to the complexity of the two dimensional convolution, i.e.Ω Andersson = N θ N ρ (log N ρ + log N θ ) + Ω L , (7.2)where Ω L is the summarized complexity of all log-polar interpolations (sinogramto log-polar and log-polar to cartesian). Here we assume that the Fourier trans-form of the kernel was pre-calculated (numerically, or analytically, like in [9])and is not taking it into account for Ω Andersson . In case of adaptive ρ selection, N ρ is being calculated using (5.7), which leads to some oversampling. For mostcommon sizes of the sinogram (i.e. 512 < N s , N x , N y < N ρ ≈ N s ), and the complexity of log-polar reconstructioncan be estimated as 8 N (log N + 1).The complexity of reconstruction, based on bst formula also depends onthe size of the final image. Much of computational complexity falls on the one-dimensional Fourier transforms, whose size is equal to the number of rays in the24 2 3 402040 k T i m e ( m s ec ) (a) bstlpnfft z T i m e ( m s ec ) (b) bstlp Figure 10: (a) Comparison of reconstruction time τ versus size of the image N = 256 × k with a constant zero-padding. (b) Reconstruction time versus zero-padding coefficient z , transforming the domain s ∈ [0 ,
1] to s ∈ [0 , z ] inpolar coordinates.sinogram, and to final 2D Fourier transform in cartesian coordinates. Hence,the complexity isΩ BST = N θ N s log ( N s ) + N x N y (log N x + log N y ) + Ω P , (7.3)where Ω P denotes the total amount of complexity for polar interpolations (po-lar to cartesian), usually Ω P = O ( N s ). Note that for clear reconstruction wealso need some oversampling on s , but not so big as in the previous case (notmore than two times). The comparison between time for reconstruction, de-pending on the size of the sinogram (for model task of Shepp-Logan phantomreconstruction) is presented on Figure 10.(a). On Figure 10.(b) we present thedependence of the reconstruction time versus the zero-padding factor, say z ,for the log-polar (adaptive) and bst . In our simulations, the zero-padding z increase the support of the sinogram from [0 ,
1] to [0 , z ] in polar coordinates.The backprojection (or reconstruction) for a slice with size 1024 × × angles) can be obtained in about 690 millisecondson a modest computer using cpu Intel(R) Core(TM) i7-3770 CPU @ 3.40GHzwith only 8 threads. Of course, the programs developed by authors can besufficiently optimized and powered for gpu , which can make the backprojectionconsiderably faster using either Andersson’s formula or bst .Figure 11 presents some benchmarks of accuracy for the developed algo-rithms with other known backprojection techniques. More precisely, we presentthe resulting mean squared error (MSE) versus the error in the input data (i.e.,the sinogram). Calculation were done for the Shepp-Logan phantom with addi-tion of Poisson noise to the analytical sinogram. From Figure 11 one can notethat for weak noises BST shows the best accuracies, while Log-Polar reconstruc-tion is very stable to strong noise.On Figure 12 we present the slice of the reconstruction of our analytical25 0 . . . . . . · − Sinogram error (%) F il t e r e db a c k p r o j ec t i o n e rr o r ( % ) (a) bstlpnfftslanbres Figure 11: (a) Mean square error (MSE) of the result, obtained using differentalgorithms, in dependence on MSE of enter data (sinogram).26 . − − . − . − . − . . . . . . . . . x B g (a) Exact BPBP with BST
Figure 12: Analytical and numerical (BST) backprojections for Example 2, acircular function presented in Section 4.example (see Section 4, Example 2). Since we obtain backprojection of the circlefunction analytically, we can compare the result of numerical backprojectionwith BST and analytical solution.
In this manuscript we have proposed a new backprojection technique ( bst ) andcompared it against two other already established algorithms, the log-polar ( lp )approach from Andersson [9] and with the use of Non-Uniform Fast FourierTransforms ( nfft ) from [18]. With the increasing size of input data, mea-sured at imaging beamlines from synchrotron facilities, the need for fast post-processing of the data becomes an immediate demand. Either conventionalreconstruction techniques like filtered backprojection or more robust iterativemethods can benefit from a fast backprojection algorithm. Finally, in order todemonstrate this possibility in practice, we have provided an application wherethe Tikhonov image is computed by means of the new bst formula in a shortamount of time, thereby enabling the practical application of interesting regu-larization schemes to very large datasets. A Integral representations
We use the following standard representation for path integrals, the proof can befound in [4]: for a continuously differentiable m : R m → R such that (cid:107)∇ m (cid:107) (cid:54) = t is true that : (cid:90) R m h ( y ) δ ( m ( y ))d y = (cid:90) C − (0) h ( y ( s )) (cid:107)∇ m ( y ( s )) (cid:107) d s ( x ) . (A.1)where d s ( x ) is the arclength measure along curve C − (0). Assuming that g ∈ V is a given sinogram and x a pixel in the reconstruction region. The backpro-jection (1.4) of g is defined as the contribution of all possible straight lines,parameterized by the angle θ , and passing through x . Using the sifting prop-erty of the Delta distribution, we have B g ( x ) = (cid:90) π g ( x · ξ θ , θ )d θ = (cid:90) π (cid:90) R g ( t, θ ) δ ( t − x · ξ θ )d t d θ (A.2)Switching the above integral from ( t, θ ) coordinates to cartesian coordinates y ∈ R we have | t | d t d θ = d y ; where B g now becomes B g ( x ) = (cid:90) R [ g ] c ( y ) δ ( m ( y )) 1 (cid:107) y (cid:107) d y (A.3)with [ g ] c ( y ) = g ( t ( y ) , θ ( y )) refering to the sinogram g in cartesian coordinates.In fact, | t | = (cid:107) y (cid:107) is the unsigned distance to the origin and θ = arctan( y y ) ∈ [0 , π ] is the angle with respect to the y -axis. Function m reads m ( y ) = t − x · ξ θ = (cid:107) y (cid:107) − x cos θ ( y ) − x sin θ ( y ) (A.4)= (cid:107) y (cid:107) − x y (cid:107) y (cid:107) − x y (cid:107) y (cid:107) = (cid:107) y (cid:107) − ( x y + x y ) (cid:107) y (cid:107) (A.5)= y · ( y − x ) (cid:107) y (cid:107) (A.6)From (A.6), (A.3) and the property δ ( au ) = | a | δ ( u ) for all u ∈ R , the backpro-jection now follows: B g ( x ) = (cid:90) R [ g ] c ( y ) δ ( κ x ( y ))d y , κ x ( y ) = y · ( y − x ) (A.7)It should be noted that, for a fixed x ∈ R , the set κ − x (0) = { y ∈ R : κ x ( y ) = 0 } is defined as a circle in the plane. Indeed, since y · ( y − x ) = y · y − y · (cid:0) x (cid:1) = (cid:13)(cid:13) y − x (cid:13)(cid:13) − (cid:13)(cid:13) x (cid:13)(cid:13) , it follows that κ − x (0) is a circle passingthrough the origin y = , centered at x and with radius (cid:107) x (cid:107) . Since κ − x (0) = { x + r ξ θ : θ ∈ [0 , π ] , r = (cid:107) x (cid:107)} is a parametric representation of the circle,the backprojection operator also reads, in an alternative form: B is a stackingoperator through circles κ − x (0) : B g ( x ) = (cid:90) κ − x (0) [ g ] c ( y ) (cid:107) y − x (cid:107) d s = 12 (cid:90) π [ g ] c (cid:18) x + 12 (cid:107) x (cid:107) ξ θ (cid:19) d θ (A.8)The above representation follows from ds = (cid:107) x (cid:107) dθ , (A.7) and (A.1) with ∇ κ x ( y ) = 2 y − x . Last equality comes from y = x + (cid:107) x (cid:107) ξ θ ∈ κ − x (0) for28ome θ . Therefore, in cartesian coordinates, the backprojection contribution fora ball { z ∈ R : (cid:107) z − x (cid:107) ≤ (cid:15) } comes from a family of circles passing throughthe ball and the origin, this fact seems to be related to the comet-tail region mentioned by [25]. B Point Source in Log-Polar Coordinates
In this section we present the details for the following log-polar representationof the backprojected image b ( e ρ ξ θ ) = (cid:90) d β (cid:90) d u δ ( a · ξ β − e u (cid:124) (cid:123)(cid:122) (cid:125) m ( u ) ) δ (1 − e ρ − u cos( θ − β ) (cid:124) (cid:123)(cid:122) (cid:125) (cid:96) ( u ) ) (B.1)= 1 (cid:112) (cos( θ − φ ) e ρ ) + ( e A − sin( θ − φ ) e ρ ) (B.2) Proof.
Using the property of the Delta distribution [20], δ ( (cid:96) ( u )) = (cid:88) k δ ( u − u k ) | (cid:96) (cid:48) ( u k ) | (B.3)where u k are roots of (cid:96) ( u ). In our case, since (cid:96) has only one zero u = ρ + ln cos( θ − β ) (B.4)and (cid:96) (cid:48) ( u ) = 1, function b reads b ( e ρ ξ θ ) = (cid:90) d β m ( u ) δ ( u − u ) (B.5)Due to the sifting property of the Delta, we obtain b ( e ρ ξ θ ) = (cid:90) d β m ( ρ + ln cos( θ − β )) (B.6)= (cid:90) d β δ ( a · ξ β − e ρ +ln cos( θ − β ) ) (B.7)= (cid:90) d β δ ( a · ξ β − e ρ cos( θ − β )) (B.8)(B.9)Using a = e A ξ φ we obtain b ( e ρ ξ θ ) = (cid:90) d β δ ( e A ξ φ · ξ β − e ρ cos( θ − β )) (B.10)= (cid:90) d β δ ( e A cos( φ − β ) − e ρ cos( θ − β ) (cid:124) (cid:123)(cid:122) (cid:125) L ( β ) ) (B.11)29ince cos( θ − β ) = cos(( φ − β ) + ( θ − φ )), function L can be rewritten as L ( β ) = e A cos( φ − β (cid:124) (cid:123)(cid:122) (cid:125) α ) − e ρ cos( φ − β ) sin( θ − φ ) (cid:124) (cid:123)(cid:122) (cid:125) S − sin( φ − β ) cos( θ − β ) (cid:124) (cid:123)(cid:122) (cid:125) C = e A cos α − e ρ S cos α + e ρ C sin α = ( e A − Se ρ ) (cid:124) (cid:123)(cid:122) (cid:125) z cos α + e ρ C (cid:124)(cid:123)(cid:122)(cid:125) y sin α Hence, the root β of L must satisfy L ( β ) = 0, i.e., tan α = − zy ≡ k andsin α = k k , cos α = k . Therefore, L (cid:48) ( β ) = z k √ k − y √ k = − z (cid:112) y + z − k (cid:112) y + z (B.12)= − (cid:112) y + z = − (cid:113) ( Ce ρ ) + ( e A − Se ρ ) (B.13)Finally, b ( e ρ ξ θ ) = (cid:90) d β δ ( β − β )) | L (cid:48) ( β ) | = 1 | L (cid:48) ( β ) | (B.14)= 1 (cid:112) ( Ce ρ ) + ( e A − Se ρ ) (B.15)= 1 (cid:112) (cos( θ − φ ) e ρ ) + ( e A − sin( θ − φ ) e ρ ) (B.16)which is the final representation of the backprojected image in log-polar coor-dinates. References [1] Stanley R Deans.
The Radon transform and some of its applications .Courier Corporation, 2007.[2] Avinash C Kak and Malcolm Slaney.
Principles of computerized tomo-graphic imaging . IEEE press, 1988.[3] F Wubbeling and F Natterer. Mathematical methods in image reconstruc-tion.
SIAM, Philadelphia , 8:16, 2001.[4] Sigurdur Helgason.
The Radon Transform on R n . Springer, 2011.[5] Heinz Pr¨ufer. Neue herleitung der sturm-liouvilleschen reihenentwicklungstetiger funktionen.
Mathematische Annalen , 95(1):499–518, 1926.306] Lawrence A Shepp and Benjamin F Logan. The fourier reconstruction of ahead section.
Nuclear Science, IEEE Transactions on , 21(3):21–43, 1974.[7] Johan Wald´en. Analysis of the direct fourier method for computer tomog-raphy.
Medical Imaging, IEEE Transactions on , 19(3):211–222, 2000.[8] F Marone and M Stampanoni. Regridding reconstruction algorithm for real-time tomographic imaging.
Journal of synchrotron radiation , 19(6):1029–1037, 2012.[9] Fredrik Andersson. Fast inversion of the radon transform using log-polarcoordinates and partial back-projections.
SIAM Journal on Applied Math-ematics , 65(3):818–837, 2005.[10] Achi Brandt, Jordan Mann, Matvei Brodski, and Meirav Galun. A fastand accurate multilevel inversion of the radon transform.
SIAM Journalon Applied Mathematics , 60(2):437–462, 2000.[11] Ashvin George and Yoram Bresler. Fast tomographic reconstruction viarotation-based hierarchical backprojection.
SIAM Journal on AppliedMathematics , 68(2):574–597, 2007.[12] Samit Basu and Yoram Bresler. O (n 2 log 2 n) filtered backprojectionreconstruction algorithm for tomography.
Image Processing, IEEE Trans-actions on , 9(10):1760–1773, 2000.[13] Eduardo Miqueles and Elias S Helou. Fast backprojection operator forsynchrotron tomographic data. In
European Conference on Mathematicsfor Industry . Springer, 2014.[14] Daniel Potts and Gabriele Steidl. New fourier reconstruction algorithmsfor computerized tomography. In
International Symposium on Optical Sci-ence and Technology , pages 13–23. International Society for Optics andPhotonics, 2000.[15] Karsten Fourmont. Non-equispaced fast fourier transforms with appli-cations to tomography.
Journal of Fourier Analysis and Applications ,9(5):431–450, 2003.[16] Frank Natterer.
The mathematics of computerized tomography , volume 32.Siam, 1986.[17] A Averbuch, RR Coifman, DL Donoho, M Israeli, and J Walden.
Fast SlantStack: A notion of Radon transform for data in a cartesian grid which israpidly computible, algebraically exact, geometrically faithful and invertible .Department of Statistics, Stanford University, 2001.[18] Stefan Kunis and Daniel Potts.
Time and memory requirements of thenonequispaced FFT . Technische Universit¨at Chemnitz. Fakult¨at f¨ur Math-ematik, 2006. 3119] James W Cooley and John W Tukey. An algorithm for the machine calcula-tion of complex fourier series.
Mathematics of computation , 19(90):297–301,1965.[20] Ron Bracewell. The fourier transform and its applications.
New York , 5,1965.[21] John D Towers. Two methods for discretizing a delta function supportedon a level set.
Journal of Computational Physics , 220(2):915–931, 2007.[22] J Kaiser and R Schafer. On the use of the i 0-sinh window for spectrumanalysis.
IEEE Transactions on Acoustics, Speech, and Signal Processing ,28(1):105–107, 1980.[23] David G Luenberger.
Optimization by vector space methods . John Wiley& Sons, 1997.[24] Gabor T Herman. Image reconstruction from projections.