A standard and linear fan-beam Fourier backprojection theorem
AA standard and linear fan-beam Fourierbackprojection theorem
Patricio Guerrero, Matheus Bernardi and Eduardo Miqueles ∗ July 2020
Abstract
We propose a theoretical formulation for the tomographic fan-beambackprojection in standard and linear geometries. The proposed formulais obtained from a recent backprojection formulation for the parallel case.Such formula is written as a Bessel-Neumann series representation in thefrequency domain of the target space in polar coordinates. A mathemat-ical proof is provided together with numerical simulations compared withconventional fan-beam backprojection representations to validate our for-mulation.
Keywords: fan-beam, backprojection, tomographic reconstruction.
Fan-beam tomographic measurements are used in different modalities of non-destructive imaging, as those obtained using x-rays. A typical tomographicdevice using a fan-beam geometry is shown in Figure 1. We assume that thedistance between the pair source-detector is high compared to the size of thesample. This is a widely used and known technique, and there are many re-construction algorithms for this configuration. After being generated with agiven aperture angle and a fixed distance source-detector, the wavefront hits thesample originating a signal (or image) on the detector. Different propagationregimes can be considered with a varying distance [9], although we will considera pure mathematical signal idealized as the Radon transform of the given ob-ject. For the parallel tomographic case by means of the classical backprojectionformulation for image reconstruction, we can use the recently backprojectionslice theorem formulation [10]. It is a formula that reduces a complete backpro-jection of a 2D sinogram from a computational cost of O ( n ) to O ( n log n ) with n × n pixels in the final reconstructed image and n tomographic rotations. In ∗ P.Guerrero, M. Bernardi and E.Miqueles are with the Scientific Computing Group fromthe Brazilian Synchrotron Light Laboratory at the Brazilian Center for Research in Energyand Materials, Brazil, Campinas-SP. Corresponding author: [email protected] a r X i v : . [ m a t h . NA ] J u l fD β (b)(a) x ∆ γ ∆ s source Figure 1: Tomographic setup for a fan-beam geometry: (a) standard detector,i.e., equispaced angular mesh with ∆ γ ; sinograms obtained here are denoted by w ( γ, β ) (b) linear detector, i.e., equispaced mesh at a perpendicular with thecentral ray and linear detector with ∆ s ; sinograms obtained here are denotedby g ( s, β ).this work, we want to take advantage of that formulation for two other popularfan-beam geometries, that are a) equispaced angles within the fan with a regularsize ∆ γ , and b) equispaced points in a linear detector with a mesh step ∆ s . Asindicated in [12], we refer to each acquisition as standard fan and linear fan ,respectively, and illustrated at Figure 1.(a) and 1.(b). The linear case is easierto be implemented at a synchrotron beamline, whereas the first is more used inthe medical case. Also, benchtop ct scanners, also based on the conventionalcone-beam geometry could take advantage of the approach presented in thiswork. For completeness, we discuss both the standard and the linear case.In [16, 18] a generalized Fourier slice theorem for the case of fan-beam ge-ometries is obtained, but presenting the same computational complexity as aconventional backprojection in the parallel case. Also, a strong rebinning pro-cess on the measured data is needed within the algorithm, which is not desirablein our case. Further relations on the frequency domain are obtained in [7] wherea rebinning is also necessary. A more elegant approach is established in [2], buta rebinning in the frequency domain is also mandatory. Further advanced rebin-ning techniques are also established in [5] using a hierarchical approach. A seriesformulation where the backprojection is presented as the first order approxima-tion for a general inversion in the standard fan-beam geometry is presented in[14].This work is organized as follows. Chapter 2 presents a review of conven-tional fan-beam backprojections algorithms for both standard and linear geome-tries. In chapter 3 we rewrite the rebinning strategy in terms of the adjoint of therebinning operator as a two-step backprojection formula. Chapter 4 presentsour fan-beam backprojection theorem and its proof as a Bessel-Neumann se-ries representation while Chapter 5 supports numerically our method using asimulated dataset. 2 Fan-beam backprojections
The two-dimensional Radon transform is defined as the linear operator R : U → V where U is the space of rapidly decreasing functions defined on R , so-called feature space; and V is the sinogram space defined on the unit cylinder R × R / π Z where R / π Z is the quotient group of reals modulo 2 π . This definitionholds when standard parallel projections are considered.Let ( t, θ ) ∈ R × R / π Z , for f ∈ U , R f is defined as the integral R f ( t, θ ) = (cid:90) R f ( x ) δ ( t − x · ξ θ ) d x , (1)where δ is the Delta distribution and ξ θ = (cos θ, sin θ ) T . We will denote f P ∈ U as a feature function acting on x = ( ρ, φ ) ∈ R + × [0 , π [ expressed in polarcoordinates. R f P is then written as R f P ( t, θ ) = (cid:90) π (cid:90) R + f P ( ρ, φ ) δ ( t − ρ cos( θ − φ )) ρ d ρ d φ. (2)Note that R f is an even function on the unit cylinder, that is, R f ( t, θ ) = R f ( − t, θ + π ) . (3) The adjoint operator.
Let T : X → Y be a bounded linear operator betweenHilbert spaces X and Y . Its adjoint operator is defined as the unique operator T ∗ : Y → X that verifies ∀ u ∈ X, ∀ v ∈ Y, (cid:104) T u, v (cid:105) Y = (cid:104) u, T ∗ v (cid:105) X . (4)The adjoint of the Radon transform [11] noted R ∗ ≡ B , closely related tothe backprojection, is the operator B : V → U defined for p ∈ V and x ∈ R by B p ( x ) = 2 (cid:90) π p ( x · ξ θ , θ ) d θ. (5)With standard L inner products on U and V , it is clear that B verifies (4),the factor 2 on the above integral comes from the symmetry property (3).Computing B p could be extremely expensive for discrete versions of thesinogram p , when x covers a domain with a large number of points (pixels inpractice) and also ( t, θ ) covers a large number of pixels and a variety of angles(according to Crowther’s criterion [3]). This is the case for synchrotron tomo-graphic projections using high-resolution detectors with more than 2048 × B was obtained in the frequency domain using polar coordinates,the Backprojection Slice Theorem ( bst ). Such theorem reads for p ∈ V and σ ∈ R , (cid:99) B p ( σ ξ θ ) = ˆ p ( σ, θ ) σ , (6)3 V V s fb p w R M s B B s Figure 2: Action diagram for operators {B , R} and {B s , R s } , with R s = M s R , B s = BM ∗ s and a generalized change of variables M s .where (cid:99) B p is the 2D Fourier transform of B p in polar frequency coordinatesand ˆ p the 1D Fourier transform of p with respect to t . The action of {R , B} ispresented in Figure 2. It is a well known fact that { p, f } are related through the Fourier slice Theorem [11], while { b, p } through the backprojection slice theorem [10]. There exist two different fan-beam parametrizations for the Radon transform[6], the first R s : U → V s is referred as the standard fan-beam transform, wherethe domain of sinograms V s lies within the interval R /π Z × R / π Z . The secondcase, R (cid:96) : U → V (cid:96) is referred as the linear fan-beam transform, with sinogramsdomain V (cid:96) varying within the interval R × R / π Z .Let us consider ( γ, β ) ∈ R /π Z × R / π Z a pair of angles and s ∈ R . Acting ona feature function f ∈ U , the linear operators R s and R (cid:96) are defined respectivelyas w ( γ, β ) = R f ( D sin γ, β + γ ) , (7)and g ( s, β ) = R f ( sD √ s + D , β + arctan sD ) , (8)where D > w ∈ V s and g ∈ V (cid:96) . The following symmetry relationships holds w ( γ, β ) = w ( − γ, β + 2 γ + π ) ,g ( s, β ) = g ( − s, β + 2 arctan sD + π ) . (9)Since both operations represent a change of variables in the classical parallelsinogram, we can use the following notation, R s f ( γ, β ) = M s R f ( γ, β ) , R (cid:96) f ( s, β ) = M (cid:96) R f ( s, β ) , (10)4here M s : V → V s and M (cid:96) : V → V (cid:96) are rebinning operators acting on aparallel sinogram p , respectively as M s p ( γ, β ) = p ( D sin γ, β + γ ) , M (cid:96) p ( s, β ) = p ( sD √ s + D , β + arctan sD ) . (11)From the experimental point of view, we assume that the origin is positionedat the center of rotation, where the sample should also be centered. In bothgeometries, for simplicity and without any loss of generality, a 1D virtual detec-tor will also be considered centered at the origin. This can be done simply byscaling numerically the detector. In the standard geometry, the angle γ from thecentral ray indicates a position in the circle centered at the beam source withradius D , this circle acts as the detector in this geometry and then it’s evidently π -periodic on γ . Whereas in the linear geometry, s is the signed height on thelinear detector perpendicular to the central ray and centered at the origin, seeFigure 1.The change of variables operation M s is depicted in Figure 2. M s and M (cid:96) are well defined operators due to β + γ and β + arctan sD both belong to R / π Z .Since is true that R s = M s R , a conventional functional relation for space U and V s provides us with the following statement, (cid:104)R s f, w (cid:105) V s = (cid:104)M s R f, w (cid:105) V s = (cid:104)R f, M ∗ s w (cid:105) V = (cid:104) f, R ∗ M ∗ s w (cid:105) U . (12)Now, since R ∗ = B , it follows that R ∗ s = BM ∗ s . In this work we propose aFourier approach for equations R ∗ s = BM ∗ s and R ∗ (cid:96) = BM ∗ (cid:96) . (13) The adjoint operator for fan-beam transforms is a weighted backprojection op-erator of the standard parallel transform [6]. Considering both fan-beam ge-ometries, we provide two main results for the adjoint transform of operators R s and R (cid:96) . We denote in the following r β the Cartesian coordinates of the source. Lemma 1.
The operator B s : V s → U , defined for w ∈ V s and x ∈ R by B s w ( x ) = (cid:90) π L β w ( γ β , β ) d β, (14)is the the adjoint of R s in the sense of (4). Here, L β = (cid:107) r β − x (cid:107) is the distanceof the source with the backprojected position x and γ β is the angle of such pointwithin the fan, i.e., cos γ β = 1 DL β r β · ( r β − x ).5 roof. Let f P ∈ U and w ∈ V s . Using the polar representation (2) of R and therelationship ρ cos( β + γ − φ ) = L β sin( γ β − γ ) [6], then (cid:104)R s f P , w (cid:105) V s = (cid:90) π (cid:90) π w ( γ, β ) R s f P ( γ, β ) d γ d β = (cid:90) π (cid:90) π w ( γ, β ) (cid:90) π (cid:90) R + f P ( ρ, φ ) δ ( L β sin( γ β − γ )) ρ d ρ d φ d γ d β, (15)which, from Fubini’s theorem, becomes (cid:104)R s f P , w (cid:105) V s = (cid:90) π (cid:90) R + f ( ρ, φ ) (cid:90) π (cid:90) π L β w ( γ, β ) δ (sin( γ β − γ )) d γ d β ρ d ρ d φ. Therefore, since δ (sin( γ β − γ )) = δ ( γ β − γ ) for γ ∈ [0 , π [ and a fixed γ β , afterperforming the above γ -integration, B s defined in (14) verifies (4). Lemma 2.
The operator B (cid:96) : V (cid:96) → U , defined for g ∈ V (cid:96) and x ∈ R by B (cid:96) g ( x ) = (cid:90) π D U β (cid:113) s β + D g ( s β , β ) d β, (16)is the adjoint of R (cid:96) in the sense of (4). U β is the ratio of the scalar projec-tion of r β − x on the central ray to the source-origin distance, and s β is thecorresponding height s for x at a source angle β , i.e., U β = D − r β · x D , s β = x · ξ β U β . (17) Proof.
As the proof of Lemma 1, it follows directly from (4) and the polarrepresentation (2).
The two-step adjoint formulations for each geometry, namely R ∗ s = BM ∗ s and R ∗ (cid:96) = BM ∗ (cid:96) derived following (12) are going to be detailed here. It is clearthat the adjoint of the rebinning operators {M s , M (cid:96) } plays an important rolein such formulations. Lemma 3 (Adjoint of a rebinning operator) . Let M : X → Y be a bijective anddifferentiable with continuous inverse rebinning operator between two Hilbertspaces X and Y . Its adjoint M ∗ : Y → X is then given by M ∗ = J M − , where J is the Jacobian determinant of the rebinning operation and M − its inverserebinning. Proof. If u ∈ X and v ∈ Y , a direct verification of (cid:104)M u, v (cid:105) Y = (cid:104) u, M ∗ v (cid:105) X through an integral representation and a change of variables will prove theresult. The bijectivity is needed to preserve the range of integration in bothdomains of X and Y and for the inverse to exist.6 emma 4. The adjoint operator of M s and M (cid:96) are respectively defined for w ∈ V s and g ∈ V (cid:96) by M ∗ s w ( t, θ ) = w ( γ ( t ) , θ − γ ( t )) J s ( t ) , (18) M ∗ (cid:96) g ( t, θ ) = g ( s ( t ) , θ − γ ( t )) J (cid:96) ( t ) , (19)where γ ( t ) = arcsin tD , s ( t ) = tD √ D − t , J s ( t ) = 1 √ D − t , J (cid:96) ( t ) = D ( D − t ) / . (20) Proof.
Being both operators bijective [12] and continuously differentiable, thisis an immediate application of Lemma 3.The formal adjoints of R s and R (cid:96) are presented in Lemmas 1 and 2. Theproblem with these formulations is the difficulty for a low-cost implementationalgorithm. To circumvent this problem, we use the fact that {R , B} are boundedoperators between Hilbert spaces U and V , here understood as the space ofrapidly decreasing functions with two-variables. Hence, we obtain the followingresult. Theorem 1.
Considering the two fan-beam geometries, formal adjoints foroperators R s and R (cid:96) are given explicitly by expressions (14) and (16), which onthe other hand, are also given by B s = BM ∗ s and B (cid:96) = BM ∗ (cid:96) respectively.Proof. This is an immediate consequence of the uniqueness of the adjoint forbounded operators on Hilbert spaces. In fact, since B and M ∗ s are bounded,the composition is also bounded. The same applies for M ∗ (cid:96) . As mentioned, the idea is a Fourier approach for equations (13) to then makeuse of the bst formula (6) for a fixed θ in the form (cid:100) B s w ( σ ξ θ ) = (cid:92) BM ∗ s w ( σ, θ ) = 1 σ (cid:92) M ∗ s w ( σ, θ ) . (21)We start announcing our Fourier-based fan-beam backprojection theoremfor standard geometry, and then we show how it is easily adapted to the linearcase. From (21), we observe the need of computing (cid:92) M ∗ s w that is written as aseries representation in the following Lemma. Lemma 5.
Given a standard fan-beam sinogram w ∈ V s , first define Z ∈ V s by Z ( γ, θ ) = w ( γ, θ − γ ) for all ( γ, θ ) ∈ R /π Z × R / π Z . As Z is 2 π -periodic7ith respect to γ we can write its Fourier series expansion with the Fouriercoefficients c n ( θ ) = 12 π (cid:90) π Z ( γ, θ ) e − inγ d γ, (22)from where we define the coefficients b n = 2 π [ c n + ( − n ¯ c n ] for n ≥ , and b = 2 πc . (23)Then we have a Bessel-Neumann series description of (cid:92) M ∗ s w as (cid:92) M ∗ s w ( σ, θ ) = ∞ (cid:88) n =0 b n ( θ ) J n ( Dσ ) , (24)where { J n } is a sequence of Bessel functions of the first kind. Proof.
Given w ∈ V s , from (18) and (20) the 1D Fourier transform of M ∗ s w ( t, θ )with respect to t is (cid:92) M ∗ s w ( σ, θ ) = (cid:90) R w ( γ ( t ) , θ − γ ( t )) e − itσ J (cid:96) ( t ) d t = (cid:90) π w ( γ, θ − γ ) e − iDσ sin γ d γ = (cid:90) π Z ( γ, θ ) e − iDσ sin γ d γ. (25)After expanding Z with c n we have (cid:92) M ∗ s w ( σ, θ ) = (cid:88) n c n ( θ ) (cid:90) π e i [ nγ − Dσ sin γ ] d γ = 2 π (cid:88) n c n ( θ ) J n ( Dσ ) , (26)from where the result claims using (23) and the property J − n ( x ) = ( − n J n ( x ).We can now state our fan-beam backprojection theorem for standard geome-tries. Theorem 2 (Standard fan-beam backprojection Theorem) . The adjoint oper-ator B s : V s → U acting on w ∈ V s can be written in the Fourier domain as theseries (cid:100) B s w ( σ ξ θ ) = 1 σ ∞ (cid:88) n =0 b n ( θ ) J n ( Dσ ) , (27) where ( σ, θ ) ∈ R + × R / π Z and the coefficients b n are computed as in Lemma5.Proof. This is an immediate consequence of (21) and Lemma 5.8
V V s V (cid:96) b g wz Z AB (cid:96) B s L ( L ∗ ) − Figure 3: Diagram for the backprojection operator B (cid:96) (without rebinning tothe space V ) through the action of ( L ∗ ) − , A and B s . Here, z = τ w ,with τ asdescribed in Theorem 3. See text for details. The adjoint operator for the linear fan-beam transform follows. Switching be-tween fan-beam geometries requires only a one-dimensional interpolation on thefirst variable. In fact, taking g ∈ V (cid:96) , the rebinning operator on the first variable L : V (cid:96) → V s defined by L g ( γ, β ) = (cid:26) g ( D tan γ, β ) , if γ (cid:54) = π/ π , otherwise ⇐⇒ R s = LR (cid:96) . (28)provides a sinogram w on the space V s . Such bijective operator, according toLemma 3, has an adjoint operator given by L ∗ w ( s, β ) = J ( s ) L − w ( s, β ) , (29)with J ( s ) = DD + s , L − w ( s, β ) = w (arctan sD , β ) . (30)The following Theorem enables us to provide a backprojection algorithm forthe linear fan-beam geometry. Theorem 3.
The adjoint operators B s and B (cid:96) are related through B (cid:96) g = B s τ L g ,with τ ( γ ) = D sec γ , for all g ∈ V (cid:96) .Proof. Starting with B s = R ∗ s and using (28) we obtain B s = ( LR (cid:96) ) ∗ = R ∗ (cid:96) L ∗ = B (cid:96) L ∗ , from where follows B (cid:96) = B s ( L ∗ ) − . Using (29) and the fact that ( L ∗ ) − = J ( s ( γ )) − L with s ( γ ) = D tan γ , the proposed equation is obtained.9ow, using the fact that B s and B are related (from Theorem 2), our resultingformulation for the linear fan-beam backprojection is obtained by the followingconstruction. Theorem 4 (Linear fan-beam backprojection Theorem) . The operator B (cid:96) : V s → U acts on g ∈ V (cid:96) as the following two-steps operation.(i) Let z ( γ, β ) = τ ( γ ) w ( γ, β ) with w = L g be his representation in the space V s . From z we compute Z = A z defined as Z ( γ, θ ) = z ( γ, θ − γ ) .(ii) The backprojection finally becomes as the Bessel-Neumann series in (27)where the coefficients b n (23) are obtained from Z .Proof. Each step is justified respectively by Theorems 3 and 2.Figure 3 illustrates the action of our method, where a backprojection b ∈ U is obtained in two-steps. The unknown backprojection B (cid:96) g is obtained in thefrequency domain from (27), using polar coordinates by two main operations toobtain sinogram Z ∈ V s through L (interpolation on the first variable) and A (interpolating at the second variable). Simulations will be concentrated on the linear case due to the standard case isone rebinning operation less than the linear, so results would be slightly betterin the last not simulated case. The phantom in Figure 4(a) will be considered.We first computed the parallel Radon transform (1) to then obtain a linear fansinogram through operation M (cid:96) (11). Discrete configuration.
The domain of the feature phantom is the unit disk (cid:107) x (cid:107) ≤ N × N uniformly sampled quadrilateral mesh domain.For the domain of U to compute (1), we have { N, N θ } points uniformly spacedwith ( t, θ ) ∈ [ − , × [0 , π ].For the linear fan-beam geometry, we consider the distance source-origin D = 10, then ¯ s = D ( D − − / is the highest detector position that measuresthe sample on the unit disk subtending an angle ¯ γ = arcsin 1 /D . We recallthat the detector is supposed to be centred, or scaled, at the origin. Due to thesymmetry property (9), ( s, β ) needs to cover the sets [ − ¯ s, ¯ s ] × [0 , π + 2¯ γ [. Forsimplicity, we sampled these sets as in the parallel geometry uniformly with N and N θ points. Finally, here we will use N = N θ = 1024. Implementation and complexity.
Implementation details of the two stepsof Theorem 4 follows given the sinogram in Figure 4(b). Step (i) is a sequenceof two 1D interpolations, L followed by A to obtain Z . L is a rebinning on thespatial variable from linear to standard fan geometries, so only needed in thelinear case, and A is a linear operation on the angular variable then not causing10 x (cid:107) ≤ β ∈ [0 , π + 2¯ γ [ s ∈ [ − ¯ s , ¯ s ] (a) (b)Figure 4: Feature phantom (a) and its linear fan-beam sinogram (b).rebinning errors with linear interpolation. Most importantly, both operationsare not expensive having complexity O ( N N θ ) and they were easily parallelizedas they are both one-dimensional.Step (ii) is more challenging, we need to compute a truncated Bessel-Neumannseries (27) until convergence. Such convergence happens due to the fact that J n ( x ) /x has a pointwise convergence to zero since | J n ( x ) | ≤ | x | n /n ! for all x ∈ R + [1, eq. 9.1.62]. Of course the number of terms in the sequence J n is thesame of b n , that itself is the number of samples of s in the sinogram g follow-ing the Nyquist sampling theorem [6]. Convergence of the series can thus bereached by first zero-padding the linear variable s of g , or equivalently, γ of Z after performing operation A .Computing coefficients b n requires N θ independent 1D Fourier transformsto obtain (22) followed by the easy operation (23). They were also parallelizedin θ having a complexity O ( N log N ) each with Fast Fourier techniques. Impor-tantly, the sequence { J n ( Dσ ) } in (27) not depending on the sample is computedonly once for a given setup and were used as a lookup table for the computationof (27), not being costly at computing time.The truncated sum (27) is finally a matrix multiplication of both matri-ces composed with elements b n ( θ ) and J n ( Dσ ) for sampled values of { θ, σ } .The frequency variable σ was sampled also according to the Nyquist samplingtheorem. Such multiplication has a complexity O ( N ) but can be reduced to O ( N . ) with fast matrix multiplication algorithms [8, 15], although we usedthe direct approach in our simulations. We recall that conventional approachespresent a cost of O ( N ). Finally, the Fourier domain of the backprojected im-age is obtained at a polar grid. In order to obtain a Cartesian grid, we simplyused bilinear interpolation because of the previously performed zero-paddingof Z in γ also helps us to reduce interpolation errors at high frequencies. Agridding strategy to pass from polar to Cartesian coordinates [12], that handles11a) (b)Figure 5: Backprojections of the fan-beam sinogram of Fig. 4(b) (a) and itsassociated filtered backprojection (b) using the Bessel-Neuman series (27) ofTheorem 4.better these interpolations errors can also be used. We finally underline thatbackprojecting an image with our method assumes a null dc component on thetarget image, due to the kernel 1 /σ . The true dc component was approximatedfrom the input sinogram as the average of a fan-beam projection for all s in thedetector at a fixed angle β .Two backprojections methods are presented, first the conventional rebinningthrough parallel geometry by operation M ∗ (cid:96) in (19) followed by the standardbackprojection (5) as presented in theorem 1. Good results are expected withthis method mainly due to our sinogram was obtained through M (cid:96) . The ob-tained results will help us as a reference for our formulation (27) in Theorem 4to be presented lastly. Backprojection of sinogram in Figure 4(b) is presented inFigure 5(a) through series (27) and its 1D profile at x = 0 .
12 in Figure 6(a) bymeans of the two mentioned methods. A nice agreement between both profilescan be appreciated.To obtain image reconstructions, we filtered the parallel sinogram with theconventional ramp filter [6] before applying the rebinning operation M (cid:96) to ob-tain a filtered fan-beam sinogram. We proceded in this way for sake of simplicitybecause the filtering processes for fan-beam sinograms, detailed in [6, 12], is notin the scope of this work. We are interested in the backprojection exclusively.Backprojection through our method of this filtered sinogram is presented inFigure 5(b) where we can appreciate a good image reconstruction. Again, theprofile x = 0 .
12 is plotted in Figure 6(b) to compare the two methods togetherwith the target phantom.Both algorithms were implemented on a multithreading architecture for eachstep. We used the python numerical library scipy [13] mainly to compute FastFourier transforms, linear interpolations and Bessel coefficients. Results are sat-isfactory however no big difference was appreciated in computing time between12 x B (cid:96) g rebinningTheo. 4 − x B (cid:96) ( g (cid:63) r a m p ) rebbiningTheo. 4phantom (a) (b)Figure 6: Comparison of profiles at x = 0 .
12 (pixel 574) of backprojections(left) and filtered backprojections (right) of the rebinning strategy and ourmethod (Theorem 4).the two methods. Although the theoretical interest of our formula (27) is themain goal of this work, computing time can also be reduced related to con-ventional backprojections because of the lower computational complexity, whenfully exploiting fast matrix multiplication and fast Fourier transform algorithmsin high performance computing facilities.
We have proposed in this work a backprojection formulation for fan-beam scan-ning with standard and linear detectors. The linear case being more interestingin synchrotron context was developed as a two-step process. The first stepconsists on a sequence of simple one-dimensional operators where a standardfan-beam sinogram is obtained. The second step, the backprojection operation,is performed by writing the Fourier transform of the backprojected image asa Bessel-Neumann series on the frequency variable σ weighted by 1 /σ . Thecoefficients of the expansion are in fact the Fourier coefficients of the sinogramobtained in the first step. This approach, based on the low cost backprojection bst formula [10] for parallel sinograms, presents a reduction on the computa-tional cost related to conventional fan-beam backprojections when fast matrixmultiplication and fast Fourier transforms algorithms are used. Iterative re-construction algorithms where the backprojection process is computed on eachiteration, as the em algorith [12] can also take advantage of this formulationto obtain robust reconstructions. A second interesting feature is the absence ofa strong rebinning process from fan to parallel beam projections as is mostly13one in other fan-backprojection algorithms. When dealing with standard fan-beam data, the first step of the algorithm is more straightforward where onlyone linear change of variables is needed without lying on interpolations errors.In addition to the formal proof, numerical results support the validity of ourmethod.Once the backprojection for linear fan sinograms is efficiently performed, acomplete inversion algorithm can be implemented for cone-beam tomographywith plane detectors. The fdk formula [4] is widely used for this aim. Briefly,the cones are considered as a set sloped fans having the same source in R andparallel line detectors placed over the plane detector. A change of variables [17]is needed to write a fan in R as a sloped fan in R and backproject it. Theconvolution filter is easily derived similar to the parallel beam case. Details ofthe full reconstruction using the fdk formula are presented in [12]. References [1] Milton Abramowitz and Irene A. Stegun (eds.),
Handbook of mathematicalfunctions with formulas, graphs, and mathematical tables , National Bureauof Standards, 1964.[2] Guang-Hong Chen, Shuai Leng, and Charles A Mistretta,
A novel exten-sion of the parallel-beam projection-slice theorem to divergent fan-beam andcone-beam projections , Medical physics (2005), no. 3, 654–665.[3] Richard Anthony Crowther, DJ DeRosier, and Aaron Klug, The reconstruc-tion of a three-dimensional structure from projections and its application toelectron microscopy , Proc. R. Soc. Lond. A (1970), no. 1530, 319–340.[4] L. A. Feldkamp, L. C. Davis, and J. W. Kress,
Practical cone-beam algo-rithm , J. Opt. Soc. Am. A (1984), no. 6, 612–619.[5] AK George and Y Bresler, A fast fan-beam backprojection algorithm basedon efficient sampling , Physics in Medicine & Biology (2013), no. 5, 1415.[6] Avinash C. Kak and Malcolm Slaney, Principles of computerized tomo-graphic imaging , Society for Industrial and Applied Mathematics, 2001.[7] Daniil Kazantsev and Valery Pickalov,
Fan-beam tomography iterative algo-rithm based on fourier transform , Nuclear Science Symposium ConferenceRecord, 2008. NSS’08. IEEE, IEEE, 2008, pp. 4138–4139.[8] Fran¸cois Le Gall,
Powers of tensors and fast matrix multiplication , Pro-ceedings of the 39th International Symposium on Symbolic and AlgebraicComputation (New York, NY, USA), ISSAC 14, Association for ComputingMachinery, 2014, p. 296303.[9] D Russell Luke, James V Burke, and Richard G Lyon,
Optical wavefrontreconstruction: Theory and numerical methods , SIAM review , no. 2.1410] Eduardo Miqueles, Nikolay Koshev, and Elias S Helou, A backprojectionslice theorem for tomographic reconstruction , IEEE Transactions on ImageProcessing (2018), no. 2, 894–906.[11] Frank Natterer, The mathematics of computerized tomography , vol. 32,Siam, 1986.[12] Frank Natterer and Frank W¨ubbeling,
Mathematical methods in image re-construction , vol. 5, Siam, 2001.[13] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland,Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, War-ren Weckesser, Jonathan Bright, St´efan J. van der Walt, Matthew Brett,Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nel-son, Eric Jones, Robert Kern, Eric Larson, CJ Carey, ˙Ilhan Polat, Yu Feng,Eric W. Moore, Jake Vand erPlas, Denis Laxalde, Josef Perktold, RobertCimrman, Ian Henriksen, E. A. Quintero, Charles R Harris, Anne M.Archibald, Antˆonio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, andSciPy 1.0 Contributors,
SciPy 1.0: Fundamental Algorithms for ScientificComputing in Python , Nature Methods (2020).[14] Yuchuan Wei, Ge Wang, and Jiang Hsieh,
Relation between the filtered back-projection algorithm and the backprojection algorithm in ct , IEEE SignalProcessing Letters (2005), no. 9, 633–636.[15] Virginia Vassilevska Williams, Multiplying matrices faster thancoppersmith-winograd , Proceedings of the forty-fourth annual ACMsymposium on Theory of computing, 2012, pp. 887–898.[16] S-R Zhao and H Halling,
A new fourier method for fan beam reconstruc-tion , Nuclear Science Symposium and Medical Imaging Conference Record,1995., 1995 IEEE, vol. 2, IEEE, 1995, pp. 1287–1291.[17] Shuang-Ren Zhao, Dazong Jiang, Kevin Yang, and Kang Yang,
Generalizedfourier slice theorem for cone-beam image reconstruction , Journal of X-rayscience and technology (2015), no. 2, 157–188.[18] Shuangren Zhao, Kang Yang, and Kevin Yang, Fan beam image recon-struction with generalized fourier slice theorem , Journal of X-ray Scienceand Technology22