0 Dominique Zosso Meritxell Bach Cuadra JeanPhilippe Thiran August 24 2007 Ecole Polytechnique F ed erale de Lausanne EPFL Signal Processing Institute LTS5 atiment ELD Station 11 CH1015 Lausanne Switzerland dominiquezosso meribach jpthiran epflch Abs ID: 23779
Download Pdf The PPT/PDF document "Direct Fourier Tomographic Reconstructio..." is the property of its rightful owner. Permission is granted to download and print the materials on this web site for personal, non-commercial use only, and to display it on your personal computer provided you do not modify the materials and that you retain all copyright notices contained in the materials. By downloading content from our website, you accept the terms of this agreement.
2 1IntroductionDirectFourierreconstruction(DFR)isoneofseveralreconstructionmethodsappliedinparallel-beamcom-putedtomography[ 6 , 1 ].Anoverviewofdifferent,transform-based–incontrasttoalgebraicoriterative–reconstructionmethodscanbefoundin[ 2 ].1.1Parallel-beamComputedtomographyInparallel-beamx-raycomputedtomography(rstgenerationCT),x-raysourceanddetectoracquireparal-lelprojections,thenrotatecircularlyaroundthescannedobjecttothenextframe.Mostoftentheacquisitionisperformedslicebyslicealongtherotationalaxis(1Dsourceanddetector),severalslicescanhoweverbescannedsimultaneouslybyusinga2Dsourceanddetector.Theacquireddatapareadiscretelysampledcylindrical(in-slice)Radontransformofthescannedvolumef,asillustratedingure 1 :p(r;q;z)=ZZf(x;y;z)d(xcos(q)+ysin(q)�r)dxdy(1)andaccordinglyhavethefollowing3dimensions: 1. Theangleofprojectionq 2. Thedistancefromthecenterr 3. TheaxialpositionzFromtheRadonprojectionstheimageofthescannedvolumecanbereconstructedusingdifferenttech-niques.1.2DirectFourierReconstructionDirectFourierReconstructionuses3majorstepstoreconstructtheimageofthescannedobject(peraxialslice): 1. 1D-FFToftheprojectiontobuildapolar2DFourierspaceusingthecentral-slicetheorem. 2. PolartoCartesianresampling. 3. Inverse2D-FFTtoobtainthereconstructedslice.Thecentral-slicetheoremstatesthattheFouriertransforminrofaRadonprojectionatagivenangleisequaltotheaxialsliceatthesameangleoftheFouriertransformoftheoriginalvolume:Pq;z(kr)=Fz(krcos(q);krsin(q))(2)wherePq;z(kr)isthe1DFouriertransforminroftheprojectionacquiredatangleqandaxialpositionz,andwhereFz(u;v)isthein-plane2DFouriertransformofthevolumeataxialpositionz.Thisrelationisillustratedingure 1 . 3 Figure1:Left:Parallel-beamx-raytomography.Foreachslicez,theprojectionsp(r;q;z)aremadeupbythelineintegralsfromsourcetodetector(A!B)throughtheobject.Right:Central-slicetheorem.The1DFouriertransformPq(kr)ofaparallelprojectionp(r;q)ofanimagef(x;y)takenatanangleqgivesasliceofthe2DFouriertransform,F(u;v),subtendingthesameangleqwiththeu-axis. ThemoststrikingadvantageofDFRisitsO(N2logN)complexityimposedbytheinverse2D-FFTinstep3,whereasthemorewidelyusedlteredbackprojection(FBP)algorithmusuallyperformsinO(N3) 1 .ThemajordrawbackofDFRliesinstep2:inappropriateinterpolationduringtheCartesianresamplingmaycauseimportantimageartifacts[ 3 ].However,thesamesourcementionssomeremediestoovercometheseartifacts:Zeropaddingoftheinputdataandoversamplingofthe2DCartesianFourierspace.Thefollowingsectionsreportontheexactalgorithmanditsparameters,thecontributedclassimplementa-tions,andanexamplewithsomesynthetictestresults.2AlgorithmAformaldescriptionofourdirectFourierreconstructionmethodisgiveninalgorithm 1 .Anillustrationofthedifferentsymbolscanbefoundingure 2 .Ingeneralthecomputationsareindex-based,ratherthaninphysicalspace,witharrayindicesbeginningat0.Asrststep,thenumberofeffectivelyrequiredinputsamplesalongtheangulardimensionisdeterminedinline1.TherelativeoversizeDofthelastangularsegment(betweenthelastangularsampleand180)isdetermined(2).Thedimensionsofthezeropaddedprojectionlineand1D-DFT,aswellastheoversampled2D-DFTarecomputed(3),andtheradiallowpasscutofffrequencygivesamaximalpolarradiusbeyondwhich2DFouriersamplesaresettozero(4).Foreachinputline,dataisshiftedtothelineendsaroundtheorigin(8)andmodulatedtoshifttheDCtothecenteroftheFourierline(9).OnceallprojectionsofaslicehavebeentransformedintoFourierspace,the2D-DFTisresampled.Theoversamplingratengshortenstheeffectivepolarradius(19),whiletheangleisdeterminedusingatan2 2 (21).Theangleisthenconvertedintoaoatingpointindex(26),usedtolinearlyinterpolatebetweenthetwoadjacentangularsamples.Priorradialinterpolationisobtainedusingb-splinesofordernb(28).Notethat,iftheupperangularindexoverows,themodiedintervalsizeDapplies,theangleistrimmedandtheradiuschangessign(30).Aftertheinverse2D-FFT,therequestedregioniscollectedfromtheedgesoftheimage(36)anddemodulatedtoundotheFFT-shift(37). 1ModiedFBPversionsexist,offeringO(N2logN)aswell2cmath(math.h):atan2(y,x)returnstheprincipalarctangentofy=x,intheinterval[�p;+p]radians.Tocomputethevalue,thefunctionusesthesignofbothargumentstodeterminethequadrant. 4 Algorithm1DirectFourierReconstruction Require: CTdataasp(r;q;z).a,bandcarethesizeoftheinputdata(alongz,qandr,respectively).aistheangularrangeindegreescoveredbytheimageseries.Theinputprojectionsarezero-paddednz-times,andthe2DDFTisng-foldoversampled.fc2[0;1]istherelativeradialcutofffrequency.x0,y0andm,naretherequestedoutputoffsetandsize.B-splineinterpolationisofordernb. Ensure: Reconstructedimagef(x;y;z) 1: b0(b180b=acCoverjustlessthan180 2: D(1+180b=a�b0Sizeoflastangularinterval 3: c0(cnz;c00(c0ngZeropadandoversample 4: rmax(fcc0=2�1Radialcutofffrequency 5: forz=0toa�1do 6: forq=0tob0�1do 7: forr=0toc�1do 8: r0((r+c0�c=2)modc0Shiftdatatolineends(origin) 9: ifr0mod2=0thenModulatetoshiftDCtoDFTcenter 10: p0(r0)(p(r;q;z) 11: else 12: p0(r0)(�p(r;q;z) 13: endif 14: endfor 15: Pq(FFT(p0)Compute1D-FFT 16: endfor 17: foru;v=0toc00�1doRe-sampleCartesian2D-DFT 18: u0(u�c00=2;v0(v�c00=2Setorigintoimagecenter 19: r00(p u02+v02=ngOversample 20: ifr00rmaxthenLowpassltering 21: Q(tan�12(v0;u0)Getangle[�p;p] 22: ifQ0thenFlipsecondhalf 23: Q(Q+p 24: r00(�r00 25: endif 26: q(bQ=pConvertangleto(oat)index 27: ifbqc+1b0thenRadiallyb-splineandangularlinearinterpolation 28: F(u;v)((1�(q�bqc))bnb[Pbqc(r00+c0=2)]+(q�bqc)bnb[Pbqc+1(r00+c0=2)] 29: elseq-overow:enlargedangularskip(D) 30: F(u;v)((1�D�1(q�bqc))bnb[Pbqc(r00+c0=2)]+D�1(q�bqc)bnb[Pbqc+1�b(c0=2�r00)] 31: endif 32: endif 33: endfor 34: f0(FFT�1(F)Inverse2D-FFT 35: forx=x0tox0+(m�1),y=y0toy0+(n�1)doCopyrequestedregion 36: x0=x+c(nzng�1=2)modc00;y0=y+c(nzng�1=2)modc00Getcorners 37: if(x0+y0)mod2=0thenDemodulate(reverseFFT-shift) 38: f(x;y;z)=f0(x0;y0) 39: else 40: f(x;y;z)=�f0(x0;y0) 41: endif 42: endfor 43: endfor 3.2itk::ComplexBSplineInterpolateImageFunction6 AlphaDirection:imagedirectionalongangles 3 AlphaRange:angularrangecoveredbythedeviceindegrees(atleast180)a ZeroPadding:projectionzeropaddingratenz Oversample:2DDFToversamplingrateng Cutoff:radialcutofffrequencyfc RadialSplineOrder:radialb-splineordernbTheclassiscurrentlynotmultithread-enabled 4 andoverridesthefollowing3image-to-imageltermethods: voidGenerateOutputInformation(), voidGenerateInputRequestedRegion()and voidGenerateData().Whilethersttwomethodsallowpropagatingregioninformationalongthepipeline(providedoutputre-gionsandnecessaryinputdata,respectively),theGenerateData()containstheactualalgorithmimple-mentation.Thereadermayrefertothesourcecodeandcommentsforfurtherimplementationdetails.3.2itk::ComplexBSplineInterpolateImageFunctionAscanbeseenintheexamplesbelow,anappropriateinterpolationschemeinFourierdomainiscrucialforgoodreconstructionquality.AscurrentlytheITKinterpolateimagefunctionsareunabletodealwithcompleximages,weprovideacomplexwrapperaround itk::BSplineInterpolateImageFunction .Theitk::ComplexBSplineInterpolateImageFunctionsplitsacomplexinputimageintoitsrealandimaginarypartsandinterpolatesthemusingtheexistingscalarb-splineinterpolator.Itde-rivesfrom itk::InterpolateImageFunction andistemplatedoverTImageType,TCoordRepandTCoefficientType.Itinstantiatesa itk::ComplexToRealImageFilter anda itk::ComplexToImaginaryImageFilter ,aswellastworespectiveb-spline-interpolatorsasmemberobjects.Uponconstruction,thesememberobjectsareinitialized. itk::BSplineInterpolateImageFunction requiresthesplineordertobesetpriortoin-putimageconnection.Inconsequencethesamerestrictionsapplytothewrapperaswell.SetSplineOrderpropagatesthesplineordertotheunderlyingreal-typeinterpolators,whileSetInputImageconnectstheinputimagetothecomplexdecompositionlters,andtheiroutputtotherespectiveinterpolators.Internally,theinterpolatorswillatthispointupdatetheirinterpolationcoefcients.TheoverriddenEvaluateAtContinuousIndexmethodbecomesthenparticularlysimple,asitonlyreturnsthemergedresultsoftwounderlyingcallstothescalarmemberinterpolators: 3Inthecode,weconsistentlyusealphainsteadofthetatodenotetheangularindex,whilethetadenotestheactualangleq.4Reconstructionofthedifferentslicesisembarrassinglyparallel,andmulti-threadingwouldbeeasytoimplementifregionsweresplitalongthez-axisonly. 7 returntypenameComplexBSplineInterpolateImageFunctionTImageType,TCoordRep,TCoefficientType&#x]TJ ; .3;ȷ ;&#x-11.;镒&#x Td[;::OutputType(&#x]TJ ; .3;ȷ ;&#x-11.;镒&#x Td[;m_RealInterpolator-EvaluateAtContinuousIndex(x),&#x]TJ ; .3;ȷ ;&#x-11.;镒&#x Td[;m_ImaginaryInterpolator-EvaluateAtContinuousIndex(x));4ExampleAnexamplethatshowshowtousetheitk::DirectFourierReconstructionImageFilterclassispro-videdinDirectFourierReconstruct.cxx.Thefollowingsectionexplainstheincorporationofthelterinasmallexampleapplication.Wethenillustrateitsperformanceonasampleimage.4.1DirectFourierreconstructionexampleapplicationFirst,includesomestandardlibrariesusedforimageIO: #include"itkImage.h"#include"itkImageFileReader.h"#include"itkImageFileWriter.h"ThenincludetheDirectFourierReconstructionImageToImageFilterheaders: #include"itkDirectFourierReconstructionImageToImageFilter.h"Wedenethepixeltype.Forthesakeofprecision,itshouldberealratherthaninteger: typedefdoublePixelType;Astheimagedimensionisxedto3,thereconstructionlterisonlytemplatedovertheinputandoutputpixeltypes: typedefitk::DirectFourierReconstructionImageToImageFilterPixelType,PixelType&#x]TJ ;IJ.;Ą ;&#x-11.;镑&#x Td[;ReconstructionFilterType;Wethenderivethecorrespondinginputandoutputimagetypes: typedefReconstructionFilterType::InputImageTypeInputImageType;typedefReconstructionFilterType::OutputImageTypeOutputImageType;anddenetheIOtypes: typedefitk::ImageFileReaderInputImageType&#x-510;ReaderType;typedefitk::ImageFileWriterOutputImageType&#x-510;WriterType;Thesinogramimageisread: ReaderType::Pointerreader=ReaderType::New();�reader-SetFileName("sinogram.hdr"); 4.1DirectFourierreconstructionexampleapplication8 Nowlet'ssetupthereconstructionlter ReconstructionFilterType::Pointerreconstruct=ReconstructionFilterType::New();andconnectthesinogramasitsinput �reconstruct-SetInput(�reader-GetOutput());Theimagedirectionsaresettothecorrespondingvalues: �reconstruct-SetRDirection(r_direction);//r�reconstruct-SetZDirection(z_direction);//slices�reconstruct-SetAlphaDirection(alpha_direction);//angleZeropaddingtheinputprojectionsandoversamplingthe2DDFTaretwomeansofreducingaliasingartifactsinthereconstructedimages(settingthesefactors,aswellastheradialinputsizetopowersof2allowsforanefcientFFTcalculation). �reconstruct-SetZeroPadding(2);�reconstruct-SetOverSampling(2);Settingtheradialcutofffrequencytovaluesbelow1:0wouldallowlow-passlteringthereconstructedimages(however,mindGibb'sphenomena): �reconstruct-SetCutoff(1.0);Whiletheangularinterpolationisstrictlylinear,thedegreeoftheradialb-splinesinterpolationcanbesetbytheuser.Degree0correspondstoanearestneighborinterpolation,degree1equalslinearinterpolation.Whilevaluesupto5arecurrentlyaccepted,degree3(cubicb-spline)ismostoftenanappropriatechoice: �reconstruct-SetRadialSplineOrder(3);Finally,theangularrangecoveredbytheprojectionshastobespecied(atleast180degree,dependingontheimagingdevice): �reconstruct-SetAlphaRange(180);Afterconnectingthepipelinetotheimagewriter,thepipelineupdatecannallybeinvoked: WriterType::Pointerwriter=WriterType::New();�writer-SetFileName("reconstructed.hdr");�writer-SetInput(�reconstruct-GetOutput());try{�writer-Update();}catch(itk::ExceptionObjecterr){std::cerr"Anerroroccurredsomewhere:"std::endl;std::cerrerrstd::endl;return1;} 4.2Shepp-LoganPhantom9 Figure3:Left:Shepp-Loganphantom.Standardslicein512512pixelresolution,providedbyMatlabR .Topright:Shepp-Logansinogram.180projectionsoftheShepp-Loganphantom,steppedat1.Bottomright:SinogramwithadditiveGaussiannoise. 4.2Shepp-LoganPhantomInputdata.Themostwidelyknownreconstruction-from-projectionstestimageistheShepp-Loganphan-tom[ 4 ].Introducedin1974itisstillincommonusetodayasareferenceimageforreconstructionalgo-rithms.Weobtaineda512512resolutionimageandcalculateditsprojectionsusingthephantomandradonfunctionsprovidedbyMatlabR .Theradonprojection(sinogram)generates180projectionsfrom0to179.Aftercroppingto512180pixels(theoriginalsinogramis729180,respectingthesizeoftheim-agediagonal),twoinstances,anoise-freeandaversionwithadditiveGaussiannoise,havebeenstackedintoa51218023Dvolumetocomplywiththelterdimensionalityrequirements.ThenativeShepp-Loganphantomandthederivedsinogramisshowningure 3 .Zeropaddingandoversampling.Figure 4 illustratesthereconstructedimagesobtainedforarangeofdif-ferentzeropaddingandDFT-oversamplingcongurationsusingradialnearestneighborinterpolationonly.Imagesreconstructedwithalowzeropaddingandoversamplingratesufferfromheavyinterpolationarti-facts,thancanbemitigatedatrelativelyhighratesonly.Interpolation.Acomparisonofimagesreconstructedusingdifferentinterpolationschemesintheradialdirectionareshowningure 5 .AlthoughtheITKb-splineinterpolateimagefunctioncurrentlyacceptssplineordersupto5,agoodreconstructionqualitycaningeneralbeobtainedusingcubicsplinesonly,incombinationwithmediumzeropaddingandoversamplingrates.Agoodinterpolationschemeimprovesreconstructionqualityandallowsreducingthenzandngrates.Angularsamplingrate.Figure 6 showsthereconstructionresultsforcubicb-splineinterpolationandrel-ativelyhighnz=ng=4rates.Althoughtheoverallreconstructionlooksveryconvincing,narrowingoftheopticalwindowrevealshigh-frequencyartifactspresentoverthewholeimage.Theyareequallyreectedintheimageintensityhistogram:thesharppeaksoftheoriginalphantomarewidenedinthereconstructedim-age.ThesestreakartifactsareduetothemainweaknesscommontodirectFourierreconstructionmethods:thelimitedhigh-frequencysampledensityinpolarFourierspacedegradestheinterpolationduetoangularundersamplingwithrespecttobandwidth[ 5 ].Ingure 7 weshowreconstructionsfromdifferentnumbersofviews,illustratingtherelationbetweenstreakartifactsandangularsamplingrate. 4.2Shepp-LoganPhantom10 nzng12481 2 4 8 Figure4:Zeropaddingandoversamplingrates.Zeropaddingratesnzincreasefromlefttoright,oversamplingratesngincreasewithlines.Congurationsusingnzng�16couldnotberunduetomemorylimitations.Asignicantimprovementinimagequalitycanbeobservedforthepriceofincreasingresourcerequirements:whilezeropaddingremovestime-domainoverlap,DFToversamplingreducesinterpolationnoise. Filtering.Theangularbandwidthcanbereducedbypre-lteringtheprojectionswithalow-passsmoothinglter.WeusedaGaussiankernelalongtheradial(sic!)directiontosmooththesinogrambeforerecon-struction.Forthepriceofreducedcontoursharpness,streakartifactscanbereducedsignicantlyinthereconstructions,asillustratedingure 8 .Incontrast,smoothingintheangulardirectionhasafarmoresevereimpactoncontoursharpness,inparticularfurtherawayfromtheimagecenter,anddoesnotresultinthedesiredartifactreduction.Noise.TheimportantadditiveGaussiannoiseissimilarlyaffectingdirectFourierreconstruction,withoutanysmoothing,andthecommonlyusedlteredbackprojection(seegure 9 ).PriorradialsmoothingwithaGaussiankernelcannotentirelyremovethenoisefromtheoutputimages. 13 5ConclusionsInthispaperwereportedonourITKimplementationofadirectFouriertomographicrecon-structionmethodforparallel-beamx-rayprojections.WepresentedtheoverallstructureofdirectFourierreconstructionanddetailedthealgorithmwedeveloped.Thisalgorithmwasimplementedasa itk::ImageToImageFilter classfor3Dimageprocessing.Inaddition,ascurrentlyno itk::InterpolateImageFunction isabletodealwithcomplexdatatypes,weprovideacomplexwrapperaroundthe itk::BSplineInterpolateImageFunction .TheuseofourDirectFourierReconstructionImageFilterclasshasbeenillustratedwithasimpletestapplication.WethenprovideextensiveresultsbasedonthestandardShepp-Loganphantomimage.Wediscussthedifferentreconstructionparametersandshowtheirrespectiveimpactonthereconstructionout-come.Weshowhowinputzeropaddingreducessignaldomainreplication(aliasingtails)and2D-DFToversamplingreducesthedishingeffectofFourierdomaininterpolation(cyclicconvolution).ChosinganappropriateradialinterpolationsplineorderforFourierdomainresamplingfurtherimprovesthereconstruc-tionquality.Wenallycomparethereconstructionresultstoastandardlteredback-projectionmethodprovidedbyMatlabR ,conrmingthehighimagequalityofourDFRimplementation.References [1] D.Gottlieb,B.Gustafsson,andP.Forss´en.OntheDirectFourierMethodforComputerTomography.IEEETransactionsonMedicalImaging,19(3):223–232,march2000. 1 [2] R.M.Lewitt.ReconstructionAlgorithms:TransformMethods.InProceedingsoftheIEEE,volume71,pages390–408,march1983. 1 [3] M.Magnusson,P.-E.Danielsson,andP.Edholm.ArtefactsandRemediesinDirectFourierTomo-graphicReconstruction.NuclearScienceSymposiumandMedicalImagingConference,2:1138–1140,october1992. 1.2 [4] L.A.SheppandB.F.Logan.TheFourierreconstructionofaheadsection.IEEETransactionsonNuclearScience,21:21–43,1974. 4.2 [5] H.Stark.Samplingtheoremsinpolarcoordinates.JournaloftheOpticalSocietyofAmerica,69(11):1519–1525,nov1979. 4.2 [6] H.Stark,J.Woods,I.Paul,andR.Hingorani.DirectFourierreconstructionincomputertomography.IEEETransactionsonAcoustics,Speech,andSignalProcessing[seealsoIEEETransactionsonSignalProcessing],29(2):237–245,April1981. 1