Variable density sampling based on physically plausible gradient waveform. Application to 3D MRI angiography.

Variable density sampling based on physically plausible gradient waveform. Application to 3D MRI angiography.


Performing -space variable density sampling is a popular way of reducing scanning time in Magnetic Resonance Imaging (MRI). Unfortunately, given a sampling trajectory, it is not clear how to traverse it using gradient waveforms. In this paper, we actually show that existing methods [1, 2] can yield large traversal time if the trajectory contains high curvature areas. Therefore, we consider here a new method for gradient waveform design which is based on the projection of unrealistic initial trajectory onto the set of hardware constraints. Next, we show on realistic simulations that this algorithm allows implementing variable density trajectories resulting from the piecewise linear solution of the Travelling Salesman Problem in a reasonable time. Finally, we demonstrate the application of this approach to 2D MRI reconstruction and 3D angiography in the mouse brain.

Variable density sampling based on physically plausible gradient waveform. Application to 3D MRI angiography.

Nicolas Chauffert, Pierre Weiss, Marianne Boucher, Sébastien Mériaux and Philippe Ciuciu
CEA/DSV/IBM NeuroSpin center, Bât. 145, F-91191 Gif-sur-Yvette, France
INRIA Saclay Ile-de-France, Parietal team, 91893 Orsay, France.
PRIMO Team, ITAV, USR 3505, Université de Toulouse.

Index Terms—  MRI, Compressive sensing, Variable density sampling, gradient waveform design, hardware constraints, angiography.

1 Introduction

Compressed Sensing (CS) provides a theoretical framework to justify the downsampling of -space (2D or 3D Fourier domain) in Magnetic Resonance Imaging (MRI). CS-MRI is usually based on independent random drawing of -space locations according to a prescribed density. From recent theoretical works [3, 4], one can derive an optimal sampling density that reduces at most the number of samples collected in MRI without degrading the image quality at the reconstruction step [5, 6]. In [7], simulations show that distributions with radial decay (see Fig. 1(a)) with full -space center acquisition perform better in numerical experiments.

However, such sampling schemes are not performed along continuous lines and thus not physically plausible in MRI because of the constraints involved on the magnetic field gradient (magnitude and slew-rate). In [8], we have proposed a new approach to design continuous sampling trajectories based on the solution of Travelling Salesman Problem (TSP), as illustrated in Fig. 1(b). The specificity of this approach is that the empirical distribution of the trajectory can approximate any prescribed distribution . Such a curve is called a -Variable Density Sampler (-VDS). Unfortunately, continuity of the sampling trajectory is not a sufficient condition in MRI and it is not clear how to design admissible gradient waveforms to traverse such a trajectory.

(a) (b)
Fig. 1: Example of 2D Variable Density Sampler. (a):  as advocated in [7]. (b): TSP-based sampling trajectory.

To the best of our knowledge, the most efficient gradient design strategies consist of finding an admissible parameterization of a given curve by using optimal control [1, 2]. However, the main drawback of this approach is that the traversal of high curvature parts of the curve is rather slow. In the extreme case of angular points such as in TSP-based trajectories, it leads to extremely large acquisition time. In [9], we have proposed a new method to design magnetic field gradients by projecting any parameterized curve onto the set of hardware constraints. This method allows one to change the curve support and thus to yield faster traversal time.

The goal of this paper is then to prove that this projection algorithm enables to implement a TSP-based VDS on MRI scanners in reasonable scanning times, in contrast to optimal reparamerization. The first part is dedicated to summarizing the projection strategy introduced in [9]. Next, we provide in Section 3 a proof-of-concept on retrospective CS simulations. In particular, we show that our algorithm yields faster sampling trajectories than the state-of-the-art for a given image reconstruction quality, and alternatively that if the scanning time is fixed, our method delivers improved reconstructions. Finally, our strategy (TSP-based sampling + projection algorithm) is applied to 3D angiography of the mouse brain before and after contrast agent injection, to demonstrate its efficiency in terms of scanning time reduction while preserving the recovery of small structures such as blood vessels.

2 Gradient waveform design using a projection algorithm

In this section, we recall the hardware constraints as generally modelled in MRI and describe methods for designing gradient waveforms in order to traverse -space sampling curves.

The gradient waveform associated with a curve is defined by , where denotes the gyro-magnetic ratio [10]. The gradient waveform is obtained by energizing gradient coils (arrangements of wire) with electric currents. Owing to obvious physical constraints, these electric currents have a bounded amplitude and cannot vary too rapidly (slew rate). Mathematically, these constraints read:

A sampling trajectory will be said admissible if it belongs to the set:

In this paper, we limit ourselves to the so-called rotation-invariant constraints. Some hardware systems, where the coils are energized independently, enable considering rotation-variant constraints. In this case, the norm is replaced by norm. The differences are discussed in [2], but the two compared methods (optimal control and projection algorithm) are able to handle both kinds of constraints.

2.1 State-of-the-art

The question of finding jointly an accurate trajectory and the admissible gradient waveform to traverse it is a difficult issue that has received special attention in [11, 12]. The most classical approaches consist of fixing a curve support [1, 2], or control points [10, 13], and finding an admissible parameterization afterwards.

In particular, the most popular approach to design an admissible curve assumes the knowledge of a parameterized curve and consists of finding its optimal reparameterization by using optimal control [1, 2]. In other words, it amounts to finding a reparameterization such that satisfies the physical constraints while minimizing the acquisition time. This problem can be cast as follows:


The resulting solution has the same support as , which might be an important feature in some applications. The problems of this approach are: i) there is no control of the sampling density, especially in the high curvature parts of where samples tend to agglutinate and ii) there is no control over the total sampling time , which can be large if the trajectory contains singular points for instance (e.g., see Fig. 1(b)). These two drawbacks are illustrated in [9].

The next part is dedicated to introducing an alternative method relaxing the constraint of the curve support.

2.2 Projection onto the set of constraints

The idea introduced in [9] is to find the projection of the input curve onto the set of admissible curves :


where . For the sake of conciseness, the theoretical grounds for the resolution and the key properties discussed below are given in [9].

Resolution. Problem 2 consists of minimizing a convex smooth function over a convex set. In [9], we have proposed a fast iterative algorithm exploiting the structure of the dual formulation of (2), which can be solved using proximal gradient methods [16].
Key properties. The two main advantages of this method are that i) the sampling density of is close to the sampling density of . In particular, if is a -VDS, the sampling density of is close to 111the closeness is quantified by Wasserstein transportation distance , see [9] for details.; and ii) the acquisition time is fixed and equal to that of the input curve . In particular, the time to traverse a curve is in general shorter than with optimal reparameterization ().

In the next part, we will emphasize that our algorithm enables to traverse VDS curves as depicted in Fig. 1(b) in a reasonable time, unlike optimal control-based reparemeterizations.

3 CS-MRI simulations

In this part, we compare the time to traverse -space along different trajectories using gradients computed either by the standard optimal control approach or by our proposed projection algorithm. For comparison between sampling schemes, we work on retrospective CS, meaning that a full dataset has been acquired, and then a posteriori downsampling is performed. We compare the reconstruction results in terms of peak signal-to-noise ratio (PSNR) with respect to the acquisition time and to the “acceleration factor”222 quantifies the reduction of the number of measurements . If the -space is a grid of pixels is commonly used in CS-MRI. .

3.1 Experimental framework

Data acquisition. The initial experimental setup aimed at observing blood vessels of living mice using an intraveinous injection of an iron oxide-based contrast agent (Magnetovibrio Blakemorei MV1). Because of natural elimination, it is necessary to speed up acquisition to improve contrast and make easier post-processing such as angiography. The experiments have been performed on a 17.2T preclinical scanner which physical rotation-invariant constraints are, for all :

A FLASH sequence (Fast Low Angle SHot) has been used to reveal the contrast induced by the injection of the contrast agent (TE/TR = 8/680 ms). The sequence was repeated 12 times to improve the signal-to-noise ratio (SNR), leading to a total acquistion time of 30 minutes to acquire the -space slice by slice. The spatial resolution achieved is .

Hypothesis. The aim of this paper is to prove that one can expect a large acquisition time reduction using partial -space measurements. The time to traverse a sampling curve is computed satisfying the gradient constraints. To achieve a fair comparison, let us mention the additional hypothesis that our acquisitions are single-shot, meaning that the partial -space is acquired after a single RF pulse. We did not take echo and repetition times into account to ensure the recovery of a -weigthed image. We only compare the time to traverse a curve using the gradients with their maximal intensity. We assume that there is no error on the -space sample locations. In practice we have to measure the three magnetic field gradients that are actually played out by the scanner to correct the trajectory and avoid distortions. We shall work on a discrete cartesian -space, and consider that a sample is measured if the sampling trajectory crosses the corresponding cell of the -space grid. Using this hypothesis, the estimated time to visit the 2D -space is 110 ms.

Strategy. We used the TPS-based sampling method [7] as input of our projection algorithm (see Fig. 2(b)), since it is a way of designing sampling trajectories that match any sampling density . The latter is central in CS-MRI since it impacts the number of required measurements [14, 15, 7]. To compare our projection method to existing reparameterization, the proposed sampling strategy is:
(i) Sample deterministically the -space center as adviced in [14, 5, 7], using an EPI sequence (see Fig. 2(a)). The scanning time can be estimated to 12 ms in 2D using optimal control.
(ii) Select a density proportional to as mentioned in [15, 7]. Draw independently points according to and join them by the shortest path to form a -VDS [7].
(iii) Parameterize the TSP path at constant speed and project this parameterization onto the set of gradient constraints, or (iii bis) Parameterize the TSP path using optimal control (the exact solution can be computed explicitely).
(iv) Form the sampling curve, define a set of the selected samples, mask the -space with , and reconstruct an image using minimization of the constrained problem. Let denote the -dimensional discrete Fourier transform and the matrix composed of the lines corresponding to . Denote also by an inverse -dimensional wavelet transform (here a Symmlet transform). Then the reconstructed image is the solution of the problem:


An approximation of is computed using Douglas-Rachford algorithm [16]. Solving the penalized form associated with (3) might be addressed by competing algorithms (ADMM, 3MG); see [17] for a recent comparison. The reconstruction results could be improved by resorting to non-Cartesian reconstruction [18], which would avoid the approximation related to the projection onto the -space grid.

3.2 Results

3.2.1 2D reconstructions

In this experiment, we considered a 2D -space () corresponding to an axial slice. We considered five sampling strategies, depicted in Fig. 2(first row): a classical EPI coverage used as reference (a); a TSP-based sampling trajectory parameterized using optimal control (b); two projected TSP-based trajectories, one with the same number of samples collected as in (b) () (c) and the other with the same scanning time as in (b) (62 ms) (d); a variable density spiral trajectory for comparison purpose in terms of time and sampling ratio (e).

As expected, the reconstruction results shown in Fig. 2(g,h) are really close, since the number of collected samples is the same, and the sampling densities are similar. However, in this comparison the gain in traversal time is significant (one half). In contrast, the longer and smoothed TSP depicted in Fig. 2(d) allows us to improve image reconstruction (1 dB gain) as illustrated by Fig. 2(i) while keeping the same acquisition time as in Fig. 2(b). For comparison purposes, we implemented spiral acquisition which consists of replacing steps (ii)-(iii) in the above mentioned sampling strategy by a spiral with density proportional to , projected onto the set of constraints. This strategy doubles the acquisition time (118 ms compared to 62 ms) whereas the acceleration factor was larger ( vs. ). In this experimental context (regridding and variable density spiral), the spiral is not appealing compared to EPI acquisition, since it is time consuming and degrades the image quality.

In each of these reconstructions, the major vessels can be recovered, although the smallest ones can only be seen for . Finally, the best compromise between acquisition time and reconstruction quality is achieved using the specific combination of TSP-based sampling and our projection algorithm onto the set of constraints shown in Fig. 2(d).

(a) (b) (c) (d) (e)

Sampling schemes

ms () ms () ms () ms () ms ()
(f) (g) (h) (i) (j)

Reconstructed slices

Reference PSNR = 25.9 dB PSNR = 25.5 dB PSNR = 26.9 dB PSNR = 26.8 dB
Fig. 2: Full -space acquisition with an EPI sequence (a) and corresponding reference image (f). Comparison between an exact parameterization of the TSP trajectory (b) and projection from TSP trajectory onto the set of constraints (c),(d). In experiments (b,c), the number of measured locations is fixed to 9% (, whereas in (b,d), the time to traverse the curve is fixed to 62 ms. (e): Spiral trajectory with full acquisition of the -space center. (g-j): Reconstructed images corresponding to sampling strategies (b-e) by solving Eq. (3).

3.2.2 3D angiography

Using the same method as in 2D, namely TSP-sampling and projection onto the set of constraints, we reconstructed volumes from 3D -space. In order to estimate the quality of the reconstructions, we compared the angiograms computed from the 3D images using Frangi filtering [19]. The results are shown in Fig. 3 for acceleration factors (Fig. 3(b,e)) and (Fig. 3(c,f)) and compared to the angiogram computed from the whole data.
Using the strategy described in Part 3.1 the time to traverse -space would be 3.53 s (full acquisition), 3.15 s () and 0.88 s (). The main drawback of TSP-based sampling schemes is that the time reduction is not directly proportional to , in contrast to classical 2D downsampling and reading out along the third dimension. Nevertheless, if the number of measurements is fixed, the TSP-based approach leads to more accurate reconstruction results since the sampling scheme may fit any density [7].
Angiograms shown in Fig. 3 illustrate that one can reduce the travel time in the -space and still observe accurate microvascular structure. If , time reduction is minor (about 10% less), but the computed angiogram is almost the same as the one obtained with a complete -space. It is interesting to notice that with a higher acceleration factor (), the acquisition time is reduced by 75%, but the computed angiogram remains of good quality. The angiogram appears a bit noisier, especially in the pre-injection setting (Fig. 3(c)), but the post-injection image allows recovering Willis polygon and most of the major vessels of the mouse brain (Fig. 3(f)).

Original (s) (s) (s)
(a) (b) (c)


PSNR=29.0 dB PSNR=26.6 dB
(d) (e) (f)


PSNR=25.5 dB PSNR=24.1 dB
Fig. 3: Angiograms computed from full -space pre-(a) and post-(d) injection data. Angiograms computed from pre-(resp., post-) injection data for decimated -space with  (b) and  (c) (resp., (e) and (f)).

4 Conclusion and Future Work

In this paper, we have shown that the projection algorithm introduced in [9] is a promising technique to design gradient waveforms. In particular, it allows us to design gradient waveforms in order to implement TSP-based VDS, while existing gradient design methods lead to extremely large acquisition time. We are currently implementing these waveforms on actual scanners to validate our method on a real CS-MRI framework.

5 Acknowledgments

This work benefited from the support of the ”FMJH Program Gaspard Monge in optimization and operation research”, and from EDF’s support to this program.


  • [1] M. Lustig, S. J. Kim, and J. M. Pauly, “A Fast Method for Designing Time-Optimal Gradient Waveforms for Arbitrary -Space Trajectories,” IEEE Trans. Med. Imag., vol. 27, no. 6, pp. 866–873, 2008.
  • [2] S. Vaziri and M. Lustig, “The fastest gradient waveforms for arbitrary and optimized k-space trajectories,” in Proc. of 10th IEEE ISBI conference, San Francisco, USA, Apr. 2013, pp. 708–711.
  • [3] E. J. Candès and Y. Plan, “A probabilistic and RIPless theory of compressed sensing,” IEEE Trans. Inf. Theory, vol. 57, no. 11, pp. 7235–7254, 2011.
  • [4] H. Rauhut, “Compressive Sensing and Structured Random Matrices,” in Theoretical Foundations and Numerical Methods for Sparse Recovery, M. Fornasier, Ed., vol. 9 of Radon Series Comp. Appl. Math., pp. 1–92. deGruyter, 2010.
  • [5] N. Chauffert, P. Ciuciu, and P. Weiss, “Variable density compressed sensing in MRI. Theoretical vs. heuristic sampling strategies,” in Proc. of 10th IEEE ISBI conference, USA, Apr. 2013, pp. 298–301.
  • [6] G. Puy, P. Vandergheynst, and Y. Wiaux, “On variable density compressive sampling,” IEEE Sig. Proc. Let., vol. 18, pp. 595–598, 2011.
  • [7] N. Chauffert, P. Ciuciu, J. Kahn, and P. Weiss, “Variable density sampling with continuous trajectories. Application to MRI,” SIAM Journal on Imaging Science, vol. 7, no. 4, pp. 1962–1992, 2014.
  • [8] N. Chauffert, P. Ciuciu, J. Kahn, and P. Weiss, “Travelling salesman-based compressive sampling,” in Proc. of 10th SampTA conference, Bremen, Germany, July 2013, pp. 509–511.
  • [9] N. Chauffert, P. Weiss, J. Kahn, and P. Ciuciu, “Gradient waveform design for variable density sampling in Magnetic Resonance Imaging,” arXiv preprint arXiv:1412.4621, 2014.
  • [10] B. A. Hargreaves, D. G. Nishimura, and S. M. Conolly, “Time-optimal multidimensional gradient waveform design for rapid imaging,” Magn. Reson. Med., vol. 51, no. 1, pp. 81–92, 2004.
  • [11] B. M. Dale, J. S. Lewin, and J. L. Duerk, “Optimal design of k-space trajectories using a multi-objective genetic algorithm,” Magn. Reson. Med., vol. 52, no. 4, pp. 831–841, 2004.
  • [12] R. Mir, A. Guesalaga, J. Spiniak, M. Guarini, and P. Irarrazaval, “Fast three-dimensional k-space trajectory design using missile guidance ideas,” Magn. Reson. Med., vol. 52, no. 2, pp. 329–336, 2004.
  • [13] M. Davids, M. Ruttorf, F. Zoellner, and L. Schad, “Fast and robust design of time-optimal -space trajectories in MRI,” Medical Imaging, IEEE Transactions on, 2014, in press.
  • [14] B. Adcock, A. Hansen, C. Poon, and B. Roman, “Breaking the coherence barrier: asymptotic incoherence and asymptotic sparsity in compressed sensing,” 2013.
  • [15] F. Krahmer and R. Ward, “Beyond incoherence: stable and robust sampling strategies for compressive imaging,” preprint, 2012.
  • [16] P. L. Combettes and J.-C Pesquet, “Proximal Splitting Methods in Signal Processing,” in Fixed-Point Algorithms for Inverse Problems in Science and Engineering, pp. 185–212. Springer, 2011.
  • [17] A. Florescu, E. Chouzenoux, J.-C. Pesquet, P. Ciuciu and S. Ciochina, “A majorize-minimize memory gradient method for complex-valued inverse problems,” in Signal Processing, vol. 103, 2014, pp. 285–295
  • [18] J. Keiner, S. Kunis, and D. Potts “Using NFFT 3 - a software library for carious nonequispaced fast Fourier transforms,” in ACM Transactions on Mathematical Software, 36(4):19 , 1998.
  • [19] A. F. Frangi, W. J. Niessen, K. L. Vincken, and M. A. Viergever, “Multiscale vessel enhancement filtering,” in MICCAI’98, pp. 130–137. Springer, 1998.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description