August 1989 7:315 研 究〉 く Comparative Studies on Various Ilerative lmage Reconstruction Algorithms for EHlission Tomography Hideo MURAYAMA*1,Eiichi TANAKA*1*2 and Norimasa NOHARA*1 ALstract: Scvcral itcrativc rccOnstruction methods apprOpriate for clnission tomography arc frcc and comparcd in rcgard to thc ratc of convcrgencc and ilnagc quality, by using thc noisc‐ noisy prttcCtiOn scts gcncratcd frbm a mathcmatical phantom using a computcr.Two modincd vcrsions Of thc silnultancous itcrativc rcconstrunction tcchniquc (SIRT) mcthOd arc intro‐ duccd as thc additivc SIRT(ASIRT)and the multiplicativc SIRT(ヽ ISIRT).ThC MSIRT mcthod is analogous to thc itcrativc spacc rcconstruction algorithln (ISRA)。 Silnulation studics prOvcd that the ヽISIRT Incthod has ttLstcr cOnvcrgcncc and bcttcr imagc quality than thc ISRA mcthod.It is suggestcd that,though thc expcctation maximization(E取 [)mcthOd givcs good cstimates for thc salnc noisy prdcctiOn set alnong all the mcthods reportcd hcrc, thc fast hcthOds such as thc iltcrcd itcrativc rcconstruction algorithm (FIRA)and thC COttugatc gradicnt method in weightcd vcrsion(CONGRWヽ り arC gOOd candidatcs for practical image rcconstruction. INTRODUCT10N Recently, iterativc rcconstruction methods havc rcccivcd considcrablc attcntion for thcir use in cmission cOmputed tomography(ECT),cSpecially in positron cmission tomography(PET) which is rapidly emerging as a vcrSatilc diagnostic tool in nuclcar medicinc. Comparcd to 1〕 analytical reconstruction methods such as the convolution backprttcCtiOn▼ (CONVO)methOd〔 and thc Fourier translorIII Incthod 〔 2〕 , iterativc mcthods havc advantages in that they arc tolcrant to incomplete sampling of prtteCtiOn data, and that it is casy to incorporatc a priori infOrmation about thc clnission distribution, such as non― pronlising methods for a stationary PET systcHェ ncgativity, ctc. ThCy are also with ine spatial rcsolution and high Sensitivity using a three dilnensiOnal detector arrangement, in which undesirable gaps among the dctcctor banks are inevitable Or the nnite distance between the attacent detcctor elemcnts 3,4〕 . resolution 〔 restricts thc spatial resolution of the image over the intrinsic detcctor It has becn alI■ ost 15 ycars sincc Gilbcrt introduced thc silnultaneous intcrative Reconロ *l Division of Physics,National lnstitutc of Radiological Scicnccs[9-1,Anagawa-4-chome, Chiba‐ shi, 260」 APAN] *2 Hamamatu Photonics K.K. (Rcccived Scp.2, 1988 and AcccPtCd OCt,7, 1988) inistry 在 This work was supportcd in part by a grant froln theヽof Education of Japan. Key words:Imagc rcconstruction,itcrativc mcthod,cmission computed tomottraphy,pOSitron cmissio tomography 7:316 EDICAL IMACING TECHNOLOGY VoI.7 No.3(1989) 1近 struction tcchniquc(SIRT)methOd〔 5〕, and COitein proposed the gradicnt(GRADY) method〔 6,7〕based On iterative least squares techniques.The cottugate gradient(CON methOd 〔 7,8〕Was alSO dCveloped in order to accelerate the rate of convergence. In the last maximization algorithm (EM)mcthOd〔9,10〕 ive ycars, thc expectationi― , introduced by /ardi〔 12-31〕 Shepp andヽ 11〕 , haS bCen studied extensively by a number of invcstigators〔 becausc it is bascd on thc rcalistic assumption that photon counts lollow a Poisson proccss. The virtues of thc EAtt method arc not only that convcrgencc to thc maxirnum likclihood ncgativity cstiIIlatc can bc prOvcn theoretically, but also that automatic inclusion of non― constraints and preservation Of tOtal image cOunts in every itcration are possible. authors〔 18,19,21,23,24〕 Since its 在mcthod have bccn dcvcloped by scvcral convergcnce ratc is iow, modined vcrsions Of the Eヽ . One silnplc mcthod is to rcducc the numbcr of vicws in thc initial itcration stcps. Another modincation is to doublc or triple the sizc of the Ettl step while still procceding in the EM direction. Lcwitt et al〔 24〕 SuggeSted accelerated algorithHls in 18〕 and KauFman〔 relaxation which thc changes to thc ilnages arc multiplied at cach iteration by an over‐ parameter. They obtained about three tilnes faster speed of convergence. approach, Tanaka〔 From a difFcrent 34〕 propOSed the altered itcrative reconstruction algorithm (FIRA) method which is based on a modined Eヽ 在 algorithm. This enhances the convergcnce spced dramatically and ilnproves thc frcquency rcsponse. method,it is well known〔 1 1,14,15,28〕that reconstructed In spite of thc virtues of thc I E脱 ilnages with the Eヽ 江 mcthod become noisy and havc large distortions ncar the edges as iterations prOcccd beyond a ccrtain point. oisc and edge In order to suppress both the ■ artifacts, ilnproved vcrsiOns of thc EMI method have recently been proposcd to include constraints with a pcnalty function 〔 23,31〕 Or tO usc a method of sicves〔 15,16,22〕 . 在1■ ethod based the Eヽ VcrklerOv and Llaccr〔 30〕intrOduccd a stopping rulc of iterations lor on statistical hypothcsis testing. As an altcrnative to the EM algorithm,the image space reconstruction algorithm(ISRA) . In terIIIs of 【 uchllchncr〔 Withcrspoon and 単 method has been proposed by Daube― 92,35〕 asymptotic theory, the ISRA mcthod provides an image which is not as good as thc maximum littclihOOd in tcrlns of precision 〔 iS advantageous since it reduces the size of largc, 33〕 , but sparsely populatcd arrays of prtteCtiOn dattt by backprttCCting them directly into a more compact imagc matrix. ising tool particularly for volume Thc ISRA method is a proH■ ilnaging, namely, three― dilnensional image reconstruction. The purposc of this papcr is to prescnt comparisons among several iterative reconstruction methods for both their algorithms and propcrties Of convergence, froIII computcr silnulation studics, and to prottidc useful infOrmation to choosc a mcthod for practical usc. mcthod as a fast algorithm is invcstigatcd, comparing with thc CONGR mcthod. modined versions of the SIRT method, which are called ASIRT and ヽ The FIRA Two 狂SIFごr, are als。 introduced in this paper. METHOD l. Reconstrucdon methods We assumc that a sct of l measurcmcnts yi(1≦ ヽ i≦ I)are availablc,wherc yi is thc nuttber August 1989 7:317 to estimate the cEliSSion prttectiOn. The basic problem iさ Of events counted in thethi‐ density xs(1≦ of」 j≦J)frOm the metturements,where j is a source pixel in a sctdements_ A sct of the J e・ liSSiOn densities xj dennes an image vector x. Let us denote the norllrlalized probability that an event4 CInitted from pixel j is assigned to prttcCtiOn i by aij. Then for each j, Σa i 5 = 1 ・ (1) i = 1 1n this paper, it is assumed that the efFects of scattering and attenuation of photons are I)are negligible, and the aij are known ttOm the geometry. A set of l elements i≦ zi(1≦ the expected values of yi, where Z i = Σa i j X 5 。 (2) 3=1 1)GRADY method〔 6,7〕 ), where lhe index n specines the This IIIcthod generates a sequcnce of sets of estilnates ・ xjく …… ・ iteratiOn,according to the following step.For each n=0,1,・ jE≦ J), , and for each j(1≦ ■ こ ), (3) )〃 1)=xj(n)十 く xj(■ β Xj(・ ・ ))名 Zxj(D=Σ i )2, Kyi一れく ノΣ(名j/σ i=l i‐ 1 th iteration, and zi(■ ) is tho expectcd where Zxj(n)is an error image at pixcl xs at thc n‐ valuc of yi, namely, る( D = Σa i j x 5 ( D (4) 1 j単 In eq.(3),β )iS a damping ttctOr which is determined so as to minimize the difFerence (・ between the measured proJectiOns yi and the estirnated proJectiOtts zi ln a least squares sensc. The function to be lninillnizcd is 2 ( x ()め ) 2 /iび 2, χ ■( め = ΣK y i 一 (5) i = 1 wherc び i is_lhe unCertainty with which yi is mcasured, and x(m)represcnts the image vector ・ ).In this method,we choseび With」 COmpOnents,xj(・ i=l for all i・ 2)CONGRコ methOd〔 7,8〕 ` This method is silnilar to the GRADY method except that the convergencc is improved by mttking an error image vcctor in each step Of iteration orthogonal to those itt the previous stepso while the arst step is taken to btt the same asin the GRADY method given by eq。 (3), the succeeding steps are given by ), )〃 ■ X j・ く x j ( lⅢ) = x j)(十 β( ・ (6) 1)ぅ ))ais/主 (D〃 x j(n二 xj(D=童 〃 i)2_γ (ai3j/σ Kyi一 ィ・ i‐ l i=1 where r(・ 'iS dCtermined so as to make all steps of iteratiOtt OrthOgonal in the sense that I J J Xj(m))/び x5(め 辞=0, Baij〃 B(透 Baij〃 )。 (透 め r 1‐l j=l j=1 て ) plays the same rolc as the d4mping ractor itt the GRADY method. (・ fOr ttl kttmj and β In thl CONGR method, we chose び i==l fOr all i. 、 デ ヽ在EDICAL IPItACING TECHNOLOGY V01.7 No,3(1989) 7:318 め CONGRW】 methOd〔 7,8〕 This lnethod is a special casc Of the weighted least squares lnethod and is run usittg the salne algorithlltl, as given by eq。 i is difFerent froln the (6), CXCept that the choice Ofび CONGR method.In the qONGRW method,び i is given by , =瞥苦携 ぃ ≧ 柱 4)EM】 method〔 11〕 The iteratiOn step for this methOd is givcn by l x j ( n + 1 ) = xΣ( j(め aijyi/4). (7) i = 1 Eqllation(7)Can be rewrittcn in thc f0110wing additive form: , x j ( n t l ) = x)j+(〃 ■x j ()・ (8) ), = Σ (yi一 ttjXj(・ 〃xs(め 孔(め )乳3Xj(D/Σ i=l j=1 where cqs。 (1)and(4)are uSed to transform eq。 (7)intO eq。 (8)。 5.FIRAコ method〔 21〕 This is a method modined fronl the EM algorithm. Each iteration step is given by n的 = C芋 ポ ;皆 ;岳 _1).β 章 壬 裾 特 )苦 戦 沖 :, Ilts, and (9) where αand βare constを C j = Σ( y i t t D /)4,( ・ ヽ i = 1 h ) s ) /( (D 孔 (yi*h)s)], U 5 = Σ[ ( y i a i(jD(*名 )*h)s)/(Zi(・ )(yi*h)s)]. As=(xjtt )Σ[(ais(Zi(・ ρ ヽ =1 二 Hcre βis a small positive constant,an4 h iSpass nlter,namely a low― h≡h(S)=(S3+2s2+3) 1(-10≦S≦10), , (10) where sis the bin number of the prtteCtiOns in the samc vicw angle ofthe prdectiOn i.Two c6nvolutions,(zi(・)*h)s and (yi*h)s, are pCrformed bn the sampling points on the s axis. The total dcnsity of thc reconstructcd imagc is nOrmalized to the total numbcr or cvcnts at every iteration. For noisy data, the correction matrix Uj fOr the high frequency component is smoothed by a nine‐ point wcighted nlter(1:2:l fOr the X and Y directions)berOrc =2 and β =0.8 in the inserting it intO (9)・ eq。 In thiS Study, we used thc paramと ters α arst itcration, but the β Was halved in the f01lowing itcrations. valuC 工〔 6)ISRAコ meth。 32〕 This method uses thc following itcratiOn step, )(め X3(n+1)=x3(・ Byャ“う/(め B4(DttD。 i=l (11) i=1 Pro」ection data are back‐proJected into inage space, and the back― proJectさ s compared d ilnagゃ ェ to an image produced by prttectiOns calculated frOm the corrcsponding itcrative reconstrucロ tion ilnage in each itcration. ヽ Thc total dcnsity Of the ilnage is norlnalizcd tO the total number of cvcnts at eヤ │ ery Augtst 1989 7:319 iteratiOn. 7)ASIRT method Wc introduce a new algorithtt givcn by ` 埼 (nt D=MaX:杓 (1)+量 (比 _名 (D)乳 ノbi,0]! (12) i=i where bi=Σaij. 1車 (13) 1 The total density oF the ilnage is normalized to the total number of events at every iteration.This algorithm is analogOus to thc additive fOrm of the Cilbert,s SIRT〔 5〕,WhiCh th ray, a4d where the number is described using components whose lettgth is that of the i“ of pixels lying along thc particular i,lh ray arc substituted for aij and bi in eq.(12). For this reason, the new method is called the additive simultancous iterative reconstruction techniquc(ASIRT). Compared with cq。 (12)iS a similar expression (8)in tho additivc form of EM, cq。 )).In ASIRT the ttctor iS proportional Zi(■ except for the factor multipling the residual(yi一 在the ttctOT is aijXj and its to ais and its nOrllflalizatiOn over all the」pittels, while in Eヽ normalizationf The rOrmer is a silnplined algebraic expression of the latter. 8)MSIRT】 method Another new algorithHl is introduced by あ 1■) 圭 (/ Dめa i j / b i 〉 x j (■ x 5 () ・ B ( y i a i j / bBi()■ i‐ l where bi is given by eq。 age is nOrlnalized to the total (13)and tOtal density Of the il■ number oF events at every iterあ tion, (14) i=1 This algorithHl is analogous to the multiplicative form Of Gllbert's SIRT〔 5〕 , and it iS thcrefore called the multiplicative simultaneous iterative reconstruction technique(MSIRT). Equation(14)iS a similar algebraic Cttpression to (11)in eq。 lsRA.At the backprdectiOn h order to correct the total de,cctiOn prObability RT uses a dividing factor, bi, ュ stage, MSI― f o r e a c h i ‐t h p r t t e C t i O n , b u t I S R A n e g l e l t s t h i S f a c t o r . This algottthIIl is derivcd by etween yi and zi as in〔 . The quantity squares distance Hleasuresむ 33〕 miniIIlizing the least‐ to be minimized with respect to the(X3)iS i乳 一 主w i ( yj=1生 jxj)2, i=l (15) matrix where the wi's are weighting ttctors.If we denote the vector(yi)by y, thc lx」 ], then this can be written as (4i3)by [A], and the diagOnal matrix (wi)by Ⅳ (y一 [A]x)TEW](yr[A]主 ), (16) where T denotes transposc. i Th1 least squares estimates,″長 of x are Obtained by solving [A]TEW][A]食 =[A]TEW]y (17) in(17)iS Written The link with the MSIRT ctt bc made by notingithat the j_th cquation a : . X r ) =yΣ i4isWi; Σa :r=lメ i(Σ i‐ i=l 1 (18) 7:320 3(1989) 1MEDICAL IMACING TECHNOLOGY Vol.7 No。 I I J l = Σ ( y i a i 5 w i () Σ / (aΣi r X r ) a i s W i ) . i=l (19) i=l r=1 ヽ在ultiplicatiOn Of bOth sides by xj provides thc following equation = yiaijwi)/(Σ 埼(ntl)=X3(D(Σ 孔(D aijWi)・ i=l (20) i=1 The weighting factor(wi)ShOuld be chosen to be about thc same sizc fOr all thc cOmponcnts in cxprcssiOn (15). If all the components of x arc equal, thcn wc should choose a対 Wi=1/Σ 。 (21) j=1 This introduces thc ヽ 在SIRT algorithm. Thc abOve assumption is reasonable when wc ch00sc a unilbrHlly distributcd source for the initial gucss of x. Thc ISRA algorithFn iS alsO derived from the samc equation,(20),for the case that wi=l for all i. 2. Simulated data The mathe甲 住atical silnulation phantom as shown in Fig。 Fig。 l B shows thc clnission dcnsity distribution. l A was gencratcd by computcr. ■is digitized in The mathematical phanto■ a64、 64 matrix and thc clnission dcnsity xs is unifrom within a square pixcl j.Wc rcstr the reconstructiOn region to a circle inscribcd in a 64× 64 Square, and chose cquispaced prtteCtiOns in which rays are prarallcl lincs with uniforIII spacing of thc salne side of a pixel frec imagc was preparcd for the silnulation study. Corrcsponding to thc A `ttruc''noisc― image vectOr T, a set of nOischfrcc prttectiOn data yl's was calculated analytically for 240 i≦I), from the following cquation.FOr cach i(1≦ cqually spaced anglcs over 180° Y i = Σt t j T j , (22) j=1 w h c r e T j i s t hleh jC― O m p o n e n t o f t h e t r u c i l n a g e v c c t o r T , J = 3 2 2 8 2a4n0d. AI =s6e8t× numbers froIIl a Of noisy prttcctiOn data was also obtaincd by generating randon■ pscudo― Poisson distribution with a mean cqual to the noisc_frec valuc of Yi. In the following silnulation study, wc uscd the threc sets of prtteCtiOn data based on the same phantoHュ : the nrst was the noisc―frce data with a tOtal count of l x 106, the Second was thc noisy data with a total cOunt of l× count of 2×105. 106, and thc third was the nOisy data with a total Logan convolution method 1〕 〔, thC reconstructed By uSing the Shepp― ilnages shown in Figs. l C, lD and E wcrc obtaincd form thc arst, sccond and third prtteCtiOn data scts, respectively. In this study, ali thc iterative reconstruction procedures were started from a unilbrln disk, namely thc initial set of estimates xj(0)is givcn by I 杓( 0 ) = Σy i / J , f O r a l l j . (23) i‐ 1 No constraints Or a priori information were used in any Of thc rcconstruction mcthOds in order to clcarly show dinお rences bctween the algorithms. Three functiOns were dcancd as measurcs of thc goodness of thc cstimatcd imagc. Thc nrst was the mean absolute error, m(X), Of thC image doaned by m(X)=童 j=1 I X j 一島 1 / 主 T j . 3=1 (24) Thc sccond was the squarc of residual nOrIIl, N2(x), Of thC prttcctiOn data deflncd by ` │ August 1989 7: 321 守( 5 ) C 2 3 讐( 0 ) 3 一 A D B Fig。 l ThC Inathclnatical phantom uscd in thc simulation(A), and itS imagc(B)・ The valucs in front of thC parcnthcscs arc diamctcrs Or distanccs in units of pixel width, and thc valucs in parcnthescs are rclativc cnlission dcnsitics. IInagcs C, D and E arc rcconstructed with thc convolutiOn backprttcctiOn oF prdectiOn data;thc noisc‐frcc data, (CONVO)mcthOd for thc three scぃ 106 and 2 x 105, and thc two noisy data scts with thc total counts of l× rcspectivcly. Thcse imagcs arc normalizcd by their total coun体 . N2(x)=Σ Kyi一 つ2. (25) i=1 The third was the log likelihood function l(x), dCancd by 1(X)=Σ (yi 10g名一名)・ (26) Thc function m(x)providcs a measure of``source convergencc",N2(x)prOVidCs a mcasure of``prttcctiOn convergence"in the sensc of lncan square error and l(x)indiCates how closely thc rcconstructcd ilnage approaches thc maxilnum likclihood cstimatce Silnulation studies allow thc cvaluation of m(X)by the direct comparison of reconstructed imagё sourcc distribution, but in practicc, thc other two functions only arc ` true imagc is usually unknown. with the truc availablc bccausc thc 3(1989) MEDICAL IMACING TECHNOLOGY Vol.7 No。 7:322 ︵ C︶X YE ” 任0任Eロ ロ,⊃ヨ0の 0 < Z < 凹 〓 ︵ 0 5 1o 15 20 40 60 oo loo NUMBER OF ITERAT10NS i n F i g . 2 P 1 0 t S o f t h c m c a n )a )b as So l au t fc u nc cr tr io or n itcration numbcr n,for thc noisc‐frcc prttectiOn sct. Thc point markcd on thc ordinate shows thc valuc of m(x)for thO m o( Fx ( t・h c CONVO imagc. RESULTS AND DISCUSS10N free data l. Convergence for ncise口 First of all,reconstruction was perforIIIcd using noisc― free prttectiOn set to delnOnstrate the dcpcndcncc of ratc of convergencc on the reconstruction method. ・2 ShOWS plots of thc F七 mcan absolute crror m(x(n))aS a function of the numbcr of interations n,for thc nrst 100 iterations with the cight difrcrent methods.The valuc of m(x)for thC CONVO image (Fig。 lC)iS l■ arkcd by ``x"on the ordinate of Fig.2. PIots of the squarc of rcsidual norェ ■N2(x(■))and the 10g likelihood l(x(n))are ShOWn in Figs. 3 and 4, respectivcly. m(X(n)) and N2(x(■ It can be sccn that, as the numbcr of iteration increascs, ))decrease monotonically, and l(x(・ ))inCrcases monofoniCally. shows thc rcconstructed ilnages at 20 itCratiOns with all mcthods. Fig. 5 As comparcd with thc CONVO image,FIRA,CONGRW and CONGR yicld superior images,while GRADY and Eヽ江provide allnOst the samc ilnage quality, and ASIRT, MSIFピ iIIlages, r and lsRA provide inferiOr We found that by the loo、 th iteration thc images of all methods are close enough to the truc image shown in F七 ・1■ ・ In ordcr tb evaluate the ratc of convergence quantitativcly, we introduced thrcc CONVO cquivalent iteration numbcrs n(m), m(N2)and n(1), WhiCh give the closest reconstruction ilnages to the CONVO ilnage by lneasurcs of the]■ norm and the log likelihOOd, respectively. ean absolute error, the square of residual Table l lists the CC)NVO equivalent iteration numbers fOr all the methods reported here. The fastest convergcnce was obtained with FIRA and the sccond ttstest with CONGRW. The Slowest is ISRA.The Eヽ than ASIRT, whilc ISRA is slower than MSIRT. ` 在 methOd is faster t Compared with GRADY, EM is infcrior in terms of both n(m)and n(N2), but Supcrior 7:323 小uguSt 1989 将︶NZ 車︵〓αOZ コ<DO︻ ︵︵c︼ の回任 ︶ 0 5 10 15 20 40 60 NUMBER OF ITERAT10NS: n orm N2(x(■ Fig.8 P10tS of the square 6f residual■ 80 100 ))as a funcdOn of the iteration number n for the noiscttfrce prdeCtiOn set.The っ fOr the point marked on the ordinate shows the vttuc of N2ぐ CONVO image. 3■5 ―free t x 10V counts(nois● 3・ ヱ 3。 ヨ叫︼ ︻コ O Oョ ︵oO一X ︶ ︵︵E X ︶ 一 い OOO工 ︻ ︶ CONVO FIRA CONGRW CONeR 3,25 0 5 10 15 20 40 60 80 100 NUMBER OF iTERAT10NS i n ))aS a function of thc iteraton Hge 4 P10tS Of the log likellh00d X(・ lぐ d On numbcr n fOr the noise"frcc pttj,cdOn set.The point ttarkゃ the ordinatc shows thc value ofo l(】 fOr the CONVO imange. in terms of n(1)s the value of■ (1)is abOut a half of n(N2)fOr EM, but n(1) iS nearly d to equal to n(N2)for GRADY. This Fenccts the fact that the EM algorithm is designも pursuc■ he車 牢 lmutt likeliho側 Ⅲo車lmぞtel‐Thatぉ ,,EM hcorporates larger,wl七 理 atchilg^rOf・ pr弱 ゃctlon,data withす N2(→ 1。 w coutttS as_comparゃ is afFected by the difFerencP tttween mPaSuredあ weight as given by cq。 htS・On the d with GRADY.The value of nd estⅢ ated pttCttiOIIS With equal (25), Whil,the matching of th1 10w count data Seriously agLcts the MEDICAL IMAGING TECHNOLOGY Vol.7 No.3(1989) 7:324 江 (E), ASIRT (C), GRADY(D), Eヽ h Rcconstructcd ilnagcs at 20 itcrations t ヽ 九 . i WC Fig. 5 FIRA(A), CONGRW (B), CONGR MSIRT(G)and ISRA(H), fOr thc nOisc― frcc prttcCtiOn sct. Table l ComparisOn of thc CC)Nヽ /O cquivalcnt itcratiOn numbcrs n(m),n(N2)and n(1),WhCrc thc rccO‐ nstructed ilnagc is closest to thc CONヽ ! /O ilnagc in Fig。 lC, by measure of thc mcan absolutc crror, the squarc of residual norm and the log 止 kelihood. CC)Nヽ ア O equivalcnt itcration number ヽIcthod n(m) n(N2) n ( 1 ) FIRA 2 3 2 CONGRW 6 6 4 CONGR 8 6 6 GRADY 25 21 20 Eヽ1 31 35 19 ASIRT 46 48 35 取ISIRT 62 57 >100 ISRA 82 86 >100 valuc of n(1)duc tO the terms of lbg zi in cq.(26). Fbrithe same r9ason, bOth MSI]RT and ISRA prOvide ttuch larger valucs Of n(1)than thc valucs of n(m)and n(N2)when compared to Eヽ在. 7=325 August 1989 。 l x 106と unts(nOisy) ISRA ASiRT MSlRT EM CRADY と_宅 ダ毛茅 ‐ : 二 10 15 20 才 40 60 6 0 ぱ ︲ 一 二伍 O Z コ< Da 中の 口 α ︶ ︵ C X ︶ Z い 田︵ 側 ︵ ︶ 6。 4・ 2・ 0 0 0 X ︶E ”EOαα回 W 卜Dコ0”口く Z< 口 Σ ︻︵E︶ l x106 counts(noに y) ISRA ASIRT MSIRT EM GRADV 80 5 NUMBER OF ITERAT10NS i n X︶一 ︵。0一X︶ ︵︵ “ aOOI︻ ヨ OOヨ ヨロ支中 C︶ 孝 80 ー ザ 戸 無 恵 0 l 10 15 20 40 60 NUMBER OF iTERAT10NS i n 5 10 15 、 20 40 60 80 100 NUMBER OF ITERATiONS:n ‐c Fig.6 P10tS of(A)thc mean absolute frror m(xく N2(x(■)), and(C)thC 10g likelihood l(x(・ n)),(B)the Square of rcsidual norm 2),rcspectively,as a functio■ oF itera‐ ムtion numbcr n for the noisy prttcctiOn sct with the total count of l x 106. 2, COnvergente of slow methOds fOr■ oisy data ln this section, we shall be concerned with the nve slow methods of GRADY, EM, ASIRT,MSIRT and ISRA only,in order to avOid complexity in the comparativゃ The others,FIRA,CONGRW ttld CONGR,will be discussed in thと discttssion. ncxt section.Inl order of convergence, the reconstruction w4s to cxalline the efFects of statistical nbisc on the ratё perforIIled fOr the nOisy prdeCtiOn set with a total count of l x 106. Ftts・ │IA,■ ,and C )) )), the square Of thl res車 立響 ■Ortt N?(x(・ show plots of the mean absolute cttor h(X(・ and the log likelihood l(X(・ )), iespectively, forithe arst 100 iterations. irhe valu● s oF In(ま ), MEDICAL Iヽ 7:‐326 在 ACING TECHNOLOGY Vol.7 Noi 3(1989) . . , 0 0 0 oisy) Σ 2 2 コ<DO中 の口伍 ︶ ¨ Xと ︵︵E︶ に 伍〇 区伍 回 ロ ト⊃コ0の 口 < Z < 口Σ X︶ Z ︵ゆ〇一X︶ ︵︵ E︶ N 2x105counts(い 5 10 15 20 40 60 NUMBER OF iTERAT10NS i n 80 5 100 10 15 20 40 60 NUMBER OF lTERAT10NS i n 80 100 2x105counts(noisy) X︶ 一 ” ︵ O一X ︶ 00 I︻ コロ支中 コ OOH o r C︶ ℃ /CONY0 ″有 5 10 15 20 40 60 80 orm Fig.7 P10tS of(A)the mean absolutl crrOr五 (x(・)),(B)thC Squarc of residual■ N2(x(■)), and (C)th1 16g likclihood l(x(n)), r6Spectively, as a:unctiott of iteraⅢ ,tiOn number n for the noisy prttectiOn set with thも total count of 2 x 105. N2(x)and l(X)fOr the CoNVO image(Fig。 lD)are marked by`牧 ''on the ordittateslof the corresponding p16tS. Although theヤ alucs or m(x(・))各 pproach zerb For noise‐ free data, Figt 6 A demonstrates altes to inlrcase after a that noisc in the proJlctiOn data prevents this attd causes thcず certain humber Of iteralions regardlも轟 tO th● recontructiott method. By 16ntrast, the ミミ ueS ))in Fig1 6 C idO not clangさ ))i■ 副ヒ・6■ an卓 =(X(・ (・ theiF血6■otonicあ bf N2(土 !behaⅣ iOr, ))4aVe almost thettame shape as r t6 the plots or Fig.3 attd 4.The plむ ts of l(x(・ 立m工今 the plots of F七 ・4, ThiS implieS that the relat市 e valdes of l(主 ))are insensitive tO (五 k l l 、 ど 、 , , , i , 1 ︰ ,・ ヽ r ト ー r ヽ 1 ユ I A l l f l ム ー ー 、4ヽ ﹂ 可 1 ,1 ヽ 1 1 11 ヽ NUMBER OF ITERAT10NS i n ― c 7:327 August 1989 ASIRT E M CRADY n=8 n=20 C2 n=100 C3 Fig.8 Rcconstructcd imagcs at 8, 20 and 100 iterations with GRADY(Al,A2,A3), Eヽ江 (Bl,B2,B3),and ASIRT(Cl,C2,C3),for thC nOisy prttcctiOn sct with the total count of 2×105. Thc top, Iliddlc and bottoln imagcs wcre obtaincd at itcration numbcrs of n=8, 20 and 100, rcspcctivcly. cxistcncc of noisc in thc prttCCtiOn data. To demonstratc a noisiCr situation, wc ran the methods for 100 iterations, using another )),N2(x(■ )) 105,Thc rcsulting plots of m(x(・ nOisy prtteCtiOn sct with the total count of 2× and l(X(・ ))are shown in Figs.7A,B and C, rcspcctively. For all the IIlethods, the ))in Fig.7 A arc obtained at a smallcr numbcr of iterations than minimum valucs of m(x(・ in Fig.6A. The monotonical bchavior of the functions N2(x(n))and l(X(n))dOCS nOt change in Figs.7B and C. Fig. 8 illustratcs thc reconstructcd ilnagcs obtaincd froHl the noisy pro」cction set with 2× 105 cOunts at 8, 20 and 100 itCrations with the GRADY,Ettl and ASIRT methods.All the eye view ilnagcs at 100 iteratiOns demonstratc ``noisc artifacts", and the reconstructed bird's― ilnages difFcr frOm thc true irnage at high iteration numbers. that both the functions l(x)and 卜 From this, it is apparcnt r2(x)are not good mcasurcs to Obtain desirablc rccons‐ truction with noisy data. Outside the boundary of thc phantom in Fig. 8, GRADY intrOduccs noisc enhanccment cvcn at low iteratiOns, duc to thc lack of non― negativity constraints, and ASIRT produccs dot-like artifacts duc to the non― lincar operation carricd out in Order to remOVe n♀ ilnage density. E脱I produces much lcss noise duc to thl automatic non‐ Onc remarkablc result from Figs. 6A and 7A is that Eヽ gatiVe negativity constraints. 狂 provides thc slnallcst valuc of m(X)among all thc mcthods for ,he same noisy prttcctiOn scte ThC minimum valuc of m(X(・ ))fOr EM is O。 129 at n=28 fOr l× 106 cOunts,and O。 188 at n=16 fOr 2× 105 cOunts. 7:328 NEIEDICAL IMAGING TECHNOLOGY V01,7 No。 l x 106 countS (n=23) 3(1989) 2x105 countS ( n = 1 6 ) E M A B Fig。9 Reconstructcd imagcs at 28 itcratiOns with Eヽ 在for the noisy prttcctiOn set with thc total count of l x lo6(A), and at 16 itcratiOns with 2× 105(B). ISRA(n=20) MSIRT(n=20) x106 counts 2x105 counts Fig。lo Reconstructcd imagcs at 20 itcrations with ISRA(Al,A2),and MSIRT(Bl, B2)。 Thc tOp and bottOm imagcs wcrc obtaincd by using thc ■ Oisy prdectiOn scts with total cOunts of l x 106 and 2 x 105, rCSpcctively. The corresponding reconstructed images arc shown in Figs。 9A a■ dB. F七 ・10 illustratcs the rccOnstructcd images at 20 iterations with ISRA and MSIRT for the noisy proJcctiOn scts with tOtal cOunts Of l× 106 and 2× 105. The ttISIFご r method has allnost the same recovery rate of spatial resOlution and the same noise propcrties as ISRA. Thc only difrerence between the two‐ methOds is that, Outside the ottcct,MSIRT has faster recovery than ISRA. This diSLrence gives an insight intO the reason why the ■ linimum valuc of m(x)fOr MISIRT is smallcr tlan that for ISRA,as shottn in both Figs.6A and 7A. Thus, we can expect that ヽ ISIFごr provides bettcr images as compared with ISRA for 7:329 August 1989 4 6 8 10 12 14 16 // 2・ 4・ い 0 0 X 〓岳 ” にO 任 伍ロ ロトコ ョOの 口 < Z< u Σ ︵︵C︶ 6・ 0 “ 硬 一 ︵︵E︺X ︶ E ” 任〇伍任 ロ ロ トD ヨ〇 切 的 く Z< 口 夏 V0 0 0 ′ c′ o′u n t s ( n o i s y ) す /′ 2 x 1 0 5 counts(noisy) 2 4 6 8 12 14 16 B A )) as a Ftlnctlon of itcratlon numbcr n for thc Plots of thc mcan absolutc crror m(x(・ and 2 x 105(B). ×1 0 6 ( A )of noisy prttcCtiOn scts with tOtal lcounts FIRA CONGRW CONGR n=2 n=4 n=6 A3 Fig. 12 18 NUMBER OF ITERAT10NS i n NUMBER OF ITERATIONS i n Fig.11 10 18 C3 B3 Reconstructcd images at 2, 4 and 6 iterations CONGRW (Bl,B2,B3)and FIRA(Cl,C2, sct with thc total count of l× 106. Thc top, obtained at itcration numbcrs of n=2, 4 and 6, with CONGR(Al,A2,A3), C3), for thC nOisy iprdCCtiOn lniddlc and botton■ imagcs wcrc rcspcctivcly. the samc noisy prttcctiOn set. 3. COnvergence of the fast methods for ncisy data are comparcd. d l n t h i s s c c t i o n , t h e t h sr te e l nf cあt h o d S O f F I R A , C O N G R t t V a n CONGR 7:330 ヽ在EDICAL I脱 CONeR IAGING TECHNOLOGY Vol.7 No.3(1989) FIRA CONGRW n=2 n=4 A2 C2 n〓 6 F 二g . 1 3 R e c o n s t r u c t c d i m a g e s a t 2 , 4 a n d 6 i t e r a t i O n s w i t h C O N G R ( A l , A 2 , A 3 ) , CONGRW (Bl,B2,B3)and FIRA(Cl,C2,C3),using the noisy prtteCtiOn sct with a total count Of 2 x 105. The tOp, Iniddlc and bottOm ilnagcs wcrc obtaincd at itcratiOn numbcrs of n=2, 4 and 6, rcspectivcly. ))fOr the nrst 18 iterations Figs.1l A and B shOw p10ts Of thc mean absolute crrOr m(x(・ using the noisy prttectiOn scts with thc tOtal cOunts of l× 106 and 2× 105, reSpectivcly. These plots demonstratc the &LSt COnvergencc of FIRA for thc noisy data as wcll as for the noise_frce data. This is bccause the frequcncy responsc of thc iterative correctiOn is signincantly ilnprovcd by incOrporating adequatc high frcquency cnhancPmCnt in the FIRA rcconstruction.A suitable number Of itcratiOns lnay be froln 2 to 4. In addition, FIRA prbved to be tt10re stable,based On the valucs of ttL(X(・ valuc had becn rcached, comparcd with CONGR and CONGRW. )),WhCn itCrating aftcr the IIlinimum This is becausc FIRA uses a selective nlter functiOn such as eq。(10), With which the cOnvergencc of each frequency component is cOntrollcd. Figs. 12 and 13 11lustrate the recOnstructed ilnages at 2, 4 and 6 iteratiOns with the threc mcthods for the twO noisy prtteCtiOn scts.It can bc sccn that FIRA prOduccs better images at a few iterations. Thc imagc Of C l with FIRA at n==2 provides a silnilar recovery of res01ution as the image OF B 3 with CONGRW at n==6, While thc fOrmer ilnagc produccs less noisc than the lattcr. This is becausc thc autOmatic non_negativity constraints involved in FIRA efFectively suppress the nOisc enhancement in areas Of 10w image density, in a silnilar way to Eヽ 在. From thc abOve cOmparison, we conclude that FIRA is a good candidate for the image reconstruction of ECT data, in practice古 From Figs.1lA and B,it is scen that CONGRW has a smallcr minimum valuc Of m ))than CONGR, while thc mhimum valucs for both methods are obtained at almost (X(・ August 1989 7:331 the same iteration number. CONGRW in Figs。 A cOmparison of the rcconstructed ilnages with CONGR and 12 and 13 ShOWS that CONGRW is superior tO CONGR due to the suppression Of the noisc in arcas outside thc ottect.ThiSis becausc CONGRW achieves better matching of prtteCtiOn data with low counts than that with high lounts by virtuc Of thc weighted lcast squares, with which thc low counts arc assigned largc weights to thc minilni‐ zing functiOn given by cq。 (5). The image quality for the fast methods is not pariticularly worse than that for the slow ))With CONGRW is o.130 fOr methods.From Fig.1l A,thc minimum valuc of m(X(・ 129)with Eヽ l× 106 cOunts, which is very c10se to thc corresponding valuc(o。 宝 from Fig. 6A.The CONGRW method secms to be practical because the minimum valuc of m(x(・ つ is achicved at n==5 iteratiOns, which is about 6 times faster compared with n==28 iteratiOns with EM.For 2× 105 cOunts, however, the minimum valuc of m(x(・ ))With CONGRW 質shown in Fig。 becomes iarger than Eヽ 在. Compared with thc imagc of E取 of CONGRW in Fig.13■ 2 iS inferior especially in areas of low image density.This is duc 9B, the imagc tO thc cxccsS enhancement of thc high― frequency component of statistical noisc, although the reconstructed ilnages are sharpened too much at early iterations. For the FIRA method, we have the frcedoIII to select the paralnetcrs aS tO ,β and theαnlter h(s)in (9),SO eq。 minimize degradation of the image quality.Further studies on FIRA will bc reported elscwhcrc〔 . 34〕 CONCLUS10N For a noisc― free prdcctiOn sct, all the methods reported here prOvide stable convegcnce; that is, the three functions, m(x), N2(x)and l(X), approach their inal valucs monotonically as the itcrations procced.Thc thrcc CONVO equivalcnt itcration numbers,n(m),n(N2)and n(1)WCrc used to evaluatc thc rate of convergcnce. From these quanfitative comparisons, the fastcst convergence can be obtained with thc FIRA mcthod. Discrepancies in the values of n(N2)and n(1)arC ObServed with some mcthodsi namely, FIRA,CONGRW,EM and ASIRT givc smaller valucs of n(1)than n(N2),while MsIRT and ISRA give larger values of n(1)・ ThiS iS beCause thc former nicthods incorporate larger weights fOr matching low count prttCCtiOns than for high count prttectiOttS,comparcd to thc latter methods. The formcr mcthods arc superior in tcrlns of noiso enhanccmcnt in the arca outsidc the ottcct, Wc intrOduced theヽ 在SIRT mcthod as a candidate for lhe recOnstruction of volume images. ed algorithm of Gilbert's SIRT and has a similar algebraic expression to ISRA. This is a modif■ It was shown that ヽ 在SIRT provides fastcr convcrgence and bcttcr imagc quality than ISRA, CSPCCially in the area outsidc the ottect. In practical usc, the FIRA mcthOd is a good choicc due to its fastcst convcrgcnce, its iow nOisc On the images with the automatic non― negativity constraints, and comparativcly stabic image quality for iterations after thc optimunl at 2-4 iteratiOns. is also useful in sOme practical cases because it provides adeq■ numbers. The CONGRW ttethod ate images at several iteratiOn 脱IEDICAL Iヽ 7:332 在AGING TECHNOLOGY Vol.7 No.3(1989) ACKNOWLEDGEMENT 在ctroロ The authors wish to thank T,Tollitani, M.Yamamoto of NIRS and H. Toyama of irokyo 、 politan lnstitutc of Gcrcntology for thcir uscful discussions. Vc予 wish to cxprcss apprcciation to K, Fukuhisa of NIRS for his tcchnical assitancc. REFERENCES 〔1〕 Shepp LA,Logan BF:Fouricr reconstruction of a head scction. IEEE Trans NucI Sci NS― 21: 21-43, 1974 〔2〕 MCrsCrcau RM:Dircct Fourier transform tcchniqucs in 3-D imagc rcconstructiono Comput Blol Med 6: 247-258, 1976 S:A mcthod FOr rcconstructing images from data obtaincd with a 〔3〕 MEcuhllchncr G, Karp」 hexagonal bar positron camera. IEEE Trans Med lmag MI-4: 134-138, 1985 〔4〕 Tanaka E,Nohara N,Tomitani T et al:Stationary positron cmission tomography and its Ied lmag MI-5: 199-206, 1986 imagc reconstruction. IEEE Trans ヽ dimcnsional rcconstruction for an ottcCt from 〔5〕 Gilbcrt P:Itcrativc m(thods for thc threc‐ prttCCtiOns.J ThCOr Bio1 86:105-117, 1972 〔6〕 Goitcinヽ I:Thrcc― dimcnsional dcnsity rccosntruction from a scrics of two‐ dimcnsional prttcC― 在cth lol: 509-518, 1972 tiOns, lWucl lnstr ヽ 〔7〕 Hucsman R,Gullbcrg GT,Grccnbcy WL,ct al: RECLBL library uscr's ma4ual; Donncr algorithms for rccOnstruction tomography. Publication PUB-214, Lawrcncc Bcrkllcy Labora― tory, 1977 〔8〕 BudingCr TF,Gullbcrg GT,Hucsman RH:Emission Computcd TomOgraphy.In: Imagc rcconstruction from prttcctiOn; Implcmcntation and application,Hcrman GT(cd.),Ncw YOrk, Spring― Verlag, pp. 147-246, 1979 〔9〕 RockmOrcメ u,M[acOvski A:A maximum likelihood approach to cmission imagc rcconstruc― tiOn frOm prdcctiOns.IEEE trans NucI Sci NS-23: 1428-1432, 1976 〔10〕 Rockmorc A」 ,M[aCOVSki A:A maximum likclihood approach to transmission image recons‐ truction from prttcctiOns,IEEE Trans NucI Sci NS-24: 1929-1935, 1977 江aximum likclihood rcconstruction for cmission tomOgraphy. IEEE 〔11〕 ShCpp LA,Vardi Y:ヽ ヽ Trans ヽ 江cd lmag MI-1: 113-122, 1982 12〕 Langc K,Carson R:EM rcconstruction algorithms for cmission and transmission tomogra― 〔 phy・」 COmput Assist Tomogr 8: 306-316, 1984 13〕 Shcpp LA,vardi Y,Ra JB ct al:` 〔 在aximum likelihood PET with rcar data.IEEE Trans NucI Sci NS-31: 910-913, 1984 〔14〕 Yardi Y,Shepp LA,KauFman L:A statistical modcl for positron emission tomography. J Amer Statist AssOc 80: 8-37, 1985 _ 15〕 SnydCr DL,ヽ 〔 fillcrヽII:Thc usc of sicvcs to stabilizc imagcs produccd with the Eヽ 江 alogr‐ ithln for c■lissiOn tomography. IEI]E Trans Ndcl Sci NS-32: 3864-3872, 1985 〔16〕 ヽ 江illCr 、在I, Snydcr DL, ヽ 在oorc Sヽ 江: An cvaluatiOn or thc usc of sicvcs ibr prOducing gorithm for PET. IEEE Trans Nucl c s t i m a t c s o f r a d i o a c t i v i t y d i s t r i b u t i o n s w i t h t h空 c aElヽ Sci NS-33: 492-495, 1986 17〕 CarsOn RE:A maximum likclihood mcthod for rcgion― 〔 graphy.」 o阜intcrcst cvaluation in cmission lomOH Comput Assit Tomogr lo: 654-663, 1986 18〕 Lcwitt RIWE, ヽ Iuchllchncr G: Accclcratcd itcrativc rcconstruction for positron cmission 〔 tomography bascd on thc Eヽ 在algorith■■for maxilnum likclihood cstimation. IEEE Trans ヽIed lmag MI-5: 16-22, 1986 :Utilization of non― negativity constraints 19〕 Tanaka E,Nohara N,Tomitani T,Yamamoto A江 〔 在cdical lmaging, in rcconstruction oF cHlission tomogralns. In: Information Processing ヽ Bacharach SL (cd。 ), Dordrccht, ヽ Iartinus Nijhor Publishcrs, pp. 379-393, 1996 20〕 Tanaka E: Rcccnt progrcss on singic Photon and PositrOn cmission tomography― 〔 from August 1989 7:333 detectors to ttgorithms――. IEEE Trans NucI Sci NS-34: 313-320, 1987 21〕 Tanaka E: A fast rcconstruction algorithm fOr stationary positron emission tomography 〔 based on a modined EM algorithm. IEEE Trans Med lmag MI-6: 98-105, 1987 22〕 SnydCr DL,Miller MI,Thomas LJ Jr:Noise and edges artittcts in maximum-1lkelihood 〔 reconstructions for emission tomography,IEEE Transヽイ cd lmag MI_6:228-238, 1987 23〕 Langc K,Bahn M,Little R:A theorctical study oF somc maximum likelihood algorithms 〔 For c■lission and translnission tomography. IEEE Trans Atted lmg MI-6: 106-114, 1987 24〕 Kaufha■ L:Implementing and accelcrating thc Eヽ f algorithm for positron emission tomou 〔 graphy. IEEE Trans Mcd lmag MI-6: 37-51, 1987 25〕 Floyd CE Jr,JaSZCZak RJ,ColCman RE:Convcrgence of thc maximum likelihood recon‐ 〔 struction algorithnl fOr cnlission computed tomOgraphy, Phys Mcd Bio1 32: 463-476, 1987 26〕 Floyd CE Jr.,」 〔 aSZCZak RJ,ColCman RE: ヽ 在aximum likclihood rcconstruction for SPECT 科たth Montc Car10 mOdclingt asymptotic bchavior. IEEE Trans NucI Sci. NS-34: 285-287, 1987 27〕 Murayama H,Tanaka E,Nohara N,ct al:A comparison oF sevcral itcrativc rccOnstruction 〔 methOds fOr ECT.JaP J Nucl Mcd 24:797-807, 1987 28〕 Llacer J,Andreac S,Vcrklcrov E et al:TOwards a practical implcmentation oF thc ヽ 〔 4LE algOrithtt for positron ellission tomography. IEEE Trans Nucl Sci NS-33: 468-477, 1986 29〕 〔 Llaccr J,VcrklCrOv E,HofFman E」 :On thC COnvcrgcnccc of thc maximum likelihOod estiロ mator method oF tomographyic image rcconstruction. LBL-21800, 1986 30〕 VerklCrOv E,Liaccr」 〔 : Stopping rule For thc MLE algorithln based on statistical hypothesis testing. IEEE Trans Med lmag MI-6: 313-319, 1986 31〕 〔 Lc宙 tan E,Hcrman GT:A mttimum oF a posteriori probability cxpcctation maximization algOrithm for image reconstruction in cmission tomography.IEEE Trans Med lmg MI-6: 185-192, 1987 32〕 Daubc― Withcrspoon ヽ IE, Muchllchncr G: An intcrativc imagc spacc reconstruction algori― 〔 thm suitablc for volumc ECT,IEEE Trans Attcd lmag MI-5: 61-66, 1986 33〕 TittCringtott DM:On the iterative imagc spacc rcconstruction algOrithm suitablc for volumc 〔 ECT,IEEE Transヽ 空 cd IInag MI-6:52-56, 1987 34〕 Tanaka E:A flltcrcd itcrativc reconstruction algorith王 ■fOr pOsitron emission tomography. 〔 In:Information Proccssing in Mcdical lmaging,dc CraaF CN,Vicrgcvcr MA(cd。 ), Ncw York, Plcnum Publishing Corp,, pp. 217-234, 1988 35〕 Daubc― WithcrsPoon` 在E, Muchllchncr G: Trcatmcnt of axial dュ 〔 PET.J Nucl Med 28: 1717-1724, 1987 ta in thrcc―dimcnsional
© Copyright 2024 ExpyDoc