Nicolas Chauffert
email: nicolas.chauffert@cea.fr

Philippe Ciuciu
email: philippe.ciuciu@cea.fr

Jonas Kahn
email: jonas.kahn@math.univ-lille1.fr

Pierre Weiss
email: pierre.armand.weiss@gmail.com

VARIABLE DENSITY SAMPLING WITH CONTINUOUS TRAJECTORIES. APPLICATION TO MRI

Keywords: Variable density sampling, compressed sensing, CS-MRI, stochastic processes, empirical measure, TSP, Markov Chains, l 1 reconstruction AMS subject classifications. 94A20, 60G20, 15A52, 94A08

Reducing acquisition time is a crucial challenge for many imaging techniques. Compressed Sensing (CS) theory offers an appealing framework to address this issue since it provides theoretical guarantees on the reconstruction of sparse signals by projection on a low dimensional linear subspace. In this paper, we focus on a setting where the imaging device allows to sense a fixed set of measurements. We first discuss the choice of an optimal sampling subspace allowing perfect reconstruction of sparse signals. Its design relies on the random drawing of independent measurements. We discuss how to select the drawing distribution and show that a mixed strategy involving partial deterministic sampling and independent drawings can help breaking the so-called coherence barrier. Unfortunately, independent random sampling is irrelevant for many acquisition devices owing to acquisition constraints. To overcome this limitation, the notion of Variable Density Samplers (VDS) is introduced and defined as a stochastic process with a prescribed limit empirical measure. It encompasses samplers based on independent measurements or continuous curves. The latter are crucial to extend CS results to actual applications. We propose two original approaches to design continuous VDS, one based on random walks over the acquisition space, and one based on Traveling Salesman Problem. Following theoretical considerations and retrospective CS simulations in magnetic resonance imaging, we intend to highlight the key properties of a VDS to ensure accurate sparse reconstructions, namely its limit empirical measure and its mixing time 1 .

1. Introduction. Variable density sampling is a technique that is extensively used in various sensing devices such as magnetic resonance imaging (MRI), in order to shorten scanning time. It consists in measuring only a small number of random projections of a signal/image on elements of a basis drawn according to a given density. For instance, in MRI where measurements consist of Fourier (or more generally k-space) coefficients, it is common to sample the Fourier plane center more densely than the high frequencies. The image is then reconstructed from this incomplete information by dedicated signal processing methods. To the best of our knowledge, variable density sampling has been proposed first in the MRI context by [START_REF] Spielman | Magnetic resonance fluoroscopy using spirals with variable sampling densities[END_REF] where spiral trajectories were pushed forward. Hereafter, it has been used in this application (see e.g. [START_REF] Tsai | Reduced aliasing artifacts using variable-density k-space sampling trajectories[END_REF][START_REF] Kim | Simple analytic variable density spiral design[END_REF][START_REF] Park | Artifact and noise suppression in GRAPPA imaging using improved k-space coil calibration and variable density sampling[END_REF] to quote a few), but also in other applications such as holography [START_REF] Rivenson | Compressive fresnel holography[END_REF][START_REF] Marim | Compressed sensing with off-axis frequency-shifting holography[END_REF]. This technique can hardly be avoided in specific imaging techniques such as radio interferometry or tomographic modalities (e.g., X-ray) where sensing is made along fixed sets of measurements [START_REF] Wiaux | Spread spectrum for imaging techniques in radio interferometry[END_REF][START_REF] Sidky | Accurate image reconstruction from few-views and limitedangle data in divergent-beam ct[END_REF].

In the early days of its development, variable density sampling was merely an efficient heuristic to shorten acquisition time. It has recently found a partial justification in the Compressed Sensing (CS) literature. Even though this theory is not yet mature enough to fully explain the practical success of variable density sampling, CS provides good hints on how to choose the measurements (i.e., the density), how the signal/image should be reconstructed and why it works. Let us now recall a typical result emanating from the CS literature for orthogonal systems. A vector x ∈ C n is said s-sparse if it contains at most s non-zero entries. Denote by a i , i ∈ {1, . . . , n} the sensing vectors and by y i = a i , x the possible measurements. Typical CS results state that if the signal (or image) x is s-sparse and if

A =    a * 1 . . . a * n  
 satisfies an incoherence property (defined in the sequel), then m = O(s log(n) α ) measurements chosen randomly among the elements of y = Ax are enough to ensure perfect reconstruction of x. The constant α > 0 depends on additional properties on x and A. The set of actual measurements is denoted Ω ⊆ {1, . . . , n} and A Ω is the matrix formed by selecting a subset of rows of A in Ω. The reconstruction of x knowing y Ω = A Ω x is guaranteed if it results from solving the following ℓ 1 minimization problem:

min z∈C n z 1 subject to A Ω z = y Ω . (1.1)
Until recent works [START_REF] Rauhut | Compressive Sensing and structured random matrices[END_REF][START_REF] Juditsky | On low rank matrix approximations with applications to synthesis problem in compressed sensing[END_REF][START_REF] Candès | A probabilistic and ripless theory of compressed sensing[END_REF], no general theory for selecting the rows was available. In the latter, the authors have proposed to construct A Ω by drawing m rows of A at random according to a discrete probability distribution or density p = (p 1 , . . . , p n ). The choice of an optimal distribution p is an active field of research (see e.g. [START_REF] Chauffert | Variable density compressed sensing in MRI. Theoretical vs. heuristic sampling strategies[END_REF][START_REF] Krahmer | Beyond incoherence: stable and robust sampling strategies for compressive imaging[END_REF][START_REF] Adcock | Breaking the coherence barrier: asymptotic incoherence and asymptotic sparsity in compressed sensing[END_REF]) that remains open in many regards.

Drawing independent rows of A is interesting from a theoretical perspective, however it has little practical relevance since standard acquisition devices come with acquisition constraints. For instance, in MRI, the coefficients are acquired along piece-wise continuous curves on the k-space. The first paper performing variable density sampling in MRI [START_REF] Spielman | Magnetic resonance fluoroscopy using spirals with variable sampling densities[END_REF] has fulfilled this constraint by considering spiral sampling trajectories. The standard reference about CS-MRI [START_REF] Lustig | Sparse MRI: The application of compressed sensing for rapid MR imaging[END_REF] has proposed to sample the MRI signal along parallel lines in the 3D k-space. Though spirals and lines can be implemented easily on a scanner, it is likely that more general trajectories could provide better reconstruction results, or save more scanning time.

The main objective of this paper is to propose new strategies to sample a signal along more general continuous curves. Although continuity is often not sufficient for practical implementation on actual scanner, we believe that it is a first important step towards more physically plausible compressed sampling paradigms. As far as we know, this research avenue is relatively new. The problem was first discussed in [START_REF] Willett | Errata: Sampling trajectories for sparse image recovery[END_REF] and some heuristics were proposed. The recent contributions [START_REF] Polak | Performance bounds for grouped incoherent measurements in compressive sensing[END_REF][START_REF] Bigot | An analysis of block sampling strategies in compressed sensing[END_REF] have provided theoretical guarantees when sampling is performed along fixed sets of measurements (e.g. straight lines in the Fourier plane), but have not addressed generic continuous sampling curves yet.

The contributions of this paper are threefold. First, we bring a well mathematically grounded definition of variable density samplers and provide various examples. Second, we discuss how the sampling density should be chosen in practice. This discussion mostly relies on variations around the theorems provided in [START_REF] Rauhut | Compressive Sensing and structured random matrices[END_REF][START_REF] Candès | A probabilistic and ripless theory of compressed sensing[END_REF]. In particular, we justify the deterministic sampling of a set of highly coherent vectors to overcome the so-called "coherence barrier". In the MRI case, this amounts to deterministically sampling the k-space center. Our third and maybe most impacting contribution is to provide practical examples of variable density samplers along continuous curves and to derive some of their theoretical properties. These samplers are defined as parametrized random curves that asymptotically fit a target distribution (e.g. the one shown in Fig. 1.1 (a)). More specifically, we first propose a local sampler based on random walks over the acquisition space (see Fig. 1.1 (b)). Second, we introduce a global sampler based on the solution of a Travelling Salesman Problem amongst randomly drawn "cities" (see Fig. 1.1 (c)). In both cases, we investigate the resulting density. To finish with, we illustrate the proposed sampling schemes on 2D and 3D MRI simulations. The reconstruction results provided by the proposed techniques show that the PSNR can be substantially improved compared to existing strategies proposed e.g. in [START_REF] Lustig | Sparse MRI: The application of compressed sensing for rapid MR imaging[END_REF]. Our theoretical results and numerical experiments on retrospective CS show that two key features of variable density samplers are the limit of their empirical measure and their mixing properties. The rest of this paper is organized as follows. First, we introduce a precise definition of a variable density sampler (VDS) and recall CS results in the special case of independent drawings. Then, we give a closed form expression for the optimal distribution depending on the sensing matrix A, and justify that a partial deterministic sampling may provide better reconstruction guarantees. Hereafter, in Sections 3 and 4, we introduce two strategies to design continuous trajectories over the acquisition space. We show that the corresponding sampling distributions converge to a target distribution when the curve length tends to infinity. Finally, we demonstrate on simulation results that our TSP-based approach is promising in the MRI context (Section 5) since it outperforms its competing alternatives either in terms of PSNR at fixed sampling rate, or in terms of acceleration factor at fixed PSNR.

Notation. The main definitions used throughout the paper are defined in Tab. 1.1.

2.

Variable density sampling and its theoretical foundations. To the best of our knowledge, there is currently no rigorous definition of variable density sampling. Hence, to fill this gap, we provide a precise definition below. Definition 2.1. Let p be a probability measure defined on a measurable space Ξ. A stochastic process X = {X i } i∈N or X = {X t } t∈R+ on state space Ξ is called a p-variable density sampler if its empirical measure (or occupation measure) weakly

p =     p1 . . . pn     , 0 pi 1, n i=1 pi = 1 R n 1 ℓ1 norm defined for z ∈ C n by z 1 = n i=1 |zi| ∞ ℓ∞ norm defined for z ∈ C n by z ∞ = max 1 i n |zi| MRI application k = kx ky or   kx ky kz   Fourier frequencies R 2 or R 3 F * n d-dimensional discrete
F A set of points ⊂ H H N C(F )
The shortest Hamiltonian path (TSP) amongst points of set F ⊂ H

T (F, H) The length of C(F ) R+ T (F, R) For any set R ⊆ H, T (F, R) := T (F ∩ R, H) R+
converges to p almost surely, that is:

1 N N i=1 f (X i ) → p(f ) a.s. or 1 T T t=0 f (X t )dt → p(f ) a.s.
for all continuous bounded f . Example 1. In the case where X = (X i ) i∈N is a discrete time stochastic process with discrete state space Ξ = {1, . . . , n}, definition 2.1 can be slightly simplified. Let

us set Z N j = 1 N N i=1
1 Xi=j . The random variable Z N j represents the proportion of points that fall on position j. Let p denote a discrete probability distribution function. Using these notations, X is a p-variable density sampler if:

lim N →+∞ Z N j = p j a.s.
In particular, if (X i ) i∈N are i.i.d. samples drawn from p, then X is a p-variable density sampler. This simple example is the most commonly encountered in the compressed sensing literature and we will review its properties in paragraph 2.1.

Example 2. More generally, drawing independent random variables according to distribution p is a VDS if the space Ξ is second countable, owing to the strong law of large numbers.

Example 3. An irreducible aperiodic Markov chain on a finite sample space is a VDS for its stationary distribution (or invariant measure); see Section 3.3.

Example 4. In the deterministic case, for a dynamical system, definition 2.1 closely corresponds to the ergodic hypothesis, that is time averages are equal to expectations over space. We discuss an example that makes use of the TSP solution in section 4.

The following proposition directly relates the VDS concept to the time spent by the process in a part of the space, as an immediate consequence of the porte-manteau lemma (see e.g. [START_REF] Billingsley | Convergence of probability measures[END_REF]).

Proposition 2.2. Let p denote a Borel measure defined on a set Ξ. Let B ⊆ Ξ be a measurable set. Let X : R + → Ξ (resp. X : N → Ξ) be a stochastic process. Let µ denote the Lebesgue measure on R.

Define µ t X (B) = 1 t µ({s ∈ [0, t], X(s) ∈ B}) (resp. µ n X (B) = 1 n n i=1 1 X(i)∈B ).
Then, the following two propositions are equivalent:

(i) X is a p-VDS (ii)
Almost surely, ∀B ⊆ Ξ a Borel set with p(∂B) = 0, lim t→+∞ µ t X (B) = p(B) a.s.

(resp.) lim n→+∞ µ n X (B) = p(B) a.s.

Remark 1. Definition 2.1 is a generic definition that encompasses both discrete and continuous time and discrete and continuous state space since Ξ can be any measurable space. In particular, the recent CS framework on orthogonal systems [START_REF] Rauhut | Compressive Sensing and structured random matrices[END_REF][START_REF] Candès | A probabilistic and ripless theory of compressed sensing[END_REF] falls within this definition.

Definition 2.1 does not encompass some useful sampling strategies. We propose a definition of a generalized VDS, which encompasses stochastic processes indexed over a bounded time set.

Definition 2.3. A sequence {{X (n) t } 0 t Tn } n∈N is a generalized p-VDS if the sequence of occupation measures converges to p almost surely, that is:

1 T n Tn t=0 f (X (n) t )dt → p(f ) a.s.
Remark 2. Let (X t ) t∈R be a VDS, and (T n ) n∈N be any positive sequence such that T n → ∞. Then the sequence defined by X (n) t = X t for 0 t T n is a generalized VDS.

Example 5. Let Ξ = R2 , and consider r : [0, 1] → R + a strictly increasing smooth function. We denote by r -1 : [r(0), r(1)] → R its inverse function and by ṙ-1 the derivative of r -1 . Consider a sequence of spiral trajectories s

N : [0, N ] → R 2 defined by s N (t) = r t N cos(2πt) sin(2πt)
. Then s N is a generalized VDS for the distribution p defined by:

p(x, y) =    ṙ-1 √ x 2 +y 2 2π r(1) ρ=r(0) ṙ-1 (ρ)ρdρ if r(0) x 2 + y 2 r(1) 0 otherwise
A simple justification is that the time spent by the spiral in the infinitesimal ring

{(x, y) ∈ R 2 , ρ x 2 + y 2 ρ + dρ} is r -1 (ρ+dρ) r -1 (ρ) dt ∝ ṙ-1 (ρ).
2.1. Theoretical foundations -Independent VDS. CS theories provide strong theoretical foundations of VDS based on independent drawings. In this paragraph, we recall a typical result that motivates independent drawing in the ℓ 1 recovery context [START_REF] Rauhut | Compressive Sensing and structured random matrices[END_REF][START_REF] Foucart | A mathematical introduction to compressive sensing[END_REF][START_REF] Candès | A probabilistic and ripless theory of compressed sensing[END_REF][START_REF] Krahmer | Beyond incoherence: stable and robust sampling strategies for compressive imaging[END_REF][START_REF] Chauffert | Variable density compressed sensing in MRI. Theoretical vs. heuristic sampling strategies[END_REF][START_REF] Bigot | An analysis of block sampling strategies in compressed sensing[END_REF][START_REF] Adcock | Breaking the coherence barrier: asymptotic incoherence and asymptotic sparsity in compressed sensing[END_REF]. Using the notation defined in the introduction, let us give a slightly modified version of [START_REF] Rauhut | Compressive Sensing and structured random matrices[END_REF]Theorem 4.2].

Theorem 2.4. Let p = (p 1 , . . . , p n ) denote a probability distribution on {1, . . . , n} and Ω ⊂ {1, . . . , n} denote a random set obtained by m independent drawings with respect to distribution p. Let S ∈ {1, . . . , n} be an arbitrary set of cardinality s. Let x be an s-sparse vector with support S such that the signs of its non-zero entries is a Rademacher or Steinhaus sequence 2 . Define:

K(A, p) := max k∈{1...n} a k 2 ∞ p k (2.1)
Assume that:

m CK(A, p)s ln 2 6n η (2.2)
where C ≈ 26.25 is a constant. Then, with probability 1η, vector x is the unique solution of the ℓ 1 minimization problem (1.1).

Remark 3. Candès and Plan have stated stronger results in the case of real matrices in [START_REF] Candès | A probabilistic and ripless theory of compressed sensing[END_REF]. Namely, the number of necessary measurements was decreased to O(s log(n)), with lower constants and without any assumption on the vector signs. Their results have been derived using the so-called "golfing scheme" proposed in [START_REF] Gross | Recovering low-rank matrices from few coefficients in any basis[END_REF]. It is likely that these results could be extended to the complex case, however it would not change the optimal distribution which is the main point of this paper. We thus decided to stick to Theorem 2.4.

The choice of an accurate distribution p is crucial since it directly impacts the number of measurements required. In the MRI community, a lot of heuristics have been proposed so far to identify the best sampling distribution. In the seminal paper on CS-MRI [START_REF] Lustig | Sparse MRI: The application of compressed sensing for rapid MR imaging[END_REF], Lustig et al have proposed to sample the k-space using a density that polynomially decays towards high frequencies. More recently, Knoll et al have generalized this approach by inferring the best exponent from MRI image databases [START_REF] Knoll | Adapted random sampling patterns for accelerated MRI[END_REF]. It is actually easy to derive the theoretically optimal distribution, i.e. the one that minimizes the right hand-side in (2.2) as shown in Proposition 2.5, introduced in [START_REF] Chauffert | Variable density compressed sensing in MRI. Theoretical vs. heuristic sampling strategies[END_REF]. (i) the optimal distribution π ∈ ∆ n that minimizes K(A, p) is:

π i = a i 2 ∞ n i=1 a i 2 ∞ (2.3) (ii) K * (A) = K(A, π) = n i=1 a i 2 ∞ . Proof. (i) Taking p = π, we get K(A, π) = n i=1 a i 2 ∞ . Now assume that q = π, since n k=1 q k = n k=1 π k = 1, ∃j ∈ {1, . . . , n} such that q j < π j . Then K(A, q) a j 2 ∞ /q j > a j 2 ∞ /π j = n i=1 a i 2 ∞ = K(A, π). So, π is the distribution that minimizes K(A, p). (ii) This equality is a consequence of π's definition.
The theoretical optimal distribution only depends on the acquisition matrix, i.e. on the acquisition and sparsifying bases. For instance, if we measure some Fourier frequencies of a sparse signal in the time domain (a sum of diracs), we should sample the frequencies according to a uniform distribution, since a i ∞ = 1/ √ n for all 1

i n. In this case, K * (F) = 1 and the number of measurements m is proportional to s, which is in accordance with the seminal paper by Candès et al. [START_REF] Candès | Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information[END_REF].

Independent drawings in MRI. In the MRI case, the images are usually assumed sparse (or at least compressible) in a wavelet basis, while the acquisition is performed in the Fourier space. In this setting, the acquisition matrix can be written as A = F * Ψ. In that case, the optimal distribution only depends on the choice of the wavelet basis. The optimal distributions in 2D and 3D are depicted in Fig. Let us mention that similar distributions have been proposed in the literature. First, an alternative to independent drawing was proposed by Puy et al. [START_REF] Puy | On variable density compressive sampling[END_REF]. Their approach consists in selecting or not a frequency by drawing a Bernoulli random variable. Its parameter is determined by minimizing a quantity that slightly differs from K(A, p). Second, Krahmer and Ward [START_REF] Krahmer | Beyond incoherence: stable and robust sampling strategies for compressive imaging[END_REF] tried to unify theoretical results and empirical observations in the MRI framework. For Haar wavelets, they have shown that a polynomial distribution on the 2D k-space which varies as 1/(k 2

x +k 2 y ) is close to the optimal solution since it verifies K(A, p) = O(log(n)). Our numerical experiments have confirmed that a decay as a power of 2 is near optimal in 2D.

In the next section, we improve the existing theories by showing that a deterministic sampling of highly coherent vectors (i.e. those satisfying a i 2 ∞ ≫ 1 n ) may decrease the total number of required measurements. In MRI, this amounts to fully sampling the low frequencies, which exactly matches what has been done heuristically hitherto.

2.2. Mixing deterministic and independent samplings. In a recent work [START_REF] Chauffert | Variable density compressed sensing in MRI. Theoretical vs. heuristic sampling strategies[END_REF], we observed and partially justified the fact that a deterministic sampling of the low frequencies in MRI could drastically improve reconstruction quality. The following theorem proven in Appendix 1 provides a theoretical justification to this approach.

Theorem 2.6. Let S ∈ {1, . . . , n} be a set of cardinality s. Let x be an s-sparse vector with support S such that the signs of its non-zero entries is a Rademacher or Steinhaus sequence. Define the acquisition set Ω ⊆ {1, . . . n} as the union of:

(i) a deterministic set Ω 1 of cardinality m 1 . (ii) a random set Ω 2 obtained by m 2 independent drawing according to distribution p defined on {1 . . . n} \ Ω 1 . Denote m = m 1 + m 2 , Ω c 1 = {1, . . . , n} \ Ω 1 and let Ω = Ω 1 ∪ Ω 2 . Assume that: m m 1 + CK(A Ω c 1 , p)s ln 2 6n η (2.4)
where C = 7/3 is a constant, and

K(A Ω c 1 , p) = max i∈{1,...,n}\Ω1 a i 2 ∞ p i
. Then, with probability 1η, vector x is the unique solution of the ℓ 1 minimization problem (1.1).

This result implies that there exists an optimal partition between deterministically and randomly selected samples, which is moreover easy to compute. For example, consider the optimal distribution

p i ∝ a i 2 ∞ , then K * (A Ω c 1 ) = i∈{1,...,n}\Ω1 a i 2 ∞ .
If the measurement matrix contains rows with large values of a i ∞ , we notice from inequality (2.4) that these frequencies should be sampled deterministically, whereas the rest of the measurements should be obtained from independent drawings. This simple idea is another way of overcoming the so-called coherence barrier [START_REF] Krahmer | Beyond incoherence: stable and robust sampling strategies for compressive imaging[END_REF][START_REF] Adcock | Breaking the coherence barrier: asymptotic incoherence and asymptotic sparsity in compressed sensing[END_REF].

A striking example raised in [START_REF] Bigot | An analysis of block sampling strategies in compressed sensing[END_REF] is the following. Assume that

A = 1 0 0 F * n-1 .
The assumed optimal independent sampling strategy would consist in independently drawing the rows with distribution p 1 = 1/2 and p k = 1/ √ n -1 for k 2. According to Theorem 2.4, the number of required measurements is 2Cs ln 2 6n η . The alternative approach proposed in Theorem 2.6 basically performs a deterministic drawing of the first row combined with an independent uniform drawing over the remaining rows. In total, this scheme requires 1 + Cs ln 2 6n η measurements and thus reduces the number of measurements by almost a factor 2. Note that the same gain would be obtained by using independent drawings with rejection.

Mixed deterministic and independent sampling in MRI. In our experiments, we will consider wavelet transforms with three decomposition levels and the Symmlet basis with 10 vanishing moments. Fig. 2.2(a)-(b) shows the modulus of A's entries with a specific reordering in (b) according to decaying values of a i ∞ . This decay is illustrated in Fig. 2.2(c). We observe that a typical acquisition matrix in MRI shows large differences between its a i ∞ values. More Precisely, there is a small number of rows with a large infinite norm, sticking perfectly to the framework of Theorem 2.6. This observation justifies the use of a partial deterministic k-space sampling, which had already been used in [START_REF] Lustig | Sparse MRI: The application of compressed sensing for rapid MR imaging[END_REF][START_REF] Chauffert | Variable density compressed sensing in MRI. Theoretical vs. heuristic sampling strategies[END_REF]. In Fig. 2.2(d), the set Ω 1 is depicted for a fixed number of deterministic samples m 1 , by selecting the rows with the largest infinite norms. Hereafter, the strategy we adopt is driven by the previous remarks. All our sampling schemes are performed according to Theorem. 2.6: a deterministic part is sampled, and a VDS is performed on the rest of the acquisition space (e.g. the high frequencies in MRI).

Variable density samplers along continuous curves.

3.1. Why independent drawing can be irrelevant. In many imaging applications, the number of samples is of secondary importance compared to the time spent to collect the samples. A typical example is MRI, where the important variable to control is the scanning time. It depends on the total length of the pathway used to visit the k-space rather than the number of collected samples. MRI is not an exception and many other acquisition devices have to meet such physical constraints amongst which are scanning probe microscopes, ultrasound imaging, ecosystem monitoring, radio-interferometry or sampling using vehicles subject to kinematic constraints [START_REF] Willett | Errata: Sampling trajectories for sparse image recovery[END_REF]. In these conditions, measuring isolated points is not relevant and existing practical CS approaches consist in designing parametrized curves performing a variable density sampling. In what follows, we first review existing variable density sampling approaches based on continuous curves. Then, we propose two original contributions and analyze some of their theoretical properties. We mostly concentrate on continuity of the trajectory which is not sufficient for implementability in many applications. For instance, in MRI the actual requirement for a trajectory to be implementable is piecewise smoothness. More realistic constraints are discussed in Section 6.

3.2.

A short review of samplers along continuous trajectories. The prototypical variable density samplers in MRI were based on spiral trajectories [START_REF] Spielman | Magnetic resonance fluoroscopy using spirals with variable sampling densities[END_REF].

Similar works investigating different shapes and densities from a heuristic point of view were proposed in [START_REF] Tsai | Reduced aliasing artifacts using variable-density k-space sampling trajectories[END_REF][START_REF] Kim | Simple analytic variable density spiral design[END_REF][START_REF] Park | Artifact and noise suppression in GRAPPA imaging using improved k-space coil calibration and variable density sampling[END_REF]. The first reference to compressed sensing appeared in the seminal paper [START_REF] Lustig | Sparse MRI: The application of compressed sensing for rapid MR imaging[END_REF]. In this work, Lustig et al have proposed to perform independent drawings in a 2D plane (defined by the partition and phase encoding directions) and sample continuously along the orthogonal direction to design piecewise continuous schemes in the 3D k-space (see Fig. 3.1). These authors have also suggested to make use of randomly perturbed spirals. The main advantage of these schemes lies in their simplicity of practical implementation since they only require minor modifications of classical MRI acquisition sequences. Recent papers [START_REF] Polak | Grouped incoherent measurements for compressive sensing[END_REF][START_REF] Bigot | An analysis of block sampling strategies in compressed sensing[END_REF][START_REF] Boyer | An algorithm for variable density sampling with blockconstrained acquisition[END_REF] have generalized CS results from independent drawing of isolated measurements to independent drawings of blocks of measurements. In these contributions, the blocks can be chosen arbitrarily and may thus represent continuous trajectories. Interestingly, these authors have provided closed form expressions for the optimal distribution on the block set. Nevertheless, this distribution is very challenging to compute in large scale problems. Moreover, the restriction to sets of admissible blocks reduces the versatility of many devices such as MRI and can therefore impact the image reconstruction quality.

In many applications the length of the sampling trajectory is more critical than the number of acquired samples, therefore, finding the shortest pathway amongst random points drawn independently has been studied as a way of designing continuous trajectories [START_REF] Willett | Errata: Sampling trajectories for sparse image recovery[END_REF][START_REF] Wang | Smoothed random-like trajectory for compressed sensing MRI[END_REF]. Since this problem is NP-hard, one usually resorts to a TSP solver to get a reasonable suboptimal trajectory. To the best of our knowledge, the only practical results obtained using the TSP were given by Wang et al [START_REF] Wang | Smoothed random-like trajectory for compressed sensing MRI[END_REF]. In this work, the authors did not investigate the relationship between the initial sample locations and the empirical measure of the TSP curve. In Section 4, it is shown that this relationship is crucial to make efficient TSP-based sampling schemes.

In what follows, we first introduce an original sampler based on random walks on the acquisition space and then analyse its asymptotic properties. Our theoretical investigations together with practical experiments allows us to show that the VDS mixing properties play a central role to control its efficiency. This then motivates the need for more global VDS schemes.

3.3. Random walks on the acquisition space. Perhaps the simplest way to transform independent random drawings into continuous random curves consists in performing random walks on the acquisition space. Here, we discuss this approach and provide a brief analysis of its practical performance in the discrete setting. Through both experimental and theoretical results, we show that this technique is doomed to fail. However, we believe that this theoretical analysis provides a deep insight on what VDS properties characterize its performance.

Let us consider a time-homogeneous Markov chain X = (X n ) n∈N on the set {1, . . . , n} and its transition matrix denoted P ∈ R n×n . If X possesses a stationary distribution, i.e. a row vector p ∈ R n such that p = pP then, by definition, X is a p-variable density sampler.

Construction of the transition matrix P.

A classical way to design a transition kernel ensuring that (i) p is the stationary distribution of the chain and (ii) the trajectory defined by the chain is continuous, is the Metropolis algorithm [START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their applications[END_REF]. For a pixel/voxel position i in the 2D/3D acquisition space, let us define by N (i) ⊆ {1, . . . , n} its neighbourhood, i.e. the set of possible measurement locations allowed when staying on position i. Let |N (i)| denote the cardinal of N (i) and define the proposal kernel P * as P * i,j = |N (i)| -1 δ j∈N (i) . The Metropolis algorithm proceeds as follows:

1. from state i, draw a state i * with respect to the distribution P * i,: . 2. accept the new state i * with probability:

(3.1) q(i, i * ) = min 1, p(i * )P * i * ,i p(i)P * i,i * .
Otherwise stay in state i. The transition matrix P can then be defined by P i,j = q(i, j)P * i,j for i = j. The diagonal is defined in a such a way that P is a stochastic matrix. It is easy to check that p is an invariant distribution for this chain 3 . It is worth noticing that if the chain is irreducible positive recurrent (which is fulfilled if the graph is connected and the distribution p positive), the ergodic theorem ensures that X is a p-VDS.

Unfortunately, trajectories designed by this technique leave huge parts of the acquisition space unexplored (see Fig. 3.2 (a)). To circumvent this problem, we may allow the chain to jump to independent locations over the acquisition space. Let P be the Markov kernel corresponding to independent drawing with respect to p, i.e. Pi,j = p j for all 1 i, j n. Define:

(3.2) P (α) = (1 -α)P + α P ∀ 0 α 1.
Then the Markov chain associated with P (0) corresponds to a continuous random walk, while the Markov chain associated with P (α) , α > 0 has a nonzero jump probability. This means that the trajectory is composed of continuous parts of average length 1/α.

3.3.2.

Example. In Fig. 3.2, we show illustrations in the 2D MRI context where the discrete k-space is of size 64 × 64. On this domain, we set a distribution p which matches distribution π in Fig. 2.1 (a). We perform a random walk on the acquisition space until 10% of the coefficients are selected. In Fig. 3.2(a), we set α = 0 whereas α = 0.1 in Fig. 3.2(b). As expected, α = 0 leads to a sampling pattern where large parts of the k-space are left unvisited. The phenomenon is partially corrected using a nonzero value of α. Remark 4. Performing N iterations of the Metropolis algorithm requires O(N ) computations leading to a fast sampling scheme design procedure. In our experiments, we iterate the algorithm until m different measurements are probed. Therefore, the number of iterations N required increases non linearly with respect to m, and can be time consuming especially when R = m/n is close to 1. This is not a tough limitation of the method since the sampling scheme is computed off-line.

Compressed sensing results.

Let us assume4 that P(X 1 = i) = p i and that X i is drawn using P as a transition matrix. The following result provides theoretical guarantees about the performance of the VDS X: Proposition 3.1 (see [START_REF] Chauffert | From variable density sampling to continuous sampling using Markov chains[END_REF]). Let Ω := X 1 , . . . , X m ⊂ {1, . . . , n} denote a set of m indexes selected using the Markov chain X.

Then, with probability 1η, if

(3.3) m 12 ǫ(P) K 2 (A, p)s 2 log(2n 2 /η),
every s-sparse signal x is the unique solution of the ℓ 1 minimization problem.

The proof of this proposition is given in Appendix 2. Before going further, some remarks may be useful to explain this theoretical result.

Remark 5. Since the constant K 2 (A, p) appears in Eq. (3.3), the optimal sampling distribution using Markov chains is also distribution π, as proven in Proposition 2.5.

Remark 6. In contrast to Theorem 2.4, Proposition 3.1 provides uniform results, i.e. results that hold for all s-sparse vectors.

Remark 7. Ineq. 3.3 suffers from the so-called quadratick bottleneck (i.e. an O(s 2 log(n)) bound). It is likely that this bound can be improved to O(s log(n)) by developing new concentration inequalities on matrix-valued Markov chains.

Remark 8. More importantly, it seems however unlikely to avoid the spectral gap O(1/ǫ(P )) using the standard mechanisms for proving compressed sensing results. Indeed, all concentration inequalities obtained so far on Markov chains (see e.g. [START_REF] Lezaud | Chernoff-type bound for finite Markov chains[END_REF][START_REF] Kargin | A large deviation inequality for vector functions on finite reversible Markov chains[END_REF][START_REF] Paulin | Concentration inequalities for markov chains by marton couplings[END_REF]) depend on 1/ǫ(P ). The spectral gap satisfies 0 < ǫ(P )

1 and corresponds to mixing properties of the chain. The closer the spectral gap to 1, the fastest ergodicity is achieved. Roughly speaking, if |i-j| > 1/ǫ(P ) then X i and X j are almost independent random variables. Unfortunately, the spectral gap usually depends on the dimension n [START_REF] Diaconis | Geometric bounds for eigenvalues of markov chains[END_REF]. In our example, it can be shown using Cheeger's inequality that

ǫ(P ) = O n -1 d
if the stationary distribution π is uniform (see Appendix 3). This basically means that the number of measurements necessary to accurately reconstruct x could be as large as O(sn 1/d log(n)), which strongly limits the interest of this CS approach. The only way to lower this number consists in frequently jumping since Weyl's theorem [START_REF] Horn | Topics in matrix analysis[END_REF] ensures that ǫ(P (α) ) > α.

To sum up, the main drawback of random walks lies in their inability to cover the acquisition space quickly since they are based on local considerations. Keeping this in mind, it makes sense to focus on more global displacement strategies that allow a faster exploration of the whole acquisition domain. In the next section, we thus introduce this global sampling alternative based on TSP-solver. Our main contribution is the derivation of the link between a prescribed a priori sampling density and the distribution of samples located on the TSP solution so as to eventually get a VDS.

Travelling salesman-based VDS.

In order to design continuous trajectories, we may think of picking points at random and join them using a travelling salesman problem (TSP) solver. Hereafter, we show how to draw the initial points in order to reach a target distribution p. In this section, the probability distribution p is assumed to be a density.

Introduction.

The naive idea would consist in drawing some points according to the distribution p and joining them using a TSP solver. Unfortunately, the trajectory which results from joining all samples does not fit the distribution p, as shown in Fig. 4.1(b)-(d). To bring evidence to this observation, we performed a Monte Carlo study, where we drew one thousand sampling schemes, each one designed by solving the TSP on a set of independent random samples. We notice in Fig. 4.1 (d) that the empirical distribution of the points along the TSP curve, hereafter termed the final distribution, departs from the original distribution p. A simple intuition can be given to explain this discrepancy between the initial and final distributions in a d-dimensional acquisition space. Consider a small subset of the acquisition space ω. In ω, the number of points is proportional to p. The typical distance between two neighbors in ω is then proportional to p -1/d . Therefore, the local length of the trajectory in ω is proportional to pp -1/d = p 1-1/d = p. In what follows, we will show that the empirical measure of the TSP solution converges to a measure proportional to p 1-1/d .

Definitions. We shall work on the hypercube H

= [0, 1] d with d 2.
In what follows, {x i } i∈N * denotes a sequence of points in the hypercube H, independently drawn from a density p : H → R + . The set of the first N points is denoted

X N = {x i } i N .
Using the definitions introduced in Tab. 1.1, we introduce γ N : [0, 1] → H the function that parameterizes C(X N ) by moving along it at constant speed T (X N , H). Then, the distribution of the TSP solution reads as follows: Definition 4.1. The distribution of the TSP solution is denoted PN and defined, for any Borelian B in H by:

PN (B) = λ [0,1] γ -1 N (B) .
Remark 9. The distribution PN is defined for fixed X N . It makes no reference to the stochastic component of X N .

Remark 10. A more intuitive definition of PN can be given if we introduce other tools. For a subset ω ⊆ H, we denote the length of C(X N ) ∩ ω as T |ω (X N , H) = T (X N , H) PN (ω). Using this definition, it follows that:

(4.1) PN (ω) = T |ω (X N , H) T (X N , H) , ∀ω.
Then PN (ω) is the relative length of the curve inside ω. The proof of the theorem is given in Appendix 4. Remark 11. The T SP solution does not define as such a VDS, since the underlying process is finite in time. Nevertheless, since PN is the occupation measure of γ N , the following result holds:

Corollary 4.3. (γ N ) N ∈N is a generalized p VDS. Remark 12. The theorem indicates that if we want to reach distribution p in 2D, we have to draw the initial points with respect to a distribution proportional to p 2 , and to p 3/2 in 3D. Akin to the previous Monte Carlo study illustrating the behavior of the naive approach in Fig. 4.1 (top row), we repeated the same procedure after having taken this result into account. The results are presented in Fig. 4.1(e)-(g), in which it is shown that the final distribution now closely matches the original one (compare Fig. 4.1(g) with Fig. 4.1(a)).

Remark 13. Contrarily to the Markov chain approach for which we derived compressed sensing results in Proposition 3.1, the TSP approach proposed here is mostly heuristic and based on the idea that the TSP solution curve covers the space rapidly. An argument supporting this idea is the fact that in 2D, the TSP curve C(X N ) does not self-intersect. This property is clearly lacking for random walks.

Remark 14. One of the drawback of this approach is the TSP's NP-hardness. We believe that this is not a real problem. Indeed, there now exist very efficient approximate solvers such as the Concorde solver [START_REF] Applegate | Concorde TSP solver[END_REF]. It finds an approximate solution with 10 5 cities from a few seconds to a few hours depending on the required accuracy of the solution. The computational time of the approximate solution is not a real limitation since the computation is done off-line from the acquisition procedure. Moreover, many solvers are actually designed in such a way that their solution also fulfil Theorem. 4.2. For example, in 2D, to reach a sampling factor of R = 5 on a 256 × 256 image, one need N ≃ 10 4 cities, and an approximate solution is obtained in 142s. In 3D, for a 256 × 256 × 256 image, N ≃ 9 10 5 and an approximate solution is obtained in about 4 hours. In each case the solutions seem to be correctly approximated. In particular they do not self-intersect in 2D.

Experimental results in MRI.

In this section, we focus on the reconstruction results by minimizing the ℓ 1 problem (1.1) with a simple MRI model: A = F * Ψ, where Ψ denote the inverse Symmlet-10 transform 5 . The solution is computed using Douglas-Rachford's algorithm [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF]. We consider an MR image of size 256 × 256 × 256 as a reference, and perform reconstruction for different discrete sampling strategies. Every sampling scheme was regridded using a nearest neighbour approach to avoid data interpolation. 6 5.1. 2D-MRI. In 2D, we focused on a single slice of the MR image and considered its discrete Fourier transform as the set of possible measurements. First, we found the best made a comparison of independent drawings with respect to various distributions in order to find heuristically the best sampling density. Then we explored the performance of the two proposed methods to design continuous schemes: random walks and Travelling Salesman Problem. We also compared our solution to classical MRI sampling schemes. In every sampling schemes, the number of measurements is 5 We focused on ℓ 1 reconstruction since it is central in the CS theory. The reconstruction quality can be improved by considering more a priori knowledge on the image. Moreover we considered a simple MRI model, but our method can be extended to parallel MRI [START_REF] Pruessmann | SENSE: sensitivity encoding for fast MRI[END_REF], or spread-spectrum techniques [START_REF] Haldar | Compressed-sensing MRI with random encoding[END_REF][START_REF] Puy | Spread spectrum magnetic resonance imaging[END_REF]. 6 We provide Matlab codes to reproduce the proposed experiments here: http://chauffertn.free.fr/codes.html the same and equals 20% of the number of pixels in the image, so that the sampling factor R is equal to 5. In cases where the sampling strategy is based on randomness (VDS, random walks, TSP...), we performed a Monte Carlo study by generating 100 sampling patterns for each variable density sampler.

5.1.1. Variable density sampling using independent drawings. Here, we assessed the impact of changing the sampling distribution using independent drawings. In all experiments, we sampled the Fourier space center deterministically as shown on Figure 5.1.

< 0.1 < 0.1 < 0.1 < 0.1 < 0.1 < 0.1 < 0.1
Table 5.1 shows that the theoretically-driven optimal distribution π is outperformed by the best heuristics. Amongst the latter, the distribution leading to the best reconstruction quality decays as 1/|k| 2 , which is the distribution used by Krahmer and Ward [START_REF] Krahmer | Beyond incoherence: stable and robust sampling strategies for compressive imaging[END_REF] as an approximation of π for Haar wavelets. The standard deviation of the PSNR is negligible compared to the mean values and for a given distribution, each reconstrucion PSNR equals its average value at the precision used in Tab. 5.1.

Continuous VDS. In this part we compared various variable density samplers:

• Random walks with a stationary distribution proportional to 1/|k| 2 and different average chain lengths of 1/α, • TSP-based sampling with distributions proportional to 1/|k| 2 and π, • Classical MRI sampling strategies such as spiral, radial and radial with random angles. The choice of the spiral follows Example 5: the spiral is parame-

terized by s : [0, T ] → R 2 , θ → r(θ/T ) cos θ sin θ where r(t) := r(0)r(1)
r(1)-t(r(1)-r(0)) , so as the spiral density decays as 1/|k| 2 . The sampling schemes are presented in Fig. 5.1 and the reconstruction results in Tab. 5.2.

(a) (b) (c) (d) (e) (f) (g) (h)
As predicted by the theory, the shorter the chains the better the reconstructions. The optimal case corresponds to chains of length 1 (α = 1) i.e. corresponding to independent VDS. When the chain is too long, large k-space areas are left unexplored, and the reconstruction quality decreases.

Besides, the use of a target distribution proportional to 1/|k| 2 instead of π for TSP-based schemes provides slightly better reconstruction results.

We also considered more classical sampling scheme. We observe that the spiral scheme and the proposed ones provide more accurate reconstruction results than radial schemes. We believe that the main reason underlying these different behaviors is closely related to the sampling rate decay from low to high frequencies, which is proportional to 1/|k| for radial schemes.

3D-MRI.

Since VDS based on Markov chains have shown rather poor reconstruction results compared to the TSP-based sampling schemes in 2D simulations, we only focus on comparing TSP-based sampling schemes to classical CS sampling schemes. Moreover, the computational load to treat 3D images being much higher than in 2D, we only perform one drawing per sampling scheme in the following experiments. Experiments in 2D suggest that the reconstruction quality is not really impacted by the realization of a particular sampling scheme, except for drawing with Markov chains or with radial with random angles, which are not considered in our 3D experiments.

5.2.1.

Variable density sampling using independent drawings. The first step of the TSP-based approach is to identify a relevant target distribution. For doing so, we consider independent drawings as already done in 2D. The results are summarized in Tab. 5.3. In this experiment, we still use a number of measurements equal to 20% of the total amount (R = 5).

The best reconstruction result is achieved with d = 2 and not the theoretically optimal distribution π. This illustrates the importance of departing from the sole sparsity hypothesis under which we constructed π. Natural signals have a much richer structure. For instance wavelet coefficients tend to become sparser as the resolution levels increase, and this feature should be accounted for to derive optimal sampling densities for natural images (see Section 6.) Let us now compare the reconstruction results using the TSP based method and the method proposed in the original CS-MRI paper [START_REF] Lustig | Sparse MRI: The application of compressed sensing for rapid MR imaging[END_REF]. These two sampling strategies are depicted in Fig. 5.2. For 2D independent drawings, we used the distribution providing the best reconstruction results in 2D, i.e. proportional to 1/|k| 2 . The TSP-based schemes were designed by drawing city locations independently with respect to a distribution proportional to p 3 2 . According to Theorem 4.2 this is the correct way to reach distribution p after joining the cities with constant speed along the TSP solution path. The experiments were performed with p = π (see Fig. 2.1 (b)), and p ∝ 1/|k| 2 , since the latter yielded the best reconstruction results in the 3D independent VDS framework. We also compared these two continous schemes to 3D independent drawings with respect to a distribution proportional to 1/|k| 2 .

Reconstruction results with an sampling rate R = 8.8 are presented in Fig. 5.4, with a zoom on the cerebellum. The reconstruction quality using the proposed sampling scheme is better than the one obtained from classical CS acquisition and contains less artifacts. In particular, the branches of the cerebellum are observable with our proposed sampling scheme only. At higher sampling rate, we still observe less artifacts with the proposed schemes, as depicted in Fig. 5.5 with a sampling rate R = 14.9. Moreover, Fig. 5.3 shows that our proposed method outperforms the method proposed in [START_REF] Lustig | Sparse MRI: The application of compressed sensing for rapid MR imaging[END_REF] by up to 2dB. If one aims at reaching a fixed PSNR, we can increase r by more than 50% using the TSP based strategy. In other words, we could expect a substantial decrease of scanning time by using more advanced sampling strategies than those proposed until to now.

The two different choices of the target density π and ∝ 1/|k| 2 provide similar results. This is a bit surprising since 3D independent VDS with these two probability distributions provide very different reconstruction results (see Tab. 5.3). A potential explanation for that behavior is that the TSP tends to "smooth out" the target distribution. An independent drawing would collect very few Fourier coefficients in the blue zones of Fig. 2.1, notably the vertical and horizontal lines crossing the Fourier plane center. Sampling these zones seems to be of utmost importance since they contain high energy coefficients. The TSP approach tends to sample these zones by crossing the lines.

Perhaps the most interesting fact is that Fig. 5.3 shows that the TSP based sampling schemes provide results that are similar to independent drawings up to important sampling rates such as 20. We thus believe that the TSP solution proposed in this paper is near optimal since it provides results similar to unconstrained acquisition schemes. The price to be paid by integrating continuity constraints is thus almost null.

6. Discussion and perspectives. In this paper, we investigated the use of variable density sampling along continuous trajectories. Our first contribution was to provide a well-grounded mathematical definition of p-variable density samplers (VDS) as stochastic processes with a prescribed limit empirical measure p. We identified PSNR (dB) to address this point in Section 2 and showed that deterministic sampling could reduce the total number of required measurements. The analysis we performed closely followed the proofs proposed in [START_REF] Rauhut | Compressive Sensing and structured random matrices[END_REF][START_REF] Candès | A probabilistic and ripless theory of compressed sensing[END_REF] and was based solely on sparsity hypotheses on the signal/image to be reconstructed. The numerical experiments we performed indicate that heuristic densities still outperform the theoretical optimal ones. This suggests that the optimality critera used so far to derive target sampling densities does not account for the whole structure of the sought signal/image. Although sparsity is a key feature that characterizes natural signals/images, we believe that introducing stronger knowledge like structured sparsity might contribute to derive a new class of optimal densities that would compete with heuristic densities.

To the best of our knowledge, the recent paper [START_REF] Adcock | Breaking the coherence barrier: asymptotic incoherence and asymptotic sparsity in compressed sensing[END_REF] is the first contribution that addresses the design of sampling schemes by accounting for a simple structured sparsity hypothesis. The latter assumes that wavelet coefficients become sparser as the resolution increases. The main conclusion of the authors is the same as that of Theorem 2.6 even though it is based on different arguments: the low frequencies of a signal should be sampled deterministically.

Finally, let us notice that the best empirical convex reconstruction techniques do not rely on the resolution of a simple ℓ 1 problem such as (1.1). They are based on regularization with redundant frames and total variation for instance [START_REF] Boyer | HYR2PICS: Hybrid Regularized Reconstruction for Combined Parallel Imaging and Compressive Sensing in MRI[END_REF]. The signal model, the target density and the reconstruction algorithm should clearly be considered simultaneously to make a substantial leap on reconstruction guarantees.

What VDS properties govern their practical efficiency? In Section 3, it was shown that the key feature characterizing random walks efficiency was the mixing properties of the associated stochastic transition matrix. In order to derive CS results using generic random sets rather than point processes or random walks, it seems important to us to find an equivalent notion of mixing properties. How to generate VDS with higher degrees of regularity? This is probably the most important question from a practical point of view. We showed that the TSP based VDS outperformed more conventional sampling strategies by substantial acceleration factors for a given PSNR value or recovers 3D images at an improved PSNR for a given acceleration factor. However, this approach may not really be appealing for many applications: continuity is actually not a sufficient condition for making acquisition sequences implementable on devices like MRI scanners or robot motion where additional kinematic constraints such as bounded first (gradients) and second (slew rate) derivatives should be taken into account. Papers such as [START_REF] Lustig | A fast method for designing time-optimal gradient waveforms for arbitrary k-space trajectories[END_REF] derive time-optimal waveforms to cross a given curve using optimal control. By using this approach, it can be shown that the angular points on the TSP trajectory have to be visited with a zero speed. This strongly impacts the scanning time and the distribution of the parametrized curve. The simplest strategy to reduce scanning time would thus consist in smoothing the TSP trajectory, however this approach dramatically changes the target distribution which was shown to be a key feature of the method. The key element to prove our TSP Theorem 4.2 was the famous Beardwood, Halton and Hammersley theorem [START_REF] Beardwood | The shortest path through many points[END_REF]. To the best of our knowledge, extending this result to smooth trajectories remains an open question 7 . Recent progresses in that direction were obtained in recent papers such as [START_REF] Le Ny | On the dubins traveling salesman problem. Automatic Control[END_REF], but they do not provide sufficient guarantees to extend Theorem 4.2. Answering this question is beyond the scope of this paper. We believe that the work [START_REF] Teuber | Dithering by differences of convex functions[END_REF] based on attraction and repulsion potentials opens an appealing research avenue for solving this issue.

Appendix 1 -proof of Theorem 2.6. For a symmetric matrix M , we denote by λ max (M ) its largest eigenvalue and by M the largest eigenvalue modulus. The crucial step to obtain Theorem 2.6 is Proposition 6.1 below. The rest of the proof is the same as the one proposed in [START_REF] Rauhut | Compressive Sensing and structured random matrices[END_REF] and we refer the interested reader to [42, Section 7.3] for further details. Proposition 6.1. Let Ω = Ω 1 ∪ Ω 2 ⊆ {1, . . . , n} be a set constructed as in Theorem 2.6. Define

ãi = a i if i ∈ Ω 1 a i / √ p i if i ∈ {1 . . . n} \ Ω 1 . and à =            ãΩ1(1) . . . ãΩ1(m1) 1 √ m2 ãΩ2(1) . . . 1 √ m2 ãΩ2(m2)            ∈ C m×n . (6.1)
Then for all δ ∈ [0, 1 2 ]:

P ÃS * ÃS -I s δ 2s exp - m 2 δ 2 CK 2 2 s
. 7 To be precise, many crucial properties of the length of the shortest path used to derive asymptotic results are lost. The most important one is subadditivity [START_REF] Steele | Subadditive euclidean functionals and nonlinear growth in geometric probability[END_REF].

where ÃS ∈ C m×s is the matrix composed of the s columns of à belonging to S. C = 7/3 is a constant.

The proof of this proposition relies heavily on the matrix Bernstein inequality below [START_REF] Tropp | User-friendly tail bounds for sums of random matrices[END_REF]. Proposition 6.2 (Matrix Bernstein inequality). Let Z k be a finite sequence of independent, random, self-adjoint matrices in C d×d . Assume that each random matrix satisfies

E(Z k ) = 0 and λ max (Z k ) R a.s. Denote σ 2 = k E(Z 2 k ) .
Then, for all t 0:

P k Z k t 2d exp - t 2 /2 σ 2 + Rt/3 .
We are now ready to prove Proposition 6.1. Proof. For any vector v ∈ C n , denote by v S ∈ C s the vector composed of the entries of v belonging to S ⊆ {1, . . . , n}. Consider the random sequence X 1 , . . . , X m2 where X i = j ∈ {1 . . . n}\Ω 1 with probability p j , and denote by Ω Xi ) 2 ) K 2 2 s(I s -M 1 ) using (ii). We can then proceed as follows using points (iii) to (iv):

E(Z 2 Xi ) = M 2 1 -2M 1 + I s + E((ã S Xi ãS * Xi ) 2 ) + 2M 1 E(ã S Xi ãS * Xi ) -2E(ã S Xi ãS * Xi ) ≤ M 2 1 -2M 1 + I s + K 2 2 s(I s -M 1 ) + 2M 1 (I s -M 1 ) -2(I s -M 1 ) = -(I s -M 1 ) 2 + K 2 2 s(I s -M 1 ) K 2 2 sI s . Then m2 i=1 E(Z 2 Xi ) m 2 K 2 2 s.

By noticing that ãS

Xi ãS * Xi -I s Z Xi ãS Xi ãS * Xi , we obtain Z Xi K 2 2 s.
Finally, by applying Bernstein inequality to the sequence of matrices Z X1 , . . . Z Xm 2 , we derive for all t 0:

P j∈Ω2 Z j t 2s exp - t 2 /2 m 2 K 2 2 s + K 2 2 st/3
.

Plugging δ := t/m 2 , and noticing that δ 1/2 ⇒ 2(1 + δ/3) 2(1 + δ/3) 7/3, the announced result is shown.

Appendix 2 -proof of Proposition 3.1. Our approach relies on the following perfect recovery condition introduced in [START_REF] Juditsky | On verifiable sufficient conditions for sparse signal recovery via ℓ 1 minimization[END_REF]:

Proposition 6.3 ([25]). If A Ω ∈ R m×n satisfies γ(A Ω ) = min Y∈R m×n I n -Y T A Ω ∞ < 1 2s ,
all s-sparse signals x ∈ R n are recovered exactly by solving the ℓ 1 minimization problem (1.1). We noted A ∞ the maximal modulus of all the entries of A. This can be seen as an alternative to the mutual coherence [START_REF] Donoho | Compressed sensing[END_REF]. We limit our proof to the real case but it could be extended to the complex case using a slightly different proof.

We aim at finding Y ∈ R m×n , such that I n -Y T A Ω ∞ < 1 2s , for a given positive integer s, where A Ω is the sensing matrix defined in Proposition 3.1. Following [START_REF] Juditsky | On low rank matrix approximations with applications to synthesis problem in compressed sensing[END_REF],

we set Θ i = aia T i pi and use the decomposition I n = A T A = n i=1 p i Θ i . We consider a realization of the Markov chain X 1 , . . . , X m , with X 1 ∼ p and X i ∼ P Xi-1,: for i > 1. Let us denote W m = 1 m m l=1 Θ X l . Then W m may be written as Y T A Ω . Lemma 6.4. ∀ 0 < t 1, (6.2) P ( I n -W m ∞ t) n(n + 1)e ǫ(P)

5 exp -mt 2 ǫ(P) 12K 2 (A, p) .

Before proving the lemma, let first recall a concentration inequality for finite-state Markov chains [START_REF] Lezaud | Chernoff-type bound for finite Markov chains[END_REF]. Proposition 6.5. Let (P, p) be an irreducible and reversible Markov chain on a finite set G of size n with transition matrix P and stationary distribution p. Let f : G → R be such that n i=1 p i f i = 0, f ∞ 1 and 0 < n i=1 f 2 i p i b 2 . Then, for any initial distribution q, any positive integer m and all 0 < t 1,

P 1 m m i=1
f (X i ) t e ǫ(P)

5 N q exp -mt 2 ǫ(P) 4b 2 (1 + g(5t/b 2 ))

where N q = ( n i=1 ( qi pi ) 2 p i ) 1/2 and g is given by g(x) = 1 2 ( √ 1 + x -(1x/2)). Now, we can prove Lemma 6.4 Proof. By applying Proposition 6.5 to a function f and then to its opposite -f , we get:

P 1 m m i=1 f (X i ) t 2e
ǫ(P)

5 N q exp -mt 2 ǫ(P) 4b 2 (1 + g(5t/b 2 )) .

Then we set f (X i ) = (I n -Θ Xi ) (a,b) /K(A, p) as real-valued function. Recall that p satisfies n i=1 p i f (X i ) = 0. Since f ∞ 1, b = 1 and t 1, we deduce 1+g(5t) < 3. Moreover, since the initial distribution is p, q i = p i , ∀i and thus N q = 1. Finally, resorting to a union bound enables us to extend our result for the (a, b)th entry to the whole infinite norm of the n × n matrix I n -W m (6.2). Finally, set s ∈ N * and η ∈ (0, 1). If m satisfies Ineq. (3.3), then To show Prop.6.7, we need to introduce the boundary TSP. For a set of points F and an area R, we denote by T B (F, R) its length on the set F ∩ R. The boundary TSP is defined as the shortest Hamiltonian tour on F ∩ R for the metric obtained from the Euclidean metric by the quotient of the boundary of R, that is d(a, b) = 0 if a, b ∈ ∂R. Informally, it matches the original TSP while being allowed to travel along the boundary for free. We refer to [START_REF] Frieze | Probabilistic analysis of the TSP[END_REF] for a complete description of this concept.

P I n -W m ∞ 1 2s < η .
We shall use a set of classical results on TSP and boundary TSP, that may be found in the survey books [START_REF] Frieze | Probabilistic analysis of the TSP[END_REF] and [START_REF] Yukich | Probability theory of classical Euclidean optimization problems[END_REF]. Lemma 6.9. Let F denote a set of n points in H.

1. The boundary TSP is superadditive, that is, if R 1 and R 2 have disjoint interiors.

T B (F, R 1 ∪ R 2 ) T B (F, R 1 ) + T B (F, R 2 ). (6.6) 2. The boundary TSP is a lower bound on the TSP, both globally and on subsets.

If R 2 ⊂ R 1 :

T (F, R) T B (F, R) (6. Let N i be the number of points of X N in ω i .

Since N i N , we may use the bound (6.9) to get: p ⊗N -a.s, ∀i.

Finally, by the law of large numbers, almost surely N i /N → p(ω i ) = ωi p(x)dx. The law of any point x j conditioned on being in ω i has density p/p(ω i ). By applying Theorem 6.8 to the hypercubes ω i and H we thus get:

Fig. 1 . 1 .

 11 Fig. 1.1. (a): Target distribution π. Continuous random trajectories reaching distribution π based on Markov chains (b) and on a TSP solution (c).

= 1 -

 1 Fourier transform on an n pixels image C n×n Ψn d-dimensional inverse discrete Wavelet transform on an image of n pixels C n×n F * n and Ψn are denoted F * and Ψ if no ambiguity VDS Ξ A measurable space which is typically {1, . . . , n} or [0, 1] d H The unit cube [0, 1] d p A probability measure defined on Ξ p(f ) = x∈Ξ f (x) dp(x), for f continuous and bounded R λ [0,1] The Lebesgue measure on the interval [0, 1] X = (Xn) n∈N * A time-homogeneous Markov chain on the state space {1, . . . , n} {1, . . . n} N * P := (Pij ) 1 i,j n the transition matrix: Pij := P(X k = j|X k-1 = i), ∀k > 1 R n×n λi(P) The ordered eigenvalues of P: 1 = λ1(λ2(P), the spectral gap of P [-1, 1]

Proposition 2 . 5 .

 25 Denote by K * (A) := min p∈∆n K(A, p).

Fig. 2 . 1 .

 21 Fig. 2.1. Optimal distribution π for a Symmlet-10 tranform in 2D (a) and a maximal projection of the optimal distribution in 3D (b).

Fig. 2 . 2 .

 22 Fig. 2.2. (a): Absolute magnitudes of A for a 2D Symmlet basis with 10 vanishing moments and 3 levels of decomposition. (b): same quantities as in (a) but sorted by decaying a i ∞ (i.e. by decreasing order). (c): decay of a i ∞. (d): Set Ω 1 depicted in the 2D k-space.

Fig. 3 . 1 .

 31 Fig. 3.1. Classical CS-MRI strategy. (a): 2D independent sampling according to a distribution π. (b): measurements performed in the orthogonal readout direction.

Fig. 3 . 2 .

 32 Fig. 3.2. Example of sampling trajectories in 2D MRI. (a) (resp. (b)): 2D sampling scheme of the k-space with α = 0 (resp. α = 0.1). Drawings are performed until 10% of the coefficients are selected (m = 0.1n).

4. 3 .

 3 Main results. Our main theoretical result introduced in[START_REF] Chauffert | Travelling salesman-based variable density sampling[END_REF] reads as follows:Theorem 4.2. Define the density p = p (d-1)/d H p (d-1)/d (x)dxwhere p is a density defined on H. Then almost surely with respect to the law p ⊗N of the random points sequence {x i } i∈N * in H, the distribution PN converges in distribution to p: PN

Fig. 4 . 1 .

 41 Fig. 4.1. Illustration of the TSP-based sampling scheme to reach distribution π. (a): distribution π. (b) (resp. (e)): independent drawing of points from distribution π (resp. ∝ π 2 ). (c) (resp (f )): solution of the TSP amongst points of (b) (resp. (e)) . (d) and (g): Monte Carlo study: average scheme over one thousand drawings of sampling schemes, with the same color scale as in (a).

Fig. 5 . 1 .

 51 Fig. 5.1. 2D continuous sampling schemes based on random walks with α = .1 (a), α = .01 (b), α = .001 (c), and based on TSP solutions with distributions proportional to π (d) and to 1/|k| 2 (e). Classical sampling schemes: spiral (f ), radial (g) and radial with random angles (h).

Fig. 5 . 3 .

 53 Fig. 5.3. Quality of 3D reconstructed images in terms of PSNR as a function of sampling rates R for various sampling strategies: independent drawings with respect to distribution ∝ 1/|k| 2 (dashed blue line), TSP-based sampling with target densities π (black line) and ∝ 1/|k| 2 (red line), and parallel lines with 2D independent drawing with respect to ∝ 1/|k| 2 distribution (green line) as depicted in Fig. 5.2[Top row].

Fig. 5 . 4 .

 54 Fig. 5.4. Reconstruction results for R = 8.8 for various sampling strategies. Top row: TSPbased sampling schemes (PSNR=42.1 dB). Bottom row: 2D random drawing and acquisitions along parallel lines (PSNR=40.1 dB). Sagital view (left) and zoom on the cerebellum (right).

Fig. 5 . 5 .

 55 Fig. 5.5. Reconstruction results for R = 14.9 for various sampling strategies. Left: TSP-based sampling schemes (PSNR=39.8 dB). Right: 2D random drawing and acquisitions along parallel lines (PSNR=38.3 dB).

2

 2 the set {X 1 , . . . X m2 }. Denote by M 1 := i∈Ω1 a S i a S i * . Consider the matrices Z j := M 1 + ãS j ãS * j -I s . According to Eq. (6.1), we get by construction: ÃS * ÃS -I s = 1 m 2 j∈Ω2 Z j . Since I s = n i=1 a S i a S i * , we notice that ∀i ∈ {1, . . . , m 2 } (i) E(Z Xi ) = 0, (ii) E(ã S Xi ãS * Xi ) = I s -M 1 . Moreover, we have (iii) 0 I s -M 1 I s and (iv) 0 M 1 I s .Using the identity (ã

Fig. 6 . 1 .

 61 Fig. 6.1. Illustration of the proof of Remark 8 in dimension 2. Left: regular grid with n = √ n × √ n nodes. Right: Graph partitioning in B and B c with p(B) = 1/2.

  7) T |R2 (F, R 1 ) T B (F, R 2 ) (6.8) 3. The boundary TSP approximates well the TSP [53, Lemma 3.7]):|T (F, H) -T B (F, H)| = O(n (d-2)/(d-1) ). (6.9) 4. The TSP in H is well-approximated by the sum of TSPs in a grid of h d congruent hypercubes [18, Eq. (33)].|T (F, H) -h d i=1 T (F, ω i )| = O(n (d-2)/(d-1) ). (6.10)We now have all the ingredients to prove the main results.Proof. [Proof of Proposition 6.7]i∈I T B (X N , ω i ) (6.6) T B (X N , H) (6.7) T (X N , H) = h d i=1 T | ωi (X N , H) (6.10) h d i=1 T (X N , ω i ) + O(N (d-1)/(d-2) )

T

  (X N , ω i ) N (d-1)/d = lim N →∞ T B (X N , ω i ) N (d-1)/d .Using the fact that there are only finitely many ω i , the following equalities hold almost surely:T |ωi (X N , H) N (d-1)/d .Since the boundary TSP is a lower bound (cf. Eqs. (6.8)-(6.7)) to both local and global TSPs, the above equality ensures that:lim N →∞ T B (X N , ω i ) N (d-1)/d = lim N →∞ T (X N , ω i ) N (d-1)/d (6.12) = lim N →∞ T |ωi (X N , H) N (d-1)/d

T

  (X N , ω i ) N (d-1)/d = β(d) ωi p(x) (d-1)/d dx p ⊗N -a.s, ∀i. and lim N →+∞ T (X N , H) N (d-1)/d = β(d) H p(x) (d-1)/d dx p ⊗N -a.s, ∀i.

Table 1 . 1

 11 General notation used in the paper.

	Notation	Definition

Table 5 . 1

 51 Quality of reconstruction results in terms of PSNR for 2D sampling with variable density independent drawings.

	π	polynomial decay: (k 2 x + k 2 y ) -d/2 d = 1 d = 2 d = 3 d = 4 d = 5 d = 6
	mean PSNR (dB) 35.6	36.4	36.4	36.3	36.0	35.5	35.2
	std dev.

Table 5 . 2

 52 Quality of reconstruction results in terms of PSNR for continuous sampling trajectories.

	Markovian drawing (α) 0.1 0.01 0.001 mean PSNR 35.7 34.6 33.5	TSP sampling ∝ π ∝ 1/|k| 2 35.6 36.1	spiral radial 35.6 34.1	radial random 33.1
	std dev.	0.1	0.3	0.6	0.1	0.1		0.4
	max value	36.0 35.1	34.8	35.9	36.2		34.0
	in Fig. 5.1:

Table 5 . 3

 53 Quality of reconstruction results in terms of PSNR for sampling schemes based on 3D variable density independent drawings, with densities ∝ 1/k d and π, and with 20% of measured samples.

	d	1	2	3	4	π
	PSNR (dB) 44.78 45.01 44.56 44.03 42.94
	5.2.2. Efficiency of the TSP sampling based strategy.

A Rademacher (resp. Steinhaus) random variable is uniformly distributed on {-1; 1} (resp. on the torus {z ∈ C; |z| = 1}).

If the neighboring system is such that the corresponding graph is connected, then the invariant distribution is unique.

By making this assumption, there is no burn-in period and the chain X converges more rapidly to its stationary distribution p.

Acknowledgments. The authors wish to thanks Yves Wiaux, Fabrice Gamboa, Jérémie Bigot, Laurent Miclo, Alexandre Vignaud and Claire Boyer for fruitful discussions and feedback. This research was supported by the Labex CIMI through a 3 months invitation of Philippe Ciuciu. This work was partially supported by ANR SPH-IM-3D (ANR-12-BSV5-0008), by the FMJH Program Gaspard Monge in optimization and operation research (MAORI project), and by the support to this program from EDF.

through both theoretical and experimental results two key features characterizing their efficiency: their empirical measure as well as their mixing properties. We showed that VDS based on random walks were doomed to fail since they were unable to quickly cover the whole acquisition space. This led us to propose a two-step alternative that consists first of drawing random points independently and then joining them using a Travelling Salesman Problem solver. In contrast to what has been proposed in the literature so far, we paid attention to the manner the points have to be drawn so as to reach a prescribed empirical measure. Strikingly, our numerical results suggest that the proposed approach yields reconstruction results that are nearly equivalent to independent drawings. This suggests that adding continuity constraints to the sampling schemes might not be so harmful to derive CS results.

We believe that the proposed work opens many perspectives as outlined in what follows.

How to select the target density? We recalled existing theoretical results

In other words, with probability at least 1η, every s-sparse signal can be recovered by ℓ 1 minimization (1.1). Remark 15. It is straightforward to derive a similar result to Theorem 2.6 and thus to justify that a partial deterministic sampling reduces the total number of measurements required for perfect recovery.

Appendix 3 -proof of Remark 8. In this part, we prove that for a random walk with uniform stationary distribution p, ǫ(P) = O(n -1 d ). We use geometric bounds known as Cheeger's inequality in [START_REF] Diaconis | Geometric bounds for eigenvalues of markov chains[END_REF] and Conductance Bounds in [START_REF] Jerrum | Approximating the permanent[END_REF][START_REF] Brémaud | Markov chains: Gibbs fields, Monte Carlo simulation, and queues[END_REF]. Let us recall a useful result concerning finite state space irreducible reversible transition matrices P.

The capacity of a set B ⊂ {1, . . . , n} is defined as p(B) := i∈B p(i) and the ergodic flow out of B is defined by F (B) := i∈B,j∈B c p(i)P i,j . The conductance of the pair (P, p) is:

Then the following result holds (see [START_REF] Jerrum | Approximating the permanent[END_REF] and [8, Theorem 4.3]): Proposition 6.6.

ϕ(P) 2 2 ǫ(P) 2ϕ(P). Now, assume that n 1/d ∈ N is even and construct a finite graph with n nodes representing a Euclidean grid of the unit hypercube of dimension d. Assume that the vertices of the graph at one grid point are the 2d nearest nodes, with periodic boundary conditions (the graph can be seen as a d-dimensional torus). Assume that the transition probability is uniform over the neighbors, thus the stationary distribution is also the uniform one. This graph is depicted in Fig. 6.1[Left], with d = 2.

Let B be the halved graph defined by the hyperplane parallel to an axis of the grid and including its center, so that p(B) = 1/2. An illustration is given in 2D in Fig. 6.1 [Right]. Since we assumed periodic boundary conditions, the number of nodes belonging to B and having a neighbor in B c is 2n (d-1)/d . Each of these nodes have 2d neighbors, but only one belonging to B c . Since the stationary distribution is equal to 1/n on each node, the ergodic flow is 2n

The following proposition is central to obtain the proof: Proposition 6.7. Almost surely, for all ω i in {ω i } 1≤i≤h d :

p ⊗N -a.s. (6.4) The strategy consists in proving that T |ωi (X N , H) tends asymptotically to T (X N , ω i ). The estimation of each term can then be obtained by applying the asymptotic result of Beardwood, Halton and Hammersley [START_REF] Beardwood | The shortest path through many points[END_REF][START_REF] Steele | Subadditive euclidean functionals and nonlinear growth in geometric probability[END_REF]: Combining this result with Eqs. (6.12) and (4.1) yields Proposition 6.7.

Proof. [Proof of Theorem 4.2] Let ε > 0 and h be an integer such that √ dh -d < ε. Then any two points in ω i are at distance less than ε.

Using Theorem 6.7 and the fact that there is a finite number of ω i , almost surely, we get:

Hence, for any N large enough, there is a coupling K of PN and p such that both corresponding random variables are in the same ω i with probability 1ε. Let A ⊆ H be a Borelian. This exactly matches the definition of convergence in the Prokhorov metric, which implies convergence in distribution.