Adaptive Autoregressive Model for Reduction of Noise in SPECT

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 (AR-OSEM-AR method). The performance of this method was compared with filtered back projection (FBP) preceded by Butterworth filtering (BW-FBP method) and the OSEM reconstruction method followed by Butterworth filtering (OSEM-BW 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 BW-FBP method showed the highest CR% values and the AR-OSEM-AR method gave the lowest CR% values for cold stacks. In the analysis of hot stacks, the BW-FBP method had higher CR% values than the OSEM-BW method. The BW-FBP method exhibited the lowest FWHM values for cold stacks and the AR-OSEM-AR method for hot stacks. In conclusion, the AR-OSEM-AR method is a feasible way to remove noise from SPECT images. It has good spatial resolution for hot objects.


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 three-dimensionally 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 (3D-MAC) [6], and it was illustrated with patient data.

AR Model.
In two-dimensional 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 low-pass filter that divides the image into two additive components, a predictable image and a prediction error image. An AR process ( 1 , 2 ) is defined by where ( 1 , 2 ) are the predictor (weighting) coefficients, indices 1 and 2 define the type of prediction region in a twodimensional array ( 1 , 2 matrix), and ( 1 , 2 ) represents prediction error, that is, the difference between the predicted value and the current value in this pixel. The predictable image pred is the image obtained by applying the AR model to the original image orig . The prediction error image err = orig − pred .
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 block-wise denoising method has been introduced also for three-dimensional ultrasound images [7]. In the AR model, a prediction region of four orthogonal neighbours of the predicted pixel with a block size of 5 × 5 pixels was used [4]. Seventy-five 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 Poisson-noise-corrupted 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 3D-MAC 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 built-in MATLAB function was used to add Poisson noise to ideal projection images of the 3D-MAC. 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 depth-independent Figure 1: Flowchart of the autoregressive denoising process. ARF: autoregressive filtering; orig : original noise-corrupted image; pred1 and pred2 : predictable images; err 1 , err 2 , err : prediction error images; err : sum of two prediction error images; final : final image. resolution was assumed in the simulated data. Thus, the only factor degrading image quality in the projection images was the Poisson noise.

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 (AR-OSEM-AR), two-dimensional Butterworth filtering before FBP reconstruction in combination with a ramp filter during Computational and Mathematical Methods in Medicine reconstruction (BW-FBP), and OSEM reconstruction followed by three-dimensional Butterworth filtering (OSEM-BW). 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). Noise-free projection images were also reconstructed using the OSEM (Ideal-OSEM) method.

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 regionof-interest (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 noise-free 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 six-pixel-thick profiles were drawn through the 10-, 20-, 40-and 60-mm-wide cylinders, respectively. The FWMH values were calculated using Hermes quality control software (version 2.0).

Patient Study.
Skeletal SPECT was performed three hours after an intravenous injection of 925 MBq of 99m Technetium-labelled 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.

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 5 × 5 pixels for the original image and a prediction region of 3 × 3 and a block size of 6 × 6 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). 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 BW-FBP method showed higher CR% values than the AR-OSEM-AR and OSEM-BW methods ( Table 3). The values were 0.003 and 0.04, respectively. The BW-FBP method had the highest CR% values for all other cold cylinders except the two smallest cylinders at the intermediate count level. The OSEM-BW method, in turn, displayed better performance than the AR-OSEM-AR method ( = 0.002). When the hot stacks were assessed, there were no statistically significant differences between the AR-OSEM-AR and BW-FBP methods nor between the AR-OSEM-AR and OSEM-BW methods, but the BW-FBP method showed higher CR% values than the OSEM-BW method ( = 0.001).
In the analysis of spatial resolution, without exception, the BW-FBP and OSEM-BW methods exhibited lower FWHM values for the cold stacks than the AR-OSEM-AR method ( = 0.002 for both comparisons), but there was no statistical difference between the BW-FBP and OSEM-BW methods (Table 4). For the hot stacks, the AR-OSEM-AR method showed lower FWHM values than the BW-FBP and OSEM-BW methods. The values were 0.01 and 0.04, respectively. Furthermore, the OSEM-BW method had lower FWHM values than the BW-FBP method ( = 0.008). It is of note that the AR-OSEM-AR 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 BW-FBP method showed lower image quality than the two other methods because of streak artefacts ( Figure 6).

Discussion
This paper presented an improved two-dimensional adaptive AR filter and introduced a three-dimensional 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 5 × 5 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 AR-OSEM-AR 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, images reconstructed from noise-free projection images using ordered subset expectation maximisation reconstruction; AR-OSEM-AR: autoregressive filtering before and after ordered subset expectation maximisation reconstruction; BW-FBP: filtered back projection reconstruction method preceded by Butterworth filtering; OSEM-BW: ordered subset expectation maximisation reconstruction followed by Butterworth filtering. Intermediate count level. Profiles were rescaled so that they had the same amount of counts as the profile of the image reconstructed from noise-free projection images using ordered subset expectation maximisation reconstruction. and by using one of the most commonly used filters in nuclear medicine, the Butterworth filter. The BW-FBP 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 built-in 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 AR-OSEM-AR 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 Poisson-distributed. 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 BW-FBP and OSEM-BW methods can be improved by increasing the cut-off frequency, but the improvement occurs at the expense of increased noise.
In our opinion, a strength of the AR-OSEM-AR method is its simplicity, but the lack of user-controlled 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 trade-off between resolution and smoothing will vary at different points within the image.
Because the AR-OSEM-AR method was only marginally better than the OSEM-BW 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 AR-OSEM-AR method has not yet been tested with positron emission tomography (PET) data, but the method should also be suitable for PET data. The signal-to-noise ratio is considerably higher in PET than in SPECT. Therefore, our model will probably provide a good fit for PET data.

Conclusions
The AR-OSEM-AR 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.