Abstract
The aim of this paper is to propose a partitioned-image filtering technique that is a new way of reducing the Gibbs phenomenon in filtered images. Methods: This technique relies on exploiting the properties of the Gibbs phenomenon and on assumptions about the structure of the image. The amplitude of the Gibbs ringing is directly proportional to the height of the image discontinuity. If the height of the discontinuity can be reduced, then the subsequent ringing will also be reduced. Separating the image into stratified layers or partitions reduces the height of the discontinuity significantly. Each partition is filtered separately and recombined nonlinearly to yield the final filtered image. This method weakens filtering of image edges that have large discontinuities, thus reducing the Gibbs phenomenon while simultaneously reducing the image noise. Results: The proposed filtering method has been applied to a simple image with only 2 intensity values to illustrate the implementation steps. The method has also been applied to 2 SPECT patient studies to show the effectiveness of the proposed filtering method, which can significantly reduce the Gibbs artifacts. Conclusion: The Gibbs phenomenon in a filtered image can be reduced by partitioning the image so that the amplitude of the discontinuity is controlled. The proposed method is efficient and simple in implementation, with fast Fourier transform.
Filtering of medical images to remove noise is performed with great care because of the risk of creating artifacts that can lead to misdiagnosis. The most common image artifact, that which arises from low-pass filtering of the image to remove high-frequency noise, is the Gibbs phenomenon, which is manifested by decaying oscillations around high-contrast edges. This paper proposes a method by which noise can be reduced while maintaining the sharp edges of image features and therefore reducing the Gibbs phenomenon.
The Gibbs phenomenon is named for Josiah Willard Gibbs, who explained it in the April 27, 1899, edition of the journal Nature. His letter to the editor was the result of a discussion in the scientific community of the “convergence of the partial sums of certain Fourier series in the neighborhood of [a signal] discontinuity.” It was noticed that for a square wave, the reconstruction showed decaying oscillations at the discontinuities (1). The amplitude of the ringing is usually in the neighborhood of 9% of the jump discontinuity (2).
The Gibbs phenomenon is usually experienced during attempts to remove noise from an image. An image contains primarily lower frequencies, with a small but significant portion of the image being of higher frequency. This high-frequency image content is responsible for the proper representation of high-contrast edges. White noise, on the other hand, is present in equal magnitudes at all frequencies; therefore, traditional noise reduction schemes involve low-pass signal filtering that passes most of the image or signal content while indiscriminately removing the high-frequency noise and image content. The removal of the upper-frequency content leads to poorly reconstructed image edges; therefore, linear low-pass filtering has a noise reduction/edge distortion tradeoff.
In Fourier domain filtering, the frequency response of an ideal low-pass filter is a rectangular window that has no transition region between the pass-band and stop-band. In the spatial domain, this kind of filtering causes maximal ringing. To reduce the spatial domain ringing, the frequency domain rectangular window is tapered to include a transition region between the pass and stop bands. This does reduce the spatial domain ringing but allows small amounts of frequencies above the frequency cutoff point to pass. This practice of tapering the rectangular window is the most common way of reducing the Gibbs phenomenon and is found in all introductory books on digital signal processing (3).
Several techniques, such as the wavelet transform, allow filtering while reducing or avoiding the Gibbs phenomenon. The application of the discrete wavelet transform to signal processing has yielded good results in avoiding the Gibbs effect. The discrete wavelet transform decomposes a signal with a variety of basis functions that have compact support in both the signal and the frequency domains. This local support allows the wavelet transform to identify not only the frequencies present in a signal but also the location of the frequencies. The high-frequency components of signal can be identified and isolated from the high-frequency noise. This provides the groundwork for the capabilities of local filtering (4,5). 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 strategies for filtering images, such as the gaussian smoothing method with a Dirichlet integral to measure the smoothness (6), the anisotropic diffusion method (7), the Rudin–Osher–Fatemi total variation method (8) and its iterated version (9), the neighborhood and the Wiener local empiric filters (10), the discrete universal denoiser (11), and unsupervised information-theoretic, adaptive filtering (12). Many of these methods have been developed for pictures comprising flat regions. MRI and CT images fall into this category, but emission nuclear medicine images do not.
Linear filtering suffers from the Gibbs phenomenon and has a noise reduction/edge distortion tradeoff. Nonlinear filtering largely avoids the Gibbs phenomenon but often suffers from telltale nonuniform (or nonlinear) noise reduction. For example, the median filter oversmooths some areas while leaving noisy areas near the edges less filtered.
We will stay with the popular and classic frequency-domain filtering methods and propose a remedy to reduce the Gibbs errors. Herein, a method is proposed that limits the effect of the Gibbs phenomenon while appropriately reducing noise. The proposed method uses an important property of the Gibbs phenomenon and exploits some prior knowledge of the structure of the image. This property of the Gibbs phenomenon is that the amplitude of the ringing is proportional to the height of the contrast edge. Our technique applies the conventional linear filtering method to nonlinearly partitioned images. We believe that our proposed method is a valuable contribution to this Gibbs phenomenon–reducing family of filtering techniques in medical imaging applications.
MATERIALS AND METHODS
A digital image is generally encoded as a matrix of gray-level values. Mathematically, an image consists of ordered pairs {i, p(i)}, where p(i) is the value at pixel i. A medical image is composed of distinct regions of similar image intensity p(i), which is referred to as an image feature in this paper. The proposed method partitions the image in an attempt to separate these distinct regions, filters each partition (or region) separately using any noise-reducing filtering method desired, and then recombines the significant portion of each partition into the final image. The degree to which the image conforms to the assumption of a piecewise smooth function determines the effectiveness of the proposed method. The interesting aspects of this partitioned-image filtering technique are the assignments of the unassigned pixels in each partition, the partition recombination scheme, and the determination of the partition threshold levels.
For simplicity, we illustrate our method by considering a simple noiseless image consisting of 2 distinct features as shown in Figure 1, where there is a large jump between these 2 features.
(A) Simple image consisting of 2 intensity values. (B) Profile drawn through image.
Step 1
Select a threshold value T = 60 between the 2 features (P = 20 and P = 100). See Figure 2.
(A) Threshold T = 60 set between intensity values of 20 and 100. (B and C) Images after image partitioning using threshold T: P1 (B) and P2 (C).
Step 2
The original image is partitioned by assigning all pixels below the threshold T to partition P1 and all pixels above to partition P2. The unassigned pixels in each image partition are set to the threshold value T. Here, P1 and P2 are 2 new images (Fig. 2). Each partition is bounded by its upper and lower threshold (here, P1 ← [0, T] and P2 ← [T, max(image intensity)]). Therefore, the unassigned pixels are set to the upper or lower bound depending on whether the original image pixels are above or below the partition pixel. This assignment mimics the continuity of the original image and achieves the desired result of reducing the height of the contrast edges. The jumps in the partitioned images are half the size of the jump in the original image.
A membership map {i, m(i)} is created for each partitioned image. In a particular partitioned image, if a pixel value p(i) is from the original image, then m(i) = 1, otherwise m(i) = 0.
Step 3
Each partitioned image is filtered separately. A simple fourth-order low-pass Butterworth filter of cutoff frequency 0.35 is used in this example. See Figure 3.
Low-pass Butterworth filter applied to partitioned images separately: filtered P1 (A) and filtered P2 (B).
Step 4
The filtered partitioned images are combined into the final image pfinal using the formula:
where m1 and m2 are the membership images for P1 and P2, respectively, as defined in Step 2. Because of the membership map weighting, this combination procedure is not simply the adding of 2 filtered images together.
This final image is compared with the direct filtering of the original image in Figure 4, where the same Butterworth filter is used for both methods. The Gibbs ringing effect is reduced with the proposed method. The overshoots and undershoots are reduced to half of those for the conventional direct filtering method.
Filtered images. (A) Filter is directly applied to original image shown in Figure 1. (B) Filter is applied to partitioned images shown in Figure 3, and filtered images are recombined according to their membership maps. (C) Close-up profiles to compare the 2 methods. Solid curve corresponds to direct filtering method, and dashed curve corresponds to proposed method.
The threshold is manually selected by visually discerning the intensity of various image features. Autodetection of image features for this method is currently under development.
To decide where to partition the image, it is helpful to study a directly filtered image without applying a partition. It is easier to determine the thresholds using the directly filtered image. The best place to partition is one that separates image features and reduces signal discontinuities. The partitioning process collects samples that lie between consecutive partitioning thresholds and places them in a partition. The undefined samples in each partition are given the upper or lower threshold depending on whether the original sample in that location is in an upper or lower partition. This defining of the undefined sample behaves somewhat similarly to the actual signal and does not introduce additional signal defects. The membership map (that recorded the partition to which each original sample went) is consulted to recombine the partitions and to collect the samples into the final signal.
RESULTS
As an example, the proposed method is applied to a medical image in Figure 5A. This is a clinical SPECT image of a human torso with an external marker used to align the images, and the marker is identified by the arrow.
(A) SPECT image of human torso with external marker (arrow). (B) Result from applying low-pass Butterworth filter to SPECT image shown in A. Large Gibbs phenomenon or decaying oscillation is seen around marker. (C) Result from applying proposed partitioned image filtering technique using same low-pass filter. Decaying oscillation around marker is reduced.
The artificial point source marker was attached to the patient's skin for coregistration purposes, and the intensity of this marker is several orders of magnitude larger than that within the patient body. When a low-pass Butterworth filter is applied to smooth the noise in the image, the Gibbs phenomenon around the marker is severe (Fig. 5B). When the proposed method is applied to the same image and the same low-pass filter is used on each of the partitions, the resulting image achieves the reduction in noise while minimizing the Gibbs phenomenon (Fig. 5C).
Another example is a patient bone SPECT image (Fig. 6). The original image, Figure 6A, was reconstructed with a filtered backprojection algorithm. A low-pass Butterworth filter was then used to obtain Figure 6B, and this is the conventional filtering method. The same Butterworth filter with the proposed portioned filtering method was used to obtain Figure 6C. From the zoom-in images, we can see that the negative undershoot Gibbs artifacts around the ribs after Butterworth filtering are significantly reduced with the partitioned filtering method.
Example of partitioned filtering of patient bone SPECT image: raw filtered backprojection reconstruction (A), direct low-pass Butterworth filtering (B), and proposed partitioned low-pass Butterworth filtering (C). D, E, and F are zoom-in subimages of A, B, and C, respectively.
DISCUSSION
It is not the aim of this paper to tell readers how to select the best postprocessing filter for their nuclear medicine images. In our hospital, the technologists prepare and process the images for the physicians to read. The technologists use a slide-bar to adjust the filtering parameters until the images on the screen appear satisfactory. Because our proposed filtering method uses the conventional fast Fourier transform–based filtering technique, the technologists can select the filtering parameters and process the images in their usual way. If the Gibbs phenomenon is severe, the suggested solution is that the image be partitioned into regions of different intensities before filtering is applied.
Whether a filter should be used or which filter should be used under a certain condition is not a topic of discussion here. This paper assumes that an appropriate filter has already been somehow selected and is helpful in producing the final image but that the side effect of the filtering—the Gibbs phenomenon—needs some attention. This paper suggests a technique to reduce the Gibbs phenomenon while the predetermined filter is still used, by applying the predetermined filter to the partitioned images. Figure 5 illustrates this point clearly. A Butterworth filter in the example can effectively suppress the noise and displays the useful features of the image. Unfortunately, the image contains an artificially attached marker with an extremely high intensity. If the Butterworth filter is directly applied to the entire image, the Gibbs phenomenon will be severe. If the partitioning technique is used before Butterworth filtering, the Gibbs phenomenon is greatly reduced.
Our proposed algorithm is closely related to the bilateral filtering method developed by Tomasi et al. (13). Bilateral filtering uses a filtering kernel that depends on both the domain distance (i.e., the distance between 2 pixels) and the range distance (i.e., the image intensity values of these 2 pixels). The main difference between our method and the bilateral filtering method is the treatment of range filtering. Our method does not filter in the range dimension but segments the image in the range dimension. Because of segmentation, our method is nonlinear.
In the bilateral filtering method, the kernel depends on the domain and range variables, and thus the common Fourier transform–based filtering methods cannot be used. On the other hand, our method filters only in the domain dimension and is the same as conventional linear filtering. Therefore, the efficient and simple fast Fourier transform can be used in our method.
Lukin further developed the bilateral filtering method of Tomasi et al. and made it into a multiresolution approach (14). In his multiresolution method, 2 images are obtained via the bilateral method: one with high-frequency noise filtered out and the other with low-frequency noise filtered out. These 2 images are then combined using mixed wavelet coefficients from the first and second images. In contrast, our method combines both high-frequency filtering and low-frequency filtering as a single conventional Fourier domain filtering procedure. The primary feature of our proposed method is its simplicity.
CONCLUSION
The main purpose of the proposed technique is to reduce the Gibbs artifact due to postfiltering. In nuclear medicine, filters are commonly applied to images to suppress noise. If no obvious Gibbs phenomenon is observed, no image partitioning is necessary. If filtering is necessary and the Gibbs phenomenon is severe, we suggest use of the proposed partitioning method to reduce the Gibbs overshoots and undershoots.
The Gibbs phenomenon is caused by the removal of high-frequency image content needed to reconstruct image edges and is identified by decaying oscillations along these edges. Our proposed filtering technique exploits the clustered image intensities of various structures in a medical image and the properties of the Gibbs phenomenon to smooth the noise while reducing the oscillations. For example, in CT images, intensities can be clustered into 3 groups: bones, soft tissues, and air. After each of these has been identified and separated, they are filtered separately and then recombined into the final image. Application of the method to SPECT medical images reduces the Gibbs phenomenon.
This partitioned-image filtering technique opens new ways of looking at signal filtering and provides exciting ideas for future work in this area. We have shown, through reasons and examples, that the Gibbs phenomenon can be reduced in an image by partitioning, independent processing, and nonlinear recombination.
Acknowledgments
We thank Dr. Chuanyong Bai of Digirad Corp. for providing us clinical images.
Footnotes
-
COPYRIGHT © 2009 by the Society of Nuclear Medicine, Inc.
References
- Received for publication December 22, 2008.
- Accepted for publication March 16, 2009.
