Abstract
This paper proposes a 2-step image reconstruction method in which the nonnegativity constraint in the iterative maximum-likelihood expectation maximization (MLEM) algorithm is used to effectively reduce Gibbs ringing artifacts. Methods: Gibbs artifacts are difficult to control during imaging reconstruction. The proposed method uses the postprocessing strategy to suppress Gibbs artifacts. In the first step, a raw image is reconstructed from projections without correction for point spread function (PSF). The attenuation correction can be performed in the first step by using, for example, the iterative MLEM or ordered-subsets expectation maximization (OS-EM) algorithm. The second step is a postprocessing procedure that corrects for the PSF blurring effect. If the target features (e.g., hot lesions) have a positive background, removing the background before application of the postprocessing filter significantly helps with target deblurring and Gibbs artifact suppression. This postprocessing filter is the image-domain MLEM algorithm. The background activity is attached back to the foreground after lesion sharpening. Results: Computer simulations and PET phantom studies were performed using the proposed 2-step method. The background removal strategy significantly reduced Gibbs artifacts. Conclusion: Gibbs ringing artifacts generated during image reconstruction are difficult to avoid if compensation for the PSF of the system is needed. The strategy of separating image reconstruction from PSF compensation has been shown effective in removal of Gibbs ringing artifacts.
- image processing
- Gibbs artifacts
- ringing artifacts
- Lucy–Richardson algorithm
- MLEM algorithm
- postprocessing
Since Lucy and Richardson's pioneering work (1,2), the Lucy–Richardson algorithm has become popular in image deblurring, because this algorithm provides satisfactory deblurred images under noisy conditions and with stationary and nonstationary point-spread-function (PSF) models. This algorithm is better known as the maximum-likelihood expectation maximization (MLEM) algorithm in the medical imaging community (3,4). Like many other deblurring algorithms, the Lucy–Richardson algorithm suffers from edge-overshoot artifacts (5). The overshoot and ringing artifacts around sharp changes in image contrast are much like the well-known Gibbs phenomenon caused by the truncation of a Fourier series of a discontinuous function. At sharp edges the contrast is high, and Gibbs artifacts are difficult to avoid when the image is being deblurred.
A typical strategy in handling Gibbs ringing artifacts is to reduce Gibbs artifacts with compromised resolution (6). The reduction of Gibbs artifacts can be achieved by blurring the input image so that the data do not contain such high-frequency components. Equivalently, one can use a narrower PSF for the kernel model in the Lucy–Richardson algorithm. Therefore, one must sacrifice some resolution to suppress the artifacts.
Different from the typical strategy, under some circumstances one can obtain a Gibbs-artifact–free image without sacrificing image resolution, taking advantage of the nonnegativity property of the Lucy–Richardson algorithm. This method was first applied by Magain et al. in deblurring astronomic images (7). Those authors pointed out that “This positivity constraint is the main inhibitor of the ringing around point sources: by forbidding the negative lobes, it automatically reduces the positive ones since the mean level must be compatible with the observed data.” However, this strategy cannot be directly applied to the iterative MLEM algorithm that is commonly used in SPECT and PET image reconstruction, because the background cannot be separated and removed from the projections. A new 2-step image reconstruction and deblurring method will be proposed in the next section.
Over the years, many investigators have devoted a lot of effort to combat Gibbs artifacts. Application of the wavelet transform has yielded good results in reducing the Gibbs effect (8). The wavelet transform is a powerful tool but can be difficult to use because of the complexity of selecting the proper basis function and the nonlinear process of filtering the wavelet coefficients. In addition to the wavelet-based filtering methods, there are many modern image filtering strategies, such as the gaussian smoothing method with a Dirichlet integral to measure the smoothness (9), the anisotropic diffusion method (10), the Rudin–Osher–Fatemi total variation method (11,12), the neighborhood and the Wiener local empiric filters (13), the discrete universal denoiser (14), and unsupervised information-theoretic, adaptive filtering (15). Most of these modern filtering techniques have been developed for nonblurred picture processing, which is difficult to incorporate with a specified PSF. These techniques are also sophisticated and have not yet attained a desirable level of applicability in routine medical image processing. This paper will apply the well-developed Lucy–Richardson (i.e., the MLEM) deblurring algorithm in the SPECT and PET reconstructed image domain and propose a new 2-step image reconstruction and deblurring method to reduce Gibbs artifacts.
MATERIALS AND METHODS
The Lucy–Richardson, or MLEM, Algorithm
It is suggested that image reconstruction in SPECT and PET be achieved in 2 steps. In the first step, the image is reconstructed from projections without PSF compensation. The attenuation correction can be performed if the user desires. The user has freedom to choose the image reconstruction algorithm, which can be analytic or iterative. Today the most popular image reconstruction algorithm is the iterative MLEM algorithm or its accelerated version—the iterative ordered-subsets expectation maximization (OS-EM) algorithm.
An image is first reconstructed in such a way that Gibbs ringing artifacts do not appear in the reconstruction. In fact, partial (not full) PSF compensation can be performed during the image reconstruction stage as long as it does not violate our requirement of “no Gibbs ringing artifacts are generated.” This partial PSF compensation can be achieved using a narrower PSF than the full PSF in the image reconstruction algorithm.
If the full PSF is used during image reconstruction, Gibbs artifacts will most likely appear in the reconstructed image. The reconstruction obtained in this first step is referred to as the raw reconstruction.
The main purpose of the first step is image reconstruction, and the main purpose of the second step is image PSF compensation. Hereafter, this paper will focus on the second step of the method.
In the second step of the proposed method, the raw reconstruction is processed by another MLEM algorithm. This MLEM algorithm is somewhat different from the one that is commonly used in SPECT and PET image reconstruction. The MLEM algorithm used in the second step of the proposed method is totally in the reconstructed image domain: both the input data and the output image are reconstructed images of the same dimension. In the rest of the paper, this reconstructed-image-domain MLEM algorithm will be referred to as the Lucy–Richardson algorithm.
The Lucy–Richardson algorithm assumes the following image blurring model:bi=∑jpijaj,Eq. 1where pij is the PSF with the constraint that ∑jpij=1; aj is a pixel in the true, original, nonblurred image; and bi is a pixel in the blurred image. The Lucy–Richardson algorithm is an iterative image-deblurring procedure (trying to recover the true image aj from a blurred image bi), which can be formulated asa(k+1)j=a(k)j∑ipijbi∑npina(k)n,Eq. 2where k is the iteration index and a(k)jis the deblurred image estimation at the kth iteration.
If an input image has a positive background and the Lucy–Richardson algorithm is used for deblurring, Gibbs artifacts are severe in the deblurred image. However, if the background is removed before the Lucy–Richardson image-deblurring procedure is applied, the artifacts are significantly reduced and the image without background is sharper. Therefore, an image-deblurring method is proposed that first removes the target background and then deblurs the target using the Lucy–Richardson algorithm. The background is different in different regions of a SPECT or PET reconstruction. Therefore, the proposed method is a local image processing method. This method works well for point-source targets. Here, a target is a region of interest to be deblurred.
If the target is not a point source but is a relatively large hot region, one can artificially determine a hard upper limit for the image in each iteration of the Lucy–Richardson algorithm. If the updated image pixel value is greater than the preset hard limit, the new image pixel value is assigned as the hard limit value.
Hot Lesion Enhancement
Our proposed method can be referred to as a local image-deblurring technique. The deblurring procedure is presented as follows. An object—say, a lesion—is first identified. Second, the background around the object is removed by threshold segmentation. Third, the iterative Lucy–Richardson algorithm is applied on the segmented region of interest with a specified PSF. The PSF is the kernel that needs to be deconvolved. If the hot lesion is relatively large, a hard upper limit is set as the maximum image value in the unfiltered large hot region. The selection of the upper limit is based on the fact that for a large, flat-top hot region, image blurring does not decrease the maximum image value and blurs only the edges. However, for a pointwise target, image blurring significantly reduces the maximum image value, and it is not proper to set an upper limit. Finally, the segmented background is attached back to the deblurred foreground object.
Cold Lesion Enhancement
The method presented above can potentially be applied to a cold-lesion enhancement task by treating the background as a large hot region and the cold lesion as the background. A cold lesion can also be deblurred by first being converted into a hot lesion, simply by the application of the inverse gray-scale representation of the image. As in hot-lesion enhancement, a region containing the lesions is first segmented. Before removal of the background, the image of the subregion is inverted; that is, the image pixel values are multiplied by −1. Now the cold lesion becomes a hot lesion. After removal of the background, the background is zero. At this point, the iterative Lucy–Richardson algorithm is applied to sharpen the hot lesion with a specified PSF. If the lesion is relatively large, a hard upper limit is set within each iteration of the Lucy–Richardson algorithm. Finally, the sharpened hot lesion is inverted back to a cold lesion and recombined with the background. In SPECT and PET, cold lesion detection is common in myocardial imaging and in some bone tumor imaging. The proposed method can be applied to estimate the size of a cold lesion. A cold lesion imaging example is illustrated in the following numeric results.
RESULTS
Example 1: Blurred Image with Hot Lesions
Figure 1 illustrates a 2-dimensional (2D) example. The phantom had 2 small hot lesions and a positive background. The image was in a 400 × 400 array. There were 2 hot lesions 30 pixels apart. The 2 lesions were identical: 1,000 units of intensity and 1 pixel in diameter. The background was a uniform constant of 1 unit of intensity. The PSF was a normalized gaussian function with an SD of 30 pixels. The number of iterations was 3,000.
Noiseless computer simulation with 2 hot lesions. First row shows images, and second row shows central horizontal profiles of corresponding images in first row. Left: input blurred noiseless image that contains 2 hot pointlike lesions and positive background. Middle: result by direct application of regular Lucy–Richardson algorithm. Right: result of proposed segmentation-filtering method.
When the background of the raw image was not removed, Gibbs artifacts appeared at a large number of iterations. After 3,000 iterations, the lesions could not be resolved. Additionally, the lesions could not be resolved by iterating even further. On the other hand, when the background was removed from the raw image first, the 2 lesions could easily be resolved and no Gibbs artifacts were generated. Because the targets were pointlike lesions, no upper limit was used in the Lucy–Richardson algorithm in this example.
No noise was added in the example shown in Figure 1. After Poisson noise was added to the raw input image, the procedure was repeated. The corresponding results are reported in Figure 2. Even with noise, the proposed method clearly outperformed the regular Lucy–Richardson deblurring method.
Noisy computer simulation with 2 hot lesions. First row shows images, and second row shows central horizontal profiles of corresponding images in first row. Left: input blurred noisy image that contains 2 hot pointlike lesions and positive background. Middle: result by direct application of regular Lucy–Richardson algorithm. Right: result of proposed segmentation-filtering method.
Example 2: Blurred Image with Cold Lesions
A 2D cold lesion example is presented in Figures 3 and 4. No noise was added in the computer simulations in Figure 3. Figure 4 is the noisy version of the scenario shown in Figure 3. As in the first example, the cold-lesion example showed the advantage of removing the background of the raw image before the Lucy–Richardson image-deblurring procedure. Because the targets were pointlike lesions, no upper limit was used in the Lucy–Richardson algorithm in this example.
Noiseless computer simulation with 2 cold lesions. First row shows images, and second row shows central horizontal profiles of corresponding images in first row. Left: input blurred noiseless image that contains 2 cold pointlike lesions and positive background. Middle: result by direct application of regular Lucy–Richardson algorithm. Right: result of proposed segmentation-filtering method.
Noisy computer simulation with 2 cold lesions. First row shows images, and second row shows central horizontal profiles of corresponding images in first row. Left: input blurred noisy image that contains 2 cold pointlike lesions and positive background. Middle: result by direct application of regular Lucy–Richardson algorithm. Right: result of proposed segmentation-filtering method.
The proposed method did not work well in a cold lesion enhancement task when the image was noisy. All ill-conditioned inverse problems are challenging, and cold lesion enhancement seems to be more challenging than hot lesion enhancement when the images are noisy. However, it is clear that the proposed method performed better than the regular Lucy–Richardson method in cold lesion deblurring within noisy environments.
Example 3: PET NEMA Phantom Experiment
A National Electrical Manufacturers Association (NEMA) torso phantom was used in an 18F study, and the projection data were acquired using a Siemens PET scanner. The phantom dimensions were 24.1 cm (height) × 30.5 cm (width) × 24.1 cm (depth). The phantom had 6 fillable spheres with inner diameters of 10, 13, 17, 22, 28, and 37 mm. At the center of the phantom was a cold cylinder 51 mm in diameter. The 4 smallest spheres were filled with a 5-fold higher concentration of 18F than the background. The 2 largest spheres contained no radioactivity.
A 3-dimensional raw image was first reconstructed using experimental PET data with 1 iteration of the ordered-subsets expectation maximization algorithm that had 21 projection data subsets. More than 1 iteration would produce a noisier reconstruction. No PSF compensation was applied during image reconstruction. A 2D slice that contained all 6 spheres was selected for postfiltering experiments.
The 2D images were in 168 × 168 arrays. The input raw reconstruction was stored in a 1-byte integer array with a maximum value of 255 and a minimum of 0. The PSF used in the iterative Lucy–Richardson algorithm was a 2D gaussian function stored in a 21 × 21 float array, and the SD of the 2D gaussian function was 4.5 pixels. A threshold of 68 was selected to segment the hot lesions. After segmentation, the minimum value of the hot lesion base was zero, and the maximum value of the hot lesions was chosen as the upper limit. The number of iterations chosen for the Lucy–Richardson algorithm was based on the trade-off of resolution recovery and noise. The lesions were fairly large compared with the point sources. Lesion images would have appeared noisier if a higher number of iterations had been used. Ten iterations of the Lucy–Richardson algorithm were used to deblur the hot lesions. Then, the sharpened hot lesions were added back to the segmented background.
To sharpen the cold lesions, a new threshold of 50 was chosen to segment the image obtained from the hot lesion deblurring steps. The upper limit was set at 50. The lower-pixel-intensity segmented image was inverted. Ten iterations of the Lucy–Richardson algorithm were applied to the inverted, lower-pixel-intensity segmented image. Finally, the filtered, inverted, lower-pixel-intensity segmented image was inverted again and added to the upper-pixel-intensity segmented image. Here, the chosen number of iterations was based on the resolution recovery and noise trade-off.
In this case, sharpening filtering was applied twice to the 2D image. The first filtering was applied to image intensities in the range of 68–255, and the second filtering was applied to the image in the intensity range of 0–50. The image in the intensity range of 50–68 was not filtered. The filtering results of this PET image are shown in the middle row of Figure 5.
PET 18F NEMA phantom experiment. First column shows images, second column shows profile along upper horizontal line (passing through cold lesion, cold central cylinder, and hot lesion), and third column shows profile along lower horizontal line (passing through 2 hot lesions). Top row: input blurred noisy raw reconstruction image. Middle row: result of proposed segmentation-filtering method. Bottom row: result by direct application of regular Lucy–Richardson algorithm.
The image in the bottom row of Figure 5 was the result of directly applying 10 iterations of the Lucy–Richardson algorithm without segmentation. Direct application of the Lucy–Richardson algorithm generated significant Gibbs ringing artifacts, whereas the proposed segmentation-filtering method was almost free of Gibbs artifacts.
All hot lesions are supposed to have the same radioactivity concentration. Because of the partial-volume effect (i.e., image blurring) in this example, hot lesions in the original blurred raw reconstruction appeared to have different concentrations, as can easily be seen by looking at the line profile drawn across the 2 hot lesions. After deblurring using the proposed method, the line profile (which is labeled as “Profile along lower line” in Fig. 5) shows that the 2 hot lesions now have the same concentration.
However, quantitative recovery of the image for the cold lesions is not as effective as that for the hot lesions. The line profiles labeled as “Profile along upper line” in Figure 5 show that the edges of the cold lesions are not effectively sharpened, and the cold lesion pixel values are above zero. The scattering effect could be making the cold lesions hotter than they should be. Photon scattering was not compensated for in this phantom study.
DISCUSSION
In this paper, the proposed 2-step Gibbs artifact–controlling method is illustrated with 2D simulated point source and PET phantom data only, with simple cases of constant background and isointensity target objects, and stationary PSF. In realistic clinical scenarios, the source distribution is 3-dimensional, the PSF is nonstationary, the background is heterogeneous, and the target objects have variable sizes. For these real-world scenarios, the challenges in applying the proposed method are in the segmentation step. Because of the heterogeneous background, it is impractical to have a universal threshold to perform segmentation. The image-sharpening step must be implemented locally. After the raw image is reconstructed in the first step, the technologist needs to interactively choose some targets to be sharpened in the second step, 1 target at a time. For each target, 2 regions must be circled by the technologist: the target and the adjacent background. This local sharpening approach can be time-consuming and is a drawback of the proposed method.
If the hot lesion and the cold lesion have the same background, according to the Poisson statistics, the hot lesion is less noisy than the background whereas the cold lesion is noisier than the background. The noise level difference makes the performance of cold lesion sharpening worse than that of hot lesion sharpening, as demonstrated in Figures 2 and 4.
As discussed by Magain et al. (7), removal of the background may destroy the Poisson property of the image noise. The Poisson noise model is not an essential requirement for the Lucy–Richardson algorithm (3).
The proposed method is not fully automatic and needs human interaction. When the background is nonuniform, the proposed method can be applied only locally, and the local region must be manually drawn before segmentation.
As shown in the middle image of Figure 5, the 2 largest hot lesions in the PET NEMA phantom study plateau at the same maximum value, which is, in fact, the preset upper limit of the image. Image sharpening is an ill-posed problem and tends to produce noisy images when the data are corrupted with noise. If an upper limit is not set, the target image value will keep increasing as the iteration number increases and the target becomes sharper. Because of noise propagation, the image of the target will get noisier at the same time that the image gets sharper. This phenomenon is often seen in a typical inverse problem. Thanks to the zero-background strategy, no Gibbs ringing artifacts are generated in the resultant image. Gibbs ringing artifacts are the common by-product of an image-deblurring procedure.
The conventional strategy for reducing Gibbs artifacts relies chiefly on compromising image resolution (16). The well-known Lucy–Richardson algorithm is able to provide Gibbs-artifact–free deblurred images if the input image is properly segmented and is iteratively filtered separately. The segmentation idea was inspired by our earlier publication (17), which does not use the Lucy–Richardson algorithm. It seems that the task of cold lesion deblurring is more difficult and less effective than hot lesion deblurring. Like all other deblurring methods, the proposed method is sensitive to image noise; but even with noise, the proposed method performs better than the regular Lucy–Richardson method as indicated by computer simulations and phantom experiments. The Lucy–Richardson algorithm has been widely used in astronomy, and the results of many astronomic studies (7,18,19) can be introduced to our nuclear medicine community.
CONCLUSION
This paper shows that if the nuclear medicine image reconstruction task is split into 2 steps, Gibbs ringing artifacts can be significantly reduced. Gibbs artifacts are better controlled in a postprocessing procedure than in the image reconstruction procedure.
Acknowledgments
I thank Dr. Roy Rowley for English editing of the manuscript and Dr. Michael Casey of Siemens Medical Solutions for providing PET 18F phantom data. This work was supported in part by the Margolis Foundation and NIH grant R21 EB006830. No other potential conflict of interest relevant to this article was reported.
Footnotes
Published online Jul. 27, 2011.
REFERENCES
- Received for publication December 8, 2010.
- Accepted for publication April 4, 2011.
