- Email: [email protected]
Simulation studies on depth of interaction effect correction using a Monte Carlo computed system matrix for brain positron emission tomography
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 820–831
journal homepage: www.intl.elsevierhealth.com/journals/cmpb
Simulation studies on depth of interaction effect correction using a Monte Carlo computed system matrix for brain positron emission tomography Nam-Yong Lee a , Yong Choi b,∗ a b
School of Computer Aided Science, Institute of Basic Sciences, Inje University, Gimhae, Gyeongnam 621-749, Republic of Korea Department of Electronic Engineering, Sogang University, 1 Sinsu-dong, Mapo-gu, Seoul 121-742, Republic of Korea
a r t i c l e
i n f o
a b s t r a c t
Article history:
The parallax errors caused by a lack of depth-of-interaction (DOI) information often degrade
Received 15 July 2011
reconstruction quality in positron emission tomography (PET). To reduce parallax errors,
Received in revised form
some PET systems employ multi-layered detectors to provide detailed DOI information,
24 April 2012
but this approach requires more complicated detector configurations and signal processing
Accepted 7 May 2012
schemes. In this paper, we conduct simulation studies on an inherited DOI effect cor-
Keywords:
maximum likelihood expectation maximization (MLEM) iterations using the Monte Carlo
rection of a brain PET with a single-layered detector. For this purpose, we compare the Positron emission tomography
computed system matrix of the single-layered detector with filtered backprojection (FBP)
Depth-of-interaction
reconstructions from projection data obtained from multi-layered detectors. We also inves-
Parallax error
tigate the benefits of multi-layered detectors in MLEM iterations using simulated data. The
Spatial resolution
quantitative comparison in this paper shows that an inherited DOI effect correction of a
Image reconstruction
Monte Carlo computed system matrix for a single-layered detector results in the associ-
Maximum likelihood expectation
ated MLEM iterations outperforming FBP reconstructions even for the case in which the
maximization
projection data for FBPs are obtained from an octuple-layered detector. It also shows that
System matrix
use of multi-layered detectors provides better results overall in MLEM reconstruction, but the improvement seems not to be substantial enough to ignore the complexity and costs required for multi-layered detectors. Based on these results, we conclude that detailed DOI information from multi-layered detectors is favorable, but unnecessary in brain PET imaging, because the inherited DOI effect correction via a Monte Carlo computed system matrix for a single-layered detector is sufficient. © 2012 Elsevier Ireland Ltd. All rights reserved.
1.
Introduction
Ideally, a PET scanner needs to have high and uniform spatial resolution across the imaging volume. A PET scanner equipped with elongated detectors or small detector rings is, however, susceptible to parallax errors because of event mispositioning caused by detection of gamma-rays incident at oblique angles on the detector. Until recently, the reduction
∗
of parallax errors has been addressed either by making the detector ring much larger than the phantom (the object to be imaged) so that the incidence angles of gamma-rays emitted from phantoms inside the field of view (FOV) remain sufficiently small or by using a detector module that can produce more detailed DOI information. More detailed DOI information can be obtained by various methods [1–5]. For example, in [1], they used a statistics-based positioning technique to obtain detailed DOI information
Corresponding author. Tel.: +82 10 9933 2624; fax: +82 2 706 4216. E-mail addresses: [email protected], [email protected] (Y. Choi). 0169-2607/$ – see front matter © 2012 Elsevier Ireland Ltd. All rights reserved. http://dx.doi.org/10.1016/j.cmpb.2012.05.007
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 820–831
821
Fig. 1 – Brain PET system consisting of 72 detector modules. Each detector module has a 4 × 4 arrays of Geigermode avalanche photodiode (GAPD) and lutetium oxyorthosilicate (LSO).
from continuous crystal detectors. In Refs. [2,3], a doublelayered detector head consisting of lutetium oxyorthosilicate/lutetium yttrium aluminum perovskite (LSO/LuYAP) phoswich crystals is used, and in [4] lutetium oxyorthosilicate/lutetium yttrium orthosilicate (LSO/LYSO) phoswich crystals are used to obtain binary DOI measurements. In [5] the authors used a quadruple-layered detector and proposed a reduction method to reduce heavy computations caused by the system matrix for the quadruple-layered detector. Obviously, the primary drawback of these approaches is that more detailed DOI information requires more complicated detector configurations and signal processing schemes. Statistical iterative methods may help to reduce the parallax errors in reconstructed images by incorporating a model of the system response into the system matrix. In the absence of noise, the imaging process can be modeled as Y = Pf, where the vector f represents the voxelized image to be reconstructed, the vector Y the measured projection data, and the system matrix P represents the relationship between f and Y. The system matrix P depends on how the projection data Y is arranged. In this work we assume that each element in Y represents the number of annihilation events detected at a crystal pair. Thus the matrix entries Pb,v can be defined as the conditional probabilities of detecting an annihilation event at a crystal pair b that was emitted from the voxel v. They are usually calculated by computing the intersection of the line of response (LOR) defined by a pair of detectors for each voxel [6,7]. More accurate approaches take into account the position of the voxel relative to the detector pair [8] or use analytical models or empirical kernels to compensate for crystal penetration effects [9,10]. It is also possible to compute the system matrix by measuring the response of the real system [11] or by using a Monte Carlo method to simulate the response of the real system [12].
Fig. 2 – Histograms of the reference phantom, which is a uniform phantom covering the entire imaging volume. Here we assume that the PET consists of 72 detector modules, and each detector module has 4 × 1 crystal array. The row and column indices in histogram figures indicate indices of crystal pairs that detect true coincidence events (here single and random events were excluded). Thus, the value at (i, j) pixel of the histogram is the number of true coincidence events detected at the pair with the ith and the jth crystals. (a) Monte Carlo computed histogram of the reference phantom. Here the number of layers is set to one. The total number of true coincidence events is 31.37 millions. (b) Histogram obtained by applying rotational symmetry to the original Monte Carlo computed histogram in (a). Here, 72 rotational averages are applied for k/36, k = 0, 1, . . ., 71, by using the fact that 72 crystal modules were located at equi-angles in a ring-shaped PET system.
The main objective of this work is to evaluate the DOI effect correction of a Monte Carlo computed system matrix of a brain PET with a single-layered detector. We conducted simulation studies to compare MLEM iterations using the system matrix generated by a Monte Carlo computation for a singlelayered detector with FBP reconstruction from projection data obtained by multi-layered detectors. We also investigated the
822
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 820–831
Fig. 3 – Procedure for FBP reconstruction for a double-layered PET detector. The images in the first column are data sets obtained by double-layered detectors. From top to bottom, (1,1), (1,2), (2,1), and (2,2) layer-paired data set are shown, where numbers inside parentheses are layer indices. The scale compensation, rebinning, and interpolation are performed for each layer-paired data set. The FBP reconstruction is applied to the sum of four sinograms obtained from four layer-paired data sets.
benefits of multi-layered detectors in MLEM iterations using simulated data.
2.
Materials and methods
2.1.
PET system
We have been developing a brain PET using 72 detector modules in a ring-shaped scanner whose inside diameter is 330 mm. Each detector module has a 4 × 4 array of Geigermode avalanche photodiode (GAPD) and LSO, where one crystal is 3 mm × 3 mm × 20 mm LSO (the crystal pitch is 3.3 mm). To increase the axial field of view, the final form of the brain PET system will have multi rings of 72 detector modules. For details, see [13] or Fig. 1. In simulation for DOI effect correction, however, we assume that the detector module has 4 × 1 (instead of 4 × 4) crystal array to restrict our attention to two-dimensional reconstruction only. Such restriction is considered to have a fair comparison between FBP and MLEM in DOI effect correction, since FBP is known to be less accurate than MLEM in three-dimensional reconstructions [14]. It is also true that twodimensional reconstruction is sufficient for comparison in DOI effect correction. Having these facts in mind, we assume that the phantom consists of one slice image of 175 × 175 voxels, where the size of the voxel is 1.7 mm × 1.7 mm × 1.7 mm. Thus, the diameter of the imaging volume (the volume of the largest phantom to be reconstructed) is 297.5 mm (175 mm × 1.7 mm) and the height of the imaging volume is 1.7 mm. We simulated detectors to have multi-layers from 1 to 8 for comparison studies. For instance, the double-layered detector had two 3 mm × 3 mm × 10 mm LSO crystals and
generated a binary DOI measurement. We assumed that the single-layered detector had a 20 mm long layer and generated a unary DOI measurement. A similar rule also holds for multi-layered detectors. For instance, each quadruple-layered detector had a layer 5 mm long and generated a quadruple DOI measurement. We assumed that there are no obstacles to photon passage between layers. Thus the multi-layered detector
Fig. 4 – FBP reconstructions for the ‘nine points’ phantom: (a) image by FBP1 , RRMS = 0.696, (b) image by FBP2 , RRMS = 0.572, (c) image by FBP3 , RRMS = 0.496, (d) image by FBP4 , RRMS = 0.452.
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 820–831
823
Fig. 5 – Center row profiles of FBP reconstructions for the ‘nine points’ phantom: (a) profile by FBP1 , (b) profile by FBP2 , (c) profile by FBP3 , (d) profile by FBP4 . Here, rectangular plots are the profiles of the center row of the original ‘nine points’ phantom. We used a normalization method that sets the maximum and minimum intensities of the image to 1 and 0, respectively.
provided more detailed DOI information than the singlelayered detector, while maintaining the same sensitivity as the single-layered detector. Our prototype system is designed for brain imaging and has a single-layered detector. The development of PET systems with multi-layered detectors is not part of our current research interests. We consider multi-layered detectors in our simulation studies to investigate the benefits of multi-layered detectors in FBP and MLEM reconstructions.
2.2.
MLEM reconstruction
The MLEM reconstruction takes the following iterations [15]: fˆvn+1 = fˆvn
P Sn Y b b,v b , Snb = b . P Pb,v fˆvn b b,v
(1)
v
Here fˆ n represents the approximation of the image f at the nth iterates, fˆvn the activity value of fˆ n at the voxel v, Yb the
observed data at the projection bin b, and Pb,v are entries of the system matrix P that determines the relationship between the unknown true activity f and the noiseless projection data P f . One can regard Y as observed data y defined by yb = v b,v v whose statistical mean is y. The projection bin b in (1) is determined by how the observed events are arranged. Rather than working directly with the raw data, it is common to rebin or interpolate the acquired data into a predetermined set of LORs that are convenient for FBP reconstruction. In our work, however, we used the raw dataset to maximize the DOI effect correction of a Monte Carlo computed system matrix. Thus, a projection bin b in (1) represents a crystal pair (c1 , c2 ). We also used a wavelet smoothing technique in MLEM iterations for noise reduction [16]. In this work, we computed the system matrix by using the method described in [12]. First, we used the GATE program [17] to simulate the system response for the reference phantom (a uniform phantom that covers the entire imaging
824
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 820–831
Fig. 7 – MLEM reconstructions for the ‘nine points’ (15)
phantom: (a) image by MLEM1 , RRMS = 0.174, (b) image by (100) MLEM2 ,
(100)
RRMS = 0.101, (c) image by MLEM3 (100)
RRMS = 0.061, (d) image by MLEM4
Fig. 6 – Plots of DOI effect correction by FBP reconstructions: (a) DOI-R and (b) DOI-L.
volume) without including single or random events. Fig. 2a shows the histogram of true coincidence events detected at each crystal pair. Throughout this work we use the term histogram to indicate the data format shown in Fig. 2a only. One pixel in Fig. 2a represents one crystal pair, and the intensity at a given pixel represents the number of true coincidence events detected at the corresponding crystal pair. For each event, we recorded the origin of the event for use in the system matrix generation. Second, we applied rotational symmetry. We used 72 detector modules in a cylindrical ring-shaped PET system. Exploiting this rotational symmetry, we rotated the origin and crystal pair that observed the event by k/36, k = 0, 1, . . ., 71, and we recorded rotated origins and rotated crystal pairs. Fig. 2b shows the rotationally symmetrized histogram of that in Fig. 2a. Finally, we formed the system matrix using the recorded information (origins of events and crystal pairs that detected events).
2.3.
FBP reconstruction
The FBP reconstruction can be explained as the inversion of the Radon transform, which is represented by line integrals.
,
, RRMS = 0.034.
Due to small gaps (approximately 1.5 mm) between modules and the inter-crystal penetration effect, however, line integrals connecting crystal pairs depend not only on the phantom distribution across the line but also on the relative position between the crystal pair and phantom. The response of crystals in the inner part of the module is very different from the response of crystals at the boundary of the module. This phenomenon explains square shaped patterns in Fig. 2 and decreased sensitivities of crystal pairs in which the angle difference in the ring-shaped PET is close to 180◦ . To compensate for such irregularities, we used the following scale compensation: Yb1 = Yb0 ·
RT b , MCb
(2)
where Y0 is the observed histogram of the phantom of interests; RT is the Radon transform of the reference phantom; MC is the Monte Carlo computed histogram of the reference phantom; Y1 is the scale-compensated histogram. We computed RT by measuring the intersection between the ray that connects the crystal pair and the phantom using the standard ray tracing technique [18]. Here the ray that connects the crystal pair is defined by the line that connects the center points of crystal faces. After scale compensation in (2), we used 160 × 281 bins to rebin the data at the histogram, where the numbers of angles and radial points are 160 and 281, respectively, and used linear interpolation for data completion. The rebinned and interpolated data will be referred to as a sinogram in this work. The standard FBP algorithm can be applied to the sinogram to complete the reconstruction. For noise reduction, various filtering methods can be applied [19]. In our simulation, we used the Hamming filter with a cut-off frequency of 0.5.
825
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 820–831
(15)
Fig. 8 – Center row profiles of MLEM reconstructions for the ‘nine points’ phantom: (a) profile by MLEM1 , (b) profile by (100) MLEM2 ,
(100) MLEM3 ,
(100) MLEM4 .
(c) profile by (d) profile by Here, rectangular plots are the profiles of the center row of the original ‘nine points’ phantom. We used a normalization method that sets the maximum and minimum intensities of the image to 1 and 0, respectively.
The FBP reconstruction can be applied to the histogram from a multi-layered detector by using scale compensation, rebinning, and interpolation for each layer-paired data. For instance, a double-layered detector produces four layer-paired data sets, e.g., (1,1), (1,2), (2,1), and (2,2) layer-paired data sets, where the (1,2) layer-paired data set is obtained by crystals at the first layer and crystals at the second layer. Fig. 3 illustrates the process of FBP reconstruction for a double-layered detector: the scale compensation, rebinning, and interpolation are performed in order to produce the layer-paired sinogram, and FBP reconstruction is applied to sum of layer-paired sinograms.
2.4.
Criteria for evaluation
For simulation studies, we used ‘nine points’, ‘uniform’, ‘hot rods’, and ‘cold rods’ phantoms. As mentioned in Section
2.1, simulation studies were performed for two-dimensional image reconstruction and phantoms formed by one slice images of 175 × 175 voxels were used in simulation. We generated histograms of the given phantoms by using the GATE program [17]. We excluded single- or random events in observation, and did not consider attenuation and scattering effects inside phantoms. We only considered true coincidence events as observed data. The numbers of true coincidence events used for the ‘nine points’, ‘uniform’, ‘hot rods’, and ‘cold rods’ phantoms were 0.028, 5.27, 1.33, and 5.47 millions, respectively. These numbers were obtained by letting GATE program to simulate the process of observation for 10,000 s. We generated histograms just one time for each phantom, and used the coordinates of the detection location to compute the multi-layered histograms. Thus, a multi-layered histogram had exactly the same number of photons as a single-layered histogram, but contained detailed DOI information.
826
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 820–831
where Ci is the center location of the ith point in the reconstructed image and pi is the center location of the ith point in the original ‘nine points’ phantom. The quantities DOI-R and DOI-L are used to measure how insufficient DOI information affects reconstruction in resolution and location. RRMS: The RRMS is defined by RRMS =
|fˆv − fv |2 v
|fv |2
,
(5)
where fˆv and fv are intensity values of the reconstructed image and the original phantom image, respectively, at the voxel v.
3.
Results (n)
We used FBPk and MLEMk to denote FBP and MLEM, respectively, where the subscript k stands for the number of detector layers used in generating the histograms and the superscript n inside the parenthesis stands for the iteration number. With this notation, what we wish to compare is the performance of MLEM1 with those of FBPk and MLEMk , k > 1.
3.1.
Fig. 9 – Plots of DOI effect correction by MLEM reconstructions: (a) DOI-R and (b) DOI-L.
We used DOI effect correction and the relative root mean square (RRMS) error in evaluating performance of the reconstruction methods. DOI: For DOI effect correction evaluation, we used the ‘nine points’ phantom. For details, see Fig. 4, which shows reconstructed images of the ‘nine points’ phantom by FBP. Each point in the ‘nine point’ phantom consists of 3 × 3 voxels, thus its support size is 5.1 mm × 5.1 mm because the voxel size is 1.7 mm × 1.7 mm × 1.7 mm. Let Fi , i = 1, 2, . . ., 9, be the full width half maximum (FWHM) of the ith point in the ‘nine points’ phantom. We define DOI-R as DOI-R =
F1 + F9 . 2F5
(3)
Note that F1 and F9 are FWHMs of points that are the farthest points from the center, and F5 is the FWHM of the point at the center. We also define DOI-L as
DOI-L =
9 i=1
|Ci − pi |2 ,
(4)
Comparison by DOI effect correction
Fig. 4a, b, c, and d shows FBP reconstructions of the ‘nine points’ phantom using single-, double-, triple-, and quadruple-layered detectors, respectively. As expected, the image quality improves as the number of layers increases. The improvement is, however, slightly noticeable in the change from the single to double-layered detector but not noticeable in other changes. Fig. 5a, b, c, and d shows profiles of the center rows of Fig. 4a, b, c, and d, respectively. Here rectangular plots are the profiles of the center row of the original ‘nine points’ phantom. In Fig. 5, we used a normalization method, which sets the maximum and minimum intensities of the image to 1 and 0, respectively. In Fig. 5, there is some improvement in resolution when the number of layers increases from 1 to 2, but other changes (such as from 2 to 3 or from 3 to 4) do not noticeably improve the resolution upon visual improvement. To examine the relationship between the number of layers and the DOI effect correction of the FBP reconstruction more clearly, we extended the FBP simulation up to an octuplelayered detector. Fig. 6 shows the DOI effect correction of the FBP reconstruction for the single- to octuple-layered detectors. As seen in Fig. 6, the improvement in DOI effect correction due to the multi-layered detectors is consistent with the DOIL criterion. The DOI-R is improved when the number of layers increases from 1 to 2, but it remains virtually the same for other changes. The FWHMs F5 of the point at the center in Fig. 5 were nearly unchanged for all multi-layered detectors. This result shows that the detailed DOI information from the multilayered detector does not improve the resolution in the center area for the FBP reconstruction. Conversely, F1 and F9 , FWHMs of the furthest points from the center, can be closer to the true ones than F5 due to the detailed DOI information from multilayered detectors. This is the reason why we had DOI-Rs that were smaller than 1 in Fig. 6a.
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 820–831
827
Fig. 10 – Plots of RRMS errors by FBP reconstructions: (a) ‘nine points’, (b) ‘uniform’, (c) ‘hot rods’, (d) ‘cold rods’.
Fig. 7a, b, c, and d shows MLEM reconstructions of the ‘nine points’ phantom using histograms by single-, double, triple-, and quadruple-layered detectors, respectively. Here, we selected images with the smallest RRMS for each simulation. The improvement in image quality attributed to detailed DOI information from multi-layered detectors is not clearly visible in Fig. 7. Fig. 8a, b, c, and d shows profiles of the center rows of Fig. 7a, b, c, and d, respectively. Here, again, as in Fig. 5, rectangular plots are the profiles of the center row of the original ‘nine points’ phantom, and the same normalization method is used. From Fig. 8 we can see that MLEM iterations using the Monte Carlo computed system matrix outperform FBP in DOI effect correction. In fact, MLEM iterations for the singlelayered detector outperform FBP reconstructions even when the projection data for FBP are obtained by the octuple-layered detector. This fact can be verified by comparing the results in Figs. 6 and 9. Fig. 9 shows the DOI effect correction of MLEM iterations for single- to quadruple-layered detectors. As seen in Fig. 9, the benefit of detailed DOI information from multi-layered detectors is not noticeable.
3.2.
Comparison by RRMS error
Fig. 10 shows RRMS errors for FBP reconstructions from multilayered detectors (the number of layers ranges from 1 to 8). As expected, the RRMS decreases as the number of layers increases. The decreasing rate is, however, noticeable only when the number of layers is less than 4. A similar phenomenon also holds for the ‘hot rod’ and ‘uniform’ phantoms. Fig. 11 shows images by FBP reconstructions of the ‘cold rods’ phantom using single-, double-, triple-, and quadruple-detectors. Some visual improvements can be seen only when fewer than four layers are used. Fig. 12 shows RRMS errors for MLEM iterations using multilayered detectors with respect to iteration numbers. Again, as expected, RRMS decreases as the number of layers increases. The improvements in RRMS due to the use of multi-layered detectors are, however, insignificant if we consider that the computational cost for one loop of MLEMn is approximately n2 times of that of MLEM1 [5]. Fig. 13 shows images from MLEM iterations of the ‘cold rods’ phantom using single-, double-, triple-, and quadrupledetectors. Here we selected images with the smallest RRMS
828
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 820–831
Fig. 11 – FBP reconstructions for ‘cold rods’ phantom: (a) image by FBP1 , RRMS = 0.144, (b) image by FBP2 , RRMS = 0.109, (c) image by FBP3 , RRMS = 0.091, (d) image by FBP4 , RRMS = 0.082.
among 100 iterations. In Fig. 13, visual improvements due to the increase of detector layers are also unclear in MLEM iterations.
4.
Discussion and conclusion
In this paper, we investigated the benefits of multi-layered detectors in a prototype brain PET using 72 detector modules, in which one detector module consists of a 4 × 4 GAPD crystal array and the crystal size is 3 mm × 3 mm × 20 mm. The prototype PET is a ring-shaped scanner with a diameter of 330 mm. The simulation results related to the FBP reconstructions (Figs. 4–6, 10, and 11) showed that the use of up to three multilayered detectors improved image quality in reconstruction, but more than three detector layers did not give noticeably better improvements. It was also true that the improvement in visual quality by the use of multi-layered detectors was hardly seen in FBP reconstructions after three detector layers (see Figs. 4 and 10). We believe that such phenomenon is partially caused by the inter-crystal scattering. It is true that the improved DOI effect correction using multi-layered detectors often does not
give the correct path of the event because of the inter-crystal scattering effect. For instance, in our simulation with the cold rods phantom, 42% of true events are affected by the intercrystal scattering effect. The inter-crystal scattering effect produces a kind of blurring on the histograms (and consequently sinograms) in Fig. 3. The FBP reconstruction used in this paper does not have a procedure to reduce such blurring, and hence cannot fully utilize the benefit of multi-layered detectors for more than three detector layers. We also believe that the fact that most events are detected at crystals in frontside layers causes the limited efficiency of multi-layered detectors. To explain this, let us consider Fig. 3, where images in the first column showed histograms obtained by double-layered detectors. Specifically, the histogram in the first row showed counts of events detected at crystals in the first layer, which will be denoted by a (1,1) histogram. Similarly, histograms in the second, third, and fourth rows showed counts of events detected by crystals in (1,2), (2,1), (2,2) layered pairs, respectively, where the numbers inside parentheses are the layer indices. As we can see in Fig. 3, counts of events from the (2,2) layer-paired crystals are much smaller than those from (1,1). To be specific, counts in (1,1), (1,2), (2,1), (2,2) layer-paired crystals are 34.2%, 24.2%, 26.5%, 15.1% of the total
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 820–831
829
Fig. 12 – Plots of RRMS errors by MLEM reconstructions: (a) ‘nine points’, (b) ‘uniform’, (c) ‘hot rods’, (d) ‘cold rods’.
count. In the double-layered detector, detectors in the second layer are still important because they are used to form (1,2) and (2,1) histograms, which contain a substantial amount of data as seen in the first column in Fig. 3. But if we use quadruplelayered detectors, the role of the fourth layer will be almost negligible, since only few events will be detected at crystals in the fourth layer. For instance, in our simulation with ‘cold rods’ phantom, the count in (4,4) layer-paired crystals is only 1.2% of the total count in the quadruple-layered detector. We believe that this phenomenon also prevents FBP reconstruction from fully utilizing the benefit of multi-layered detectors for the case when there are more than three detector layers. In MLEM iterations using a Monte Carlo computed system matrix, the use of multi-layered detectors provides smaller RRMS as the number of layers increases, as seen in Fig. 12. But improvements seem to be insufficient to overcome the increased memory and computation requirements for the use of the system matrix associated with multi-layered detectors. Note that the size of the system matrix is proportional to the square of the number of detectors. Moreover, in DOI comparisons (Fig. 9), the benefit of detailed DOI information from multi-layered detectors is diminished. Thus, for MLEM
iterations, a single-layered detector seems to be a practical choice for our prototype brain PET. In our simulations, MLEM iterations using a single-layered detector outperformed FBP reconstructions up to at least an octuple-layered detector. We believe that the superior performance of MLEM iterations must be supported by the accuracy of the system matrix. In our simulation, we used the GATE Monte Carlo program [17] to compute the system matrix accurately. By doing so, we incorporated the blurring effect from inter-crystal scattering into the system matrix; hence, MLEM iterations using the system matrix can remove artifacts from the blurring gradually as the iteration proceeds. Roughly speaking, an inherited DOI effect correction of a Monte Carlo computed system matrix for a single-layered detector is even better than the real DOI effect correction due to the octuplelayered detector as long as the former DOI information is used in MLEM iterations and the latter DOI information is used in FBP reconstructions. The fact that MLEM iterations using a Monte Carlo computed system matrix outperforms FBP reconstructions or MLEM iterations using a ray tracing method based system matrix in PET imaging is not new. There have been
830
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 820–831
(100)
Fig. 13 – MLEM reconstructions for the ‘cold rods’ phantom: (a) image by MLEM1 RRMS = 0.049, (c) image by
(100) MLEM3 ,
RRMS = 0.047, (d) image by
several reports that show the benefit of correctly computed system matrices, whether they are computed by measuring the response of real PET systems [11] or by Monte Carlo simulations [4], in parallax error reduction. So far, however, comparisons studies have been performed between reconstruction methods under the same number of detector layers. A novel result from this work is the quantitative evaluation of the DOI effect correction of a Monte Carlo computed system matrix by comparison studies of FBP and MLEM with a different number of detector layers. The simulation results in this paper shows that in our prototype brain PET, the single-layered detector is sufficient since multi-layered detectors provide rather small improvement. For instance, in MLEM, the improvement by the multi-layered detectors is not sufficient to ignore increased computational requirements, which are proportional to the square of the number of layers [5]. It also shows that MLEM iterations with a single-layered detector outperform FBP reconstructions with up to octuple-layered detectors. Based on these results, we conclude that the detailed DOI information from multi-layered detectors is favorable but not necessary in brain PET imaging, since the inherited DOI effect correction of a Monte Carlo computed system matrix for a single-layered detector is sufficient.
(100) MLEM4 ,
(100)
, RRMS = 0.054, (b) image by MLEM2
,
RRMS = 0.045.
Acknowledgements This study was supported in part by the Converging Research Center Program (2011K000715) and by the Nuclear R&D Program (2010-0026096) of the Ministry of Education, Science and Technology, Republic of Korea.
references
[1] T. Ling, T.K. Lewellen, R.S. Miyaoka, Depth of interaction decoding of a continuous crystal detector module, Physics in Medicine and Biology (2007) 2213–2228. [2] Y.H. Chung, Y. Choi, G. Cho, Y.S. Choe, K.H. Lee, B.T. Kim, Optimization of dual layer phoswich detector consisting of LSO and LuYAP for small animal PET, IEEE Transactions on Nuclear Science (2005) 217–221. [3] J.H. Jung, Y. Choi, Y.H. Chung, O. Devroede, P. Bruyndonckx, S. Tavernier, Optimization of LSO/LuYAP phoswich detector for small animal PET, Nuclear Instrumentation and Methods in Physics Research A (2007) 669–675. [4] C.-M. Kao, Y. Dong, Q. Xie, C.-T. Chen, Image reconstruction of a dual-head small-animal PET system by using
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 820–831
[5]
[6]
[7]
[8]
[9]
[10]
Monte-Carlo computed system response matrix, in: 2007 International Meeting on Fully Three-Dimensional Image Reconstruction in Radiology and Nuclear Medicine, 2007, pp. 398–401. T. Yamaya, N. Hagiwara, T. Obi, M. Yamaguchi, K. Kita, N. Ohyyama, K. Kitamura, T. Hasegawa, H. Haneishi, H. Murayama, DOI-PET image reconstruction with accurate system modeling that reduces redundancy of the imaging system, IEEE Transactions on Nuclear Science (2003) 1404–1409. G.T. Gullberg, R.H. Huseman, J.A. Malko, N.J. Pelc, T.F. Budinger, An attenuated projector-backprojector for iterative SPECT reconstruction, Physics in Medicine and Biology (1985) 799–816. G.T. Herman, L.B. Meyer, Algebraic reconstruction techniques can be made computationally efficient, IEEE Transactions on Medical Imaging (1993) 600–609. A. Terstegge, S. Weber, H. Herzog, H. Halling, High resolution and better quantification by tube of response modeling in 3D PET reconstruction, in: IEEE Nuclear Science Symposium, 1997, pp. 1603–1607. V. Selivanov, Y. Picard, J. Cadorette, S. Rodrigue, R. Lecomte, Detector response models for statistical iterative image reconstruction in high resolution PET, IEEE Transactions on Nuclear Science (2000) 1168–1175. T. Frese, N.C. Rouze, C.A. Bouman, K. Sauer, G.D. Hutchins, Quantitative comparison of FBP, EM and bayesian reconstruction algorithms for the IndyPET scanner, IEEE Transactions on Medical Imaging (2003) 258–276.
831
[11] V.Y. Panin, F. Kehren, C. Michel, M. Casey, Fully 3-D PET reconstruction with system matrix derived from point source measurements, IEEE Transactions on Medical Imaging (2006) 907–921. [12] M. Rafecas, B. Mosler, M. Dietz, M. Pögl, A. Stamatakis, M.P. McElroy, S.I. Ziegler, Use of a Monte-Carlo based probability matrix for 3D iterative reconstruction of MADPET-II data, IEEE Transactions on Nuclear Science (2004) 2597–2605. [13] K.J. Hong, Y. Choi, J.H. Jung, J.H. Kang, W. Hu, B.J. Min, Y.S. Huh, S.S. Kim, J.W. Jung, K.B. Kim, M.S. Song, H.W. Park, A MR insertable brain PET using tileable GAPD arrays, in: 2010 IEEE Nuclear Science Symposium and Medical Imaging Conference, M03-5, 2010. [14] C.A. Johnson, J. Seidel, R.E. Carson, W.R. Gandler, A.A. Sofer, M.V. Green, M.E. Daube-Witherspoon, Evaluation of 3D reconstruction algorithms for a small animal PET camera, IEEE Transactions on Nuclear Science (1997) 1303–1308. [15] L.A. Shepp, Y. Vardi, Maximum likelihood reconstruction for emission tomography, IEEE Transactions on Medical Imaging (1982) 113–122. [16] N.-Y. Lee, Y. Choi, A modified OSEM algorithm for PET reconstruction using wavelet processing, Computer Methods and Programs in Biomedicine (2005) 236–245. [17] S. Jan, et al., GATE: a simulation toolkit for PET and SPECT, Physics in Medicine and Biology (2004) 4543–4561. [18] R.L. Siddon, Fast calculation of the exact radiological path for a three dimensional CT array, Medical Physics (1985) 252–255. [19] A.C. Kak, M. Slaney, Principles of Computerized Tomographic Imaging, IEEE Press, 1988.
journal homepage: www.intl.elsevierhealth.com/journals/cmpb
Simulation studies on depth of interaction effect correction using a Monte Carlo computed system matrix for brain positron emission tomography Nam-Yong Lee a , Yong Choi b,∗ a b
School of Computer Aided Science, Institute of Basic Sciences, Inje University, Gimhae, Gyeongnam 621-749, Republic of Korea Department of Electronic Engineering, Sogang University, 1 Sinsu-dong, Mapo-gu, Seoul 121-742, Republic of Korea
a r t i c l e
i n f o
a b s t r a c t
Article history:
The parallax errors caused by a lack of depth-of-interaction (DOI) information often degrade
Received 15 July 2011
reconstruction quality in positron emission tomography (PET). To reduce parallax errors,
Received in revised form
some PET systems employ multi-layered detectors to provide detailed DOI information,
24 April 2012
but this approach requires more complicated detector configurations and signal processing
Accepted 7 May 2012
schemes. In this paper, we conduct simulation studies on an inherited DOI effect cor-
Keywords:
maximum likelihood expectation maximization (MLEM) iterations using the Monte Carlo
rection of a brain PET with a single-layered detector. For this purpose, we compare the Positron emission tomography
computed system matrix of the single-layered detector with filtered backprojection (FBP)
Depth-of-interaction
reconstructions from projection data obtained from multi-layered detectors. We also inves-
Parallax error
tigate the benefits of multi-layered detectors in MLEM iterations using simulated data. The
Spatial resolution
quantitative comparison in this paper shows that an inherited DOI effect correction of a
Image reconstruction
Monte Carlo computed system matrix for a single-layered detector results in the associ-
Maximum likelihood expectation
ated MLEM iterations outperforming FBP reconstructions even for the case in which the
maximization
projection data for FBPs are obtained from an octuple-layered detector. It also shows that
System matrix
use of multi-layered detectors provides better results overall in MLEM reconstruction, but the improvement seems not to be substantial enough to ignore the complexity and costs required for multi-layered detectors. Based on these results, we conclude that detailed DOI information from multi-layered detectors is favorable, but unnecessary in brain PET imaging, because the inherited DOI effect correction via a Monte Carlo computed system matrix for a single-layered detector is sufficient. © 2012 Elsevier Ireland Ltd. All rights reserved.
1.
Introduction
Ideally, a PET scanner needs to have high and uniform spatial resolution across the imaging volume. A PET scanner equipped with elongated detectors or small detector rings is, however, susceptible to parallax errors because of event mispositioning caused by detection of gamma-rays incident at oblique angles on the detector. Until recently, the reduction
∗
of parallax errors has been addressed either by making the detector ring much larger than the phantom (the object to be imaged) so that the incidence angles of gamma-rays emitted from phantoms inside the field of view (FOV) remain sufficiently small or by using a detector module that can produce more detailed DOI information. More detailed DOI information can be obtained by various methods [1–5]. For example, in [1], they used a statistics-based positioning technique to obtain detailed DOI information
Corresponding author. Tel.: +82 10 9933 2624; fax: +82 2 706 4216. E-mail addresses: [email protected], [email protected] (Y. Choi). 0169-2607/$ – see front matter © 2012 Elsevier Ireland Ltd. All rights reserved. http://dx.doi.org/10.1016/j.cmpb.2012.05.007
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 820–831
821
Fig. 1 – Brain PET system consisting of 72 detector modules. Each detector module has a 4 × 4 arrays of Geigermode avalanche photodiode (GAPD) and lutetium oxyorthosilicate (LSO).
from continuous crystal detectors. In Refs. [2,3], a doublelayered detector head consisting of lutetium oxyorthosilicate/lutetium yttrium aluminum perovskite (LSO/LuYAP) phoswich crystals is used, and in [4] lutetium oxyorthosilicate/lutetium yttrium orthosilicate (LSO/LYSO) phoswich crystals are used to obtain binary DOI measurements. In [5] the authors used a quadruple-layered detector and proposed a reduction method to reduce heavy computations caused by the system matrix for the quadruple-layered detector. Obviously, the primary drawback of these approaches is that more detailed DOI information requires more complicated detector configurations and signal processing schemes. Statistical iterative methods may help to reduce the parallax errors in reconstructed images by incorporating a model of the system response into the system matrix. In the absence of noise, the imaging process can be modeled as Y = Pf, where the vector f represents the voxelized image to be reconstructed, the vector Y the measured projection data, and the system matrix P represents the relationship between f and Y. The system matrix P depends on how the projection data Y is arranged. In this work we assume that each element in Y represents the number of annihilation events detected at a crystal pair. Thus the matrix entries Pb,v can be defined as the conditional probabilities of detecting an annihilation event at a crystal pair b that was emitted from the voxel v. They are usually calculated by computing the intersection of the line of response (LOR) defined by a pair of detectors for each voxel [6,7]. More accurate approaches take into account the position of the voxel relative to the detector pair [8] or use analytical models or empirical kernels to compensate for crystal penetration effects [9,10]. It is also possible to compute the system matrix by measuring the response of the real system [11] or by using a Monte Carlo method to simulate the response of the real system [12].
Fig. 2 – Histograms of the reference phantom, which is a uniform phantom covering the entire imaging volume. Here we assume that the PET consists of 72 detector modules, and each detector module has 4 × 1 crystal array. The row and column indices in histogram figures indicate indices of crystal pairs that detect true coincidence events (here single and random events were excluded). Thus, the value at (i, j) pixel of the histogram is the number of true coincidence events detected at the pair with the ith and the jth crystals. (a) Monte Carlo computed histogram of the reference phantom. Here the number of layers is set to one. The total number of true coincidence events is 31.37 millions. (b) Histogram obtained by applying rotational symmetry to the original Monte Carlo computed histogram in (a). Here, 72 rotational averages are applied for k/36, k = 0, 1, . . ., 71, by using the fact that 72 crystal modules were located at equi-angles in a ring-shaped PET system.
The main objective of this work is to evaluate the DOI effect correction of a Monte Carlo computed system matrix of a brain PET with a single-layered detector. We conducted simulation studies to compare MLEM iterations using the system matrix generated by a Monte Carlo computation for a singlelayered detector with FBP reconstruction from projection data obtained by multi-layered detectors. We also investigated the
822
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 820–831
Fig. 3 – Procedure for FBP reconstruction for a double-layered PET detector. The images in the first column are data sets obtained by double-layered detectors. From top to bottom, (1,1), (1,2), (2,1), and (2,2) layer-paired data set are shown, where numbers inside parentheses are layer indices. The scale compensation, rebinning, and interpolation are performed for each layer-paired data set. The FBP reconstruction is applied to the sum of four sinograms obtained from four layer-paired data sets.
benefits of multi-layered detectors in MLEM iterations using simulated data.
2.
Materials and methods
2.1.
PET system
We have been developing a brain PET using 72 detector modules in a ring-shaped scanner whose inside diameter is 330 mm. Each detector module has a 4 × 4 array of Geigermode avalanche photodiode (GAPD) and LSO, where one crystal is 3 mm × 3 mm × 20 mm LSO (the crystal pitch is 3.3 mm). To increase the axial field of view, the final form of the brain PET system will have multi rings of 72 detector modules. For details, see [13] or Fig. 1. In simulation for DOI effect correction, however, we assume that the detector module has 4 × 1 (instead of 4 × 4) crystal array to restrict our attention to two-dimensional reconstruction only. Such restriction is considered to have a fair comparison between FBP and MLEM in DOI effect correction, since FBP is known to be less accurate than MLEM in three-dimensional reconstructions [14]. It is also true that twodimensional reconstruction is sufficient for comparison in DOI effect correction. Having these facts in mind, we assume that the phantom consists of one slice image of 175 × 175 voxels, where the size of the voxel is 1.7 mm × 1.7 mm × 1.7 mm. Thus, the diameter of the imaging volume (the volume of the largest phantom to be reconstructed) is 297.5 mm (175 mm × 1.7 mm) and the height of the imaging volume is 1.7 mm. We simulated detectors to have multi-layers from 1 to 8 for comparison studies. For instance, the double-layered detector had two 3 mm × 3 mm × 10 mm LSO crystals and
generated a binary DOI measurement. We assumed that the single-layered detector had a 20 mm long layer and generated a unary DOI measurement. A similar rule also holds for multi-layered detectors. For instance, each quadruple-layered detector had a layer 5 mm long and generated a quadruple DOI measurement. We assumed that there are no obstacles to photon passage between layers. Thus the multi-layered detector
Fig. 4 – FBP reconstructions for the ‘nine points’ phantom: (a) image by FBP1 , RRMS = 0.696, (b) image by FBP2 , RRMS = 0.572, (c) image by FBP3 , RRMS = 0.496, (d) image by FBP4 , RRMS = 0.452.
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 820–831
823
Fig. 5 – Center row profiles of FBP reconstructions for the ‘nine points’ phantom: (a) profile by FBP1 , (b) profile by FBP2 , (c) profile by FBP3 , (d) profile by FBP4 . Here, rectangular plots are the profiles of the center row of the original ‘nine points’ phantom. We used a normalization method that sets the maximum and minimum intensities of the image to 1 and 0, respectively.
provided more detailed DOI information than the singlelayered detector, while maintaining the same sensitivity as the single-layered detector. Our prototype system is designed for brain imaging and has a single-layered detector. The development of PET systems with multi-layered detectors is not part of our current research interests. We consider multi-layered detectors in our simulation studies to investigate the benefits of multi-layered detectors in FBP and MLEM reconstructions.
2.2.
MLEM reconstruction
The MLEM reconstruction takes the following iterations [15]: fˆvn+1 = fˆvn
P Sn Y b b,v b , Snb = b . P Pb,v fˆvn b b,v
(1)
v
Here fˆ n represents the approximation of the image f at the nth iterates, fˆvn the activity value of fˆ n at the voxel v, Yb the
observed data at the projection bin b, and Pb,v are entries of the system matrix P that determines the relationship between the unknown true activity f and the noiseless projection data P f . One can regard Y as observed data y defined by yb = v b,v v whose statistical mean is y. The projection bin b in (1) is determined by how the observed events are arranged. Rather than working directly with the raw data, it is common to rebin or interpolate the acquired data into a predetermined set of LORs that are convenient for FBP reconstruction. In our work, however, we used the raw dataset to maximize the DOI effect correction of a Monte Carlo computed system matrix. Thus, a projection bin b in (1) represents a crystal pair (c1 , c2 ). We also used a wavelet smoothing technique in MLEM iterations for noise reduction [16]. In this work, we computed the system matrix by using the method described in [12]. First, we used the GATE program [17] to simulate the system response for the reference phantom (a uniform phantom that covers the entire imaging
824
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 820–831
Fig. 7 – MLEM reconstructions for the ‘nine points’ (15)
phantom: (a) image by MLEM1 , RRMS = 0.174, (b) image by (100) MLEM2 ,
(100)
RRMS = 0.101, (c) image by MLEM3 (100)
RRMS = 0.061, (d) image by MLEM4
Fig. 6 – Plots of DOI effect correction by FBP reconstructions: (a) DOI-R and (b) DOI-L.
volume) without including single or random events. Fig. 2a shows the histogram of true coincidence events detected at each crystal pair. Throughout this work we use the term histogram to indicate the data format shown in Fig. 2a only. One pixel in Fig. 2a represents one crystal pair, and the intensity at a given pixel represents the number of true coincidence events detected at the corresponding crystal pair. For each event, we recorded the origin of the event for use in the system matrix generation. Second, we applied rotational symmetry. We used 72 detector modules in a cylindrical ring-shaped PET system. Exploiting this rotational symmetry, we rotated the origin and crystal pair that observed the event by k/36, k = 0, 1, . . ., 71, and we recorded rotated origins and rotated crystal pairs. Fig. 2b shows the rotationally symmetrized histogram of that in Fig. 2a. Finally, we formed the system matrix using the recorded information (origins of events and crystal pairs that detected events).
2.3.
FBP reconstruction
The FBP reconstruction can be explained as the inversion of the Radon transform, which is represented by line integrals.
,
, RRMS = 0.034.
Due to small gaps (approximately 1.5 mm) between modules and the inter-crystal penetration effect, however, line integrals connecting crystal pairs depend not only on the phantom distribution across the line but also on the relative position between the crystal pair and phantom. The response of crystals in the inner part of the module is very different from the response of crystals at the boundary of the module. This phenomenon explains square shaped patterns in Fig. 2 and decreased sensitivities of crystal pairs in which the angle difference in the ring-shaped PET is close to 180◦ . To compensate for such irregularities, we used the following scale compensation: Yb1 = Yb0 ·
RT b , MCb
(2)
where Y0 is the observed histogram of the phantom of interests; RT is the Radon transform of the reference phantom; MC is the Monte Carlo computed histogram of the reference phantom; Y1 is the scale-compensated histogram. We computed RT by measuring the intersection between the ray that connects the crystal pair and the phantom using the standard ray tracing technique [18]. Here the ray that connects the crystal pair is defined by the line that connects the center points of crystal faces. After scale compensation in (2), we used 160 × 281 bins to rebin the data at the histogram, where the numbers of angles and radial points are 160 and 281, respectively, and used linear interpolation for data completion. The rebinned and interpolated data will be referred to as a sinogram in this work. The standard FBP algorithm can be applied to the sinogram to complete the reconstruction. For noise reduction, various filtering methods can be applied [19]. In our simulation, we used the Hamming filter with a cut-off frequency of 0.5.
825
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 820–831
(15)
Fig. 8 – Center row profiles of MLEM reconstructions for the ‘nine points’ phantom: (a) profile by MLEM1 , (b) profile by (100) MLEM2 ,
(100) MLEM3 ,
(100) MLEM4 .
(c) profile by (d) profile by Here, rectangular plots are the profiles of the center row of the original ‘nine points’ phantom. We used a normalization method that sets the maximum and minimum intensities of the image to 1 and 0, respectively.
The FBP reconstruction can be applied to the histogram from a multi-layered detector by using scale compensation, rebinning, and interpolation for each layer-paired data. For instance, a double-layered detector produces four layer-paired data sets, e.g., (1,1), (1,2), (2,1), and (2,2) layer-paired data sets, where the (1,2) layer-paired data set is obtained by crystals at the first layer and crystals at the second layer. Fig. 3 illustrates the process of FBP reconstruction for a double-layered detector: the scale compensation, rebinning, and interpolation are performed in order to produce the layer-paired sinogram, and FBP reconstruction is applied to sum of layer-paired sinograms.
2.4.
Criteria for evaluation
For simulation studies, we used ‘nine points’, ‘uniform’, ‘hot rods’, and ‘cold rods’ phantoms. As mentioned in Section
2.1, simulation studies were performed for two-dimensional image reconstruction and phantoms formed by one slice images of 175 × 175 voxels were used in simulation. We generated histograms of the given phantoms by using the GATE program [17]. We excluded single- or random events in observation, and did not consider attenuation and scattering effects inside phantoms. We only considered true coincidence events as observed data. The numbers of true coincidence events used for the ‘nine points’, ‘uniform’, ‘hot rods’, and ‘cold rods’ phantoms were 0.028, 5.27, 1.33, and 5.47 millions, respectively. These numbers were obtained by letting GATE program to simulate the process of observation for 10,000 s. We generated histograms just one time for each phantom, and used the coordinates of the detection location to compute the multi-layered histograms. Thus, a multi-layered histogram had exactly the same number of photons as a single-layered histogram, but contained detailed DOI information.
826
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 820–831
where Ci is the center location of the ith point in the reconstructed image and pi is the center location of the ith point in the original ‘nine points’ phantom. The quantities DOI-R and DOI-L are used to measure how insufficient DOI information affects reconstruction in resolution and location. RRMS: The RRMS is defined by RRMS =
|fˆv − fv |2 v
|fv |2
,
(5)
where fˆv and fv are intensity values of the reconstructed image and the original phantom image, respectively, at the voxel v.
3.
Results (n)
We used FBPk and MLEMk to denote FBP and MLEM, respectively, where the subscript k stands for the number of detector layers used in generating the histograms and the superscript n inside the parenthesis stands for the iteration number. With this notation, what we wish to compare is the performance of MLEM1 with those of FBPk and MLEMk , k > 1.
3.1.
Fig. 9 – Plots of DOI effect correction by MLEM reconstructions: (a) DOI-R and (b) DOI-L.
We used DOI effect correction and the relative root mean square (RRMS) error in evaluating performance of the reconstruction methods. DOI: For DOI effect correction evaluation, we used the ‘nine points’ phantom. For details, see Fig. 4, which shows reconstructed images of the ‘nine points’ phantom by FBP. Each point in the ‘nine point’ phantom consists of 3 × 3 voxels, thus its support size is 5.1 mm × 5.1 mm because the voxel size is 1.7 mm × 1.7 mm × 1.7 mm. Let Fi , i = 1, 2, . . ., 9, be the full width half maximum (FWHM) of the ith point in the ‘nine points’ phantom. We define DOI-R as DOI-R =
F1 + F9 . 2F5
(3)
Note that F1 and F9 are FWHMs of points that are the farthest points from the center, and F5 is the FWHM of the point at the center. We also define DOI-L as
DOI-L =
9 i=1
|Ci − pi |2 ,
(4)
Comparison by DOI effect correction
Fig. 4a, b, c, and d shows FBP reconstructions of the ‘nine points’ phantom using single-, double-, triple-, and quadruple-layered detectors, respectively. As expected, the image quality improves as the number of layers increases. The improvement is, however, slightly noticeable in the change from the single to double-layered detector but not noticeable in other changes. Fig. 5a, b, c, and d shows profiles of the center rows of Fig. 4a, b, c, and d, respectively. Here rectangular plots are the profiles of the center row of the original ‘nine points’ phantom. In Fig. 5, we used a normalization method, which sets the maximum and minimum intensities of the image to 1 and 0, respectively. In Fig. 5, there is some improvement in resolution when the number of layers increases from 1 to 2, but other changes (such as from 2 to 3 or from 3 to 4) do not noticeably improve the resolution upon visual improvement. To examine the relationship between the number of layers and the DOI effect correction of the FBP reconstruction more clearly, we extended the FBP simulation up to an octuplelayered detector. Fig. 6 shows the DOI effect correction of the FBP reconstruction for the single- to octuple-layered detectors. As seen in Fig. 6, the improvement in DOI effect correction due to the multi-layered detectors is consistent with the DOIL criterion. The DOI-R is improved when the number of layers increases from 1 to 2, but it remains virtually the same for other changes. The FWHMs F5 of the point at the center in Fig. 5 were nearly unchanged for all multi-layered detectors. This result shows that the detailed DOI information from the multilayered detector does not improve the resolution in the center area for the FBP reconstruction. Conversely, F1 and F9 , FWHMs of the furthest points from the center, can be closer to the true ones than F5 due to the detailed DOI information from multilayered detectors. This is the reason why we had DOI-Rs that were smaller than 1 in Fig. 6a.
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 820–831
827
Fig. 10 – Plots of RRMS errors by FBP reconstructions: (a) ‘nine points’, (b) ‘uniform’, (c) ‘hot rods’, (d) ‘cold rods’.
Fig. 7a, b, c, and d shows MLEM reconstructions of the ‘nine points’ phantom using histograms by single-, double, triple-, and quadruple-layered detectors, respectively. Here, we selected images with the smallest RRMS for each simulation. The improvement in image quality attributed to detailed DOI information from multi-layered detectors is not clearly visible in Fig. 7. Fig. 8a, b, c, and d shows profiles of the center rows of Fig. 7a, b, c, and d, respectively. Here, again, as in Fig. 5, rectangular plots are the profiles of the center row of the original ‘nine points’ phantom, and the same normalization method is used. From Fig. 8 we can see that MLEM iterations using the Monte Carlo computed system matrix outperform FBP in DOI effect correction. In fact, MLEM iterations for the singlelayered detector outperform FBP reconstructions even when the projection data for FBP are obtained by the octuple-layered detector. This fact can be verified by comparing the results in Figs. 6 and 9. Fig. 9 shows the DOI effect correction of MLEM iterations for single- to quadruple-layered detectors. As seen in Fig. 9, the benefit of detailed DOI information from multi-layered detectors is not noticeable.
3.2.
Comparison by RRMS error
Fig. 10 shows RRMS errors for FBP reconstructions from multilayered detectors (the number of layers ranges from 1 to 8). As expected, the RRMS decreases as the number of layers increases. The decreasing rate is, however, noticeable only when the number of layers is less than 4. A similar phenomenon also holds for the ‘hot rod’ and ‘uniform’ phantoms. Fig. 11 shows images by FBP reconstructions of the ‘cold rods’ phantom using single-, double-, triple-, and quadruple-detectors. Some visual improvements can be seen only when fewer than four layers are used. Fig. 12 shows RRMS errors for MLEM iterations using multilayered detectors with respect to iteration numbers. Again, as expected, RRMS decreases as the number of layers increases. The improvements in RRMS due to the use of multi-layered detectors are, however, insignificant if we consider that the computational cost for one loop of MLEMn is approximately n2 times of that of MLEM1 [5]. Fig. 13 shows images from MLEM iterations of the ‘cold rods’ phantom using single-, double-, triple-, and quadrupledetectors. Here we selected images with the smallest RRMS
828
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 820–831
Fig. 11 – FBP reconstructions for ‘cold rods’ phantom: (a) image by FBP1 , RRMS = 0.144, (b) image by FBP2 , RRMS = 0.109, (c) image by FBP3 , RRMS = 0.091, (d) image by FBP4 , RRMS = 0.082.
among 100 iterations. In Fig. 13, visual improvements due to the increase of detector layers are also unclear in MLEM iterations.
4.
Discussion and conclusion
In this paper, we investigated the benefits of multi-layered detectors in a prototype brain PET using 72 detector modules, in which one detector module consists of a 4 × 4 GAPD crystal array and the crystal size is 3 mm × 3 mm × 20 mm. The prototype PET is a ring-shaped scanner with a diameter of 330 mm. The simulation results related to the FBP reconstructions (Figs. 4–6, 10, and 11) showed that the use of up to three multilayered detectors improved image quality in reconstruction, but more than three detector layers did not give noticeably better improvements. It was also true that the improvement in visual quality by the use of multi-layered detectors was hardly seen in FBP reconstructions after three detector layers (see Figs. 4 and 10). We believe that such phenomenon is partially caused by the inter-crystal scattering. It is true that the improved DOI effect correction using multi-layered detectors often does not
give the correct path of the event because of the inter-crystal scattering effect. For instance, in our simulation with the cold rods phantom, 42% of true events are affected by the intercrystal scattering effect. The inter-crystal scattering effect produces a kind of blurring on the histograms (and consequently sinograms) in Fig. 3. The FBP reconstruction used in this paper does not have a procedure to reduce such blurring, and hence cannot fully utilize the benefit of multi-layered detectors for more than three detector layers. We also believe that the fact that most events are detected at crystals in frontside layers causes the limited efficiency of multi-layered detectors. To explain this, let us consider Fig. 3, where images in the first column showed histograms obtained by double-layered detectors. Specifically, the histogram in the first row showed counts of events detected at crystals in the first layer, which will be denoted by a (1,1) histogram. Similarly, histograms in the second, third, and fourth rows showed counts of events detected by crystals in (1,2), (2,1), (2,2) layered pairs, respectively, where the numbers inside parentheses are the layer indices. As we can see in Fig. 3, counts of events from the (2,2) layer-paired crystals are much smaller than those from (1,1). To be specific, counts in (1,1), (1,2), (2,1), (2,2) layer-paired crystals are 34.2%, 24.2%, 26.5%, 15.1% of the total
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 820–831
829
Fig. 12 – Plots of RRMS errors by MLEM reconstructions: (a) ‘nine points’, (b) ‘uniform’, (c) ‘hot rods’, (d) ‘cold rods’.
count. In the double-layered detector, detectors in the second layer are still important because they are used to form (1,2) and (2,1) histograms, which contain a substantial amount of data as seen in the first column in Fig. 3. But if we use quadruplelayered detectors, the role of the fourth layer will be almost negligible, since only few events will be detected at crystals in the fourth layer. For instance, in our simulation with ‘cold rods’ phantom, the count in (4,4) layer-paired crystals is only 1.2% of the total count in the quadruple-layered detector. We believe that this phenomenon also prevents FBP reconstruction from fully utilizing the benefit of multi-layered detectors for the case when there are more than three detector layers. In MLEM iterations using a Monte Carlo computed system matrix, the use of multi-layered detectors provides smaller RRMS as the number of layers increases, as seen in Fig. 12. But improvements seem to be insufficient to overcome the increased memory and computation requirements for the use of the system matrix associated with multi-layered detectors. Note that the size of the system matrix is proportional to the square of the number of detectors. Moreover, in DOI comparisons (Fig. 9), the benefit of detailed DOI information from multi-layered detectors is diminished. Thus, for MLEM
iterations, a single-layered detector seems to be a practical choice for our prototype brain PET. In our simulations, MLEM iterations using a single-layered detector outperformed FBP reconstructions up to at least an octuple-layered detector. We believe that the superior performance of MLEM iterations must be supported by the accuracy of the system matrix. In our simulation, we used the GATE Monte Carlo program [17] to compute the system matrix accurately. By doing so, we incorporated the blurring effect from inter-crystal scattering into the system matrix; hence, MLEM iterations using the system matrix can remove artifacts from the blurring gradually as the iteration proceeds. Roughly speaking, an inherited DOI effect correction of a Monte Carlo computed system matrix for a single-layered detector is even better than the real DOI effect correction due to the octuplelayered detector as long as the former DOI information is used in MLEM iterations and the latter DOI information is used in FBP reconstructions. The fact that MLEM iterations using a Monte Carlo computed system matrix outperforms FBP reconstructions or MLEM iterations using a ray tracing method based system matrix in PET imaging is not new. There have been
830
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 820–831
(100)
Fig. 13 – MLEM reconstructions for the ‘cold rods’ phantom: (a) image by MLEM1 RRMS = 0.049, (c) image by
(100) MLEM3 ,
RRMS = 0.047, (d) image by
several reports that show the benefit of correctly computed system matrices, whether they are computed by measuring the response of real PET systems [11] or by Monte Carlo simulations [4], in parallax error reduction. So far, however, comparisons studies have been performed between reconstruction methods under the same number of detector layers. A novel result from this work is the quantitative evaluation of the DOI effect correction of a Monte Carlo computed system matrix by comparison studies of FBP and MLEM with a different number of detector layers. The simulation results in this paper shows that in our prototype brain PET, the single-layered detector is sufficient since multi-layered detectors provide rather small improvement. For instance, in MLEM, the improvement by the multi-layered detectors is not sufficient to ignore increased computational requirements, which are proportional to the square of the number of layers [5]. It also shows that MLEM iterations with a single-layered detector outperform FBP reconstructions with up to octuple-layered detectors. Based on these results, we conclude that the detailed DOI information from multi-layered detectors is favorable but not necessary in brain PET imaging, since the inherited DOI effect correction of a Monte Carlo computed system matrix for a single-layered detector is sufficient.
(100) MLEM4 ,
(100)
, RRMS = 0.054, (b) image by MLEM2
,
RRMS = 0.045.
Acknowledgements This study was supported in part by the Converging Research Center Program (2011K000715) and by the Nuclear R&D Program (2010-0026096) of the Ministry of Education, Science and Technology, Republic of Korea.
references
[1] T. Ling, T.K. Lewellen, R.S. Miyaoka, Depth of interaction decoding of a continuous crystal detector module, Physics in Medicine and Biology (2007) 2213–2228. [2] Y.H. Chung, Y. Choi, G. Cho, Y.S. Choe, K.H. Lee, B.T. Kim, Optimization of dual layer phoswich detector consisting of LSO and LuYAP for small animal PET, IEEE Transactions on Nuclear Science (2005) 217–221. [3] J.H. Jung, Y. Choi, Y.H. Chung, O. Devroede, P. Bruyndonckx, S. Tavernier, Optimization of LSO/LuYAP phoswich detector for small animal PET, Nuclear Instrumentation and Methods in Physics Research A (2007) 669–675. [4] C.-M. Kao, Y. Dong, Q. Xie, C.-T. Chen, Image reconstruction of a dual-head small-animal PET system by using
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 1 0 8 ( 2 0 1 2 ) 820–831
[5]
[6]
[7]
[8]
[9]
[10]
Monte-Carlo computed system response matrix, in: 2007 International Meeting on Fully Three-Dimensional Image Reconstruction in Radiology and Nuclear Medicine, 2007, pp. 398–401. T. Yamaya, N. Hagiwara, T. Obi, M. Yamaguchi, K. Kita, N. Ohyyama, K. Kitamura, T. Hasegawa, H. Haneishi, H. Murayama, DOI-PET image reconstruction with accurate system modeling that reduces redundancy of the imaging system, IEEE Transactions on Nuclear Science (2003) 1404–1409. G.T. Gullberg, R.H. Huseman, J.A. Malko, N.J. Pelc, T.F. Budinger, An attenuated projector-backprojector for iterative SPECT reconstruction, Physics in Medicine and Biology (1985) 799–816. G.T. Herman, L.B. Meyer, Algebraic reconstruction techniques can be made computationally efficient, IEEE Transactions on Medical Imaging (1993) 600–609. A. Terstegge, S. Weber, H. Herzog, H. Halling, High resolution and better quantification by tube of response modeling in 3D PET reconstruction, in: IEEE Nuclear Science Symposium, 1997, pp. 1603–1607. V. Selivanov, Y. Picard, J. Cadorette, S. Rodrigue, R. Lecomte, Detector response models for statistical iterative image reconstruction in high resolution PET, IEEE Transactions on Nuclear Science (2000) 1168–1175. T. Frese, N.C. Rouze, C.A. Bouman, K. Sauer, G.D. Hutchins, Quantitative comparison of FBP, EM and bayesian reconstruction algorithms for the IndyPET scanner, IEEE Transactions on Medical Imaging (2003) 258–276.
831
[11] V.Y. Panin, F. Kehren, C. Michel, M. Casey, Fully 3-D PET reconstruction with system matrix derived from point source measurements, IEEE Transactions on Medical Imaging (2006) 907–921. [12] M. Rafecas, B. Mosler, M. Dietz, M. Pögl, A. Stamatakis, M.P. McElroy, S.I. Ziegler, Use of a Monte-Carlo based probability matrix for 3D iterative reconstruction of MADPET-II data, IEEE Transactions on Nuclear Science (2004) 2597–2605. [13] K.J. Hong, Y. Choi, J.H. Jung, J.H. Kang, W. Hu, B.J. Min, Y.S. Huh, S.S. Kim, J.W. Jung, K.B. Kim, M.S. Song, H.W. Park, A MR insertable brain PET using tileable GAPD arrays, in: 2010 IEEE Nuclear Science Symposium and Medical Imaging Conference, M03-5, 2010. [14] C.A. Johnson, J. Seidel, R.E. Carson, W.R. Gandler, A.A. Sofer, M.V. Green, M.E. Daube-Witherspoon, Evaluation of 3D reconstruction algorithms for a small animal PET camera, IEEE Transactions on Nuclear Science (1997) 1303–1308. [15] L.A. Shepp, Y. Vardi, Maximum likelihood reconstruction for emission tomography, IEEE Transactions on Medical Imaging (1982) 113–122. [16] N.-Y. Lee, Y. Choi, A modified OSEM algorithm for PET reconstruction using wavelet processing, Computer Methods and Programs in Biomedicine (2005) 236–245. [17] S. Jan, et al., GATE: a simulation toolkit for PET and SPECT, Physics in Medicine and Biology (2004) 4543–4561. [18] R.L. Siddon, Fast calculation of the exact radiological path for a three dimensional CT array, Medical Physics (1985) 252–255. [19] A.C. Kak, M. Slaney, Principles of Computerized Tomographic Imaging, IEEE Press, 1988.