Denis Mariano-Goulart

Jean-François Crouzet

Comptes Rendus Physique

/phs thats the 2Dribution.olationution.ndriques.. Il utiliseformée desionnelr autoriserC. R. Physique 6 (2005) 133–137http://france.elsevier.com/direct/COMRENPhysics/Mathematical physics, theoretical physicsA new reconstruction method in the Fourier domainfor 3D positron emission tomographyDenis Mariano-Goularta, Jean-François Crouzetba Department of nuclear medicine, Montpellier Medical University, 371, avenue du Doyen Giraud, 34295 Montpellier, Franceb Department of mathematics, Montpellier Science University, place Eugène Bataillon, 34095 Montpellier, FranceReceived 21 June 2004; accepted after revision 22 November 2004Available online 28 December 2004Presented by Maurice GoldmanAbstractThis article describes a new exact 3D reconstruction algorithm dedicated to cylindrical positron emission tomogradoes not require an estimation of missing projection data nor complex 3D interpolation procedures. This algorithm useFast Fourier Transform of non-transaxial projections to place suitable voxel values in the 3D FFT of the radioactive distThis leads to a direct fully 3D reconstruction algorithm with a limited amount of computation that requires only 1D interpprocedures and benefits from redundant projection data to improve the signal to noise ratio in the radioactive distribTocite this article: D. Mariano-Goulart, J.-F. Crouzet, C. R. Physique 6 (2005). 2004 Académie des sciences. Published by Elsevier SAS. All rights reserved.RésuméUne nouvelle méthode de reconstruction 3D dans le domaine de Fourier en tomographie par émission de positons. Cetarticle décrit un nouvel algorithme exact de reconstruction 3D dédié aux tomographes par émission de positons cyliIl ne nécessite pas l’estimation des données de projection manquantes ni d’interpolations délicates dans l’espace 3Dla transformée de Fourier 2D de signaux construits à partir des projections obliques acquises pour calculer la transFourier 3D de la distribution de radioactivité à déterminer. On obtient ainsi un algorithme de reconstruction tridimenrapide qui ne nécessite que des interpolations 1D et tire profit de la redondance des données de projection 3D pouun meilleur rapport signal sur bruit dans l’objet reconstruit.Pour citer cet article : D. Mariano-Goulart, J.-F. Crouzet, C. R.Physique 6 (2005). 2004 Académie des sciences. Published by Elsevier SAS. All rights reserved.Keywords: Positron emission tomography; 3D reconstruction; Noise; Fourier transformMots-clés : Tomographie par émission de positons ; Reconstruction 3D ; Bruit ; Transformation de FourierE-mail addresses: d-mariano_goulart@chu-montpellier.fr (D. Mariano-Goulart), crouzet@math.univ-montp2.fr (J.-F. Crouzet).1631-0705/$ – see front matter 2004 Académie des sciences. Published by Elsevier SAS. All rights reserved.doi:10.1016/j.crhy.2004.11.011134 D. Mariano-Goulart, J.-F. Crouzet / C. R. Physique 6 (2005) 133–137se of 3Doincidentive tracerred witholid angleaximume filteredtisfied fordioactiveequires arojectionationale directreover,r statisticfor eachnstructionuction orical PETithm usesioactive.edmeltor,effectiveralaryates2,1. IntroductionThe recent development of Positron Emission Tomography (PET) in functional imaging has made possible the uprojection data in the tomographic reconstruction process. As collimators are no longer necessary when recording cevents in a cylindrical PET scanner, this imaging technique provides measurements of line integrals of the radioactdistribution (LOR) that are not restricted to lie within a plane perpendicular to the axis of the detector surface. Compathe usual 2D reconstruction procedures used in Single Photon Computed Tomography, the LOR measured in a large ssurrounding the radioactive distribution are redundant so that the reconstruction algorithms that benefit from the mnumber of oblique projections are likely to produce a 3D reconstruction with a minimum amount of stochastic noise.The generalization to 3D of usual 2D reconstruction algorithms based on the Fourier Slice Theorem (such as thbackprojection algorithm) requires measurements of projection data in a complete hemisphere. This condition is not sacommercially available cylindrical PET scanners so that not all the projections that are necessary to reconstruct the radistribution are recorded. As a consequence, the generalization to 3D of the usual backprojection algorithm [1–3] rfirst estimate of the unmeasured projection data [3–6]. This estimate is generally computed thanks to the forward pof a first volumetric distribution reconstructed using 2D filtered backprojection (2D FBP) [6]. To date, the computburden of this fully 3D reconstruction algorithm, as well as the complex 3D interpolation procedures required by thimplementation of the Fourier slice theorem make these methods hardly compatible with clinical applications [7]. Mothe truncated projections are computed from a first noisy volumetric estimate of the radioactive distribution so that theicharacteristics are decreased compared to the projections that are actually recorded.A different approach uses approximations that sort the 3D projections into a 2D dataset containing one sinogramtransaxial slice to be reconstructed. The resulting sinograms are processed using conventional 2D tomographic recoalgorithms. However, these rebinning algorithms are based on approximations that limit the resolution of the reconstrbreak down when the axial aperture of the scanner increases [7].The purpose of this theoretical paper is to describe a new exact 3D reconstruction algorithm dedicated to cylindrscanners that does not require estimation of missing projection data or complex 3D interpolation procedures. This algorthe 2D Fast Fourier Transform (FFT) of non-transaxial projections to place suitable voxel values in the 3D FFT of the raddistribution. This leads to a direct fully 3D reconstruction algorithm that requires only a limited amount of computation2. Posing the problemLet us consider a cylindrical PET scanner of radiusR and lengthH , with its horizontal axis along thex3-axis (Fig. 1). Anyplane orthogonal to the axial directionx3 will be called a transaxial plane. Any pointa on the detector surface can be definby a = (R cosθa,R sinθa, za), with 0 θa < 2π and 0 za  H . A couple of photons emerging from an elementary voluof activity corresponds to a certain LOR intersecting the detector at two opposite pointsa a da′, and having an unit directionavectorω = (cosθ sinϕ,sinθ sinϕ,cosϕ), with the Euler’s angles ranging 0 θ < 2π and 0 ϕ  π .For whole body measurements, the heightH of the cylinder is much smaller than the heightL of distribution of activity inthe patient. This is the reason why it is necessary to operate a sequence of translations of the PET scanner along thex3-axis. L tus denoteλH (typically with λ = 1/4) the length of the superposition of two consecutive translations of the PET detecγandγ ′ the distances covered by the translations of the PET scanner before and after the radioactive distribution. Thefield of view of the PET scanner will be the cylindrical volumeV = {(x1, x2, x3) ∈ R3, x21 + x22  R2, 0  x3  γ + L + γ ′}(Fig. 2).If one considers a couple of opposite points(a, a′) and all the vectorsω of the unit sphere, one can make a first naturestriction on the angleϕ: 0  ϕ  π2 . Let us consider the directionω of a given coincidence ray coming from any elementdistribution of radioactivity located in the cylindrical volumeV . The invariance by translations of the 3D projection data stthat for all pointa on the detector surface, the line integral corresponding to(a,ω) is actually recorded. As shown in Fig.Fig. 1. Transverse and longitudinal projection views of a cylindrical PET scanner.D. Mariano-Goulart, J.-F. Crouzet / C. R. Physique 6 (2005) 133–137 135e 3Ded inued-ityal PETFig. 2. Three translations of a cylindrical PET scanner illustrating the definitions ofϕc , γ andγ ′.this condition leads to tanϕ  (2R)/(λH). We will demonstrate in the next paragraph that invariance by translations of thprojection data as well as a last condition stating thatϕ = π2 are necessary and sufficient to construct the algorithm proposthis article. This leads to the final conditions:arctan(2RλH)= ϕc  ϕ < π2, (1)with γ > λH andγ ′ > λH .Remark 1. In the case of focal acquisitions, no translation of the cylindrical PET is necessary and one gets the valϕc =max(arctan(2Rγ ),arctan(2Rγ ′ )).3. The Fourier-based 3D reconstruction methodLet f ∈ L2(R3) with support in{(x1, x2, x3) ∈ R3, x21 +x22 < R2, γ < x3 < L+γ } be the finite energy signal corresponing to the distribution of radioactivity. The positions of the detectorsa which record signals from the distribution of radioactivcorrespond toza ∈ [0,L + γ ] along thex3-axis. Consider the line integralPf (a,ω) measured by the pair of detectors(a, a′)in the directionω. For allθa ∈ [0,2π [, za ∈ [0,L + γ ], θ ∈ [0,2π [, ϕ ∈ [ϕc, π2 [ :Pf (a,ω) =∫Rf (a + tω)dt =∫Rf (R cosθa + t cosθ sinϕ,R sinθa + t sinθ sinϕ, za + t cosϕ)dt. (2)The change of variableu = za + t cosϕ, with cosϕ > 0, leads toPf (a,ω) =∫Rf(R cosθa + u − zacosϕcosθ sinϕ,R sinθa + u − zacosϕsinθ sinϕ,u)ducosϕ. (3)If one denotes(X1,X2) = (R cosθa − za tanϕ cosθ, R sinθa − za tanϕ sinθ), (4)Eq. (3) becomesPf (a,ω) = 1cosϕ∫ L+γγ f (X1 + u tanϕ cosθ,X2 + u tanϕ sinθ,u)du. Let us define,∀(x1, x2) ∈ R2, ∀θ ∈[0,2π [, ∀ϕ ∈ [ϕc, π2 [pθ,ϕ(x1, x2) = 1cosϕL+γ∫γf (x1 + u tanϕ cosθ, x2 + u tanϕ sinθ,u)du. (5)This leads toPf (a,ω) = pθ,ϕ(X1,X2). (6)We now look for a necessary and sufficient condition under which the set of projections recorded by the cylindricscanner makes possible the computation ofpθ,ϕ over its whole support.Let us define the line segmentLθ,ϕ = {(−za tanϕ cosθ,−za tanϕ sinθ), 0  za  L + γ }, for any given values forθandϕ (Fig. 3). A LOR can be recorded only ifθ − θa ∈]π2 ,3π2 [. Thus forθa ∈]θ + π2 , θ + 3π2 [, and forza ∈ [0,L + γ ],136 D. Mariano-Goulart, J.-F. Crouzet / C. R. Physique 6 (2005) 133–137tyrtandlputationm inFig. 3. The domainDθ,ϕ .(X1,X2) describes a reunion of semi-circles with radiusR, centered onLθ,ϕ . If one assumes that tanϕ  (2R)/(λH), thelength(L+γ ) tanϕ of Lθ,ϕ is greater than the diameter 2R of the semi-circles. Then the points(X1,X2) fill the whole domainDθ,ϕ = {x = (x1, x2) ∈ R2, d(x,Lθ,ϕ)  R}, where d(x,Lθ,ϕ) is the Euclidian distance from a pointx to the setLθ,ϕ .Consider any pointx = (x1, x2) /∈ Dθ,ϕ verifying d(x,Lθ,ϕ) > R.Thus,∀u ∈ [0,L + γ ], (x1 + u tanϕ cosθ)2 + (x2 + u tanϕ sinθ)2 > R2. As the support of the distribution of radioactivif is included in a cylinder of radiusR, whose horizontal axis is thex3-axis, the term in the integral of Eq. (5) is zero.Then(x1, x2) /∈ Dθ,ϕ ⇒ pθ,ϕ(x1, x2) = 0. Thus the values ofpθ,ϕ are known on a domainDθ,ϕ that contains the suppoof pθ,ϕ .Reciprocally, let us assume that tanϕ < (2R)/(λH), that isλH < 2R/(tanϕ). For any pair of detectorsa anda′, say forany ray(a,ω), one easily shows thatza′ − za = −(2R cos(θ − θa))/ tanϕ  (2R)/(tanϕ). As a consequence, there existsinfinite set of couples of detectorsa anda′ for which λH < za′ − za  2R/ tanϕ whose LORPf (a,ω) cannot be recordeby the PET (Fig. 1). AsX1 andX2 in Eq. (4) are continuous with respect to the angular variableθa and to the longitudinavariableza , there exists subsets ofDθ,ϕ of non-zero measure inR2 that are not reached by any(X1,X2).This completes the proof thatϕ  ϕc = arctan((2R)/(λH)) is necessary and sufficient to calculatepθ,ϕ over its wholesupport from the projection data measured in a cylindrical PET scanner. Thus the condition (1) makes possible the comof its 2D Fourier transform̂pθ,ϕ .Let us denotef̃ the Fourier transform of along the two first variables. The inversion formula of the Fourier transforL2 and the definition (5) leads topθ,ϕ(X1,X2) = 1cosϕL+γ∫γ(∫ ∫R2f̃ (ξ1, ξ2, u)e2iπ((X1+u tanϕ cosθ)ξ1+(X2+u tanϕ sinθ)ξ2) dξ1 dξ2)du. (7)Using a truncation of the integrals onR2, Fubini’s theorem and the definition of the Fourier transform inL2, one can showthat Eq. (7) leads topθ,ϕ(X1,X2) = 1cosϕ∫ ∫R2( L+γ∫γf̃ (ξ1, ξ2, u)e2iπ(u tanϕ(ξ1 cosθ+ξ2 sinθ)) du)e2iπ(X1ξ1+X2ξ2) dξ1 dξ2 (8)= 1cosϕ∫ ∫R2( +∞∫−∞f̃ (ξ1, ξ2, u)e−2iπ(−u tanϕ(ξ1 cosθ+ξ2 sinθ)) du)e2iπ(X1ξ1+X2ξ2) dξ1 dξ2 (9)= 1cosϕ∫ ∫2f̂(ξ1, ξ2,− tanϕ(ξ1 cosθ + ξ2 sinθ))e2iπ(X1ξ1+X2ξ2) dξ1 dξ2. (10)RD. Mariano-Goulart, J.-F. Crouzet / C. R. Physique 6 (2005) 133–137 137com-).ibed in athe 2Da given-isy, this-resultsbe able0 (1976)6 (1987)87–313.3D PETLet us callp̂θ,ϕ the 2D Fourier transform inL2 of the functionpθ,ϕ . As we have shown thatpθ,ϕ(X1,X2) provides acomplete knowledge ofpθ,ϕ , the uniqueness of the inverse Fourier transform leads to∀(ξ1, ξ2) ∈ R2, p̂θ,ϕ(ξ1, ξ2) = 1cosϕ f̂(ξ1, ξ2,− tanϕ(ξ1 cosθ + ξ2 sinθ)). (11)Last, let us show that for all realξ3 and(ξ1, ξ2) = (0,0), there exists a couple(θ,ϕ), with 0  θ < 2π andϕc  ϕ < π2 ,verifying ξ3 = − tanϕ(ξ1 cosθ + ξ2 sinθ).Sinceξ3 describes the interval[− tanϕ√ξ21 + ξ22 , tanϕ√ξ21 + ξ22 ] and limϕ→ π2 tanϕ = +∞, thenξ3 describesR. There-fore (ξ1, ξ2, ξ3) describesR3 except the two half-lines defined by{(ξ1, ξ2, ξ3) ∈ R3, |ξ3| > 0, ξ1 = ξ2 = 0}.As a consequence, the Fourier transformf̂ is known almost everywhere so that its inverse Fourier transform can beputed.This completes the proof thatf can be determined from the projection dataPf (a,ω) by means of Eqs. (6) and (11The practical conditions for performing the various FFTs are not discussed in the present article. They will be descrforthcoming publication.4. ConclusionIn this theoretical article, we proved that a volumetric distribution of radioactivity can be reconstructed usingFast Fourier Transform of data derived from non-transaxial projections and 1D interpolation procedures. As for(ξ1, ξ2) ∈ R2∗, the function(θ,ϕ) → − tanϕ(ξ1 cosθ + ξ2 sinθ) is surjective but not injective onR, Eq. (11) allows redundant evaluation ofˆ almost everywhere in the 3D space. As the projections acquired in an actual PET scanner are noredundancy can be averaged to improve the signal to noise ratio in the evaluation off̂ . This will provide a full 3D reconstruction algorithm in which the use of oblique projections will allow a better control of stochastic noise compared with theachieved with usual 2D reconstruction algorithms.Moreover, the use of only one 3D FFT in the reconstruction of the whole distribution ensures that this algorithm willto provide reconstructions more rapidly than other usual 3D reconstruction algorithms.References[1] S.S. Orlov, Theory of three dimensional reconstruction. I. Conditions for a complete set of projections, Sov. Phys. Crystal. 2312–314.[2] S.S. Orlov, Theory of three dimensional reconstruction. II. The recovery operato, Sov. Phys. Crystal. 20 (1976) 429–433.[3] J.G. Colsher, Fully three-dimensional PET, Phys. Med. Biol. 25 (1980) 103–115.[4] J.G. Rogers, R. Harrop, P.E. Kinahan, The theory of three-dimensional image reconstruction for PET, IEEE T. Med. Imaging239–243.[5] M. Defrise, D.W. Townsend, R. Clack, Image reconstruction from truncated two-dimensional projections, Inverse Probl. 11 (1995) 2[6] P.E. Kinahan, J.G. Rogers, Analytic 3D image reconstruction using all detected events, IEEE T. Nucl. Sci. 36 (1989) 964–968.[7] M. Defrise, P.E. Kinahan, W. Townsend, C. Michel, M. Sibomana, D.F. Newport, Exact and approximate rebinning algorithms fordata, IEEE T. Med. Imaging 16 (1997) 145–158.

A new reconstruction method in the Fourier domain for 3D positron emission tomography

https://comptes-rendus.academie-sciences.fr/physique/item/10.1016/j.crhy.2004.11.011.pdf

A new reconstruction method in the Fourier domain for 3D positron emission tomography

Abstract

Similar works

Full text

Available Versions

Comptes Rendus Physique