ORIGINAL ARTICLE 
Annals of Nuclear Medicine Vol.10, No.1, 33-4O, 1996 

Iterative correction method for shift-variant blurring caused 
by collimator aperture in SPECT 

Koichi OGAWA and Haruto KATSU 
Department of Electronic Informatics, College of Engineering, Hosei University 

A collimation system in single photon computed tomography (SPECT) induces blurring on reconstructed images. The blurring varies with the collimator aperture which is determined by the shape of the hole (its diameter and length), and with the distance between the collimator surface and the object. The blurring has shift-variant properties. This paper presents a new iterative method for correcting the shift-variant blurring. The method estimates the ratio of "ideal projection value" to "measured projection value" at each sample point. The term "ideal projection value" means the number of photons which enter the hole perpendicular to the collimator surface, and the term "measured projection value" means the number of photons which enter the hole at acute angles to the collimator aperture axis. If the estimation is accurate, ideal projection value can be obtained as the product of the measured projection value and the estimated ratio. The accuracy of the estimation is improved iteratively by comparing the measured projection vale with a weighted summation of several estimated projection value. The simulation results showed that spatial resolution was improved without amplification of artifacts due to statistical noise. 

Key words : processing single photon emission CT, iterative correction method, shift-variant blurring, image 

INTRODUCTION 

IN SlNGLE PHOTON EMISSION COMPUTED TOMOGRAPHY (SPECT), it is important to restore images distorted by attenuation and scattering of photons, and blurring due to the collimator. This paper describes a method for correcting the blurring which varies with the collimator aperture determined by the shape of the hole (its diameter and length), and with the distance between the collimator surface and the object,1 The blurring system is called 'shift-variant.' Figure I shows a parallel-hole collimator which consists of holes of the same shape and orientation. The collimator is aimed to absorb the needless photons which enter the scintillation crystal out of perpendicular to its surface. The collimator is so designed as to have an array of holes separated by thin septa which are made of a dense, high-atomic number material such as lead or tungsten. Although, the collimator hole has a finite diameter and length, it is difficult to prevent photons entering the hole at acute angles to the collimator aperture axis. This degrades spatial resolution of reconstructed images, it is therefore important to estimate the photons which enter the hole perpendicular to the collimator surface. 

In the past, many techniques were investigated to improve the spatial resolution of nuclear medicine images. These are roughly divided into two categories. The first is based on the assumption that the blurring is shift-invariant. On this assumption the distortion model is simplified and filtering techniques based on the theory of the linear system are used. These filtering techniques involve the Metz-Wiener filter2-5, Gaussian filter6 and Butterworth-Wiener filter.7 The second approach, including deconvolution techniques8-11 and an iterative method, 12 considers the shift-variant property. In the deconvolution technique,8-11 the shift-variant property is assumed to be shift-invariant in the Fourier domain and the deconvolution is performed to restore the high frequency components which are reduced by the collimator aperture. The second approach yields higher resolution images than the methods which do not consider the shift-variant property, and our study therefore focuses on the second approach. These deconvolution techniques,8-ll however, have the disadvantage of being greatly affected by statistical noise.

For this reason, we developed an iterative method which was not so affected b statistical noise. 
In this method, the term "measured projection value" is defined as the number of photons which enter the collimator aperture at acute angles, and the term "ideal projection value" is defined as the number of photons which enter the hole at right angles to the collimator surface. Considering the geometry shown in Figure 2, the measured projection value can be assumed as a Gaussian-weighted summation of several ideal projection value of acute angles to the collimator aperture. Using the weighting function, the method estimates the ratio of the ideal projection value to the measured projection value. If the estimated ratio is accurate at each point, the ideal projection value can be obtained as the product of the estimated ratio and the measured projection value. The accuracy of the estimation is improved iteratively by comparing the measured projection value with the weighted summation of several estimated projection value. A corrected image is then reconstructed from a set of the estimated projection value by the filtered backprojection method. The simulation results, in which we considered a two-dimensional collimation system for simplification, showed that spatial resolution was improved without amplification of artifacts due to statistical noise. 

METHODS 

Figure 2 shows the relation between the measured projection value Pm(s,THETA) and the ideal projection value P(s,THETA). The solid line l denotes the collimator surface at the rotation angle THETA, and the broken line l'denotes the collimator surface at the rotation angle THETA+ PHY. D is the distance between the collimator surface and the rotation center of the detector. The collimator aperture at the sample point (s,THETA) is shown as the shaded area, and the maximum angle of the aperture is PHY0. The measured projection value Pm(s.THETA) is the number of photons which enter the hole from the shaded area. The ideal projection value P(s,THETA) is the number of photons which enter the hole from the narrow beam area n which is perpendicular to the line l, and P(s',THETA+PHY) is the number of photons which enter the hole from the narrow beam area n' which is perpendicular to line l'. It is easy to see that the shaded area can be divided into several narrow beam areas. The contribution of photons which enter the hole from these narrow beam areas can be approximated to a Gaussian function of PHY. From the relation Pm(s.THETA) can be determined as the weighted summation of P(s',THETA+PHY), that is; 
(1)
where W(PHY) is a Gaussian function as follows: 
(2) 
In Eq. (2), b is a parameter which is determined by the shape of the hole. 
On the other hand, the ideal projection value P(s,THETA) is assumed to be the product of a factor A(s,THETA) and the measured projection value Pm(s,THETA), as follows: 
(3) 
The physical meaning of A(s,THETA) is the ratio of the photons which enter the hole from the narrow beam area n to the photons which enter the hole from the shaded area, and we named it the "elimination factor." The factor is rewritten using equations (1) and (3) as follows: 
(4)
Our task is to decide factor A(s,THETA). If the factor is accurate at each point, the ideal projection value can be obtained as the product of the factor and the measured projection value. To solve this problem, we estimated the factor using the following iterative process. The flowchart which describes the procedure in discrete form is shown in Figure 3. 
Step 1: 
To decide the initial elimination factor A [0](s, THETA), we replace the ideal projection value P(s,THETA) by the mea-sured projection value Pm(s,THETA) in Eq. (4), as follows: (5)
Step 2:
The initial estimated projection value Pn[1](s,THETA) is obtained as the product of the elimination factor A[0](s,THETA) and the measured projection value Pm(s,THETA). 
Step 3: 
The k-th elimination factor A[k](s,THETA) is re-placed by the new one which is the ratio of the k-th estimated projection value Pn[k](s.THETA) to the measured projection value Pm(s,THETA). 
Step 4:
The projection value Pa[k](s,THETA) is calculated as the weighted summation of several estimated projection value Pn[k](s,THETA). 
Step 5:
If the calculated projection value Pa[k](s,THETA) is equal to the measured projection value Pm(s,THETA) at each sample point, the accurate estimation is established. To evaluate the accuracy of the estimation, we calculate the error projection P[k](s,THETA) which is the difference between pa[k](s, THETA) and Pm(s,THETA). Moreover, we evaluate the total error Q[k] which is a summation of the absolute values for Pe[k](s,THETA) of all sample points. 
Step 6:
When the total error Q[k] becomes a smaller value than the preset one, this iterative correction process Is finished and the corrected image is reconstructed from a set of the last estimated projection value using the filtered backprojection method. Otherwise go to Step 7. 
Step 7:
The k + 1 th estimated projection value Pn[k+1](s, THETA) can be obtained by subtracting A[k](s, THETA) Pe[k](s,THETA) from Pn[k](s,THETA). Then, k is increased by one and go to Step 3. 

SIMULATION 

In order to evaluate the validity of the proposed method, we carried out simulations with a numerical phantom. The phantom is shown in Figure 4. The matrix size is 128 x 128 pixels and the pixel size is 2.34 mm/pixel. In the simulations we evaluated the following two types of collimator shape: a high resolution (H.R.) type with PHY0 = 1deg. and a high sensitivity (H.S.) one with PHY0 = 5deg.. The spatial resolution (FWHM) is denoted by y*tan PHY0, where y is the distance between the collimator surface and the source, and PHY0 is the maximum aperture angle (see Fig. 1 ). The contribution of photons which enter the hole at acute angles to its axis in the collimator aperture can be approximated to the Gaussian function as follows: 
(6) 
 In both collimator types (H.R. and H.S.), the projection value was calculated by means of the blurring function Eq. (6) for 720 angular samples with an interval of 0.5 degrees in parallel beam geometry. The number of views was 360 and the number of bins was 128. The distance D between the collimator surface and the rotation center of the collimator was 15cm. The value used for the convergence test (Step 6) in METHODS was 0.01 . The resulting images were reconstructed by the filtered-backprojection method with the Shepp & Logan filter. For simplicity, we didn't consider attenuation and scattering of photons. 
 Practically the measured value is influenced by statistical noise. The noise varies with the aperture angle and data acquisition time. To investigate noise effects in the correction process, the following Gaussiannoise was added to the measured projection data: 
(7) 
where T is the time factor and N[O,1] denotes the normal distribution in which the average and the standard deviation are zero and one respectively. 
 Moreover, in order to evaluate the reconstructed values quantitatively, we set three regions of interest (ROIs) shown in Figure 5. The size of these ROIs was 10 x 10 pixels. We then calculated the average, the standard deviation, and the percent fractional standard deviation (% FSD) for these ROIs. 

RESULTS AND DISCUSSION 

For noise free cases, resulting images for both collimator types are shown in Figure 6; (a) is an image degraded by the H.R. collimator, (b) is a corrected image (converged) of (a), (c) is an image degraded by the H.S. collimator, and (d) is a corrected image (converged) of (c). Figure 7 shows profiles on the y-axis of the images, i,e., the ideal image described below, the degraded image by the H.S. collimator, and the corrected image (converged). The ideal image, the profile of which was shown as the dotted line, was reconstructed from the ideal projection data which wasn't influenced by the collimator aperture. It can be seen that spatial resolution is much improved. 
 For noise considering cases, resulting images for T= l are shown in Figure 8; (a) is a reconstructed image considering the narrow beam model which isn't influenced by the collimator aperture, (b) is an image degraded by the H.R. collimator, (c) is a corrected image (converged) of (b), (d) is an image degraded by the H.S. collimator, and (e) is a corrected image (converged) of (d). The total count for each image was approximately 5,000. Figure 9 shows resulting images for T = O.1 . The notation of Figure 9 is the same as that of Figure 8. From Figures 8 and 9, it can be seen that the blurring caused by the H.S. collimator is eliminated. 
 The results of quantitative evaluation for these reconstructed images are shown in Table 1. In each ROI, fluctuation of the corrected image for each time factor is almost the same as that of the reconstructed image considering the narrow beam model. 
 Followings are some discussions of these two models. 
 A. Noise Free Model. The proposed method was applied to restore the blurred images due to the collimator aperture. The examined collimator types were the H.R. collimator (PHY0 = 1deg.) and the H.S. collimator (PHY0 = 5deg.). In the case of the H.R. collimator, it is difficult to see differences between the degraded image (a) and the corrected image (b) in Figure 6, because the H.R. collimator hardly degrades the spatial resolution. On the other hand, in the case of the H.S. collimator, as shown in Figure 6-(c), it is easy to see the blurring caused by the collimator aperture. In Figure 6-(d) the blurring is much reduced by the proposed method. To evaluate the accuracy of the correction, profiles of the images are compared in Figure 7. It is observed that the smoothed edges are clearly improved. The variations in the total errors Q for the noise free and the noise models are shown in Figure 10, where (a) is the results for the H.R. collimator, and (b) is those for the H.S. collimator. In each case, the total error is con-verged into a small value. These results show that the elimination of the shift-variant blurring can be achieved by the proposed method. 
B. Noise Model. To evaluate the effects of noise in the iterative correction process, the statistical noise was added to the blurred projection data as shown in Eq. (7). The resulting images for the ti,me factor T= 1 are shown in Figure 8. In case of the H.R. collimator, there is hardly seen any difference between the blurred image (b) and the corrected image (c). In the case of the H.S. collimator, there is less statistical fluctuation in blurred image (d) than in blurred image (b). The reason is that more counts can be acquired by the H.S, collimator than the H.R, collimator, and, therefore, the S/N ratio for the reconstructed image using the H.S. collimator is higher than that using the H.R. collimator. The blurring of the image Figure 8-(d) is eliminated in the corrected image Figure 8-(e). It seems that there is no artifact caused by the iterative method but more conspicuous fluctuation occurred in the corrected image (e) than in the blurred image (d). The enhancement of fluctuation is also seen in the corrected image Figure 9-(e) to which projection data was added more noise (T=0.1) than in the case of Figure 8 (T = 1). Table 1 shows a comparison of the fluctuation in both acquisition times T= 1 and T= 0.1 among the reconstructed images using the narrow beam model, the degraded images, and the corrected images. The average value is approximately 140 in ROI 1 , 100 in ROI 2, and 50 in ROI 3. In each ROI, it can be seen that % FSD of the corrected image for each time factor is almost the same as that of the narrow beam model. The variations in total errors Q are shown in Figure 10 for both H.R. and H.S. collimator types, and the total error was converged into a small value in each case. These results show that spatial resolution degraded by the collimator aperture can be restored without amplification of artifacts due to the statistical noise. 

CONCLUSIONS 

In this paper we proposed the iterative correction method to eliminate the shift-variant blurring caused by the collimator aperture. The results in Figures 6, 8 and 9 show that the correction is performed well. The advantages of the method are: 1) The method can improve spatial resolution degraded by the collimator aperture, 2) The correction process is not greatly influenced by the statistical noise, and 3) The projection value which is not influenced by the collimator aperture is uniquely obtained from the measured projection value, but the accuracy depends on the aperture angle, the source distribution, and the data acquisition time. Further investigation of this method should include improvement of accuracy and restraint of the statistical fluctuation. 

REFERENCES 

l. Payne JT, Williams LE, Ponto RA, Goldberg ME, Loken MK. Comparison and performance of anger cameras. Radiology 109: 381-386, 1973. 
2. King MA. Doherty PW, Schwinger RB. A Wiener filter for nuclear medicine images. Med Phys 10: 876-880, 1983. 
3. King MA, Doherty PW, Schwinger RB, Jacobs DA. Kidder RE. Miller TR. Fast count-dependent digital filtering of nuclear medicine images: Concise communication. J Nucl Med 24: 1039-l045, 1983. 
4. King MA, Schwinger RB. Doherty PW, Penny BC. Two-dimensional filtering of SPECT images using the Metz and Wiener filters. J Nucl Med 25: 1231-1240, 1984. 
5. Hon TC, Rangayyan RM. Hahn LJ. Kloiber R. Restoration of gamma camera-based nuclear medicine images. IEEE Trans Med Imag 8 (4): 354-363, 1989. 
6. Madsen MT, Park CH. Enhancement of SPECT images by Fourier filtering the projection image set. J Nucl Med 26 (4): 395-402, 1985. 
7. Honda N, Machida K, Tsukada J. Kaizu H, Hosoba M. 
Optimal preprocessing Butterworth-Wiener filter for Tl-201 myocardial SPECT. Eur J Nucl Med 13: 404-407, l987. 
8. Ogawa K, Paek S, Nakajima M, Yuta S. Kubo A. Hashimoto S. Correction of collimator aperture using a shift-variant deconvolution filter in gamma camera emission CT. Proceedings of the SPIE. Med Imag 11 914: 699-706, 1988. 
9. Edholm PR, Lewitt RM, Lindholm B. Novel properties of the Fourier decomposition of the sinogram, Proceedings of the SPIE, International Workshop on Physics and Engineering of Computerized Multidimensional Imaging and Processing, vol. 671 , pp. 8-18, 1986. 
10. Lewitt RM, Edholm PR, Xia W. Fourier method for correction of depth-dependent collimator blurring, Proceedings of the SPIE. Med Imag 111 1092: 232-243, 1989. 
l1. Xia W, Lewitt RM, Edholm PR. Resolution improvement in SPECT by spatial filtering. Conference Record of the IEEE Med and Biol Soci, vol. 12, pp. 379-380, 1990. 
12. Xia W, Lewitt RM. Iterative correction for spatial collimator blurring in SPECT. Conference Record of the IEEE Nucl Sci Symposium 2, pp. 1158-1162, 1990.