# Suppression of MRI Truncation Artifacts Using  Total Variation Constrained Data Extrapolation

Authors: Kai Tobias Block; Martin Uecker; Jens Frahm
Journal: International Journal of Biomedical Imaging
 (2008)

---

## Abstract

The finite sampling ofk-space in MRI causes spurious image artifacts, known as Gibbs ringing, which result from signal truncation at the border of k-space. The effect is especially visible for acquisitions at low resolution and commonly reduced by filtering at the expense of image blurring. The present work demonstrates that the simple assumption of a piecewise-constant object can be exploited to extrapolate the data in k-space beyond the measured part. The method allows for a significant reduction of truncation artifacts without compromising resolution. The assumption translates into a total variation minimization problem, which can be solved with a nonlinear optimization algorithm. In the presence of substantial noise, a modified approach offers edge-preserving denoising by allowing for slight deviations from the measured data in addition to supplementing data. The effectiveness of these methods is demonstrated with simulations as well as experimental data for a phantom and human brain in vivo.


---

## Body

## 1. Introduction

In MRI, spatial information is obtained from the
object using magnetic field gradients, which link the Larmor frequency of the
excited spins to their spatial location. Thus, the received signal is the
continuous Fourier transform of the object's proton density(1)S(k→)=∫ρ(x→)e−ik→⋅x→dx→,where the k-space position k→ can be calculated from the time course of the
applied gradients. In practice, the proton density is further modulated by spin
relaxation, off-resonance effects, and other mechanisms all neglected here.It is well known that objects with compact support
have a Fourier transform with nonlimited support. For example, the Fourier
transform of a rectangle is composed of sinc functions in each dimension.
Because only a single location of the Fourier space can be measured at a time,
it is impossible to fully sample such Fourier transform by travelling the MRIk-space with magnetic field gradients. Hence, there are two experimental
restrictions for MRI. First, the continuous Fourier transform is sampled
discretely, which can be seen as a multiplication with a comb-function in frequency
space. In image space, this corresponds to a convolution with a reciprocally
spaced comb-function and leads to periodic object copies with a spacing inverse
to the sample distance in k-space. Second, the Fourier transform is sampled
only within a finite region around the k-space center with all other
information missing.In the conventional case, a discrete Fourier
transformation of the finitely measured data is performed to reconstruct an
image. This strategy implicitly assumes that the Fourier transform is zero
everywhere outside the sampled region. It is clear that the assumption is not
very appropriate for finite objects, although the corresponding reconstruction
totally complies with all data measured. In fact, any solution that coincides
at the sampling positions is a valid reconstruction, because the finite
sampling pattern opens degrees of freedom from the null space of the projection
evoked by finite sampling. Setting this null space to zero is a simple and
convenient solution. Unfortunately, however, the procedure corresponds to a
multiplication of the true object's Fourier transform with a rect-function (in
case of Cartesian sampling) which, in image space, results in a convolution of
the true object with a sinc-function. This effect is well known as truncation
artifact or Gibbs ringing and mainly appears
as an oscillating overshoot of the image intensity near
discontinuities [1, 2]. Although the problem may be
reduced by increasing the measured k-space, many practical applications still
rely on acquisitions with a relatively low-matrix resolution in at least one
image dimension, and therefore suffer from respective artifacts.So far, various methods have been developed to
ameliorate image disturbances due to finite sampling [2–5]. However, in the majority of
MRI applications and, in particular, for most commercially available MRI
systems, only a simple data filtering is routinely employed. In this case,
visual reduction of the ringing artifacts is achieved by a smearing of the
intensity oscillations, which leads to an undesired loss of image resolution.
Alternative methods attempt to extrapolate the measured data and thereby avoid
a sharp cut-off in k-space [6–9]. A key difference to the filtering approach is that
the actually measured data is not changed but supplemented with synthetic data—a reasonable strategy as the measured data is not incorrect but only
incomplete. This can be achieved by exploiting a priori knowledge about the
true object and, consequently, all extrapolation techniques rely on certain
assumptions, where the existing methods follow different strategies. In this
regard, the present work demonstrates that also the very unspecific assumption
of a piecewise-constant object can be utilized to successfully extrapolate data
in k-space and concomitantly reduce the ringing artifacts without compromising
image resolution.


## 2. Theory

Figure1 compares the one-dimensional profile of a
rectangle reconstructed by Fourier transformation from only 96 Fourier samples
to that of the original function. It clearly illustrates severe ringing
artifacts, although the true function is piecewise constant and free of any
oscillations. Such oscillations can be quantified using the total variation
(TV), which sums the modulus of jumps between all neighboring pixels of a
reconstructed image I(x,y)(2)TV(I)=∑y=0N∑x=0N|I(x,y)−I(x−1,y)|+|I(x,y)−I(x,y−1)|.The TV concept was initially
introduced to image processing by Rudin et al. [10] for denoising applications
because noise patterns create a high TV value relative to that of a noise-free
image, and they become particularly reduced when modifying the image in such a
way that the TV value is minimized. As a specific property, edges are preserved
during this procedure, and thus TV minimization emerged as one the most popular
denoising techniques. In recent years, the TV concept is attracting strong
interest in the field of compressed sensing [11] because for specific sampling techniques, the TV
value can be utilized to identify and to remove artifacts from undersampling,
offering a remarkable reduction of the measurement time [12]. In a similar manner,
truncation artifacts lead to an increased TV value relative to that of the true
object, so that the TV may also be taken as a measure of the artifact strength
for finite k-space sampling, which has been recognized by Landi et al. as
well [13]. Therefore,
the proposed idea is to add a set of synthetic frequencies v→ to the measured data y→,
which is specifically chosen such that the TV value of the image reconstructed
from the combination of the measured and synthetic data is minimized(3)v→=argminv→TV(ℱ{v→⊕y→}),where ℱ denotes the discrete Fourier transformation.
Interestingly, by searching for the set of synthetic frequencies v→,
the unmeasured k-space data is recovered to a certain degree if the assumption
of a piecewise-constant object is appropriate.Figure 1
(Top) One-dimensional profile of a rectangle
reconstructed by Fourier transformation from 96 Fourier samples (solid line) in
comparison to the true function (dotted). While the true function is piecewise
constant, the Fourier reconstruction exhibits severe ringing artifacts due to
truncation of the Fourier coefficients, which causes an increased total
variation (TV) value. (Bottom) Magnified view.Estimation of the synthetic data can be achieved by
minimizing3 with a nonlinear numerical optimization technique. The present
proof-of-principle implementation used the CG-Descent algorithm [14], which is a recent variant
of the nonlinear conjugate gradient method that allows to rather efficiently
solve large-scale problems. The algorithm can be used in a black-box manner,
requiring only the evaluation of a cost function and its gradient for given
estimate vectors v→.
The cost function is needed to quantify the goodness of a given estimate (i.e.,
it is small for a good estimate and large otherwise), and for the problem
defined in 3 it simply has the form(4)Φ(v→)=TV(ℱ{v→⊕y→}).The gradient of the cost
function corresponds to the derivative of this function with respect to all
components of the estimate vector v→.
Because the discrete Fourier transformation is a unitary operation, it can be
evaluated conveniently by calculating the gradient of the TV term in the image
domain (i.e., estimating a vector that describes how the TV value changes for
modifications of the individual pixels), followed by an inverse Fourier
transformation to the frequency domain.


### 2.1. Extended TV Formulation

Calculation of the TV value according to2 uses only
the first-order derivative of the image with respect to its x- and
y-directions. This value is minimized if an image consists of areas with
constant signal intensity, so that the extrapolation procedure yields a
solution primarily with constant areas. While desirable for truly flat objects
like numerical phantoms, it tends to create images with a slightly blocky or
patchy appearance for real-world objects. Therefore, it is advisable to
additionally include second-order derivatives into the TV term, which then
allows for intensity gradients in the images and yields more naturally looking
solutions(5)TV2(I)=∑y=0N∑x=0Nσ⋅(|I(x,y)−I(x−1,y)|+|I(x,y)−I(x,y−1)|)⁢+(1−σ)⋅(|I(x−1,y)−2⋅I(x,y)+I(x+1,y)|+|I(x,y−1)−2⋅I(x,y)+I(x,y+1)|+|I(x,y)−I(x−1,y)−I(x,y−1)+I(x−1,y−1)|).Here, σ∈[01] is a weighting factor which can be used to
tune the images between a slightly more blocky looking and a slightly smoother
appearance. For the reconstructions presented, it was set to σ=0.77 based on the considerations by Geman and Yang  [15].


### 2.2. Edge-Preserving Denoising

In practice, experimental MRI data can be
significantly contaminated by Gaussian noise. While the aforementioned approach
is still able to reduce visible truncation artifacts under these circumstances,
it does not reduce image noise because the measuredk-space data remains
unchanged. On the other hand, an additional denoising may be achieved by
loosing the fixed bound on the measured data, that is by introducing a data
fitting term. In this case, the algorithm not only adds synthetic frequencies
to obtain a TV minimization, but is also allowed to find a solution that
slightly diverges from the measured data, which yields an effective
edge-preserving denoising. Therefore, the estimate vector v→ has to be extended such that it contains both
synthesized frequencies as well as frequencies from the measured part of
k-space, which is indicated by writing v→d instead.In the denoising case, the cost function takes the
form(6)Φ(v→d)=λ⋅∥v→d⊖y→∥22+TV(ℱ{v→d}),where ⊖ denotes an operation that calculates the
residual between the measured values and the corresponding entries of the
estimate, which are now contained in the vector v→d.
Further, λ is a weighting factor that allows to select
the desired denoising strength. While a low weight permits considerable
divergences from the measured values and, thus, leads to an effective removal
of noise, it can also cause a loss of object detail if selected too low.
Therefore, the weight has to be adjusted with respect to the signal-to-noise
ratio of the measurement sequence, where a reasonable strategy is to estimate a
fixed value once for each protocol by computing a set of test images with
different λ values and selecting the value yielding the
desired degree of denoising.


### 2.3. Phase Variations

Although the basic physical quantity measured by MRI,
that is, the spin-density modulated by relaxation or saturation effects, should
be real-valued and nonnegative in theory, inherent experimental phase
variations usually cause the observed object to be complex-valued. Moreover,
modern MRI systems often use multiple receive coils with complex-valued
sensitivity profiles, yielding differently modulated views of the object. As a
consequence, spatially varying transitions between the real and imaginary
component occur as well as local intensity changes, which conflict with the
assumption of a piecewise-constant quantity and prevent a direct application of
the TV constraint. Therefore, some mechanism is required to cope with the phase
variations and the multicoil scenario.In this proof-of-principle study, phase variations
were removed in a preprocessing step by performing a Fourier transformation of
the data from each coil and calculating the sum-of-squares of all channels in
the image domain. Subsequently, an inverse Fourier transformation of the sum-of-squares
data was performed to obtain a combined data set with real-valued and
nonnegative values in the image domain, which enables a calculation of the TV
value using only the real part of the image. While this simple technique turned
out to be sufficient for demonstrating a removal of truncation artifacts by
TV-constrained data extrapolation, routine applications will probably require a
more sophisticated procedure, in particular when combined with advanced
techniques such as parallel imaging and when using complex coil configurations
with more localized sensitivities of the individual receiver elements.


## 2.1. Extended TV Formulation

Calculation of the TV value according to2 uses only
the first-order derivative of the image with respect to its x- and
y-directions. This value is minimized if an image consists of areas with
constant signal intensity, so that the extrapolation procedure yields a
solution primarily with constant areas. While desirable for truly flat objects
like numerical phantoms, it tends to create images with a slightly blocky or
patchy appearance for real-world objects. Therefore, it is advisable to
additionally include second-order derivatives into the TV term, which then
allows for intensity gradients in the images and yields more naturally looking
solutions(5)TV2(I)=∑y=0N∑x=0Nσ⋅(|I(x,y)−I(x−1,y)|+|I(x,y)−I(x,y−1)|)⁢+(1−σ)⋅(|I(x−1,y)−2⋅I(x,y)+I(x+1,y)|+|I(x,y−1)−2⋅I(x,y)+I(x,y+1)|+|I(x,y)−I(x−1,y)−I(x,y−1)+I(x−1,y−1)|).Here, σ∈[01] is a weighting factor which can be used to
tune the images between a slightly more blocky looking and a slightly smoother
appearance. For the reconstructions presented, it was set to σ=0.77 based on the considerations by Geman and Yang  [15].


## 2.2. Edge-Preserving Denoising

In practice, experimental MRI data can be
significantly contaminated by Gaussian noise. While the aforementioned approach
is still able to reduce visible truncation artifacts under these circumstances,
it does not reduce image noise because the measuredk-space data remains
unchanged. On the other hand, an additional denoising may be achieved by
loosing the fixed bound on the measured data, that is by introducing a data
fitting term. In this case, the algorithm not only adds synthetic frequencies
to obtain a TV minimization, but is also allowed to find a solution that
slightly diverges from the measured data, which yields an effective
edge-preserving denoising. Therefore, the estimate vector v→ has to be extended such that it contains both
synthesized frequencies as well as frequencies from the measured part of
k-space, which is indicated by writing v→d instead.In the denoising case, the cost function takes the
form(6)Φ(v→d)=λ⋅∥v→d⊖y→∥22+TV(ℱ{v→d}),where ⊖ denotes an operation that calculates the
residual between the measured values and the corresponding entries of the
estimate, which are now contained in the vector v→d.
Further, λ is a weighting factor that allows to select
the desired denoising strength. While a low weight permits considerable
divergences from the measured values and, thus, leads to an effective removal
of noise, it can also cause a loss of object detail if selected too low.
Therefore, the weight has to be adjusted with respect to the signal-to-noise
ratio of the measurement sequence, where a reasonable strategy is to estimate a
fixed value once for each protocol by computing a set of test images with
different λ values and selecting the value yielding the
desired degree of denoising.


## 2.3. Phase Variations

Although the basic physical quantity measured by MRI,
that is, the spin-density modulated by relaxation or saturation effects, should
be real-valued and nonnegative in theory, inherent experimental phase
variations usually cause the observed object to be complex-valued. Moreover,
modern MRI systems often use multiple receive coils with complex-valued
sensitivity profiles, yielding differently modulated views of the object. As a
consequence, spatially varying transitions between the real and imaginary
component occur as well as local intensity changes, which conflict with the
assumption of a piecewise-constant quantity and prevent a direct application of
the TV constraint. Therefore, some mechanism is required to cope with the phase
variations and the multicoil scenario.In this proof-of-principle study, phase variations
were removed in a preprocessing step by performing a Fourier transformation of
the data from each coil and calculating the sum-of-squares of all channels in
the image domain. Subsequently, an inverse Fourier transformation of the sum-of-squares
data was performed to obtain a combined data set with real-valued and
nonnegative values in the image domain, which enables a calculation of the TV
value using only the real part of the image. While this simple technique turned
out to be sufficient for demonstrating a removal of truncation artifacts by
TV-constrained data extrapolation, routine applications will probably require a
more sophisticated procedure, in particular when combined with advanced
techniques such as parallel imaging and when using complex coil configurations
with more localized sensitivities of the individual receiver elements.


## 3. Methods

Simulations were performed with the Shepp-Logan
phantom, which is composed of several ellipses. Because the Fourier transform
of a single ellipse is given by a Bessel function, an analytical Fourier
transform of the phantom is obtained by a superposition of respective Bessel
functions. Truncation artifacts can be studied by evaluating the noncompact
analytical transform at the sampling positions along the trajectory, here
yielding a matrix of96×96 Fourier samples. All simulations and
processing of experimental data were done offline using an in-house software
package written in C/C++.MRI experiments were conducted at 2.9 T (Siemens
Magnetom TIM Trio, Erlangen, Germany) with use of a receive only 12-channel
head coil equipped with hardware signal combiners, yielding four receiver
channels with different combinations of the coil elements. Measurements were
performed for a water phantom and human brain in vivo, where written informed
consent was obtained from all subjects prior to each examination. For
demonstration purposes, the image acquisitions were done with a simple
slice-selective spin-echo sequence at a200×200 mm2 field of view, covered by a 96×96 acquisition matrix. Different sequence
settings were used to obtain data sets with low and high level of noise, where
the latter was achieved by reducing the flip angle and slice thickness while
increasing the receiver bandwidth. Further, one data set was acquired with a
slice-selective gradient-echo sequence, which allowed for the rapid measurement
of a full 288×288 acquisition matrix.All images were reconstructed on a288×288 matrix corresponding to an extrapolation
factor of 3. The proposed algorithm was run for a fixed number of 3000
iterations, which takes about 2-3 minutes on a standard microprocessor. In
cases where an additional data fitting term was used, the weighting factor λ was adjusted manually to yield a reasonable
solution as judged by visual inspection. Zero-padded solutions with and without
filtering were calculated for comparison. Here, a simple Lanczos sigma filter,
that is, multiplication with a sinc-function, was applied, where the window
width was selected such that the sinc-function's first null coincides with the
border of the measured k-space. Although other filters might perform better, it
serves to demonstrate the general problem related to data filtering. In
addition, a two-dimensional version of the extrapolation method described by
Constable and Henkelman [6] was implemented with a window width of P=2 for the edge-preserving sigma filter. In our
implementation, the filter parameter Δ was selected according to Δ=c⁢⋅I,
where I denotes the intensity of the pixel to be
filtered and c is a global coefficient that was set to c=0.1 based on visual inspection. Finally, all
images were magnified and cropped to improve the visibility of the artifacts.


## 4. Results

Figure2 shows different reconstructions of the
Shepp-Logan phantom (left column) together with the respective Fourier
transforms (right column). It is clearly visible that the zero-padded solution
(zero) suffers from severe ringing artifacts around all edges of the phantom.
The extent of the measured k-space can be seen in its Fourier transform. Most
ringing artifacts disappear after filtering (filter), however, at the expense
of a significant loss of image resolution. In contrast, the image reconstructed
with the proposed method (TV) is neither affected by ringing artifacts nor by
blurring, and it presents with considerably
sharper edges relative to the zero-padded solution. Its Fourier transform
reveals that the measured data has been properly extrapolated into the
uncovered areas of k-space. For comparison, a full reconstruction from 288×288 samples is shown in the bottom row (full).(Left) Images of the numerical Shepp-Logan
phantom (96×96 samples, 288×288 reconstruction matrix) and (right)
corresponding k-space representations reconstructed using zero-padding (zero),
filtered zero-padding (filter), and the proposed extrapolation method (TV). For
comparison, a data set with a fully sampled 288×288 matrix is shown (full). Arrow  =  truncation
artifact.Figure3 demonstrates the application of the method to
experimental data obtained for a phantom (left column) and a human brain in
vivo (right column) in comparison to zero-padded (zero) and filtered
zero-padded solutions (filter). Again, the ringing artifacts obtained for zero
padding (indicated by arrows) are significantly reduced when using
TV-constrained data extrapolation with only first-order (TV) or additionally
second-order derivatives (TV2). The blocky appearance of the TV reconstruction
becomes much more smoother for the TV2 approach, although both solutions (TV
and TV2) look somewhat more blocky than the zero-padding solution.Spin-echo images (96×96 samples, 288×288 reconstruction matrix) of (left) a phantom
(TR/TE  =  4000/8 ms, BW 243 Hz/pixel, FA 70∘,
3 mm slice) and (right) a human brain in vivo (TR/TE  =  4000/25 ms, BW 180 Hz/pixel, FA 90∘,
2 mm slice) using zero-padding (zero), filtered zero-padding (filter), the
proposed extrapolation method with first-order (TV), and additionally
second-order derivatives (TV2). Arrows  =  truncation artifacts.In Figure4, the proposed approach (TV2) is compared
to a 2D version of the extrapolation method by Constable and Henkelman (comp)
for the Shepp-Logan phantom (left column) and an experimental study of the
human brain in vivo (right column). It can be seen that the performance of the
proposed extrapolation method is slightly better for the simulated data, while
both approaches yield an effective suppression of ringing artifacts for the
experimental data (note that the alternative method leads to a slight denoising
of the image). However, in our implementation, the algorithm by Constable and
Henkelman tends to be sensitive to the parameter selection for the initial
filter that is used to detect true edges of the object, whereas a selection of
respective parameters is not needed in the TV-based approach.(Left) Images of the Shepp-Logan
phantom (96×96 samples, 288×288 reconstruction matrix) and (right)
gradient-echo images of the human brain in vivo (TR/TE  =  500/8 ms, BW 80 Hz/pixel, FA 30∘,
2 mm slice) reconstructed using zero-padding (zero), the extrapolation method
by Constable and Henkelman (comp), and the proposed method with first- and
second-order derivatives (TV2). For comparison, reconstructions from fully
sampled 288×288 matrices are shown (full). Arrows  =  truncation
artifacts.Figure5 shows reconstructions of the Shepp-Logan
phantom from noisy data using zero-padding (zero), the proposed extrapolation
approach (TV), and its combination with denoising (TVdns). While the basic
extrapolation approach leads to a reduction of truncation artifacts also for
noisy data, it does not reduce the noise patterns. However, when extending the
TV penalty to the measured data, the method effectively flattens noise patterns
in addition to the suppression of ringing artifacts.Figure 5
Images of the numerical Shepp-Logan phantom
reconstructed from noisy data (96×96 samples, 288×288 reconstruction matrix) using zero-padding
(zero), the proposed extrapolation method (TV), and the proposed method
combined with denoising (TVdns). Arrow  =  truncation artifact.Finally, corresponding reconstructions from
experimental data with a high degree of noise are shown in Figure6. Here, a
combination of first- and second-order derivatives was used for the TV
calculation. As in the simulations, the proposed method leads to a reduction of
truncation artifacts (TV2), while the extension to data fitting yields an
additional edge-preserving denoising (TV2dns).Spin-echo images of (top) a phantom (TR/TE  =  4000/100 ms, BW 789 Hz/pixel, FA50∘,
1 mm slice) and (bottom) a human brain in vivo (TR/TE  =  4000/15 ms, BW 401 Hz/pixel, FA 70∘,
1 mm slice) reconstructed from noisy data (96×96 samples, 288×288 reconstruction matrix) using zero-padding
(zero), the proposed extrapolation method (TV2), and the proposed method
combined with denoising (TV2dns). Arrows  =  truncation artifacts.


## 5. Discussion

### 5.1. Accuracy and Limitations

Both simulations and experiments demonstrate that
TV-constrained data extrapolation effectively reduces truncation artifacts due
to finitely sampled MRI acquisitions. Usually, the images exhibit a more blocky
appearance compared to zero-padding. However, it should be noted that the
smoothness observed for zero-padding originates to a significant degree from
the convolution with the sinc-function. As a consequence, a sharp edge of the
object is mapped as a rather smooth pattern, which might appear more familiar
to the viewer than a blocky image, but strictly represents an image artifact.
Hence, the extrapolation technique may even lead to a slight gain of resolution
due to a sharpening of the point-spread function, following from the
reciprocity property of the Fourier transformation. This effect can be best
seen in Figure6 when comparing the borders of the dark brain vessels obtained
for zero padding (zero) with the proposed method (TV2).Residual image artifacts are explained by multiple
reasons. First, the method is based on the assumption that the true object is
piecewise constant, which is only approximately valid for real-world objects.
In the presence of additional experimental effects like flow artifacts, the
assumption might be even less appropriate. Hence, the extrapolation performance
depends on the object's conformance with the assumption that it is
piecewise-constant. Moreover, if the true object contains strongly varying
patterns, the algorithm may erroneously soften such patterns by supplementing
respective high frequencies. On the other hand, in the majority of cases, the
assumption of a piecewise-constant object seems to be more appropriate than
that of all conventional reconstructions, namely, a Fourier transform of the
object that is zero outside the sampledk-space area.Second, the proposed method synthesizes only a finite
number of additional frequencies, whereas an infinite number ofk-space samples
would be required to completely eliminate all truncation effects. In practice,
however, it turned out that there is no perceivable benefit of extrapolating by
a factor of higher than three. The reason is that the method yields an implicit
filtering of the extrapolated data: assuming that the extrapolation procedure
would recover the unmeasured k-space samples exactly, then a new truncation
effect would arise at the extended border and again lead to ringing artifacts
in image space (though with a higher oscillation frequency). Because this would
increment the TV value, the method automatically lowers outer frequencies during
the extrapolation procedure to prevent the upcoming of new ringing artifacts.
Hence, the extrapolated values diverge categorically from the true frequencies
which, in this case, is a rather desirable feature as the prime target is to
reduce visually annoying ringing artifacts rather than to gain
super-resolution.Third, if a completely artifact-free reconstruction of
the object would be available, then respective frequency samples could be
calculated with a discrete Fourier transformation of the given image.
Interestingly, these samples would diverge from the experimentally measured
frequencies, because image pixels are discrete and, thus, the Fourier transform
of the image is periodic such that outer frequencies from neighboring copies
(of the true object's noncompact Fourier transform) overlap. This is different
from the experimental situation where the object is continuous and the outer
frequencies are missing instead of overlapping. Consequently, an artifact-free
discrete reconstruction can only be obtained if the samples used for the
reconstruction specifically diverge from the measured frequencies. A complete
artifact removal, therefore, requires to alter the measured frequencies instead
of keeping them unchanged. Unfortunately, the
information  how the samples have to
be adjusted is not available, so that in practice a data fitting term might be
the best solution when a complete removal of ringing artifacts is needed.
However, this might cause a loss of object detail as described in the theory
section.


### 5.2. Implementation Issues

The modulus function in the TV formula2 has a
fundamental role for the success of TV-based image processing. Because
subtraction of neighboring pixels—performed before taking the modulus—can
be seen as applying a difference operator to the estimate, TV minimization
yields a solution with minimum L1-norm in the difference basis. Due to the
specific character of the modulus function, this solution tends to be sparse in
the difference basis: it has few large jumps and most differences between
neighboring pixels are near zero, which directly translates into a
piecewise-constant image (and explains the edge preserving character of
TV-based denoising). If the modulus would be replaced by a square function,
then the optimizer would try to find a minimum L2-norm solution with minimal
jumps between all neighboring pixels. This corresponds to a globally smooth
image, which is usually undesired due to a loss of sharp edges. While it is
rather simple to obtain a minimum L2-norm solution as its cost function is
strictly convex, finding a minimum L1-norm solution is much more challenging;
and many optimization algorithms fail if directly applied to the TV problem.
One major reason is that the derivative of the modulus function is just ±1,
which does not help to guess a reasonable step size toward the function's
minimum. However, it turned out that the CG-Descent algorithm is capable to
handle the problem as it comprises a powerful line-search procedure, but it is
probably not the optimal method for finding the solution. In particular, the
convergence tends to be somewhat sensitive to the scaling of the data. In order
to ensure convergence, it was, therefore, necessary to introduce a scaling
factor that limits the modification strength for each iteration and to run the
algorithm in turn for a high number of iterations (e.g., 3000 iterations as
arbitrarily chosen here). Nevertheless, this issue should not be seen as a
drawback of the proposed extrapolation approach itself, but rather as a
technical aspect of the optimization method utilized in this proof-of-principle
study. Employing a dedicated algorithm for TV minimization should render a
scaling factor unnecessary and significantly improve the convergence rate.
Although such enterprise promises reconstructions in a fraction of the current
processing time, it is outside the scope of the present study.


## 5.1. Accuracy and Limitations

Both simulations and experiments demonstrate that
TV-constrained data extrapolation effectively reduces truncation artifacts due
to finitely sampled MRI acquisitions. Usually, the images exhibit a more blocky
appearance compared to zero-padding. However, it should be noted that the
smoothness observed for zero-padding originates to a significant degree from
the convolution with the sinc-function. As a consequence, a sharp edge of the
object is mapped as a rather smooth pattern, which might appear more familiar
to the viewer than a blocky image, but strictly represents an image artifact.
Hence, the extrapolation technique may even lead to a slight gain of resolution
due to a sharpening of the point-spread function, following from the
reciprocity property of the Fourier transformation. This effect can be best
seen in Figure6 when comparing the borders of the dark brain vessels obtained
for zero padding (zero) with the proposed method (TV2).Residual image artifacts are explained by multiple
reasons. First, the method is based on the assumption that the true object is
piecewise constant, which is only approximately valid for real-world objects.
In the presence of additional experimental effects like flow artifacts, the
assumption might be even less appropriate. Hence, the extrapolation performance
depends on the object's conformance with the assumption that it is
piecewise-constant. Moreover, if the true object contains strongly varying
patterns, the algorithm may erroneously soften such patterns by supplementing
respective high frequencies. On the other hand, in the majority of cases, the
assumption of a piecewise-constant object seems to be more appropriate than
that of all conventional reconstructions, namely, a Fourier transform of the
object that is zero outside the sampledk-space area.Second, the proposed method synthesizes only a finite
number of additional frequencies, whereas an infinite number ofk-space samples
would be required to completely eliminate all truncation effects. In practice,
however, it turned out that there is no perceivable benefit of extrapolating by
a factor of higher than three. The reason is that the method yields an implicit
filtering of the extrapolated data: assuming that the extrapolation procedure
would recover the unmeasured k-space samples exactly, then a new truncation
effect would arise at the extended border and again lead to ringing artifacts
in image space (though with a higher oscillation frequency). Because this would
increment the TV value, the method automatically lowers outer frequencies during
the extrapolation procedure to prevent the upcoming of new ringing artifacts.
Hence, the extrapolated values diverge categorically from the true frequencies
which, in this case, is a rather desirable feature as the prime target is to
reduce visually annoying ringing artifacts rather than to gain
super-resolution.Third, if a completely artifact-free reconstruction of
the object would be available, then respective frequency samples could be
calculated with a discrete Fourier transformation of the given image.
Interestingly, these samples would diverge from the experimentally measured
frequencies, because image pixels are discrete and, thus, the Fourier transform
of the image is periodic such that outer frequencies from neighboring copies
(of the true object's noncompact Fourier transform) overlap. This is different
from the experimental situation where the object is continuous and the outer
frequencies are missing instead of overlapping. Consequently, an artifact-free
discrete reconstruction can only be obtained if the samples used for the
reconstruction specifically diverge from the measured frequencies. A complete
artifact removal, therefore, requires to alter the measured frequencies instead
of keeping them unchanged. Unfortunately, the
information  how the samples have to
be adjusted is not available, so that in practice a data fitting term might be
the best solution when a complete removal of ringing artifacts is needed.
However, this might cause a loss of object detail as described in the theory
section.


## 5.2. Implementation Issues

The modulus function in the TV formula2 has a
fundamental role for the success of TV-based image processing. Because
subtraction of neighboring pixels—performed before taking the modulus—can
be seen as applying a difference operator to the estimate, TV minimization
yields a solution with minimum L1-norm in the difference basis. Due to the
specific character of the modulus function, this solution tends to be sparse in
the difference basis: it has few large jumps and most differences between
neighboring pixels are near zero, which directly translates into a
piecewise-constant image (and explains the edge preserving character of
TV-based denoising). If the modulus would be replaced by a square function,
then the optimizer would try to find a minimum L2-norm solution with minimal
jumps between all neighboring pixels. This corresponds to a globally smooth
image, which is usually undesired due to a loss of sharp edges. While it is
rather simple to obtain a minimum L2-norm solution as its cost function is
strictly convex, finding a minimum L1-norm solution is much more challenging;
and many optimization algorithms fail if directly applied to the TV problem.
One major reason is that the derivative of the modulus function is just ±1,
which does not help to guess a reasonable step size toward the function's
minimum. However, it turned out that the CG-Descent algorithm is capable to
handle the problem as it comprises a powerful line-search procedure, but it is
probably not the optimal method for finding the solution. In particular, the
convergence tends to be somewhat sensitive to the scaling of the data. In order
to ensure convergence, it was, therefore, necessary to introduce a scaling
factor that limits the modification strength for each iteration and to run the
algorithm in turn for a high number of iterations (e.g., 3000 iterations as
arbitrarily chosen here). Nevertheless, this issue should not be seen as a
drawback of the proposed extrapolation approach itself, but rather as a
technical aspect of the optimization method utilized in this proof-of-principle
study. Employing a dedicated algorithm for TV minimization should render a
scaling factor unnecessary and significantly improve the convergence rate.
Although such enterprise promises reconstructions in a fraction of the current
processing time, it is outside the scope of the present study.


## 6. Conclusion

The present work demonstrates that the simple
assumption of a piecewise-constant object can be exploited to extrapolate
measured data ink-space. This allows for a significant reduction of ringing
artifacts that arise from data truncation in k-space. In contrast to commonly
used filtering approaches, the method does not degrade the spatial resolution
of the reconstructed image and rather leads to a mild resolution enhancement
due to sharpening of the point-spread function. If the measured data is
seriously contaminated by noise, an extended approach offers edge-preserving
denoising by slightly altering the measured data in addition to supplementing
synthetic data. Both variants can be implemented as a pure postprocessing
procedure and are also applicable for partial Fourier acquisitions. Therefore,
no modification of the MRI sequence is required. While the current
implementation suffers from a relatively high computational load, the use of a
dedicated TV optimization algorithm promises a processing speed suitable for
routine applications.



---

*Source: 184123-2008-09-07.xml*