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
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
Retrospective artifact correction (RAC) improves image quality post
acqu...
read it
Random tessellations are well suited for the probabilistic modeling of
t...
read it
Combining fast MR acquisition sequences and high resolution imaging is a...
read it
Recent developments in the context of minimum residual finite element me...
read it
Image reconstruction is an inverse problem that solves for a computation...
read it
Highly-directional image artifacts such as ion mill curtaining, mechanic...
read it
The mapped bases or Fake Nodes Approach (FNA), introduced in [10], allow...
read it
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
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
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
We are grateful to Valerij G. Kiselev for very helpful discussions.
Comments
There are no comments yet.