D'Asseler, Yves

Lemahieu, Ignace

Staelens, Steven

Van Holen, Roel

Vandenberghe, Stefaan

English

Ghent University Academic Bibliography

Fast 3D image reconstruction for rotating slatcollimated gamma camerasRoel Van Holen, Stefaan Vandenberghe, Steven Staelens, Yves D’Asseler and Ignace LemahieuAbstract—As an alternative to the use of traditional parallelhole collimators, SPECT imaging can be performed using rotat-ing slat collimators. A gain in image quality could be expectedfrom the higher photon collection efficiency, however, the planeintegral data measured by this type of collimator require a recon-struction to invert the three-dimensional (3D) Radon transform.The use of iterative methods to do fully 3D reconstruction iscomputationally much more expensive compared to a classicalSPECT reconstruction. A computationally attractive alternativeto direct 3D reconstruction could be to use two subsequent in-version steps, where the first step would consist of reconstructingthe plane integral data to conventional SPECT sinograms andthe second step would be a classical SPECT reconstruction. Inthis work, an algorithm of the Maximum Likelihood ExpectationMaximization (MLEM) type which integrates these two steps intoone algorithm is proposed as an alternative to two subsequentMLEM reconstructions. A simulation study validates the newapproach with respect to a two step MLEM reconstruction.Results indicate that the integration of two MLEM iteration loopsinto one can provide improved image quality regarding contrastrecovery and noise. Moreover it avoids determining a stoppingrule for the first step in the 2-step approach.Index Terms—SPECT, image reconstruction.I. INTRODUCTIONIN SPECT imaging, image quality is limited by the intrinsicspatial resolution versus sensitivity trade-off resulting fromthe geometric properties of parallel hole collimators, whichare traditionally used to obtain photons from the directionperpendicular to the collimator. While maintaining the sameField Of View (FOV), rotating slat collimators (RSC) providea better spatial resolution versus sensitivity compromise. Thisis due to the much higher photon collection efficiency of theplane integrals, measured by this device. Therefore, RSC canpotentially provide a gain in image quality.The collection of plane integrals of a 3D activity distributionby spinning the detector around its own axis at each regularSPECT angle can be seen as the 3D Radon transform whenignoring the effects of attenuation, depth dependent blurringand depth dependent sensitivity. For image reconstruction,a 3D algorithm which inverts the 3D Radon transform ismandatory. A Filtered Back-Projection (FBP) implementationof this 3D Radon inversion has been previously developed byLodge [1] and by Zeng [2], an iterative RBI-EM approachhas been proposed by Wang [3]. The iterative approachoffers better noise characteristics, since it models the PoissonR. Van Holen, S. Vandenberghe, S. Staelens, Y. D’Asseler and I. Lemahieuare with the Department of Electronics and Information Systems, GhentUniversity, B-9000 Ghent, Belgium e-mail: (Roel.VanHolen@UGent.be).This work was supported in part by the Institute for the Promotion ofInnovation by Science and Technology, Flanders, IWT, Belgium.statistics of the data while it can also incorporate the effects ofdepth dependent blur and sensitivity into the system matrix. Adrawback of the iterative approach is that it is computationallyvery expensive which makes it less attractive for practical use.A computationally less demanding alternative to the 3D Radoninversion is a 2-step algorithm which first reconstructs theplane integral data to conventional sinograms before recon-structing the 3D object from these sinograms using a classicalSPECT reconstruction. This approach has been proved to befeasible [4], but is less accurate than a fully 3D reconstruction.In this work we propose a new technique for reconstructingplane integral data based on an integrated iterative reconstruc-tion algorithm. The newly developed method differs from the2-step approach in the sense that both, previously separateiteration loops, are now integrated in 1 iteration loop asillustrated in Fig. 1. First, the current image estimate is forwardprojected to sinogram data, which are then forward projectedto plane integral data. Next, the ratio of the measured planeintegral data and the forward projected plane integral datais back projected to sinogram space. At this point, a backprojection to 3D image space is performed, yielding the nextimage estimate.We performed a simulation study to obtain contrast to noiseproperties of both a 2-step MLEM (2S-MLEM) reconstructionand the newly developed Integrated MLEM (I-MLEM). Thecomputation time will be discussed and compared to a fully3D algorithm.II. METHODSA. Phantom simulationThe data used in this study were generated using a MonteCarlo simulations. Once the origin and direction of a photonwas calculated using Monte Carlo techniques, the photonwas further tracked geometrically. Both the detector and thecollimator were modeled to be a perfect absorber. Thus, everyphoton striking the collimator was discarded while everyphoton arriving at the detector was detected. No attenuationor scatter was modeled, but the data are Poisson distributed.The plane integral data were simulated to contain a totalof 50 million events. A numerical image quality phantomwas simulated and consisted of a warm cylinder (diameter:180 mm, height 216 mm), containing 4 hot spheres (diameters:5.4 mm, 9 mm, 12.6 mm, 16.2 mm) and two cold spheres(diameters:19.8 mm and 23.4 mm). The sphere-to-backgroundactivity ratio was set 8:1 for the hot spots. In the simulation,the detector rotated around the phantom in 120 discrete angles,equally spread over 360◦. At each of these SPECT angles,2(a)(b)Fig. 1. In the 2S-MLEM reconstruction (a), sequence 1-2 (step A) is the first MLEM algorithm which reconstructs the sinogram data. Once the sinogramdata are reconstructed, the second MLEM algorithm is performed by iterating loop 3-4 (step B). This yields the final 3D image. In the I-MLEM algorithm(b), 1 iteration involves sequence 1-2-3-4.the detector rotated around its own axis in 120 discrete spinangles, again equally spread over 360◦. The detector itselfconsisted of 192 by 192 detection elements, each measuring1.8 mm by 1.8 mm. Parallel to the detector elements in onedirection, a slat with finite width of 0.3 mm was simulated inbetween every two neighboring detector pixel rows. The heightof the slats was set to 40 mm, while their length was equalto the length of the detector, being 345.6 mm. This collimatorconfiguration results in a collimator spatial resolution of about5 mm at 10 cm collimator distance. The SPECT rotation radiusof the detector was 18 cm.B. Image reconstructionWe used 2 different image reconstruction algorithmsbased on the Maximum Likelihood Expectation Maximization(MLEM) algorithm, namely the 2S-MLEM and the I-MLEMalgorithm. These reconstruction methods are illustrated inFig. 1. The reason for choosing MLEM is to avoid interferenceof any extra reconstruction parameters (such as the number ofsubsets when using OSEM).1) 2S-MLEM: The first step of the algorithm (step A)involves a calculation of the sinogram data, independently forevery fixed SPECT angle θ using MLEM:st+1i =sti∑j Aij∑jAijgjg¯tj, (1)withg¯tj =∑nAnjstn, (2)and where sti is the estimate at iteration t of the sinogram dataat bin i, gj are the observed plane integral data at projectionbin j and Aij is the system matrix that models the step togo from sinogram data to plane integral data. We want everySPECT angle to be treated separately, so Aij is zero if i and jdo not correspond to the same SPECT angle θ (i.e. θi 6= θj).In system matrix A, we do not include any sensitivity model-ing. This can be justified since the sensitivity in a plane parallelto the detector is quite uniform, in contrasts to the cosαdependency when a strip detector is used [4]. Moreover, theresolution which is only dependent on the source to collimatordistance, will not be modeled in this step.After performing T iterations of this first step, the obtainedsinogram is used as an input for the second step. This step (stepB) is analogous to the first except for the system matrix Bkiwhich includes both a model for depth dependent sensitivityand depth dependent blur. This matrix models the contributionfrom image voxels to sinogram elements. This is similar to aclassical 3D SPECT reconstruction, except for the sensitivitymodeling. The algorithm can be written as:f t+1k =f tk∑i Bki∑iBkisTis¯ti, (3)3withs¯ti =∑nBniftn, (4)where f tk is the object estimate at iteration t of voxel k, sTi arethe sinogram data at bin i, obtained by iterating the previousstep T times.The main problem with this method is the determination of T ,or in other words, to define a stopping rule for the first MLEMalgorithm. If T is chosen too small, we can never reach thedesired contrast in the images produced by step B. If, on theother hand, T is chosen too high, the output of step B wouldimmediately be too noisy. In this study, the stopping rule willbe derived empirically for the phantom under consideration.Step A will be stopped at different iteration points before stepB reconstruction. On the resulting images, a contrast to noiseanalysis will point us to the optimal number T .2) Integrated MLEM: The integration of both steps in oneiteration loop results in an algorithm of the formf t+1k =f tk∑i Bki∑iBkist−1i∑jAij∑j Aijgjg¯tjs¯ti, (5)withg¯tj =∑lAlj s¯tl (6)=∑lAlj{∑nBnlftn}. (7)The main difference with the previous method is that whenintegrating both iteration loops, the estimate of the sinogramdata sj at SPECT angle θ is dependent not only on gj atthat angle θ, but on the complete plane integral data. Withthis method, the data are kept more consistent in the sensethat the all data now are used to provide information on eachSPECT angle through the iteration process. In step A of the2S-MLEM, only the data measured at SPECT angle θ provideinformation on SPECT angle θ of sj .III. RESULTSA. Determination of TTo determine the number of iterations for step A in the2S-MLEM reconstruction, step A was stopped at 4 differentiteration points: for T equal to 25, 50, 75 and 100. All foursinograms were reconstructed using step B. A contrast to noiseanalysis on the largest hot lesion in the images resulting fromstep B (Fig. 2) shows that for T > 50, the noise increaseswithout substantially improving the contrast.On the other hand, we can see that the contrast for T = 25is limited because step A did not converge yet. Since thesefindings can be extrapolated to the other lesions, T = 50 ischosen for further comparison with the I-MLEM algorithm.Fig. 2. Contrast to noise plot for different values of T (for 16.2 mm diameterhot lesion). It can be seen from this plot that the best contrast noise trade-offis found for T = 50.B. Comparison of 2S-MLEM and I-MLEMThe 2S-MLEM (50 iterations of the first step) and I-MLEMreconstruction techniques are compared by means of contrastto noise plots for both the 16.2 mm diameter hot lesion andthe 23.4 mm diameter cold lesion. Results are shown in Fig. 3and Fig. 4 respectively. Every fifth iteration is plotted as adiscrete point, starting from iteration number 5. It can be seenthat for a noise level of about 40%, we need about 20 step Biterations for the 2S-MLEM while we need about 40 iterationsof the I-MLEM algorithm. At this point, the hot spot andcold spot contrast reached are respectively aboutby 8% and13% higher for the I-MLEM compared to 2S-MLEM. Similarimprovements were found for the other lesion sizes. In Fig. 5,transversal slices through the reconstructed images are showntogether with a central profiles. The noise in both imageswas matched at 40%. The slight increase in contrast recoveryappears as a lower background activity in the I-MLEM slice.The profile show that especially for the cold lesion, the contrastis better for the integrated reconstruction method.Fig. 3. Contrast noise plot for 16.2 mm diameter hot lesion. The contrastreached by the I-MLEM is about 8% higher than the contrast reached by the2S-MLEM method.4Fig. 4. Contrast noise plot for for 23.4 mm diameter cold lesion. The contrastimproves with about 13% for the integrated method.(a) (b)(c)Fig. 5. (a) and (b) respectively show a transversal slice through the 2S-MLEM reconstructed image and through the I-MLEM reconstructed image.Profiles profile drawn across the central horizontal line are shown in (c).C. Computation timeThe computation time required for one iteration of I-MLEMis about 18 minutes. This is equal to the time required forone iteration of both steps in 2S-MLEM (Fig. 6). The reasonfor this is that both algorithms use the same projector/backprojector pairs. For comparison, the reconstruction time forone iteration of our fully 3D MLEM algorithm is also plotted.It can be seen that the speed up by using a two step algorithminstead of a real 3D algorithm reduces the computation timewith two orders of magnitude.IV. DISCUSSIONReconstructing plane integral data based on a 2-step ap-proach is a computationally attractive alternative to fully 3Dreconstruction. The problems concerning the stopping rule ofFig. 6. The computation time for one iteration of 2S-MLEM and I-MLEMis shown together with the time required for one fully 3D iteration.the first step when using iterative reconstruction methods how-ever limit the flexibility of this approach since this parameter isobject and activity dependent. The integration of both MLEMloops, replacing the two separate reconstructions, cancels outthis problem. Furthermore, due to the data that are kept moreconsistent, an improved contrast versus noise balance could beobtained by integrating both steps. Since matrix Bki includesboth sensitivity and resolution modeling, the computation timeof the projections and back projections associated with itwill be dominating the computation time. Therefore, for thephantom used in this study, where we need more iterationsof the I-MLEM than for the 2S-MLEM reconstruction, the I-MLEM will take more computation time (about the double)than the 2S-MLEM method. However, looking at the hugegain in speed compared to a fully 3D algorithm, this increasein computation time for I-MLEM is small.V. CONCLUSIONIn this paper we show that a 2-step MLEM algorithm canbe successfully replaced by an algorithm which integrates bothMLEM reconstructions. Not only did we see an increase inimage quality arising from the integration, we also avoided theproblems concerning the stopping rule for the first step in a 2-step MLEM approach. The reduction of the computational bur-den associated with a fully 3D MLEM reconstruction, whichwas the motivation for this work, seemed to be successful. Areduction in computation time by two orders of magnitudecould be obtained. Future work will show how the imagequality of I-MLEM compares to that of fully 3D MLEM.REFERENCES[1] M.A. Lodge, S. Webb, M.A. Flower, en D.M. Binnie, “A prototyperotating slat collimator for single photon emission computed tomography,”IEEE Transactions on Medical Imaging, vol. 15, no. 4, pp. 500–511, 1996.[2] G. L. Zeng, D. Gagnon, C. Matthews, J. Kolthammer, J. Radachy,en W. Hawkins, “Image reconstruction algorithm for a rotating slatcollimator,” Medical Physics, vol. 29, no. 7, pp. 1406–1412, 2002.[3] W. Wang, W. Hawkins, en D. Gagnon, “3D RBI-EM reconstruction withspherically-symmetric basis function for SPECT rotating slat collimator,”Physics in Medicine and Biology, vol. 49, no. 11, pp. 2273–2292, 2004.[4] G. L. Zeng, “Detector blurring and detector sensitivity compensation fora spinning slat collimator,” IEEE Transactions on Nuclear Science, vol.53, no. 5, pp. 2628–2634, 2006.

Fast 3D image reconstruction for rotating slat collimated gamma cameras

https://biblio.ugent.be/publication/416974/file/448148.pdf

Fast 3D image reconstruction for rotating slat collimated gamma cameras

Abstract

Similar works

Full text

Available Versions

Ghent University Academic Bibliography