Gibbs-Ringing Artifact Removal Based on Local Subvoxel-shifts

01/30/2015 ∙ by Elias Kellner, et al. ∙ 0

Gibbs-ringing is a well known artifact which manifests itself as spurious oscillations in the vicinity of sharp image transients, e.g. at tissue boundaries. The origin can be seen in the truncation of k-space during MRI data-acquisition. Consequently, correction techniques like Gegenbauer reconstruction or extrapolation methods aim at recovering these missing data. Here, we present a simple and robust method which exploits a different view on the Gibbs-phenomena. The truncation in k-space can be interpreted as a convolution with a sinc-function in image space. Hence, the severity of the artifacts depends on how the sinc-function is sampled. We propose to re-interpolate the image based on local, subvoxel shifts to sample the ringing pattern at the zero-crossings of the oscillating sinc-function. With this, the artifact can effectively and robustly be removed with a minimal amount of smoothing.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 5

page 6

page 7

page 8

Code Repositories

gibbs_removal

Javascript implementation of Gibbs ringing removal from image volumes based on the method described by Kellner et. al. at https://arxiv.org/abs/1501.07758


view repo
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

Abstract

Gibbs-ringing is a well known artifact which manifests itself as spurious oscillations in the vicinity of sharp image transients, e.g. at tissue boundaries. The origin can be seen in the truncation of k-space during MRI data-acquisition. Consequently, correction techniques like Gegenbauer reconstruction or extrapolation methods aim at recovering these missing data. Here, we present a simple and robust method which exploits a different view on the Gibbs-phenomena. The truncation in k-space can be interpreted as a convolution with a sinc-function in image space. Hence, the severity of the artifacts depends on how the sinc-function is sampled. We propose to re-interpolate the image based on local, subvoxel shifts to sample the ringing pattern at the zero-crossings of the oscillating sinc-function. With this, the artifact can effectively and robustly be removed with a minimal amount of smoothing.

Key words: Gibbs-Ringing — Ringing-Artifact — Truncation-Artifact

Introduction

In MRI, images are not gained directly, but with a reconstruction from acquisitions of corresponding Fourier expansion coefficients in k-space,

[1]

where denotes the image value at voxel index and are the expansion coefficients. In practice, only a finite number of expansion coefficients can be acquired. This truncation of Fourier space introduces artifacts if the expansion coefficients do not decay fast enough with increasing [1]. This is the case for sharp image transitions, where all higher frequency components are required to properly reconstruct the edge. The artifact can nicely be illustrated when reconstructing the image on a fine grid using zero-filling by setting the missing high frequency components for to zero. This operation corresponds to a multiplication of the ‘true’ k-space with a rectangular window, which is equivalent to a convolution with a sinc function in the image domain. The side lobes of the sinc result in oscillations (‘ringing’) in the neighborhood of sharp edges. The issue is illustrated in Fig. 1 for the image of a single edge.

Several approaches have been proposed to reduce ringing artifacts. Most straightforwardly, image filters, e.g. the Lanczos -approximation [2, 3] or a median filter can be used to smooth the oscillations. However, as filtering implies global blurring, the spatial resolution is effectively reduced, leading to a loss of fine image details. More advanced methods have been developed based on piecewise re-reconstruction of smooth regions using Gegenbauer-Polynomials [2, 4]. A drawback of these methods is the requirement of an edge detection and potential instabilities in some applications [5] due to the involved choice of parameters. Another modern approach consists of combined de-noising and ringing removal using using total variation constrained data exploration [6].

The method we propose in this work is based on a different view on the effect: In fact, due to the finite number of expansion coefficients, the image is not reconstructed on a continuous domain, but on a discrete grid. Obviously, the strength of the ringing in the reconstructed image depends on the precise location of the edge relative to the sampling grid, i.e. how the sinc-function is sampled. If the side lobes of the sinc-pattern are sampled at its extrema, the ringing amplitude becomes maximal, whereas it disappears when sampled at the zero crossings (see Fig. 1). This feature has been demonstrated in experimental data in the context of the dark rim artifact in cardiovascular imaging [7]. Finding the optimal subvoxel-shift for pixels in the neighborhood of sharp edges in the image can therefore minimize the oscillations. However, as there are multiple edges present in an image, the correction cannot be achieved by a single, global shift, but must be performed on a local basis. In the next section, we describe a non-iterative method which conducts this task.

Figure 1: An image with a single discontinuity (black edge) is reconstructed from truncated k-space data. The resulting image (blue dots) is discrete and exhibits ringing artifacts. The amplitude of the ringing depends on whether the underlying sinc pattern arising from the windowing in k-space (red curve) is sampled at its extrema (A), or at the zero-crossings (B).

Methods

One-dimensional Case

Let be the original, discrete image, and its Fourier expansion coefficients. From this, a set of images , where , with subsequent subvoxel-shifts is created by multiplication with phase-ramps in Fourier-space:

[2]

From this dataset, for each pixel , the optimal shift which minimizes potential oscillations in the neighborhood is found. The corresponding measure can be calculated with any oscillation-sensitive kernel, e.g. the absolute differences in the neighborhood. It seems advisable to measure the ringing to both sides of the pixel individually, and select the smaller value. This way, in the neighborhood of an edge, the ringing is always measured in the direction opposite to the edge (see Fig. 2. Thus, the step itself does not contribute, and interferences from closely located edges are minimized. We decided to use a simple absolute differences approach to measure the oscillations.

where defines the window size which is used for measuring the amount of oscillations. There are reasons to exclude the central point itself, using a window of the form . This leads to more stability for points directly on the edge, and minimizes increase in noise correlation, as the oscillation measurement itself is not correlated with the actual point we are interested in. The choice of the window constitutes the only parameters of the method. The results will indicate that the proposed algorithm is relatively robust against this choice.

Now, for each point the optimal shift is determined by finding the shift such that the total variation is minimal. First, independently for the right (+) and left (-) side:

[3]
[4]

Finally, we decide whether the overall minimum comes from left or right and this shift is determined to be the optimal shift . Now, we know the optimal shift at each image location, which minimizes TV and hence also the ringing artifact. But, finally we have to go back to original grid, i.e. evaluating at the non-integer position . In order to do so we can use any kind of alternative interpolation which is not a sinc-interpolation. Formally, we have

We decided to do this final ’back interpolation’ by a simple linear interpolation. In fact, one could also use alternative high-order interpolation schemes (like polynomial interpolation), but they usually lead to ringing artifacts again.

Figure 2: 1-D image with two edges. A set of images with increasing subvoxel-shifts is created. The optimal shift is different for the two edges. For the right edge, image A) is optimal whereas it is C) for the left edge. Consequently, the optimization must be performed on a local basis, i.e. for each pixel individually. The optimization criteria is given by the amount of ringing, which can be characterized by calculating the absolute differences in a certain range for each pixel. In C), this is exemplary illustrated for the red pixel. The total variation is calculated separately to both, the left and right side of the center pixel. This ensures that the ringing is measured away from the edge, and the edge itself is not contributing. Further, the center pixel itself is excluded in order to minimize the introduced correlations between neighboring pixels. Finally, the pixel value is recalculated by interpolating the shifted image D) at the original pixel position, which in this example is 0.5 pixels to the left side (see red dot at bottom). Mathematical details are given in the text.

Two-dimensional Case

The extension to the two dimensional case is not straightforward as diagonal edges produce checkerboard-like ringing patterns, as can be observed in the phantom image in Fig. 4). Hence, it is not possible to find the optimal shift in both dimensions simultaneously. As a solution, we correct the image in both dimensions separately resulting in two one-dimensionally corrected images and . These are then combined in Fourier space to the final image via

[5]

where

denotes the Fourier transform. We propose to use the ‘weighting functions’

and with a saddle-like structure in Fourier Space of the form

[6]

This way, the high frequency components along the direction of the correction are enhanced, while they are dampened along the non-corrected direction. Due to the normalization , artifact-free images, where , are left unchanged by Eq. 6. This ensures that minimal smoothing is introduced.

Figure 3: Weighting functions (left) and (right) for multiplication with the 1D-corrected images in Fourier Space. High frequencies are enhanced along the direction of the correction and suppressed along the other direction.

Reference methods

We compare the proposed method against other, non-iterative filtering methods. In order to minimize unwanted smoothing, low-pass filter with a high cutoff-frequency seem preferable. This can be achieved with the Lanczos sigma-approximation [2, 3], where the filtered image is given by

[7]

where the parameter controls the filter strength. In this study, we use two different settings for . On one hand, we apply the standard choice of . This induces, however, a rather strong smoothing of the image. Hence, for a sound comparison with the proposed method, we further set such that the increase in the correlation of the noise are the same for both, the proposed method and the filtering approach. The corresponding reference noise correlation was measured in a pure-noise region on the mean-free image via

[8]

Another popular, non-linear filtering approach which preserves edges is given by the median filter. For comparison, we also applied a median filter, where the ‘width’ of the filter was fixed to a 2x2 neighborhood.

Numerical Phantoms

The method was applied to two numerical phantoms (Fig. 4). For the first phantom, a polygonal shape with some stripes and small structures was simulated. Starting with a high-resolution image, the artifact was simulated by reconstructing the image from a truncated k-space with a reduction factor of 20. Second, for a more realistic brain phantom, data from a -weighted post-contrast MRI measurement of the brain was used. The artifact was artificially enhanced by re-reconstructing the image from a smaller k-space with a reduction factor of 4. In both cases, a ‘ground truth’ image without artifact, but with the same decreased spatial resolution was generated by convolving the high-resolution image with a boxcar function, and sampling the result on the corresponding low-resolution grid. Gaussian noise with a signal-to-noise ratio of 100 was added.

MRI Measurements

The method was applied to diffusion-weighted-images (DWI) with 70 directions, , using a gradient echo EPI sequence with , matrix size 104x104, resolution , performed on a 3T scanner (Siemens TIM TRIO, Siemens, Erlangen, Germany). No distortion correction was applied, as the involved correction methods already lead to significant filtering of the artifact, especially in phase direction. Due to the different image contrast for different -values, the artifact might even be amplified during post-processing of the diffusion parameters. Therefore, we also calculated diffusion maps, without and with artifact correction. We further applied the method to a -weighted image acquired with a turbo-spin-echo sequence, , resolution 1x1x5 on a 1.5T scanner (Siemens SONATA, Siemens, Erlangen, Germany). Both dataset were acquired in the context of clinical routine, written consent was obtained to use the data for scientific use.

Results

Numerical Phantoms

The results of the phantom simulations are shown in Fig. 4. We show results using the median filter, the Lanczos approximation with , results obtained with the proposed method using different parameters, and results using the Lanczos approximation with filter parameters adapted to yield an equal noise correlation as the proposed method.

Obviously, the median filter preserves the edges better than the Lanczos-approximation with , at a smaller increase in the average smoothing, indicated by the smaller increase in noise correlation. However, it shows a stronger residual of the artifact, and fine image details like small, peak-like structures are destroyed. With the proposed method, on the other hand, the artifact can effectively be removed with minimal smoothing of edges. The method is rather robust against the choice of the kernel parameter. A larger neighborhood results in less smoothing, but comes at the price of slightly reduced artifact removal. A kernel size of seems to be a appropriate compromise between artifact removal and noise correlation. This setting is used for the application to the MRI images.

These findings are basically the same for the phantom constructed from the -weighted image in Fig. 5. Also here, the artifact can most effectively be removed using the proposed method, while preserving fine image details.

MRI images

The results for DWI measurements are shown for one slice in Fig. 6. Apparently, the -images exhibit strong ringing artifacts, which is even more emphasized after diffusion calculation. The artifact can be reduced with both, the median filter and the Lanczos approximation with , however, at the cost of strong smoothing. With the proposed method on the other hand, the artifact can virtually completely be removed with minimal filtering. Results from the -weighted image given in Fig. 7. The findings are basically the same as for the DWI measurement.

Figure 4: Phantom image with multiple edges and noise. Inserts correspond to a scaled area with pure noise, correlation. The amount of noise correlation (given by r) reflects the strength of the smoothing introduced by the methods and their parameters. Both, median filter and Lanczos approximation with lead to rather strong filtering. The median filter additionally introduces artifacts on point-like structures (see e.g. the black dots in the images). Results using the proposed method are given in the second row, the Lanczos approximation with corresponding equal noise correlation in the bottom row. With the proposed method, the artifact can virtually completely be removed with preservation of the edges and without destroying fine image details, whereas the Lanczos approximation with equal noise correlation still suffers from significant artifacts.
Figure 5: weighted post-contrast MRI images with artificially enhanced artifact. Here, the same characteristics as in Fig. 4 can be observed. Both, median filter and Lanczos approximation with lead to blurring of image details (compare e.g. areas at red arrows in ground truth). Again, in the Lanczos approximation with adapted parameters, there is still residual ringing visible (see e.g. blue arrows), even though the difference to the proposed method is here less pronounced compared to the phantom in Fig. 4.

b=0 images                                          Radial diffusivity
 

Figure 6: DWI measurement. The original image exhibits strong ringing artifacts, which are even more amplified in the calculated radial diffusivity maps (right maps). The artifact can most effectively be removed with the proposed method at a minimal image smoothing. The Lanczos approximation shows residual artifacts.
Figure 7: -weighted image with obvious ringing artifacts. With the proposed method, the artifact can effectively be removed. The Lanczos approximation shows residual artifacts (see e.g. blue arrow).

Discussion

Even though the Gibbs-ringing artifact is omnipresent in MRI, most vendors do not include removal techniques in the standard image reconstructing pipeline to date. One reason for this might be that standard filtering approaches inevitably reduce the effective image resolution, and more advanced methods like Gegenbauer re-reconstruction [4, 2, 8, 9] or data extrapolation methods [10, 11] are practically difficult to handle due to their complexity and, the requirement of an edge detection, and potential instabilities induced by the dependency on the choice of the parameters.

Another approach consists of optimizations based on total variation [6]. These methods have proven to effectively remove the artifact, however, they treat noise and artifact equally, in contrast to the proposed method, which explicitly aims at separating both contributions. Further, the outcome of total variation applications strongly depends on the strength of the filter, which must be adapted to the respective application.

The proposed method is rather robust to the choice of its parameter, the kernel width . Further, this parameter is independent of the image size, as oscillation pattern of the ringing occurs always in the distance of one voxel, and hence scales with the matrix size. The method can therefore applied to any image with a universal value for . We found that constitutes a good compromise between artifact removal and smoothing.

Conclusions

In this work, we presented a non-iterative method for removal of ringing artifacts based on re-sampling the image such that the source of the ringing pattern, the sinc-function, is sampled at its zero crossings. We demonstrated that the method effectively removes the artifact while introducing minimal smoothing. The method has a low computational cost, consists of a rather simple mathematical framework and is very stable against the choice of its few parameters. This robustness suggests it as a suitable candidate for a robust implementation in the standard image processing pipeline in clinical routine. Even though designed in the context of MRI, the method might also prove its applicability in other areas such as Fourier-based data compression algorithms.

Acknowledgements

We are grateful to Valerij G. Kiselev for very helpful discussions.

References

  • 1. Gottlieb D, Shu CW, Solomonoff A, Vandeven H. On the Gibbs phenomenon I: recovering exponential accuracy from the Fourier partial sum of a nonperiodic analytic function. Journal of Computational and Applied Mathematics. 1992;43(1):81–98.
  • 2. Gottlieb D, Shu CW. On the Gibbs phenomenon and its resolution. SIAM review. 1997;39(4):644–668.
  • 3. Jerri AJ. Lanczos-like -factors for reducing the Gibbs phenomenon in general orthogonal expansions and other representations. Journal of Computational Analysis and Applications. 2000;2(2):111–127.
  • 4. Archibald R, Gelb A. A method to reduce the Gibbs ringing artifact in MRI scans while keeping tissue boundary integrity. Medical Imaging, IEEE Transactions on. 2002;21(4):305–319.
  • 5. Boyd JP. Trouble with Gegenbauer reconstruction for defeating Gibbs’ phenomenon: Runge phenomenon in the diagonal limit of Gegenbauer polynomial approximations. Journal of Computational Physics. 2005;204(1):253–264.
  • 6. Block KT, Uecker M, Frahm J. Suppression of MRI truncation artifacts using total variation constrained data extrapolation. International journal of biomedical imaging. 2008;2008.
  • 7. Ferreira P, Gatehouse P, Kellman P, Bucciarelli-Ducci C, Firmin D. Variability of myocardial perfusion dark rim Gibbs artifacts due to sub-pixel shifts. Journal of Cardiovascular Magnetic Resonance. 2009;11(1):1–10.
  • 8. Shizgal BD, Jung JH. Towards the resolution of the Gibbs phenomena. Journal of Computational and Applied Mathematics. 2003;161(1):41–65.
  • 9. FENG Qj, HUANG X, FENG Yq. A Fast Algorithm to Reduce Gibbs Ringing Artifact Based on the Chebyshev Polynomials. Journal of Image and Graphics. 2006;8:013.
  • 10. Amartur S, Haacke EM. Modified iterative model based on data extrapolation method to reduce Gibbs ringing. Journal of Magnetic Resonance Imaging. 1991;1(3):307–317.
  • 11. Constable R, Henkelman R. Data extrapolation for truncation artifact removal. Magnetic resonance in medicine. 1991;17(1):108–118.