Research Article  Open Access
Adaptive Autoregressive Model for Reduction of Noise in SPECT
Abstract
This paper presents improved autoregressive modelling (AR) to reduce noise in SPECT images. An AR filter was applied to prefilter projection images and postfilter ordered subset expectation maximisation (OSEM) reconstruction images (AROSEMAR method). The performance of this method was compared with filtered back projection (FBP) preceded by Butterworth filtering (BWFBP method) and the OSEM reconstruction method followed by Butterworth filtering (OSEMBW method). A mathematical cylinder phantom was used for the study. It consisted of hot and cold objects. The tests were performed using three simulated SPECT datasets. Image quality was assessed by means of the percentage contrast resolution (CR%) and the full width at half maximum (FWHM) of the line spread functions of the cylinders. The BWFBP method showed the highest CR% values and the AROSEMAR method gave the lowest CR% values for cold stacks. In the analysis of hot stacks, the BWFBP method had higher CR% values than the OSEMBW method. The BWFBP method exhibited the lowest FWHM values for cold stacks and the AROSEMAR method for hot stacks. In conclusion, the AROSEMAR method is a feasible way to remove noise from SPECT images. It has good spatial resolution for hot objects.
1. Introduction
Numerous methods for removing noise from SPECT images have been proposed [1, 2]. This indicates the difficulty of the task. Noise removal can be performed before reconstruction (prefiltering), during reconstructions or after reconstruction (postfiltering). In modern iterative methods, collimator correction denoises images during reconstruction, but the reconstructed images may still require postfiltering [3]. Earlier we introduced an adaptive autoregressive (AR) filter to reduce noise in scintigraphic planar images or projection images of a SPECT study [4]. In the present work, the AR filter was further improved to reduce noise from the projection images and also from threedimensionally reconstructed data. It is important to apply the best AR filter to the projection data of SPECT, because a small change in the projection data may cause a large change in the estimated transaxial image [5]. Our method was compared with two established methods for improving image quality in SPECT. The methodical comparison was carried out using a threedimensional mathematical cylinder phantom (3DMAC) [6], and it was illustrated with patient data.
2. Methods
2.1. AR Model
In twodimensional AR modelling, each value of an image is regressed on its neighbourhood pixel values, called the prediction region. An AR model can be regarded as a lowpass filter that divides the image into two additive components, a predictable image and a prediction error image. An AR process is defined by where are the predictor (weighting) coefficients, indices and define the type of prediction region in a twodimensional array (, matrix), and represents prediction error, that is, the difference between the predicted value and the current value in this pixel. The predictable image is the image obtained by applying the AR model to the original image . The prediction error image .
In a typical scintigraphic image, there are large local spatial variations in the count number of the image. Therefore, the same model cannot be applied to the entire image, but the model must be adapted to the variations. In this adaptive method, the image area is divided into smaller blocks and the AR model is then fitted into each block separately by using MATLAB subroutines. Recently, a blockwise denoising method has been introduced also for threedimensional ultrasound images [7]. In the AR model, a prediction region of four orthogonal neighbours of the predicted pixel with a block size of pixels was used [4]. Seventyfive percent overlap of the image blocks in combination with one iteration of the filtering procedure was used. The two error term images were summed up and subjected to AR filtering, and the resulting image was then added to the iteratively filtered image (Figure 1). In the present study, we tested the effect of using another AR model for the summed error term images than for the original image. We used the same transaxial slice of the Zubal phantom [8] and the same simulation conditions as in our previous work [4], and image quality was assessed by means of the mean squared error (MSE) of the image. It is of note that the Poissonnoisecorrupted slice of the phantom actually represented an artificial scintigraphic planar image or a projection image of a SPECT study. The AR model with the lowest MSE was then used to prefilter the SPECT projection images and also to postfilter iteratively reconstructed data. The filter was applied to each set of orthogonal plane images separately. The software was based on MATLAB subroutines (The MathWorks, Inc.).
2.2. Phantom
Data were simulated using a 3DMAC phantom [6]. The phantom measured 200 mm in both diameter and length. It comprised three imbedded objects: two hot objects and a cold one. Each object consisted of five stacked cylinders. The cylinders had diameters of 4, 10, 20, 40, and 60 mm and a length of 30 mm. The smallest cylinder was not utilised in the present study because its dimensions were beyond the resolution of the simulated SPECT system used. Relative activities were 1, 0, 2, and 4 for the background, a cold stack, and two hot stacks, respectively. The tests were performed using three SPECT datasets with different image statistics. Total counts of the projection images were approximately 50000 (low level), 100000 (intermediate level), and 150000 (high level) per projection. A builtin MATLAB function was used to add Poisson noise to ideal projection images of the 3DMAC. The mean counts of a pixel in the projection images were 12, 24, and 37, respectively, and the range of pixel values was 0–52, 0–104, and 0–156, respectively. The matrix size was 64 × 64 pixels, pixel size was 4 mm, and the number of projections was 120. There was no scatter or attenuation component and perfect depthindependent resolution was assumed in the simulated data. Thus, the only factor degrading image quality in the projection images was the Poisson noise.
2.3. Reconstruction Methods
Transaxial slices were reconstructed using either the filtered back projection method (FBP) [9] or an iterative ordered subset expectation maximisation (OSEM) algorithm [10]. The reconstruction methods were implemented on the Hermes SPECT (G) reconstruction software (version 3.8) and reconstruction engine of Hermes HybridRecon (Hermes Medical Solutions, Stockholm, Sweden), respectively. Three methods were compared: AR filtering before and after ordered OSEM reconstruction (AROSEMAR), twodimensional Butterworth filtering before FBP reconstruction in combination with a ramp filter during reconstruction (BWFBP), and OSEM reconstruction followed by threedimensional Butterworth filtering (OSEMBW). The Butterworth filter was originally designed for onedimensional data [11]. In the OSEM method, the number of subsets was set to 8 and the number of iterations to 10. Postfiltering was performed using Multimodality software (Hermes Medical Solutions, Stockholm, Sweden). Noisefree projection images were also reconstructed using the OSEM (IdealOSEM) method.
2.4. Assessment of Image Quality
To obtain a fair comparison of the methods, the same amount of filtering was applied in each method. This was done by drawing a circular regionofinterest (ROI) 150 mm in diameter in the uniform part of the phantom and calculating the percentage coefficient of variation (CoV%) in the ROI, that is, the ratio of the standard deviation to the mean multiplied by 100. This kind of presentation ensures that filtering between each method is equal.
Percentage contrast resolution (CR%) values for the activity in each cylinder and uniform activity were calculated. CR% can be expressed by the following formula [12]: where is the count value of uniform activity and is the count value in each cylinder. Activity in each cylinder was analysed using a circular ROI with the same diameter as the cylinder. The ROIs were drawn on noisefree transaxial slices and were copied to each set of reconstructed data, so their position and area were equal in every image. The ROIs were drawn using Multimodality software. The CR% values were obtained using the average counts in the ROIs.
Spatial resolution was estimated by the full width at half maximum (FWHM) of the line spread functions of the cylinders. One, two, four and sixpixelthick profiles were drawn through the 10, 20, 40 and 60mmwide cylinders, respectively. The FWMH values were calculated using Hermes quality control software (version 2.0).
2.5. Patient Study
Skeletal SPECT was performed three hours after an intravenous injection of 925 MBq of methylene diphosphonate. The images were obtained over a 360° arc, using 64 projections at 20 sec per projection. The images were acquired into a 128 × 128 matrix with a pixel size of 4.8 mm. Total counts of the projection images were 41006–66830 counts per projection.
2.6. Statistical Methods
The data were analysed using WinSTAT for Excel (version 2007.1; R. Fitch Software, Staufen, Germany). Pairwise comparisons were performed with the nonparametric Wilcoxon’s ranksum test. Comparisons were made between the AROSEMAR and BWFBP methods, the AROSEMAR and OSEMBW methods, and the BWFBP and OSEMBW methods. Data from the cold stacks and the pooled hot stacks were analysed separately. For each cylinder, the paired difference between the values of a variable was computed. The values of the differences were sorted to get a rank order. Finally, the mean rank of negative differences was compared with that of positive differences. Wilcoxon’s ranksum test determines to what extent the difference in mean rank is significant. A value of less than 0.05 was considered significant.
3. Results
The MSE of the images improved when a different AR model was used for the summed error term image rather than for the original image. A prediction region of four orthogonal neighbours with a block size of pixels for the original image and a prediction region of and a block size of for the summed error term image produced the lowest MSE, although the differences were small (Table 1). Part of the counts at the edges of the image could be returned to the filtered image to reduce blurring of the image (Figure 2).

(a)
(b)
(c)
(d)
Butterworth filtering was chosen so that the methods had the same amount of statistical fluctuation in the uniform part of the phantom, as confirmed by the CoV% values (Table 2). For the cold stacks, the BWFBP method showed higher CR% values than the AROSEMAR and OSEMBW methods (Table 3). The values were 0.003 and 0.04, respectively. The BWFBP method had the highest CR% values for all other cold cylinders except the two smallest cylinders at the intermediate count level. The OSEMBW method, in turn, displayed better performance than the AROSEMAR method (). When the hot stacks were assessed, there were no statistically significant differences between the AROSEMAR and BWFBP methods nor between the AROSEMAR and OSEMBW methods, but the BWFBP method showed higher CR% values than the OSEMBW method ().
 
AROSEMAR: autoregressive filtering before and after ordered subset expectation maximisation algorithm; BWFBP: Butterworth prefiltering and filtered back projection; OSEMBW: ordered subset expectation maximisation algorithm and Butterworth postfiltering; RelAct: activity relative to background activity of 1; CF: cutoff frequency. The order of the filter was 2; —: not definable. 
 
A: low count level; B: intermediate count level; C: high count level; AROSEMAR: autoregressive filtering before and after ordered subset expectation maximisation algorithm; BWFBP: Butterworth prefiltering and filtered back projection; OSEMBW: ordered subset expectation maximisation algorithm and Butterworth postfiltering; RelAct: activity relative to background activity of 1; Ø: diameter. 
In the analysis of spatial resolution, without exception, the BWFBP and OSEMBW methods exhibited lower FWHM values for the cold stacks than the AROSEMAR method ( for both comparisons), but there was no statistical difference between the BWFBP and OSEMBW methods (Table 4). For the hot stacks, the AROSEMAR method showed lower FWHM values than the BWFBP and OSEMBW methods. The values were 0.01 and 0.04, respectively.
 
A: low count level; B: intermediate count level; C: high count level; AROSEMAR: autoregressive filtering before and after ordered subset expectation maximisation algorithm; BWFBP: Butterworth prefiltering and filtered back projection; OSEMBW: ordered subset expectation maximisation algorithm and Butterworth postfiltering; RelAct: activity relative to background activity of 1; Ø: diameter. 
Furthermore, the OSEMBW method had lower FWHM values than the BWFBP method (). It is of note that the AROSEMAR method showed better resolution than the other two methods in the analysis of the two smallest hot stacks, with two exceptions in the analysis of cylinders with two times the background activity and a diameter of 10 mm (Table 4).
Visually, the differences between the images produced by the three methods were small (Figures 3, 4, and 5). When comparing the skeletal SPECT data, the BWFBP method showed lower image quality than the two other methods because of streak artefacts (Figure 6).
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(a)
(b)
(c)
4. Discussion
This paper presented an improved twodimensional adaptive AR filter and introduced a threedimensional adaptive AR model for reduction of noise in SPECT images. We demonstrated that the quality of scintigraphic images can be improved when the same AR procedure is not applied to the original image and the summed error term image. We have previously shown that if a prediction region of four orthogonal neighbours of the predicted pixel with a block size of pixels is used for both the original image and the summed error term image in the same simulation conditions, then the mean squared errors for the three different images with Poisson statistics are 0.85, 2.23, and 7.12 [4]; that is, this combination exhibits lower performance than any of those presented in Table 1.
The goal of filtering in SPECT is to suppress statistical noise and simultaneously preserve contrast and spatial resolution [1]. In the present study we showed that the AROSEMAR method simultaneously provides both efficient noise rejection and good spatial resolution for hot objects. The methodological comparison was done using the wellknown de facto reconstruction standards, FBP and OSEM, and by using one of the most commonly used filters in nuclear medicine, the Butterworth filter.
The BWFBP method produced better performance than the two other methods in the analysis of the cold stacks. FBP’s good performance with cold features has been noticed before [5, 13]. OSEM’s builtin nonnegativity constraint explains its poor contrast in cold regions. In the analysis of the hot stacks, the magnitude of the differences between the three methods proved to be small, but the AROSEMAR method had statistically the best performance. This is obviously due to the fact that part of the counts at the edges of the error term images could be returned back to the filtered image. Adding postfiltering to the method produced efficient noise reduction without compromising contrast or spatial resolution significantly.
The FBP method consists of filtering of the projection data and back projection of the filtered data [5, 10]. Prefiltering is generally not applied in the OSEM method because it lowers spatial resolution. Secondly, OSEM assumes that projection pixel values are independent and the number of counts is Poissondistributed. Filtering might hamper these assumptions. OSEM reconstructed images are usually postfiltered because the images become noisier as the iterations proceed.
The disadvantage of FBP is that it can produce radial streak artefacts because filtered noisy projection profiles do not cancel each other out in back projection. In the present study, the phenomenon was seen in the clinical data. Iterative reconstruction algorithms also provide some other advantages over FBP. They permit the use of several important corrections, such as scatter, attenuation, and collimator response corrections, which can be included in the image reconstruction procedure. Incorporation of anatomical information derived from magnetic resonance imaging or computerized tomography is possible as well [14]. For the abovementioned reasons, FBP has in recent years been progressively replaced with iterative reconstruction algorithms.
The Butterworth filter is defined by two parameters: cutoff frequency and order [2]. In the present study, order was set to 2 because a ringing artefact is imperceptible to Butterworth filters of order 2 but can become a significant factor in filters of a higher order [15]. Filter orders much higher than 2 are often seen in clinical praxis. Edge sharpness in images produced by the BWFBP and OSEMBW methods can be improved by increasing the cutoff frequency, but the improvement occurs at the expense of increased noise.
In our opinion, a strength of the AROSEMAR method is its simplicity, but the lack of usercontrolled variables can also be regarded as a limitation. Sometimes, adjustable parameters are needed. No particular filter can emerge as the best filter for any organ system. However, filtering should be performed locally in the spatial domain, not globally in the frequency domain, because the correct tradeoff between resolution and smoothing will vary at different points within the image.
Because the AROSEMAR method was only marginally better than the OSEMBW method, an additional study is needed to find out whether image quality will be even better if the AR method is applied to the intermediate results in between the iterations. Secondly, the AROSEMAR method has not yet been tested with positron emission tomography (PET) data, but the method should also be suitable for PET data. The signaltonoise ratio is considerably higher in PET than in SPECT. Therefore, our model will probably provide a good fit for PET data.
5. Conclusions
The AROSEMAR method is a feasible denoising method in SPECT. It has good spatial resolution for hot features and it is simple to use. It does not have any adjustable parameters.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
References
 M. Lyra and A. Ploussi, “Filtering in SPECT image reconstruction,” International Journal of Biomedical Imaging, vol. 2011, Article ID 693795, 14 pages, 2011. View at: Publisher Site  Google Scholar
 P. Zanzonico, “Technical requirements for SPECT: instrumentation, data acquisition and processing, and quality control,” in Clinical SPECT Imaging, E. L. Kramer and J. J. Sanger, Eds., pp. 7–41, Raven Press, New York, NY, USA, 1995. View at: Google Scholar
 T. S. Kangasmaa, J. T. Kuikka, E. J. Vanninen, H. M. Mussalo, T. P. Laitinen, and A. O. Sohlberg, “Halftime myocardial perfusion SPECT imaging with attenuation and Monte Carlobased scatter correction,” Nuclear Medicine Communications, vol. 32, no. 11, pp. 1040–1045, 2011. View at: Publisher Site  Google Scholar
 R. Takalo, H. Hytti, and H. Ihalainen, “Adaptive autoregressive model for reduction of poisson noise in scintigraphic images,” Journal of Nuclear Medicine Technology, vol. 39, no. 1, pp. 19–26, 2011. View at: Publisher Site  Google Scholar
 J. Biemond, R. L. Lagendijk, and R. M. Mersereau, “Iterative methods for image deblurring,” Proceedings of the IEEE, vol. 78, no. 5, pp. 856–883, 1990. View at: Publisher Site  Google Scholar
 H. Onishi, N. Motomura, M. Takahashi, M. Yanagisawa, and K. Ogawa, “A 3dimensional mathematic cylinder phantom for the evaluation of the fundamental performance of SPECT,” Journal of Nuclear Medicine Technology, vol. 38, no. 1, pp. 42–48, 2010. View at: Publisher Site  Google Scholar
 L. Li, W. Hou, X. Zhang, and M. Ding, “GPUbased blockwise nonlocal means denoising for 3D ultrasound images,” Computational and Mathematical Methods in Medicine, vol. 2013, Article ID 921303, 10 pages, 2013. View at: Publisher Site  Google Scholar
 I. G. Zubal, C. R. Harrell, E. O. Smith, Z. Rattner, G. Gindi, and P. B. Hoffer, “Computerized threedimensional segmented human anatomy,” Medical Physics, vol. 21, no. 2, pp. 299–302, 1994. View at: Publisher Site  Google Scholar
 R. A. Brooks and D. G. Chiro, “Principles of computer assisted tomography (CAT) in radiographic and radioisotopic imaging,” Physics in Medicine and Biology, vol. 21, no. 5, pp. 689–732, 1976. View at: Google Scholar
 H. M. Hudson and R. S. Larkin, “Accelerated image reconstruction using ordered subsets of projection data,” IEEE Transactions on Medical Imaging, vol. 13, no. 4, pp. 601–609, 1994. View at: Publisher Site  Google Scholar
 S. Butterworth, “On the theory of filter amplifiers,” Wireless Engineering, vol. 7, pp. 536–541, 1930. View at: Google Scholar
 W. A. Edelstein, P. A. Bottomley, H. R. Hart, and L. S. Smith, “Signal, noise, and contrast in nuclear magnetic resonance (NMR) imaging,” Journal of Computer Assisted Tomography, vol. 7, no. 3, pp. 391–401, 1983. View at: Publisher Site  Google Scholar
 S. D. Wollenweber and K. L. Gould, “Investigation of cold contrast recovery as a function of acquisition and reconstruction parameters for 2D cardiac PET,” in Proceedings of the IEEE Nuclear Science Symposium Conference Record (NSS/MIC '05), vol. 5, pp. 2552–2556, October 2005. View at: Publisher Site  Google Scholar
 G. Gindi, M. Lee, A. Rangarajan, and I. G. Zubal, “Bayesian reconstruction of functional images using anatomical information as priors,” IEEE Transactions on Medical Imaging, vol. 12, no. 4, pp. 670–680, 1993. View at: Publisher Site  Google Scholar
 R. C. Gonzales and R. E. Woods, Digital Image Processing, Prentice Hall, Upper Saddle River, NJ, USA, 2002.
Copyright
Copyright © 2015 Reijo Takalo et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.