Keywords:

Set of atoms learned from the CamCan database [Taylor b] where

Ciuciu for his continuous support during these three years. I really appreciated his patience, the trust and the freedom he gave me to conduct my work. In addition to the people directly involved in my work, I would like to thank my fellows at NeuroSpin, who contributed making the time spent here memorable. To Caroline Le Ster, I could never thank you enough for spending those periods of doubts with me, you were always here and ready to give me some precious advice, I hope that our friendship will last for years. To Arnaud Gloaguen thank you for your help and for all the useful conversations we had on optimization and proximity operators. To Bruno Pinho Menses thank you for always laughing at my jokes and for the great conversations on picanha. To my former colleagues Lisa Leroi, Raphaël Tomi-Tricot, Emilie Poirion you were models. Not forgetting my formers interns Hamza Cherkaoui, Hanaé Carrié, hoping you will enjoy your research life in Singapore and Benoit Sarthou, I have learned more than I was expecting. I also wanted to thank the aperitivo team: Nicolas Guigui, Arnaud Gloagen and Cathy Philippe, Venice wouldn't have been the same without you. To Léonie Borne, Slim Karkar and Arnaud Gloaguen you were the perfect iii running mates.

I also want to thank all my Parietal colleagues, it has been a great pleasure to work with you, you were a true inspiration to me. To Olivier Grisel, thank you for your kindness and for the great conversations we had, it has been a pleasure having you as professor for Deep Learning classes.

Finally I also want to thank my relatives and friends for supporting me during the PhD thesis and in my life in general. To Maha, Kenza and Zineb thank you for always being here for your kindness and support. To Salma, Oussama and Narjiss thank you for your presence during the defense, it meant a lot to me. To my parents, my brother and my sister this thesis is the result of the many sacrifices you have made and the guidance with which you have provided me for all these years. A very warm thank you to my father-in-law Michel for always supporting me and being the perfect father in law. And of course thank you Pierre, I could never find words to express how grateful I am to have you in my life.

List of Tables

σ electrical conductivity kg -1 • m -2 • s 3 • A 2 (S • m -1 )
ϕ radiofrequency phase degree (°), rad χ magnetic susceptibility -

General mathematical notations

The mathematical notations, in this manuscripts, follows general typographical conventions.

Symbol Description

N

The set of Natural number. This PhD thesis has been actually dedicated to bring the SPARKLING acquisition schemes to the clinical routine. The primary goal of this PhD thesis was to develop an online reconstruction framework in order to provide MR images by the end of acquisition.

To do so, we need first to accelerate reconstruction algorithms but also to interleave image acquisition and reconstruction in order to feed the reconstruction part with new data. The main advantage of the proposed approach is to provide some feedback to the physician on the MR image quality during data acquisition. In this thesis, we will illustrate all the developments on T * 2 imaging at Tesla but the corresponding algorithmic contributions will be soon available on the Tesla scanner at NeuroSpin as SPARKLING has been deployed on this MR system as well.

Contributions

The fist goal of this thesis was to improve the reconstruction accuracy in the high-resolution highly accelerated setupi.e. using multi-channel phased array acquisition and non-Cartesian encoding schemes. The latter means that the data are no longer collected over a grid or lattice of samples but instead that the measurements lie on the continuous Fourier domain, also called k-space in MRI.

While various reliable reconstruction strategies have been successfully developed for

Cartesian acquisitions using different model-based approaches (e.g. parallel imaging, regularization for inverse problems in imaging, etc.) and different kinds of prior knowledge (e.g. sparsity, low-rank structure), non-Cartesian reconstruction is often addressed using a gridding step, which impedes the final image quality. Based on an inverse problem formulation, we take into account non-Cartesian multi-channel data with and propose a dedicated reconstruction algorithm that leverages structured sparsity information across the multiple channels. We also get rid of the complicated and painful calibration step of the sensitivity maps associated with the multiple channels, which makes our CS reconstruction algorithm a single-step approach. This approach remains in the convex optimization setting which guarantees convergence to a stable solution. Our results demonstrate an enhanced image quality achieved in a more time-efficient manner.

While the scan time has been drastically reduced using non-Cartesian trajectories, the application of compressed sensing to daily-routine exams remains limited owing to the long iterative reconstruction process. Therefore the second objective of this thesis was to propose a method to speed up this reconstruction. Taking into account the segmented acquisition properties of k-space data, a novel online framework with interleaved data acquisition and image reconstruction has been proposed. Mathematically, we used a minibatch formulation to adapt existing optimization algorithms to this online reconstruction setting. The online reconstruction allows one to give partial feedbacks during the scan time making Compressed Sensing and high-resolution MRI available in clinical routines.

A third significant contribution of this thesis was software development as all pieces of work have been implemented in an open source package called PySAP (Python Sparse data Analaysis Package), which is developed at CEA between the NeuroSpin and CosmoStat teams in the context of the DRF Impulsion COSMIC project led by my main supervisor, Philippe Ciuciu. On top of that, we interfaced the algorithmic developments done in PySAP with the client/slave Gadgetron interface system in order to push the online reconstruction to a remote workstation. Asynchronous communication between the MR system and the computer dedicated to image reconstruction was thus implemented.

Thesis Outline

This thesis is organized as follows:

Chapter : Background in Magnetic Resonance Imaging

In the first chapter, the principles of MR physics that will be needed in the rest of the manuscript will be summarized. This allows to make the document self-contained.

In particular, we discuss the principles of signal detection, sampling and image formation using various encoding schemes as they are at the heart of accelerated acquisition strategies. As the NMR signal is collected on the k-spacei.e. Fourier conjugate space -we briefly describe tools for image reconstruction when the samples do not fall into a Cartesian grid. Furthermore, we present specific reception coils that first allow to increase the signal-to-noise ratio (SNR) and second empower scan acceleration with dedicated acquisition and reconstruction methods, i.e. parallel-imaging. However, the accelerations achieved by parallel imaging are limited and the goal is to go beyond this acceleration.

Compressed Sensing theory is the right framework to drastically shorten scan time and will be naturally discussed in the next chapter.

Chapter : Compressed Sensing in MRI

In this chapter, we briefly introduce Compressed Sensing theory and its application to MRI. Therefore, its three main ingredients are discussed in the following order: first, the way to collect k-space data as incoherent as possible using variable density sampling, first along Cartesian lines and then using more sophisticated non-Cartesian trajectories.

Second, CS theory relies on the sparsity of the sought object (i.e. MR image) in an appropriate domain. Although MR images are not sparse in the image domain, they are at least compressible in a transform domaini.e. they can be efficiently represented or approximated by a few non-zeros coefficients in this domain. We review different classes of sparsifying transforms, whether they are orthogonal/complete, overcomplete embedding redundant information, fixed or data driven. Third, we review the different classes of optimization problems associated with different assumputions on sparsity (analysis vs synthesis prior) and present the efficient optimization algorithms. Last but not least, as this thesis is dedicated to image reconstruction, we will present and choose a suitable image quality metric to assess the performances of the algorithms proposed along this manuscript.

Chapter : Multi-channel MRI reconstruction

Achieving high-resolution MR imaging requires high-input SNR, which is often obtained using specific phased array coil to collect the data. Phased array coils are composed of small elementary coils that are combined together to boost the SNR. However this combination at the reconstruction stage is non-trivial, therefore this chapter will cover the way image reconstruction is performed in the setting of multi-channel array coil acquisition.

In particular, we present a first contribution to this field, namely a simple yet efficient two-step approach for coil calibration and image reconstruction in the SENSitivity Encoding (SENSE) framework. Also, we will pay attention to the compatibility of the different algorithms to on-line reconstruction constraints, which are at the heart of this PhD thesis.

Chapter : Calibrationless multi-channel reconstruction

Due to the lack of compatibility of two-step SENSE-based approaches with online constraints and the limited applicability of alternative methods (e.g. LORAKS) to non-Cartesian sampling schemes, this chapter will explore a different track, called calibrationless reconstruction methods. To avoid the calibration step, these approaches aim to reconstruct an image per channel. As multiple images of the same object, i.e. organ, are reconstructed, they share information therefore and in order to handle the redundant information shared by the coil images, structured sparsity prior knowledge is imposed.

In this field, one key contribution of this thesis is the application of new convex regularization schemes based on mixed norms (e.g. OSCAR, k-support norm) that promote structured sparsity in a more adapted way than conventional group-LASSO (which requires the explicit knowledge of the group-structure). The set of MR images is thus given by the minimization of a global nonsmooth objective function and primal-dual algorithm is used to perform the optimization. The validation -on magnitude ad phase imagesis conducted on prospectively highly accelerated T * 2 d MR imaging from ex-vivo human and baboon brain with in-plane resolution of µm.

Background in Magnetic Resonance Imaging

M

agnetic resonance imaging (MRI) is a powerful technique to non-invasively probe the human body with a vast variety of contrasts. As the underlying phenomenon of magnetic resonance may be thought of rather complex, the goal of the first part of this chapter is to give a brief overview of the basic physical principles underlying MRI .

Importantly, information contained in the pixels of a given MR image is not directly probed. Instead, the measurements are collected in the k-space, which in its simplest and ideal (i.e. artifact-free) form corresponds to the spatial frequency domain or Fourier space.

In Section . , we will thus review the spatial encoding principle of MR images in k-space.

In a second step we recall the principle behind the recording of the magnetic signal. As these two processes of spatial encoding and data recording are at the heart of the complex nature of MR images, a particular attention will be paid to the demodulation. Next, as MRI is hampered by the long scan time, accelerating techniques allowing to shorten MRI Background in Magnetic Resonance Imaging exams have been developed for the last two decades. It has actually been a major field of research with outcomes nowadays available in clinical routine. Based on the use of multi-channel receiver coils, the parallel imaging method and the acceleration it provides will be described in the last section.

. Nuclear Magnetic Resonance phenomenon

Magnetic resonance imaging principles are based on Nuclear Magnetic Resonance (NMR) phenomenon, that was first described in by Isidor Isaac Rabi [Rabi ], work for which he won the Nobel Pize in . The applications to molecular structure exploration quickly followed with the independent works of Felix Bloch [Bloch ] and Edward

Mills Purcell [Purcell ], for which they were jointly awarded the Physics Nobel Prize in . The principles described by these scientists are briefly explained in Sections . . The NMR phenomenon can be observed in any atom that has a non-zero spin, a quantum-mechanical property conveying an intrinsic form of angular momentum. Atoms composed of an odd number of nucleons all have a non-zero spin. This is the case for the naturally abundant isotopes of hydrogen, fluorine, sodium and phosphorus ( 1 H, 19 F, 23 Na and 31 P, respectively) as well as for rarer isotopes of carbon ( 13 C), sodium ( 15 N) or oxygen ( 17 0).

As human body tissues are mostly made of water and fat compartments, they are extremely rich in hydrogen atoms 1 H, whose nucleus consists of a unique proton with a 1 2 -spin. For this reason, most clinical applications of MRI focus on proton imaging. Nonproton -also called "exotic" or "X"--nuclei spectroscopy and imaging is a wide and rich area of research, and their imaging is a dynamic and emerging field with a promising future for clinical research [Coste ], especially at ultra-high magnetic field ( Tesla).

Indeed, their lower abundance in the human body associated with a low signal-to-noise ratio (SNR) may be compensated by a strong static magnetic field. However, their analysis remains beyond the scope of this thesis, which focuses solely on 1 H imaging.

. . Magnetization

The fundamental concept behind nuclear magnetic resonance imaging is the interaction of a proton spin with an external magnetic field B 0 . Nuclei with a non-zero spin can be characterized by their magnetic moment, and behave like microscopic magnets. When the nuclei is exposed to a static magnetic field B 0 , its moment aligns with the static field, either pointing in the same direction (called parallel) or in the opposite direction (called

-: Classical representation of macroscopic net magnetization M 0 (hollow vector) in an isochromat consisting of several adjacent spins immersed in an external static magnetic field (B 0 , grey vector). It is induced by the slight excess of magnetic moments (black vectors) parallel to B 0 , compared to anti-parallel ones. ω 0 is the Larmor angular frequency at which each spin precesses around B 0 . From [Tomi-Tricot , ch. ].

anti-parallel), and start spinning like a gyroscope. This gyroscopic movement is called precession. The angular frequency of this gyroscopic spinning, also referred as angular precession frequency or Larmor frequency, is described by the Larmor equation given by:

ω 0 = γB 0 ( . )
with γ being a nuclei-specific constant called gyromagnetic ratio. In the case of proton 1 H, the constant is equal to γ H = 267.5 × 10 6 rad/s/T which corresponds to 1 2π γ H = 42.57 MHz/T. At the thermal equilibrium, the spins are distributed over the two existing parallel and anti-parallel states corresponding respectively to low and high-level energy, the proportion of which being given by the Boltzmann statistics. The latter depends on the Boltzmann constant, the temperature, the gyromagnetic ratio and the field strength.

Due to thermal energy the difference between spin populations at low-and high-energy level is extremely small. This difference, called the spin excess, is about over at 1.5 T. Even though the spin excess is a million time smaller than the number of protons, the NMR experiments consider a mesoscopic scale called isochromat, which consists of million of billions of protons. In average this isochromat creates a mesoscopic magnetic moment called net magnetization M also denoted at the equilibrium in a static field B 0 by M 0 . This is particularly convenient as this net magnetization operates under the laws of classical physics.

At thermal equilibrium, the distribution of low and high-energy spins within the isochromat follows a Maxwell-Boltzmann law. If we denote by N the number of spins considered, T the temperature, the reduced Planck constant and k B the Boltzmann constant (see 'Symbols Used' section on page xiii), and if we consider that the energy difference is negligible compared to k B T , then the macroscopic magnetization at equilibrium pictured

Background in Magnetic Resonance Imaging

in Figure . -can be written as:

M 0 ≈ N γ 2 2 4k B T B 0 ( . )
During NMR experiments the observable variable is the net magnetization M 0 . Interestingly, Eq. ( . ) shows that M 0 depends on the strength of B 0 . This constitutes the major reason for the development of ultra-high field MR systems as the stronger B 0 , the larger the magnetization and so the potential NMR signal.

. . Excitation and Relaxation

Even in a macroscopic body, spin excess does not guarantee a detectable signal, the magnetization vector must be tipped away from the static field direction B 0 in order to set it into precession. This step is performed by applying a Radio Frequency (RF) magnetic field, B + 1 generated by specific RF coils, for a short period of time also referred as pulse, the latter step is called "Excitation". In order to interact with the spins the generated B + 1 magnetic field must oscillate at the Larmor frequency. On top of this frequency, the generated field is orthogonal to B 0 which creates a transverse component of the net magnetization as depicted in Figure . -. The magnetization vector can then be decomposed into two components, namely the transverse and longitudinal magnetizations, denoted respectively by M xy and M z and defined as follows:

M = M xy + M z with: M xy =     M x M y 0     and M z =     0 0 M z     ( . )
Once the magnetization has been tilted, the B + 1 field is stopped and the system starts returning to its initial equilibrium state. This phenomena is called "Relaxation", the rate for which the system returns to its equilibrium state defines the relaxation time. On one hand, the longitudinal rate is characterized by the time needed for the M z component to retrieve % of its initial value, this time is called T 1 . Since the T 1 is related to the spin interactions with their surrounding, it is therefore referred to as the spin-lattice relaxation. On the other hand, the transverse decay is characterized by the time T 2 needed for M xy to decay to % of its initial value. The phenomena emanates from the loss of phase coherence between spins and is therefore referred to as the spin-spin relaxation (see Figure . -). In practice, the B 0 field is not uniform and the transverse decay is faster than expected, hence the accessible time is not directly the T 2 but rather T * 2 that is linked through the following relation:

1 T * 2 = 1 T 2 + 1 T 2i ( . )
where 1/T 2i = γ∆B i is the relaxation rate contribution attributable to field inhomogeneities (∆B i ) across a voxel. Both T 1 and T 2 are intrinsic tissue properties, however only T 1 depends on the B 0 field strength. Figure . -: Longitudinal (grey) and transverse (black) relaxation of magnetization for a given T 1 and T 2 . From [McRobbie ].

While magnetization relaxes to its equilibrium state, it is possible to measure the evolution of its transverse component, using either the same coil that was used for transmission (Tx) of B + 1 , or a dedicated receive coil (Rx). Indeed, the rotation of M xy in the transverse plane induces a measurable voltage in the coil [Hoult ]. The recorded magnetic field -also called B - 1 -generates a signal, called free induction decay (FID) [Hahn ]; as shown in Figure . -, it oscillates at ω 0 , with an exponential attenuation corresponding to T * 2 . Note that the amplitude of FID depends on the initial longitudinal magnetization M 0 that is tipped into the transverse plane: it is therefore related to Eq. ( . ).

The behavior of magnetization exposed to a magnetic field B was formalized by Felix Bloch [Bloch ] and written in the laboratory frame (xyz) as:

dM dt lab = γM × B - 1 T 2 M xy - 1 T 1 (M z -M 0 ) ( . )
As pictured in Figure . -, it is convenient to describe the evolution of magnetization in a frame (xyz) rotating according to the rotation vector Ω:

Ω = 0 0 ω T ( . )
Equation ( . ) then becomes:

dM dt rot = γM × B - Ω γ -     M x /T 2 M y /T 2 (M z -M 0 )/T 1     ( . )
For instance, Figure . -b is obtained, ignoring relaxation, with B = B 0 + B + 1 and ω = ω 0 = ω rf : in (xyz), B 0 disappears from Eq. ( . ), and B + 1 is fixed. While the Bloch equations describe the spins dynamics for a large voxel in a homogeneous sample, they do not describe how signals generated from different voxels are discriminated. Hence, the following section will explain how these signals can be distinguished using spatial encoding.

.

Image formation

Spatial encoding in MRI is not based on a direct measurement of the pixels values. The reason for that is rather simple, the maximum achievable spatial resolution is governed by the Rayleigh criterion. In MRI, the wavelength is in the order of dozens of centimeters making the direct encoding of the body infeasible in practice. It is only in that the works of Paul Lauterbur [Lauterbur ] and Sir Peter Mansfield [Mansfield a] allowed to encode the NMR signal spatially in a clever way. The goal of this section is to provide tools to understand the concepts underlying spatial encoding. . .

Spatial Encoding

The For example, when G x is applied, the magnetic field B(x) will vary with respect to the position according to B(x) = B 0 + G x x. This variation causes the precession angular frequency to vary linearly in space as follows:

ω(x) = γ(B 0 + G x x) ( . )
Since the precession frequency linearly depends on the position, therefore positive positions (with respect to the iso-center located at the origin) will get higher frequency than negative ones.

2D Imaging

In d MRI, the slice to be imaged is first selectively excited with a slice-selecting gradient G Slice (usually in the z direction with G z ). The data corresponding to this slice is then encoded in terms of spatial frequencies along the two other dimensions. Two techniques are used together to fill the d k-space, line by line:

• Applying a frequency-encoding (or readout) gradient G Read changes the precession frequency of isochromats while signal is acquired, thus allowing to discriminate their position in that direction. We collect as such a line in the acquisition space (also called k-space as defined in Section . . ).

• For now on, we have only acquired one line at the center of k-space. To fill the whole k-space we need to encode spatial frequencies along the remaining dimension, therefore to cover multiple lines. We thus apply a third gradient, G Phase , called phase-encoding gradient, orthogonal to G Read . It is applied between excitation and readout and induces a variation of phase along its direction, which will hold during the signal acquisition (readout) as described above.

Background in Magnetic Resonance Imaging

To acquire a full d k-space, it is necessary to repeat the G Phase → G Read pattern as many times as there are lines in k-space, changing the intensity of G Phase every time, so as to explore all needed spatial frequencies in this direction. The relation between a gradient G -over any combination of axes -applied for a certain duration t and the corresponding encoded spatial frequency k (in m -1 ) is:

k(t) = γ t 0 G(τ ) dτ ( . )
At the end of the d imaging process, we have encoded k-space data for one slice, and we are able to reconstruct the corresponding image (using a d transform). We will study the reconstruction process in detail in Section . . . In d imaging, one can repeat the process described above several times, selecting different slices, to acquire a three-dimensional Field Of View (FOV).

However, resolution is limited in the third dimension (slice thickness), because of three main factors. First, selecting a thinner slice requires higher intensity gradients, which eventually run into hardware limitations. Secondly, it requires increasingly sharp excitation profiles to avoid overlap between adjacent slices (at a lower slice resolution, we usually leave a gap between slices to prevent this overlap). Last but not least, when a slice is excited, the acquired signal comes from all isochromats embodied in it: a lower signal thus comes from a thinner slice resulting in a loss of image quality. Therefore to overcome those issues D imaging have to be considered.

3D Imaging

All the mentioned factors make high isotropic resolution hard -even impossible -to achieve with d acquisitions. Fortunately, it is possible to excite the whole FOV at once and perform d imaging by acquiring a three-dimensional k-space. To do, so we apply the same steps as described earlier and we add a second phase-encoding gradient, G Part , along which the slice -also called partition -direction is encoded by inducing phase variations through the partition-encoding gradient. The latter is generally played simultaneously with The next section, which explains how the NMR signal is detected and will clarify the origin of Fourier encoding.

. . Reception: Detection, Conversion & Demodulation

This step is at heart of MRI as it justifies the complex nature of MR images. In NMR experiments, we seek to measure the time-varying magnetization of the sample. For this purpose we rely on Faraday's law of electromagnetic induction, where a variation of the flux magnetic field creates an electromotive force which is converted to current in a reception coil .

The signal acquired by the coil is first amplified with a low-noise amplifier (LNA). If we assume that the scanned sample is made up of a single material, i.e. the resonance frequency is unique, and we apply a gradient varying such that the magnetic field varies linearly with the position, then the signal acquired by the coil can be expressed as follows:

s(t) ∝ - d dt sample M(r, t) • B - 1 (r) ( . )
with B - 1 the magnetic field per unit current that would be produced by the coil at a point r, the latter field is an example of the reciprocity principle, and M the magnetization produced by the spin precession and defined as follows [Haacke ]: r) .

M(r, t) =     e -t/T 2 M ⊥ (r) cos(γB 0 t + ϕ 0 (r)) e -t/T 2 M ⊥ (r) sin(γB 0 t + ϕ 0 (r)) e -t/T 1 M z (r, 0) + (1 -e -t/T 1 )M 0     , ( . ) with M xy (r, 0) = M ⊥ (r)e -ıϕ 0 (
We recall that depending on the tissue γ can be either positive or negative and that the phase ϕ 0 and the magnitude M ⊥ are determined by the initial RF pulse conditions.

Hence after interchanging the integral and time derivative operators and neglecting the While it is more common to have dedicated coils for transmission (Tx) (excitation) and reception (Rx), a single coil could ensure both roles, namely excitation and reception (Rx/Tx).

derivative of the e -t/T 1 and e -t/T 2 factors compared to that of the e -ω 0 t factor, the signal gets proportional to:

s(t) ∝ ω 0 Vs e -t/T 2 (r) |M xy (r)| |B xy (r)| sin (Φ M (r, t) -Φ B (r)) dr ( . )
where ω 0 = γB 0 is the Larmor frequency in radians per second;

B xy (r) = B ⊥ (r)e -ıΦ B (r)
is the transverse component of the receive coil B - 1 field; ; B ⊥ (r) is the magnitude and Φ B (r) the phase of B xy in the laboratory frame Φ M (r, t) is the phase of M xy , i.e. its angle with the x-axis.

The spectrum of the signal in Eq. ( . ) induced in the reception coil includes two frequency bands around ω 0 . For proton imaging (i.e. 1 H), these frequencies are located around and MHz at 3 T and 7 T, respectively. As working with such high frequency level can be problematic, the first step in the reception chain consists in removing this fast oscillation -the demodulation step -by multiplying the signal with a cosinusoid and a sinusoid whose frequency is around the transmit one.

. . Reconstruction

Fourier Transform

After the demodulation step the sampled signal s dem (t) is recombined into a complex expression that can be read as follows:

. . Image formation

s dem (t) ∝ sample ρ(r)e -ır•γ t 0 G(τ )dτ dr ( . )
with ρ being the spin density proportional to ω 0 , the transverse magnetization M xy and the transverse receptive field B - 1 ⊥ . The integral γ t 0 G(τ ) dτ is often denoted by:

k(t) = γ t 0 G(τ ) dτ ( . )
Eq. ( .) is nothing but the Fourier adjoint variable of the spatial location r. By combining Eqs. ( . ) and ( . ), we get

s dem (k(t)) ∝ sample ρ(r)e -ır•k(t) dr = F[ρ(r)] ( . )
which tells us that the k-space measurements actually correspond to the spatial Fourier transform of the spin density ρ(r). The acquisition space is thus called the k-space. Hence the gradients G x , G y and G z are driving the frequency samples to be collected in the k-space. The simplest way to acquire an image in MRI is to design gradient such that the k-space samples are falling into a Cartesian grid (as depicted in Section . . ). The image is then reconstructed by applying an inverse Fourier transform that can be computed using the Fast Fourier Transform (FFT) algorithm originally proposed by [Cooley ].

Non-Cartesian Fourier Transform

We have previously seen that the gradients are driving the k-space trajectory. Although it is more common to fill a Cartesian grid and then to apply an inverse Fourier Transform to recover the image, in some cases the k-space trajectory leaves the grid. In such scenario, the sampling is called non-Cartesian and there are several possibilities to deal with non-Cartesian acquisitions. First, the data can be interpolated and resampled on a

Cartesian grid before applying a conventional inverse FFT. This technique is referred to as "gridding" [Pauly ], and it is a much faster than computing the Discrete Fourier Transform (DFT). Alternatively, one can deal with continuous measurements out of the grid using the Nonuniform or Nonequispaced FFT (NFFT) [Fessler , Keiner ].

Hereafter, we present first two categories of gridding methods and then provide details about the NFFT.

Grid-driven.

The idea behind grid-driven approach is to estimate the value at each grid point based on the neighboring data (i.e. k-space samples). Figure . -illustrates the principle of grid-driven methods. Although they are easy to implement, they are suboptimal as they do not exploit the whole input data to perform interpolation, especially in heavily sampled k-space regions, loosing part of the available SNR in those area. To overcome this issue, the target Cartesian grid can be oversampled leading to a finer grid, and then the data is interpolated with a simple kernel such as a bilinear interpolator. Using the finer grid, grid-driven interpolation can yield high-fidelity reconstructions especially when the sampling density is varying or when the non-Cartesian samples do not fall far from the finer Cartesian grid.

In an ideal artifact-free acquisition scenario. Figure . -: (a) Data samples (o) lie on diameters in k-space, each data point is considered convolved with a small kernel and the value of that convolution is added to the adjacent kspace grid point (+). (b) Each sampled point (dashed line) is convolved with a gridding kernel and that convolution is evaluated at the adjacent grid points.

In practice, this approach is seldom if ever used, as the MRI community converges on the use of data-driven interpolations, which are discussed next.

Data-driven.

The idea behind data-driven approach is to add the contribution of each data point to the neighboring grid points. Figure . -illustrates the data-driven gridding operation. Consequently, each input sample is considered convolved with a gridding kernel [Jackson , O'sullivan ]. The latter is chosen wide enough to cover the neighboring grid points. In this way, each data sample is distributed over adjacent grid locations. The striking difference between this approach and the grid-driven one lies on the use of the whole data set in data-driven methods.

However to mitigate the fact that some k-space region (e.g. the center) might be over

:= x j ∈ T d : j = 0, . . . , M -1 is drawn from the d-dimensional taurus T d ∼ = -1 2 , 1 2 d
with the number of nodes equal to |χ| = M . For each dimension t = 0, . . . , d -1, the bandwidth N t ∈ 2N is defined as a fixed even number. Let us define the multi-index set I N which is a representation of all possible frequencies in a transform:

I N := Z d ∩ d-1 t=0 - N t 2 , N t 2 = k = (k t ) t=0,...,d-1 ∈ Z d : - N t 2 ≤ k t < N t 2 , t = 0, . . . , d -1
Given Fourier coefficients fk ∈ C, k ∈ I N as input, the nonequispaced Discrete Fourier Transform (NDFT) is defined as the evaluation of the corresponding trigonometric polynomial f ∈ T N at the set of M arbitrary nodes χ, i.e., the calculation of the sums:

f j = k∈I N f k e -2ıπk•x j , (∀j = 0, . . . , M -1) ( . )
The NDFT can also be represented as a matrix-vector product defined as follows:

f = A f ( . )
with the columnwise vectors f := (f j ) j=0,...,M -1 , f := f k k∈I N , and the nonequispaced

Fourier matrix A := e -2πikx j k∈I N j=0,...,M -1
, typically A is not square. The adjoint operator (or matrix) is then defined by:

h = A H f ( . )
which is equivalent to the sums:

∀k ∈ I N , h k = M -1 j=0 f j e 2ıπk•x j . ( . )
The NFFT C-library is a fast approximation algorithm that computes the sums in

Nyquist criterion in k-space sampling

In the field of digital signal processing, the sampling theorem is fundamental as it establishes a sufficient condition on the sampling rate that allows a discrete sequence of samples to capture all information contained in a band-limited signal.

K max x = 1 2∆ x and K max y = 1 2∆ y ( . b)
where ∆k x (resp. ∆k y ) is the frequency interval in the x (resp. y) k-space direction, and

K max x (resp. K max y
) is the maximum frequency along the x (resp. y) k-space direction.

Due to the short lifespan of the NMR signal (few tens of ms) and the requirements in ( . ) related to the sampling theorem, the acquisition time gets often too long in MRI.

Consequently, reducing this time has become a major issue especially for high resolution imaging, as the latter means larger K max x and K max y . Therefore, over the last two decades, many contributions have tried to dramatically reduce the scan time. In the next section, we will review one of the most widely used techniques to shorten the acquisition time, which is called parallel imaging or parallel-MRI . The latter is based on phased array coils. Consequently, we will first explain the physics principles underlying the way these coils work and second we will highlight their pros and cons. Last, we will summarize the acceleration and reconstruction methods in the parallel imaging field.

which is accelerated for computing FFT for integers that do not read as 2 n with n ∈ N

Parallel-MRI might be confusing since it can refer to both, the use of multi-channel receiver coil to acquire the signal and the acceleration.

.

Multi-channel receiver coil and Parallel Imaging

In Section . . , we only considered a coil with a quadrature design to record the NMR signal. Since the magnitude of the electric voltage induced in the receiver coil is pretty small, namely in the order of millivolts, the most important properties of an RF receiver coil are the maximal achievable signal-to noise ratio (SNR) and a large coverage of the RF response over the imaged volume. Hereafter, we discuss both aspects.

. . Signal-to-Noise Ratio increase

Many factors determine the SNR available in a nuclear magnetic resonance experiment.

In the early days of MRI, two somewhat parallel paths towards more efficient MRI experiments were drawn: the first consisted of improving the gradient technology and pulse sequence design to increase spatial resolution and decrease imaging time, and the second relied on the development of the hardware RF coil technology. Regardless of the field strength, the key requirement of any receiver coil is achieving the maximum SNR in order to obtain the best possible image quality. This section provides tools to understand how phased array coils maximize the input-SNR.

Noise origin

As it is not straightforward to exhibit a simple relation between the sample scanned and the noise associated with the RF coil, the discussion hereafter will suppose perfectly homogenized fields. The expression for the noise in the coil is based upon thermodynamic principle and is given by [Haacke ]:

noise ∝ 4kT coil ∆f R ef f ( . )
where T coil refers to the temperature of the coil, ∆f the bandwidth used for the experiment, and R ef f is the effective resistance, which includes contributions from the coil R coil , the electronics R electronics and the sample being imaged R sample [Haacke ], such that:

R ef f = R coil + R electronics + R sample ( . )
Although R coil ∝ V coil , with V coil being the coil volume, the latter could be neglected owing to the recent improvement in electronics. Hence, it is seen that

R ef f ≈ R sample ,
which is proportional to the volume of the region of the body (V sensitivity ) from which the signal arises. From the signal expression described in Eq. ( . ) and the noise definition given in Eq. ( . ), the SNR associated with a particular coil reads as follows :

SN R ∝ ω 7/4 B ⊥ T coil ∆f V sensitivity ( . )
From Eq. ( . ) one can deduce that optimal coils are the ones that only cover the imaging region of interest.

Background in Magnetic Resonance Imaging

Phased array coil

In MRI, a phased array coil generally refers to a set of receiver coils whose signals are combined to obtain a uniform image over a region larger than any individual coil could cover while taking advantage of the high SNR available from the smaller individual coils, as shown in Eq. ( . ). Figure . -shows the single channel birdcage coil available at

Neurospin and an example of multi-channel array head coil, where the small elements provide a higher input-SNR compared to the single channel coil.

In order to obtain optimal SNR from phased array coil, it is necessary to make sure that the noise from coil to coil is largely uncorrelated. On top of dedicated electronic circuit for each coil, minimal electromagnetic interaction should exists between the coils to make sure that the noise samples remain uncorrelated. The SNR is not the only parameter to take into account for the design of phased array coils, the penetration depth is as important as the SNR. The penetration depth refers to the depth at which the coils sensitivity drops to % of that at the coil center. As a rule of thumb, the penetration depth of a circular, sample noise-dominated loop-coil is approximately equal to its diameter [Haase ].

Figure . -illustrates the trade-off between penetration depth and sensitivity. While a large number of channels, e.g. , with small diameter provides a significant SNR gain on the cortical surface, the central SNR shows % SNR loss compared to the -channel coil for root-sum-of-square combination (cf Section . . for more details on coil combination).

Channel combination: retrieving magnitude and phase information

When multiple receivers are used, a complex-valued image is reconstructed from each receiver coil using the inverse FFT. However the individual coil images are not really useful as each image is sensitive to a particular brain region. Instead, one prefers to combine all of them in order to cover the whole FOV. In the next part, we address the problem of coil combination for both the phase and magnitude images. The quest of optimal combination of coil specific images (I k ) 1≤k≤L has been a longstanding issue. It can be optimally tackled by injecting the knowledge of the receiver field (B - 1 ) 1≤ ≤L associated with the th -coil in the following relation:

I(x, y) = L k, =1 B - 1 * (x, y) Σ -1 k, I k (x, y) L k, =1 B - 1 * (x, y) Σ -1 k, B - 1 (x, y) ( . )
where Σ is the noise correlation matrix associated with the coils, and each entry Σ k, is the correlation between the k th and the th coil given by:

Σ k, = E [N * N k ] -E [N * ] E [N k ] ( . )
with E [•] being the expectation operator, taken over the noise measurement N during a short acquisition without gradient or RF pulse and for which the configuration of the readout bandwidth must be the same as the one used for the "real" scan. The receiver field B - 1 is also known as the sensitivity profile of the coil. Figure . -illustrates the sensitivity profile of a -channel reception coil. Since the coil sensitivity profiles depend on the loading of the receiver coils (cf Eq. ( . )), they need to be estimated for each scan (Section . will discuss methods that estimate those maps). In many situations, to avoid the estimation of the coil sensitivity profile, and potential errors due to this estimation, a method that combines the data without detailed knowledge of the receiver fields and that preserves,at the same time, a high SNR for the phased array, is really appealing. In [Roemer ], the authors have proposed an approximation for the coil combination in Eq. ( .) by proposing to end up with a single image I(x, y) computed as the sum-of squares (SOS) of coil-specific images I k (x, y), ∀k = 1, . . . , L:

Background in Magnetic Resonance Imaging

I(x, y) = L =1 I * (x, y)I (x, y) . ( . )
Due to the inconvenience of estimating the receiver fields B - 1 , the sum-of-squares approximation is more commonly used today. the authors of [Roemer ] have also shown that the SNR loss resulting from this SOS approximation usually ranges around a few percent only. Moreover, it should be noted that the optimal coil combination in (Eq. ( . ))

provides a complex-valued MR image, whereas the SOS method only yields a magnitude image. This makes its application problematic with methods that requires phase information, in applications where phase brings relevant knowledge such as flow/velocity, temperature mapping or shimming. Combining the phase information is therefore an critical aspect for the mentioned techniques.

Based on the slowly varying phase assumption, [Parker ] have proposed to create a virtual coil reference to which the measurements of each receiver coil are aligned and then combined to obtain an optimal phase distribution estimate. The method is called virtual coil reconstruction and defines a virtual coil as follows: Figure . -: Top: Coil combination methods on an ex-vivo human brain acquired at 7 T with a -channel head coil: (a) magnitude image associated with the complex summation, (b) approximate combination proposed by [Roemer ] using the square root of the SOS. Bottom: (c) phase image associated with the complex summation, (d) phase obtained with the virtual coil method as proposed by [Parker ].

I virtual = L =1 w I = L =1 |w | e ıϕ ref |I
averaged phase can be considered instead:

ϕ ref = -angle x,y p (x, y) x,y |p (x, y)| ( . )
where x, y are D positions over the entire image. This can be easily extended to D imaging of course. Using the virtual coil image, coil-specific phase images are subtracted to the virtual phase and low-pass filtered to remove the noise phase. Hence, the virtual coil is used as phase reference.

Figure . -shows the coil-combination of an ex-vivo human brain acquired at 7 T with a -channel head coil, using either a simple complex sum of the coil-specific images or the SOS [Roemer ] to combine the magnitude images and the virtual coil [Parker ] to retrieve phase information. The magnitude images depicted in Figures . -a and. b demonstrate the good performances of the SOS method while the simple complexvalued sum yields a strong variation of the SNR across pixels, hence poorer performances.

Also, the virtual coil method shows the benefits for phase combination as the resulting phase image gets smoother. It should be also noted that the phase is wrapped therefore variations form -π to π are actually smooth.

Background in Magnetic Resonance Imaging

. . Parallel MRI: first MR accelerations

In the early 's, it has been noticed that phased array coils might not only be dedicated to improve the SNR, but also to accelerate the acquisition [Pruessmann ], yielding to the known parallel imaging techniques. The basis of these methods is that the scan time is proportional to the number of phase encoding lines in a Cartesian acquisition. Increasing the distance between phase-encoding lines in k-space by a factor of R, while keeping the maximal spatial resolution fixed, reduces the scan time by the same factor. However increasing the distance between phase-encoding lines also decreases the FOV (as shown in Eq. ( . )), leading to aliasing artifacts. In parallel imaging, dedicated reconstruction techniques have been proposed to prevent aliasing artifact to occur. These methods are based on the spatial dependence of the B - 1 field on the receiver coils (also called sensitivity maps). Generalised autocalibrating partially parallel acquisition (GRAPPA) [Griswold

, Blaimer ] and SENsitivity Encoding (SENSE) [Pruessmann ] methods are the most popular parallel imaging reconstruction techniques that will be described hereafter . In short, SENSE methods perform image reconstruction in the image space whereas GRAPPA methods operate in the native k-space domain.

Sensitivity Encoding methods for fast MRI

With parallel imaging the scan time is reduced by a factor of R if we collect one line every R lines in the phase-encoded direction. The FOV of the reconstructed image gets therefore reduced by the same factor. Since Nyquist criterion is no longer satisfied, some pixels will be aliased. For each spatial location of reduced FOV images, the pixel value

Other complementary acceleration methods exist such as Controlled Aliasing In Parallel Imaging Results In Higher Acceleration (CAIPIRINHA) [Breuer ], however those are only specific to D acquisitions.

is actually the superposition of the original object and R shifted replicates, the distance between two replicates being equal to FOV y /R if we assume k y to be the phase encoding direction.

The key idea to separate the signals lies in the fact that in each single-coil image I , the signal superposition occurs with different weights according to local coil sensitivities S , with ∈ {1, . . . , L}. The forward model in the image domain relating the unfolded image I to the folded ones (I ) L =1 thus reads:

∀ = 1, . . . , L, I = S I (+N ) ( . )
where N stands for some additional noise related to receiver coil . As matrices S are rectangular of dimensions N x N y /R × N x N y , the SENSE solution I sense , introduced by [Pruessmann ], is thus given by the least squares solution:

I sense = L =1 S H S -1 S H I . ( . )
Of course, it is worth mentioning that this method relies on the extraction and prior knowledge of sensitivity profiles (S ) L =1 . This topic will be discussed in Section . .

Generalized Auto-calibrating Partially Parallel Acquisition

In order to avoid sensitivity profile estimation, methods like GRAPPA [Griswold ] formulate the parallel imaging reconstruction as an interpolation problem in k-space. This method relies on the acquisition of auto-calibration central lines (ACL) that are used to calibrate an interpolation kernel in k-space. It is worth mentioning that the acquisition of those extra-lines slows down the overall scan time. To shortly describe the GRAPPA method, let us define a set of block operators.

The operator R r represents the selection a block of k-space samples in the neighborhood of the position r over all coils. GRAPPA was originally developed for Cartesian under-sampling, so the k-space neighborhood was defined over the grid. Also, the operator P r represents the local sampling pattern that selects the collected k-space samples from a given block. Let y be the multi-coil k-space samples concatenated into a columnwise vector in which unseen data are zero filled, and let y be the reconstructed k-space. The product P r R r y thus defines a vector containing only the the measurements located in the k-space neighborhood around position r. The recovery of the missing values in the th -coil at position r is given by:

y (r) = (P r R r y) g r, ( . )
where g r, is the GRAPPA kernel, i.e. a set of reconstruction weights. The full grid is reconstructed by the evaluation of Eq. ( . ) for all coils at each missing k-space position, while the GRAPPA kernel g r, are obtained by solving the same relation at different position where y i is known. Typically this is done from the ACL region in the center of k-space. Comparison SENSE vs GRAPPA up MR acquisition even further at the cost of more complex image reconstruction but with limited impact on image quality. T his section introduces some theoretical and practical concepts underlying Compressed sensing (CS) for MRI.

.

Compressed Sensing theory

In this thesis, we consider the case of discrete images, i.e. defined on a grid with finite number of pixels denoted by x ∈ C n with n being the total number of pixels. Let y ∈ C m be the observed signal and E ∈ C m×n the encoding matrix which encodes the physical principles underlying MRI acquisition, either Fourier transform in the easiest case treated in this thesis or more complex models including off-resonance effects, chemical shift, etc, see for instance [Doneva ] for an extensive review on computational models in MRI.

The relation between x, y and E -also called the forward model -reads as follows:

y = Ex ( . )
If E is invertible (in particular for m = n) then the reconstruction is performed through the inversion of matrix E. However, the Compressed Sensing (CS) theory supposes that m n. This assumption will hold in the rest of this thesis.

. . Sparsity and incoherence

In the CS context (i.e. m n) the number of solutions of Eq. ( .) is infinite and the problem is considered ill-posed. To solve Eq. ( . ), prior information must be injected to guarantee the uniqueness of the solution. CS theory works with sparsity assumptions on x. We recall here what this means. We define a s-sparse vector on C n as a vector with only s n non-zero entries. If the image x is not naturally sparse then one can find an orthonormal transform Ψ :

C n → C p which transforms x → z = Ψx such that
z is s-sparse (or at least compressible), and Ψ H being its adjoint i.e. transpose conjugate which in the case of orthonormal transforms corresponds to its inverse. The choice of the sparsifying transform will be discussed in Section . . . Then the following optimization tells us how to reconstruct the image:

z = arg min z∈C p z 0 subject to y = Az. ( . )
where 0 stands for the pseudo-norm corresponding to the number of non-zeros coefficients and A = EΨ H . Once a solution z of Eq. ( .) is found the image is reconstructed as follows: x = Ψ H z. Due to the combinatorial properties of the 0 pseudo-norm and its non-convexity, problem Eq. ( .) is NP-complete meaning that, any solution of Eq. ( .) can be tested quickly in polynomial time. However the existence of an algorithm that solves problem ( . ) in a efficient manner (in polynomial time too) has not been demon- for all s-sparse z vector, it exists δ s ≥ 0 such that:

(

1 -δ s ) z 2 ≤ Az 2 ≤ (1 + δ s ) z 2 . ( . )
Assuming the noiseless case, [Candès ] proved the following results:

Theorem . . . if δ 2s ≤ √ 2 -1, then the solution z of Eq. ( .
) is unique.

However, due to the restricted assumptions of the RIP condition only a small subset of encoding matrices -such as random matrices with identically independently distributed (i.i.d.) entries -fulfill the RIP conditions and the RIP constant δ s cannot be calculated explicitly.

Over the same period of time, [Rauhut , Candès ] have extended the theory to random linear projections from orthogonal bases. They considered a sensing matrix A = (a * i ) 1≤i≤n ∈ C m×n to be constructed by randomly drawing rows of an orthogonal matrix A 0 ∈ C n×n given by:

A 0 =      a * 1 . . . a * n      ( . )
Then, matrix A can be constructed by randomly drawing rows of A 0 as follows:

A =      a * J i . . . a * Jm      ( . )
The"P=NP" problem is one of the seven Millennium Prize Problems, for which the Clay Mathematics Institute offered a $ million prize for the first correct solution.

where {J 1 , • • • , J m } are drawn from an independent and identically distributed manner over the set {1, . . . , n}: this means they are copies of a uniform random variable X such that P(X = i) = π i = 1/n, ∀1 ≤ i ≤ n. The coherence κ of matrix A 0 can be defined by:

κ(A 0 ) = n • max 1≤i≤n a i 2 ∞ ( . )
Then the following theorem holds: Theorem . . . Let z be an s-sparse vector. If the number of measurements m satisfies:

m C • s • κ(A 0 ) • log n η ( . )
where C is a universal constant, then z is the unique minimizer of Eq. ( . ) with a probability exceeding 1 -η.

The latter theorem is particularly interesting since it bridges the gap between two fundamental principles underlying CS theory, namely sparsity and incoherence. First, from Theorem . . it can be noticed that the number of measurements required is proportional to the sparsity level s. Second, the other important value is the coherence κ(A 0 ) which varies from to n for orthogonal matrices. The coherence measures how a vector represented in the sparse space spreads in the acquisition space. In the case of discrete Fourier transform F the coherence is minimal since κ(F ) = 1. However, in MRI (as for most applications), the sensing matrix . .

A 0 = F Ψ H is coherent (κ(A 0 ) = O(n)),

Breaking the coherence barrier

The above mentioned result does not assume any structure -apart from sparsity -on the signal of interest. Recovering arbitrary sparse vectors is a very demanding property that precludes the use of CS in many practical settings. If some structure of the sought signal is assumed, then it can be recovered with fewer measurements than those previously uncovered, however dedicated sensing matrices such as variable density sampling should be considered [Chauffert , Adcock ]. In that case, the coherence of the matrix A 0 is defined as follows:

κ(A 0 , π) := max 1≤i≤n a i 2 ∞ π i ( . ) where π = (π 1 , • • • , π n ) denotes the probability distribution on {1, • • • , n}. It has been
shown that the probability density π minimizing κ(A 0 , π) (in Eq. ( .)) and thus that minimizes the number of measurements m in Eq. ( . ) reads as follows:

π i = a i 2 ∞ n k=1 a k 2 ∞ ( . )
The density π defined above, samples more frequently the measurements (i.e. k-space positions in MRI) associated with larger a i 2 ∞ values, which happens in the low frequency regime. Consequently, the higher frequencies are sparsely sampled compared to the lower ones. Thus, variable density sampling (VDS) is achieved by drawing m raws (J 1 , . . . , J m ) from A 0 according to π to build up A. The next theorem thus gives us tighter bounds for exact recovery [Chauffert ]:

Theorem . . . Let z an s-sparse vector. If the number of measurement m satisfies:

m C • s • n k=1 a k 2 ∞ • log n η ( . )
where C is a universal constant, then z is the unique minimizer of Eq. ( . ) with a probability exceeding 1 -η.

This theorem is actually more appealing for MRI as one can show that n k=1 a k assumed. Then vector z can then be recovered by the 1 minimization with probability greater than 1 -ε.

Although the latest works by [Boyer , Adcock ] derived theorems that guarantee exact reconstruction using block sampling sensing matrices in the 1 -minimization and two major properties, i.e. extra and intra block coherence, few major drawbacks can be identified: i) the evaluation of the extra and intra coherence can be hard;

ii) only blocks of parallel lines are considered;

iii) the design of the optimal sampling strategy is still an open issue, especially in non-Cartesian acquisitions; iv) only the 1 minimization with orthonormal basis have been proposed, whereas we could expect an improvement of the reconstruction assuming more sophisticated norms or transforms.

Overall the existing theoretical results only provide some guidelines on the design of sampling schemes (e.g. variable density with a locally uniform coverage) and regarding reconstruction (nonlinear reconstructions promoting sparsity). In the following sections we focus on the application of CS to MRI by first considering the acquisition of k-space following the aforementioned considerations, before handling the reconstruction.

.

MR acquisition consideration

Compressed sensing (CS) theory has been applied to MRI by Lustig [Lustig ] shortly after its invention by Candès group [Candès a] and Donoho [Donoho ]. Since the design of the sensing matrix is at the heart of CS theory, it is crucial to understand the constraints that govern it. Therefore, the first part of this section will be dedicated to the hardware constraints involved on magnetic field gradients for data acquisition. Next, we will review traditional methods to under-sample k-space according to CS recipes, before introducing the SPARKLING algorithm to generate physically plausible and maximally efficient k-space trajectories [Lazarus , Lazarus a].

. . Gradients constraints

The k-space acquisition is governed by different gradients as it has been described in Section . . . Due to the hardware and physiological constraints -such as peripheral nerve stimulation -the coils used for the generation of those magnetic fields are bounded by their maximum current intensity in Eq. ( . a) and their first-order derivative in Eq. ( . b), also called slew-rate which are defined ∀t ∈ [0, T obs ] as follows:

G(t) ≤ G max ( . a) dG(t) dt ≤ S max ( . b)
where T obs is the readout time. The norm used is either the 2 norm or the ∞ depending on the electronic configuration of the coil, i.e. if for each coil the current is generated either by the same generator or not. Since the trajectory is driven by the gradients following the relation defined in Eq. ( .), the k-space constraints are derived ∀t ∈ [0, T obs ] as follows:

dk(t) dt ≤ γ G max ( . a) d 2 k(t) dt 2 ≤ γ S max ( . b)
Eq. ( . a) defines the maximum sampling speed which is given by γ G max while Eq. ( . b)

defines the maximum acceleration bounded by γ S max . Since these constraints are gradientspecific, it might explain why most of the widely used trajectories are based on simple geometrical patterns such as lines (Cartesian or not) or spirals. Also, the more complex the trajectories, the more likely the actual curves played by the scanner may deviate from the originally prescribed ones (see [Vannesjo , Vannesjo ]). In addition, the rapid decay of the MR signal ( ms) usually prevents the measurement of all the needed data at once. For these reasons, k-space trajectories are generally composed of multiple seg- It is worth mentioning that single shot trajectories exist too -such as echo planar imaging (EPI) -however, those trajectories may induce distortions. Therefore their usage is limited to situations where fast scan is mandatory such as in functional or diffusionweighted MRI. Consequently single shot EPI is never used for high-resolution anatomical imaging.

. . State-of-the art on accelerated trajectories

As explained before, owing to the physical constraints over the magnetic field gradients, it is easier to use simple geometrical patterns such as straight lines than more sophisticated curves. For that reason, first acceleration methods reduce the number of shots n s . Hence for D imaging, D acceleration was achieved by subsampling the collected shots along the phase encoding direction, for instance using a D variable density centered around the central line (k y = 0). In D imaging, shots were drawn on a D plane -defined by the phase and partition encoding directions -following a variable density according to the Poisson disk-sampling principle , while the readout (frequency encoding) was performed along the third orthogonal direction. However, the two major drawbacks of those strategy are:

i) the sampling efficiencyi.e. the k-space portion covered per unit of time -is not maximal;

ii) the optimization of the sampling lies on a lower-dimensional subspace, i.e. the D cross-section for D imaging. Hence, variable density is not implemented along the third dimension and consequently redundant information is collected.

Therefore to maximize sampling efficiency, non-Cartesian trajectories need to be considered. D variable density sampling can either be achieved thanks to radial or spiral patterns, with the maximum sampling efficiency achieved for spiral imaging.

The next paragraph summarizes the most popular approach to perform k-space undersampling, originally introduced in [Lustig ] in the MRI field, namely Poisson-disk sampling.

D Poisson disk sampling

Poisson disk under-sampling basically consists first in drawing independent random points from a D sampling density typically with a radially decaying profile in order to perform variable density and promote low frequencies. Second, the k-space measurements are collected over orthogonal lines (readout) whose cross-section is actually given by the drawn iii) it avoids to play with complex sampling trajectories and thus to manipulate nonuniform or non-equispaced Fourier transforms.

However, Poisson-disk sampling is barely tested in D acquisition setups but rather in D simulations, i.e. on retrospective under-sampling scenarios. Instead, its performances are most often demonstrated in d imaging regardless of the MR acquisition specification.

It is actually particular case of SPARKLING trajectories where we no longer activate the projection constraints over the magnetic field gradients.

D variable density lines

d acquisitions can also fulfill the requirements of CS theory, by adopting similar idea to those used in parallel imaging: reducing the number of shotsi.e. TR-therefore considering fewer lines in the phase encoding direction. However, in contrast to parallel MRI the lines are not drawn deterministically but rather randomly with a variable density distribution along the phase encoding direction, i.e. in that case only D VDS is implemented.

An example of the reconstruction using variable density lines is presented

Spiral and Radial trajectories

Before the advent of Compressed Sensing theory and besides parallel-imaging, few trajectories were used to under-sample the k-space and thus accelerate scans. Among all k-space trajectories, radial patterns are the most widely and successfully used in MRI and bridge the gap with tomographic imaging. Radial sampling has a particular appeal for compressed sensing since the k-space center is naturally oversampled compared to higher frequencies. Moreover its repeated sampling of the k-space center results in a signal averaging in the image space making this trajectory robust to motion artifacts. For dynamic imaging purposes, radial trajectories have been stacked together to form a D stack of stars, which is known to be efficient for dynamic application (e.g. cardiac MRI) as it has been proved in [Chandarana ] for free-breathing abdominopelvic exams.

Under-sampled spiral trajectories were also used to accelerate the scan as they present

. . Design of the SPARKLING trajectories

Recently a new method to design sampling patterns, called SPARKLING, has been proposed [Lazarus a]. SPARKLING stands for Spreading Projection Algorithm for Rapid K-space samplING. Its objective is to yield multi-shot trajectories to efficiently undersample the k-space while being compliant with the gradient constraints, implementing d variable density in D imaging, performing locally uniform k-space coverage and eventually manipulating contrast affine constraints by specifying for instance the echo time TE, i.e. the k-space center crossing time.

The SPARKLING method relies on optimization to automatically generate k-space trajectories k(t) under the aforementioned hardware constraints on the gradients as defined in Eq. ( .) -which can be expressed in terms of trajectory constraints as defined in Eq. ( . ) -by minimizing a tailored distance between the sample distribution and any prescribed density π. The objective becomes fitting the density π over the trajector-

In D we get k(t) = [kx(t), ky(t)].
min k∈Qp dist(π, v(k)) = min k∈Qp 1 2 h (v(k) -π) 2 2 ( . )
where h is an interpolation kernel, v(k) is the probability of measure supported by the trajectory k and Q p is the set of admissible curves respecting to the constraints defined in Eq. ( .). An 2 distance has been chosen for the sake of simplicity. The discretization of Eq. ( . ) leads to the following formulation:

min k∈Qp 1 p 2 1≤i,j≤p k[i] -k[j] 2 2 Fr(k) - 1 p p i=1 Ω x -k[i] 2 2 π(x)dx Fa(k) ( . )
The latter formulation can be interpreted as the minimization of a potential energy, composed of a repulsive term F r (to avoid clusters) and an attractive term F a (to concentrate the samples according to the target density π). Due to the non-convexity of Eq. ( . ), this optimization is sensitive to the initialization as there as potentially many local minimizers. The SPARKLING trajectories have also been implemented at 7 T for prospective acquisition within a gradient recalled echo (GRE) Fast Low Angle SHot (FLASH) sequence.

The method has been compared to variable density centered-in/out radial -used as initialization -and spiral [Lee ] imaging. The results of this comparison are presented in Figure . -. They indicate that SPARKLING trajectories maintain high image quality even for large acceleration factors, i.e. high similarity to the fully sampled Cartesian reference . Due to the good performances of SPARKLING trajectories, the experimental setup used in the rest of this thesis will be the radially initialized SPARKLING trajectories.

. Reconstruction

Due to the under-determination of the MR image reconstruction inverse problem from highly under-sampled data, there exists an infinite number of solutions. CS theory has provided the right framework to compute a unique solution to the inverse problem by imposing sparsity in the sought image. Originally, in the noise-free context, the image solution was computed by minimizing the 1 norm subject to the data consistency problem (see Eq. ( .)). As we the data is noisy and due to its robustness to perturbation, the regularized version of reconstruction is often preferred, it reads as follows:

z = argmin z∈C N Ψ 1 2 y -F Ω Ψ s z 2 2 F (z) + λ z 1 R(z) ( . )
where y is a vectorized version of the k-space measurements, z is the sparse representation in the dictionary Ψ s of the image to be estimated, F Ω is the under-sampled Fourier operator over the binary mask or support Ω that picks up the collected measurements: it can be a single FFT for Cartesian under-sampling or the NFFT operator for non-Cartesian acquisitions. Ψ s : C N Ψ → C N is the synthesis version of the sparsifying transform.

Therefore, the problem formulated in Eq. ( .) is also called the synthesis formulation as the estimated image is given by x = Ψ s z. However, one can also introduce an analysis formulation that computes the solution directly in the image space as follows:

x = argmin x∈C N 1 2 y -F Ω x 2 2 F (x) + λ Ψ a x 1 R(x) ( . )
where x corresponds to the image and Ψ a : C N → C N Ψ the analysis version of the sparsifying transform. The difference between the synthesis and the analysis formulations (Eq. ( . ) i.e. vs Eq. ( . )) will be discussed in the following section. Moreover it should be mentioned that their adjoint operator is defined as the conjugate-transpose Ψ H s for the synthesis transform and Ψ H a for the analysis one. Whatever the retained formulation, the reconstruction problem relies on two main ingredients that will be discussed in depth in this part: i) the sparsifying transform Ψ, which should be chosen as incoherent as possible with respect to the acquisition basis;

An extensive comparison including several image contrasts can be found in [Lazarus a].

ii) the optimizer -solver -that affects the convergence speed, hence the computation time.

In what follows, these two points will be discussed in detail while the last subsection will be dedicated to the definition of quantitative image quality metrics. Indeed, as this work focuses on image reconstruction, it gets critical to choose an appropriate quantitative metric to assess the performances of the reconstruction algorithms proposed along this work.

. . Sparse decomposition

One of the key tools in CS theory is the incoherence between the sparse representation and the acquisition basis. However, the natural MR images are not sparse in the canonical basis and need therefore to be sparsified through a decomposition that can either have an analytical fixed formulation -such as Wavelet Transform -or learned from the data itself such as in the Transform or Dictionary Learning (TL vs DL) [Ravishankar ,

Ravishankar

, Wen ] framework. In what follows, we will study the impact of the synthesis or analysis formulation on fixed dictionaries, as well as the choice of complete vs over-complete transforms. In a second step we explore the performances of the DL on the reconstructed images. shearlets [Guo ] or undecimated Bi-orthogonal WT [Starck ] combined with the analysis formulation outperform the standard orthogonal wavelet transform. However owing to the redundant information manipulated in that case -N < N Ψ -optimization is performed in the image space to save memory usage and limit the computational load, especially in the high-resolutions context.

Analytical fixed sparsifying dictionaries

Adaptive data driven representations

While analytical over-complete dictionaries already provide good results, one could expect an improvement on the image quality using learned dictionaries. DL methods [Ravishankar , Ravishankar ] extract and learn atoms which represent the image to be reconstructed with few non-zero coefficients. Two approaches have been developed so far: either the dictionary is learned off-line and then injected as a sparsifying transform or it can be directly inferred in a blind bilinear reconstruction problem where both the reconstruction and the dictionary estimation steps are interleaved [Ravishankar , Ravishankar ].

The two step procedure can be described as follows: first, the dictionary learning step can be formulated as in Eq. ( . a) then the reconstruction step Eq. ( . b) is defined using the previously learned dictionary.

( D, ẑ) = argmin

D∈C N ×nc ,z∈C nc ∀k≤nc, D k 2 =1 1 2 x -Dz 2 2 + λ z 1 ( . a) x = argmin x∈C N 1 2 y -F Ω x 2 2 + λ Dx 1 ( . b)
where n c is the number of components used to represent the vectorized image x ∈ C n [Ravishankar , Mairal ], the learned dictionary D is composed of atoms that are illustrated in Figure . -. However these approaches either require a large database of images for learning atoms or are computationally slow. Another drawback concerns the reconstruction time [Carrié ], since at each iteration, the image needs to be decomposed in the dictionary, hence for non-Cartesian acquisitions like SPARKLINGi.e. using NFFT operator in Eq. ( . b) -and single channel acquisition setting, the reconstruction takes approximately min against min with the wavelet decomposition for N = 512 × 512.

. . Proximal optimization method

In general CS-based reconstruction -either based on the analysis Eq. ( . ) or synthesis Eq. ( . ) model -is composed of a sum of data-consistency term -denoted by F (•) -and regularization -denoted by R(•) -which promotes sparsity through nonsmooth penalty. In this thesis we mainly considered convex methods as it guarantees convergence to the global minimizer of a convex but nonsmooth objective function

(F + R)(•) [Combettes ].
Hereafter we define the notation and the class of algorithms we have used to solve the reconstruction problem.

Definitions

We denote by C N the usual N -dimensional complex Hilbert space. Standard definitions and notations from convex analysis will be assumed [Rockafellar ] . The domain of a function f :

C N → ]-∞, +∞
] is denoted by domf and defined such as domf =

x ∈ C N |f (x) < +∞ . A convex function is defined as follows:

Definition . . . A function f : C N → R is convex if domf is a convex set and if ∀x, y, ∈ domf and θ ∈ R | 0 ≤ θ ≤ 1, we have f (θ x + (1 -θ)y) ≤ θf (x) + (1 -θ)f (y) .
Let Γ 0 (C N ) be the class of lower semi-continuous convex functions from

C N → ]-∞, +∞].
The convex conjugate of the function f denoted by f ∈ Γ 0 (C N ), is defined as follows:

f : C N → ] -∞, +∞ ] u → sup x∈C N { x, u -f (x)} ( . )
The Moreau subdifferential is defined as follows:

Definition . . . The Moreau subdifferential of a function f :

C N → ] -∞, +∞ [ at x ∈ C N is defined as : ∂f : C N → 2 C N x → u ∈ C N | (∀y ∈ C N ) f (y) f (x) + y -x, u Any vector u in ∂f (x) is called subgradient of f at x.
We define the proximity operator as follows:

Definition . . . Let f ∈ Γ 0 (C) N .
For every x ∈ C N , the minimization problem:

min y∈C N f (y) + 1 2 x -y 2 ( . )
admits a unique solution called proximity operator of f denoted by prox f (x) :

C N → C N .
On top of the notations defined page xv.

Compressed Sensing in MRI

Moreover the proximity operator of f can be characterized by the following inclusion:

∀(x, p) ∈ C N × C N p = prox f (x) ⇔ x -p ∈ ∂f (p) ( . )
We also define a β-Lipschitz gradient of a function f as:

∀(x, y) ∈ C N ×N ∇f (x) -∇f (y) 2 ≤ β x -y 2 ( . )
Since the data-fidelity term (F (•) as defined in Eqs. ( . ) and ( .)) is smooth and has a β-Lipschitz gradient and the regularization term (R(•) as defined in Eqs. ( .) and ( . )) is proper (i.e. domR = ∅) closed convex function, the reconstruction problem can be efficiently solved using proximal splitting methods that are described in the next part.

Proximal majorize minimize method: Forward-Backward Splitting method

A classical approach to solve Eq. ( .) is to use the Forward-Backward Splitting algorithm (FBS) [Levitin , Combettes ]. It can be seen from the Taylor expansion of F (y)

and simple linear algebra that the iterations sequence to solve Eq. ( . ) takes the form:

z k+1 := prox γ k R backward step (z k -γ k ∇F (z k )) forward step , with γ k ∈ ] 0, 2/β [ ( . )
with z 0 ∈ C N . In the particular case where the 1 -norm is used for R the algorithm correspond to the Iterative Soft-Thresholding Algorithm (ISTA) [Daubechies ]. Moreover, it is well established that FBS iterates converge as long as γ k ∈ ] 0, 2/β [. However the convergence speed is quite slow (with a rate in the order of O (1/k)). In the literature, two complementary paths were developed to speed up convergence, namely adaptive restart [O'donoghue ] and the use of a momentum strategies [Nesterov ], both being discussed hereafter. In order to speed the convergence, many variants of Eq. ( .)

have been proposed in the literature, most of them rely on the setting of a momentum term [Nesterov ] that can be summarized as follows: Until now, in Eq. ( . ) we have only considered the synthesis formulation of the reconstruction problem as defined in Eq. ( .). The analysis formulation can be handled similarly when Ψ a is an orthogonal sparsifying transform as the proximity operator could be calculated explicitly [Combettes , Table . -x], for more general decompositions this is no longer the case. Hence, the proximity operator has to be computed iteratively which slows down the reconstruction again. In these specific cases, to avoid the iterative calculation of proximity operators, primal-dual methods have been developed.

Algorithm : Accelerated Forward Backward Set z 0 ∈ C N Ψ , w 0 = z 0 , θ 0 = 1, γ k ∈ ] 0, β ] and q ∈ [0, 1]; for k = 0, . . . do θ k+1 solves θ 2 k+1 = (1 -θ k+1 ) θ 2 k + qθ k+1 ; α k+1 = θ k (1 -θ k ) / θ 2 k + θ k+1 ; w k+1 = z k+1 + α k+1 (z k+1 -z k ); z k+1 := prox γ k R (w k+1 -γ k ∇F (w k+1

Proximal optimization: Primal-Dual algorithm

We have previously seen that the use of over-complete dictionaries could noticeably improve image reconstruction accuracy. However, solving it with one of the FBS variants is time consuming. Hopefully, the analysis problem Eq. ( . ) can be solved by minimizing both the primal and dual formulation as suggested by [Komodakis ] . The dual problem can be formulated as follows:

z ∈ argmin z∈C N Ψ F (-Ψ H z) + R (z) ( . )
where F and R are the convex conjugate of F and G, respectively. It should be noted that the primal and dual problems can be reformulated so as to be solved together, hence the Lagrangian saddle point is found [Bauschke , Chap. ]:

( x, z) ∈ argmin x∈C N argmax z∈C N Ψ F (x) -R (z) + Ψx, z ( . )
Moreover if we suppose that : i) Eq. ( . ) admits at least one solution ii) ri(dom R) ∩ Ψ dom F = ∅, with ri the restricted interior;

then the dual problem admits a solution and at convergence, the dual gap is null meaning that:

In this case we only considered the sum of smooth and non-smooth functions however the primal-dual approaches are more general and could be applied to any problem defined as follow:

F (x) + G(x) + R(Ψx),
where F , R and Ψ are defined as before and G is is proper closed convex function.

F ( x) + R(Ψ x) = -F (-Ψ H z) -R ( z) ( . )
The first approach developed to solve the primal-dual problem was the Alternating Direction Method of Multiplier (ADMM) [Komodakis , Boyd , Gabay ].

Although it has been widely used in various fields, the setting of its parameters may be tricky. It should also be mentioned that a variant of ADMM in the dual space has been proposed in [Eckstein ].

Another primal-dual algorithm has been developed simultaneously by [Condat ] and [Vũ ]. The sequence of its iterates for solving Eq. ( .) can be summarized as follows:

Algorithm : Primal-Dual Condat-Vũ Algorithm Set z 0 ∈ C N , z 0 ∈ C N Ψ , (σ, τ ) ∈ ] 0, +∞ [ 2 , and ∀k ∈ N, λ k ∈ ] 0, +∞ [ ; for k = 0, . . . do u k := x k -τ (∇F (x k ) + Ψ H z k + e F,k ); w k := prox σR (z k + σΨ(2u k -x k )) + e R,k ; x k+1 := λ k u k + (1 -λ k ) x k ; z k+1 := λ k w k + (1 -λ k ) z k ;
end A symmetric version of the Condat-Vũ algorithm has also been proposed and consists in switching Lines . and . in Algorithm . Moreover, the convergence proof of the iterates has been established by the following theorem: , Theorem . ] Let γ > 0, τ > 0 and the sequences (λ k ) ∀k∈N , (e F,k ) ∀k∈N and (e R,k ) ∀k∈N be the parameters of Algorithm . Let β > 0 be the Lipschitz constant of the gradient of F , and suppose that the following holds:

Theorem . . . [Condat
i) 1 τ -σ Ψ 2 ≥ β 2 ii) ∀k ∈ N, λ k ∈] 0, δ [ , where we set δ := 2 -β 2 1 τ -σ Ψ 2 -1 ∈ [1, 2] iii) k∈N λ k (δ -λ k ) = +∞ iv) k∈N λ k e F,k < +∞ and k∈N λ k e R,k < +∞ Then there exists a pair ( x, z) ∈ C N × C N Ψ solution

. . Metric for Image comparison

As the quantification of the MR image quality is important for the evaluation of imaging techniques both at the acquisition and reconstruction levels, this part shortly reviews some standard metrics to assess image quality when a ground truth or reference image is available.

Mean Square Error (MSE) and Structural Similarity Index (SSIM) are the most frequently used indices for assessing image quality between a reference and a reconstructed image. The implicit assumption is that they correlate well with the physicians assessment.

In [Zanforlin ], the authors have compared the SSIM score with the mean observation error (MOS). Interestingly they drew a comparison between the similarity score and the MOS and proposed a table that summarizes the correlation between SSIM and the visual perception, see Table . , However the results stated by [Zanforlin ] were considered for natural images and not for medical ones. Recently, [Mason ] has assessed the correlation between various quality metrics and physician' views. To do so, the experimental protocol was defined as follows: nine brain scans selected from hospital images and corrupted versions

Compressed Sensing in MRI to the SSIM score on the regime where it fits the physician perception i.e. SSIM > 0.9.

. . Conclusion

As seen in the previous chapters, MRI is a powerful non-invasive technique that suffers from long scan time. Although methods such as parallel-imaging permit to shorten this time, they acceleration power remains too limited for high-resolution MRI. The very reason for that is twofold: one one hand, under-sampling is performed in a deterministic manner with a target uniform sampling density. On the other hand, most often the reconstruction itself is purely linear. To go one step further, the application of CS theory to MRI has revolutionized the way data can be collected and images are reconstructed. CS provides higher acceleration factors compared to parallel imaging or partial Fourier techniques but at the cost of more computationally demanding reconstruction. However, to become even more efficient regarding the way SNR is preserved, CS has also been combined with multi-channel acquisition over phased array coils. Then a natural question has emerged on the optimal manner to combine CS and multi-channel acquisition for making image reconstruction more efficient and accurate. These issues will be discussed in the next chapters. T his chapter introduces and discusses the MR image reconstruction methods based on multi-channel array coil acquisition. A particular attention will be paid to their compatibility with non-Cartesian sampling schemes.

. Introduction

Compressed Sensing provides a way to drastically reduce the scan time, however the maximum under-sampling factor depends on the available input-SNR and the image resolu-

tion . An empirical study on the maximum under-sampling factor for T * 2 MR images [Lazarus c] illustrates several important results, for a targeted image quality:

• at a given resolution, the under-sampling factor can be increased as the input-SNR improves until a stationary regime is reached for high input-SNR.

• high resolution images offers the possibility to reach high under-sampling factors as long as the input-SNR remains sufficient, while the under-sampling ratio is limited for low-resolution images.

• at each resolution, an optimal point, which combines the minimum requisite input SNR with the maximum under-sampling factor allowing to maintain the targeted image quality (e.g. SSIM score of . ), can be determined.

Figure . -presents those relations in a slightly different way: the D map provides the evolution of the maximum under-sampling factor achievable for a range of image resolution and input-SNR which preserves a fair image quality, corresponding to an SSIM score above

. according to the mean opinion score (Table . ). The highlighted red line corresponds to the available input-SNR on the Magnetom 7 T scanner (Siemens-Healthineers, Erlangen, Germany) using the available single-channel coil.

To achieve high acceleration factors in the high resolution setting, the input SNR has to be boosted. This can be reached using a multi-channel coil array, as demonstrated in [Roemer ] and presented Section . . . However such reception coil will increase the reconstruction complexity and raise questions of Compressed Sensing theory in such a configuration. This chapter will review CS reconstruction methods for multi-channel acquisitions. We will focus on the non-Cartesian case since the latter permits higher acor image dimension as we assume the FOV fixed. .

Parallel MR image reconstruction based on sensitivity maps

In the Compressed Sensing (CS) context, sparse multi-coil image reconstruction has first been applied using a variant of the SENSE formulation (Section . . ) called sparse-SENSE and originally proposed by [Liang ]. The reconstruction problem reads as follows:

x = arg min

x∈C N 1 2 L =1 F Ω S x -y Ω, 2 2 + λ Ψx 1 , λ > 0. ( . )
with Ψ being a sparsifying transform. Note however that a prior work in the parallel imaging context -with deterministic under-sampling along the phase encoding direction -has been proposed by [Chaâri ] and further extended in [Chaâri ].

In the SENSE formulation, the calibration of the spatial sensitivity profiles associated with the multiple channels gets critical as they heavily affect the final reconstructed image. Originally those sensitivity profiles were estimated using a dedicated pre-scan (lowresolution acquisition). However, this method suffers from two limitations: first it increases the acquisition time and second it is sensitive to the patient's motion between the calibration and accelerated scans. To overcome those issues, many methods have been proposed

to estimate the sensitivity profiles. One sensible approach consists of selecting a calibration region in the k-space domain where the signal has been sampled (at least) at the Nyquist rate. As these maps are spatially smooth, one can actually extract the center of k-space (low frequency domain) and then apply the inverse Fourier transform to recover an estimate. In what follows, we describe the different strategies in detail. Then, we propose a simple yet efficient extraction method well suited for non-Cartesian sampling patterns. they could be considered as inherently calibrated as it has been mentioned by [Yeh ].

. . Coil sensitivity extraction

Direct extraction

Using a gridding operation, low-resolution coil-specific images are then reconstructed using a zero-filled inverse Fast Fourier transform and used as sensitivity profiles for non-Cartesian sampling.

In the same context, we have proposed a similar approach to extract the sensitivity information. First, akin to [Yeh ], we selected only a portion θ of the central surface of k-space and set the other values to zero. Second, instead of using a gridding step, the NFFT [Keiner ] adjoint operator F H was applied to the data ỹ Ω, = y LF Ω, in order to get the n × n low-resolution images: 0] ỹ Ω, where LR stands for low resolution. Third, the square root of the Sum of Squares (sSOS) was computed: x LR = L =1 x LR 2 . Finally, the sensitivity maps were given by the pixel-wise ratio between the coil-specific images and the sSOS:

x LR = F H [Ω |θ% ,
s LR ij = diag [S ] ij = x LR ij / x LR ij , ∀ ∈ 1, .
. . , L and ∀i, j ∈ 1, . . . , N . Because of this sSOS operation, our method is less dependent on the threshold θ compared to [Yeh ], who directly exploited the x LR images as sensitivity maps. In order to attenuate the noisy background, one could mask the estimated sensitivity maps by a binary mask computed by thresholding x LR , where the actual threshold value is given by a -cluster k-means algorithm. The binary mask is eventually defined as the largest connected component.

Indirect extraction

While the sensitivity maps estimation in the previous works relied on a direct k-space extraction, [Uecker ] proposed an indirect approach, based on eigenvector decomposition called Eigenvalue Self-consistent Parallel Imaging Reconstruction using iTerative algorithm (ESPIRiT). Because the latter is considered as the state-of-the art in parallel imaging reconstruction we will explain it hereafter. In the original paper, the proposed method rests upon a calibration matrix A to extract the sensitivity information (such as SENSE) but it also enforces self-consistency with respect to the calibrated data (such as GRAPPA). Hence the authors claim to bridge the gap between SENSE and GRAPPA.

Construction

Extraction of the coil sensitivity profiles. With highly correlated rows, matrix

A is supposed to be low-rank, hence only a small set of vectors span the k-space. To analyze the calibration matrix, the authors perform its singular value decomposition, which decomposes A into U , Σ and V H such that A = U ΣV H , where the columns of V are the basis of the rows of A and therefore a basis for the overlapping blocks in the calibration data. The authors also show that V can be split into V and V ⊥ that respectively span for the row-space and the null-space of A. Similarly to SENSE, the AC stands for Auto Calibrated. authors propose to explicitly identify the sub-space spanned by the coil sensitivities from the calibration matrix. Given the operator R r , which extracts and flattens the patches around the r th k-space position, and ŶΩ = y Ω,1 , . . . , y Ω,L the unknown k-space data, one can re-write the null space condition thanks to the decomposition of subspace V :

V V H + V ⊥ V H ⊥ = Id.
The null space condition reads as follows:

r R H r V ⊥ V H ⊥ R r Ŷ = 0 r R H r R r -1 r R H r V V H R H r W Ŷ = Ŷ
Assume that image x is not null and satisfies the forward model, i.e. FSx = Ŷ , it follows that the sensitivity matrices are represented by the eigenvectors of F H WF associated with the eigenvalue 1. Those eigenvectors could be used to compute the sensitivity profiles, however as reported by the authors this leads to a high computation cost.

. . Blind bi-linear reconstruction

Although for "calibrated"-Cartesian sampling schemes the extraction of sensitivity matrices can be efficiently performed using one of the above mentioned methods, their adaptation to non-Cartesian, random or uncalibrated Cartesian sampling patterns is far less straightforward. In these cases, sensitivity matrices cannot be extracted prior to reconstruction.

However to extend the SENSE formulation, a new category of methods, which aims at recovering both the coil sensitivity profiles and the image solution, has emerged. Such methods are often referred to as blind bi-linear inverse approaches.

Hereafter we define the general regularized blind bi-linear reconstruction problem as follows:

x, S 1 , . . . , S L = argmin

x∈C N ∀ ∈{1,...,L}, S ∈C N L =1 1 2 F Ω S x -y 2 2 + R sens (S ) + R im (x) ( . )
with R sens being the penalty term on the coil sensitivity profiles (S ) , which promotes their spatial smoothness, while R im encodes the regularization term on the image x, which can enforce sparsity in the CS framework.

The regularized blind bi-linear reconstruction was first introduced by [Uecker ].

They originally proposed a Tikhonov regularization [Tikhonov , Tikhonov ] as a prior on x and the Sobolev regularization [Sobolev ] for (S ) ∈{1,...,L} . The Sobolev penalty term is a special case of Tikhonov regularization since it penalizes the weighted 2 -norm of the Fourier coefficients instead of the 2 -norm in the original domain (e.g. image space). The regularization terms proposed by [Uecker ] were defined as follows:

∀x ∈ C N , R im (x) = α 2 x 2 2 ( . a) ∀ ∈ {1, . . . , L}, S ∈ C N ×N , R sens (S ) = β 2 1 + k 2 2 γ/2 FS 2 2 ( . b)
with k 2 2 being the distance to the origin of k-space, F the fully sampled Fourier transform and α, β and γ (the polynomial degree) being positive hyper-parameters that need to be set.

In with ∇x being the spatial gradient of x i.e. (∇x) i,j = (x i+1,j -x i,j , x i,j+1 -x i,j ), hence ∇v is linked to the second order derivative of x. While the overall criterion in Eq. ( .)

is not jointly convex in x, (S ) , the authors of both works [Uecker , Knoll ] initialize x and (S ) to zero-valued vector and matrices and claim that the method was not sensitive to the initialization.

Besides the

∀x ∈ C N , R im (x) = λ Ψx 1 + α ∇x 1 ( . a) ∀ ∈ {1, . . . , L}, S ∈ C N ×N , R sens (S ) = β ∇S 1 ( . b) Contrary to [Knoll , Uecker ],
[She ] set up the image x with zeros and the sensitivity profiles (S ) with low resolution coil images. They also studied the impact of hyper-parameters setting and showed the benefits of 1 -regularization.

One of the main limitation of blind-bi-linear methods concerns their computational cost. For example, in the benchmark presented in [She ], the authors considered D

Cartesian data collected over a eight-channel coil and reconstructed with a target resolution of 256 × 256, for a FOV of 20.4 × 20.4cm 2 . In this scenario, the computation time for blind bi-linear methods (i.e. J-SENSE, Sparse-BLIP and IRGN-TV) was between . and times slower than Sparse-SENSE (see Table . ). This constitutes the major drawback of blind bi-linear methods. In highly accelerated non-Cartesian sampling schemes, we can anticipate even more demanding computational load and memory usage as: i) the target base resolution is often larger than the one proposed by [Shin ] for the same FOV;

ii) the FFT is replaced by the NUFFT wich has a higher numerical complexity (Section . . );

iii) a larger number of receivers is used in multi-channel coil acquisition to boost the SNR. While most recent blind methods rely on the regularization of sensitivity profiles, their dependency on the initial guess, the hard setting of hyper-parameters and the slow conver- ii) K-space sampling is locally uniform and the parallel MRI reconstruction problem is in fact over-determined due to the redundancy between channels.

In particular, [Zhang ] were the first to exploit Assumption i) and proposed a reconstruction method called Parallel Reconstruction Using Null-Operations (PRUNO) that exploits a k-space low-rank constraint. We explain it in detail hereafter.

Let us consider a Cartesian acquisition scenario. If we denote by y a the k-space measurements concatenated over the channels and y m the missing ones over all channels, then y = y a y m corresponds to the concatenation of the measured and missing samples.

The missing k-space samples are supposed to be estimated using an interpolation kernel C that is calibrated from the calibration region in k-space (usually around its center).

This convolution can be written as follows:

y m = Cy a ( . ) Id -C y m y a = 0 ( . )
The optimization algorithm alternates between a conjugate gradient step to solve Eq. ( .)

and
subject to W = H H A, ΩW = Y ( . )
with W the estimated k-space and Ω the sampling mask, hence, ΩW = Y is a consistency term with respect to the collected k-space measurements. To solve Eq. ( .), the authors rely on a projected gradient descent algorithm. Recently, [Lobos ] has studied the optimal window shape (green frame in Figure . -) and proved that circular shapes maximize the reconstruction quality. Although these methods can be easily applied to any

Cartesian sampling (phase-encoded variable density sampling or D Cartesian Poisson disk sampling etc...), their application to non-Cartesian readouts is less trivial, mainly because of the H operator, which is more complex to calculate in that case. This limitation could be relaxed using a gridding procedure such as GROG [Seiberlich ], however the latter may dramatically impair the reconstruction.

The two methods are identical, yet they have been derived independently and are based on different theoretical assumptions. While SAKE [Shin ] is based on multi-channel k-space relationships, p-LORAKS [Haldar ] was originally motivated by the k-space support modeling and phase assumptions of the images.

. . Domain based methods

While k-space based methods are efficient for calibrationless image reconstruction from Cartesian data, they are prone to artifacts for non-Cartesian readouts as they require a gridding step. Hence, for non-Cartesian acquisitions an inverse problem formulation is beneficial for maximizing image quality. However as seen in Section . . , those reconstruction methods require the explicit knowledge of coil sensitivity profiles. In some cases, the extraction of the latter can be either computationally demanding, unstable or inaccurate (depending on the initialization), leading to unexpected artifacts on the magnitude or phase images. Some alternatives have proposed to reconstruct channel-by-channel images using an inverse problem formulation, which basically consists of minimizing a variational form that regularize the reconstruction of multiple images through a joint sparsity penalization over the channels. Consequently, we focus here on methods that formulate the problem as follows: non-overlapping patches), however to avoid the well-know patching artifact this work can be easily extended to overlapping patches where s < p.

X = arg min X∈C N ×L 1 2 L =1 F Ω x -y Ω, 2 2 + R(Z) . ( . ) with X = x 1 , . . . , x L , Z = Ψx 1 , . . . Ψx L and
This method has two major drawbacks: the large number of hyper-parameters to be tuned (λ, s and p) and the computational load that depends on the stride setting. While small strides prevent patching artifacts they increase the computational cost.

The group-sparsity regularization has also been used for MR reconstruction by [Majumdar a]. In the latter work, the authors have used a constraint optimization method leading to solving the following problem:

X = argmin X∈C N ×L ΨX 2,1 , such that F Ω X -Y 2 ≤ ε
Using similar ideas, we have relaxed the constraint to adopt a penalized formulation and derive more sophisticated norms that will be extensively described in Chapter for shots were generated all together using the algorithm proposed in [Boyer ] to draw samples according to a variable density with a polynomial decay of (i.e., h(k) = 1/ k 2 with k = (k x , k y )). A Cartesian reference scan composed of 512 shots of 512 samples each has also been acquired using the same coil. For comparison purposes, we used the square root of the sum-of-squares. The acceleration factor in time AF , corresponding to the time ratio between the accelerated acquisition and the Cartesian reference, the slice thickness and the object under investigation are reported in Table . . While exvivo acquisitions are relevant for quantitative comparisons, quantitative comparisons for in-vivo acquisitions should be considered carefully due to possible motion artifacts or misregistration between the Cartesian reference scan and the accelerated one. Hence, the

. . Experiments & Results

latter will be only compared visually for illustration purposes. The reader is referred to Appendix A for in-vivo comparison. ( shots)

Reconstruction

In this benchmark we considered:

. Sparse-SENSE formulation with the proposed extraction [El Gueddari b];

. 1 -ESPIRiT [Uecker ] implemented within BART toolbox [Uecker ];

.

p-LORAKS [Haldar ] implementation provided by the authors of [Kim b]

using a gridding step before the reconstruction;

. CLEAR method [Trzasko ] with various patch shape with overlaps.

We omitted blind-bi-linear reconstruction since the reported computation time (see Table . )

was considered excessively long for relatively simple cases, i.e. Cartesian sampling,channel coil, and 256 × 256 image size. The CALM method was also omitted since Chapter proposes an extensive comparison with this method and its variants. All the regularization parameters were tuned using a grid-search procedure to maximize the Structural Similarity Index (SSIM) score [Wang ]. For channel-by-channel MR reconstruction the coil combination was performed using the SOS. We considered the Symmlet basis as sparsifying transform when needed. All experiments were run on a machine with GB of RAM and an -core ( . GHz) Intel Xeon E -v Processor.

. . Results

To study the impact of different reconstruction methods on prospectively accelerated non-Cartesian acquisitions, we proceeded in two steps. First, we considered the SENSE-like reconstruction to evaluate the impact of sensitivity profiles on the reconstructed image. As no ground-truth exists for these sensitivity maps, in a second step we validated the results on the reconstructed images where we compared the SSIM scores for the different techniques (i.e., Sparse-SENSE with the developed sensitivity profile extraction, 1 -ESPIRiT, p-LORAKS and CLEAR methods).

Sensitivity matrix impact on the reconstructed images

First, we studied the impact of the sensitivity profiles over the reconstructed image using a SENSE-like prior, with the extraction method proposed [El Gueddari b]. Using our method, we investigated the impact of the k-space portion used to extract sensitivity information. The results are presented in Figure . -, where we reported the SSIM score of the reconstructed image for different values of θ ∈ {1%, 5%, 25%, 50%, 100%}. As regards the setting of parameter λ, we performed a grid-search to find its optimal value. This figure suggests that for an optimal λ value -5 × 10 -6 -small k-space portionsi.e.

θ < 25% -The SSIM score is maximized. However, images reconstructed with sensitivity profiles extracted from small k-space portion are more sensitive to the tuning of λ. In fact when parameter λ is underestimated (e.g. λ < 10 -6 ), medium to large k-space coverageθ > 25% -provide better reconstruction. In contrast, a slight overestimation of λ impacts less image quality after reconstruction.

Figure . -: Impact of sensitivity profiles estimation on the final reconstructed image quality.

The last two settings (θ ∈ {0.5, 1}) yield superimposed traces indicating that there is no gain in collecting high-frequency information for the recovery of sensitivity profiles.

Visually, examples of the extracted sensitivity profiles, estimated with 1 -ESPIRiT and our approach are presented in Figure . -for the mm slice-thickness baboon brain acquisition (cf Table . ( .)). The striking difference between the two approaches lies in the patterns of S images (see Figure . -): because of the SVD decomposition involved in 1 -ESPIRiT, the corresponding sensitivity profiles are less structured and smoother than the maps yielded by our approach. The latter appear sharper and clearly show the sensitivity area of each channel. Moreover, our approach is much more efficient in terms of computing time since the average cost is about min using a Matlab R a-based implementation as compared to min for 1 -ESPIRiT on the same hardware configuration.

Proposed extraction

1 -ESPIRiT th sensitivity map th sensitivity map th sensitivity map th sensitivity map

Figure .

-: out of sensitivity maps extracted using our method (left) and 1 -ESPIRiT (right).

Reconstructed images

Table . presents the quantitative comparison with respect to the ground truth based on the SSIM score. According to the SSIM score, the Sparse-SENSE reconstruction provides the best results for the three data sets under study. A visual illustration is provided in Figure . -on the mm ex-vivo baboon brain which confirms these quantitative results. we could only consider a sliding window of size and (Figure . -green frame), since we were running out of memory on a GB RAM workstation for larger window sizes.

Last, CLEAR results presents very similar results compared to 1 -ESPIRIT in term of SSIM score, however the reconstruction is more affected by the noise. Figure . -presents the hyper-parameter setting for CLEAR. The optimal parameter values were obtained with high overlapping factor of = 4 -corresponding to a stride s of -and a patch size p of 32 × 32. Also, large values of p provide a significant gain in robustness with respect to the setting of parameter λ irrespective of the overlapping factor: a high SSIM score is maintained overs almost five orders of magnitudes (λ ∈ 10 -6 , 10 -2 ). Importantly, while large overlapping factors prevent from patching artifact, they are more computationally demanding although they can be parallelized.

Moreover when we visually investigated the results of CLEAR method, we found some

. Conclusion

For high-resolution MR images, the input-SNR is a crucial issue as it allows to increase both the spatial resolution and the acceleration factor. As getting higher input-SNR often means collecting the signal with a multi-channel receiver coil, the input data set dimension increases and the reconstruction problem becomes harder. Multi-channel reconstruction algorithms can be split in two categories: either they deliver a single image from all the coils or one image per coil.

The first category relies on the knowledge of the coil sensitivity profiles that need to be estimated for each scan since they are scan-dependent. The estimation of coil sensitivity profiles is often based on a heavily sampled k-space region (also termed a calibration region) where much data have been acquired. Direct calibration methods use the coil-specific images as sensitivity maps (see e.g. [El Gueddari b]) whereas indirect approaches estimate the sensitivity maps using eigenvalue decomposition (cf ESPIRiT [Uecker ] and its 1 extension), still the direct approaches are more efficient in terms of computational time and memory footprint. As an alternative, the blind bi-linear method also relies on SENSE reconstruction but jointly reconstructs the coil sensitivity profiles and a combined image. However, due to the non-convexity of the overall objective function and the underlying numerical complexity, the computation time often remains the major bottleneck.

The second category performs coil-by-coil reconstruction and implements regularization in the k-space or image space. While k-space regularization provides good performances for Cartesian sampling, its application to non-Cartesian acquisitions often requires a gridding procedure which degrades the final image quality as shown in this chapter. Alternative image space coil-by-coil reconstruction (such as CLEAR or CALM), relies on joint sparsity promotion over the multiple channels. Image based coil-coil reconstruction presents three main advantages: i) Efficiency for non-Cartesian k-space acquisition as image reconstruction is solved using an inverse problem formulation.

ii) Independence with respect to the coil sensitivity profiles, no need to be extracted/estimated prior to reconstruction.

iii) Reliability as this formulation relies on the minimization of a convex criterion: any reached local minimum is global.

In terms of contribution, in this chapter we presented a new method to extract sensitivity maps and used a Sparse-SENSE formulation leading to a new self-calibrated method for non-Cartesian variable density sampling schemes. Compared to the state-of-the-art, our approach is much more efficient and the sensitivity profiles are easier to interpret.

We demonstrated on several anatomical T * 2 data collected at 7 T that our global selfcalibrated method is both more accurate and efficient than the competing alternatives such as 1 -ESPIRiT, p-LORAKS and CLEAR. The proposed method has also been exten-

. Outlook

In T he purposes of this chapter is to assess the impact of regularization on the image quality in a calibration-less reconstruction framework. We emphasize the importance of promoting structured sparsity across the multiple receivers to improve the overall image quality.

. Introduction

Multi-channel image reconstruction becomes challenging when it has to comply with timing constraints such as the ones involved for online reconstruction: basically, the SENSE formulation is no longer the most convenient as the sensitivity maps cannot be extracted from incomplete data, which does not cover a wide portion of the center of k-space. Hence, the self-calibrating methods presented in Chapter are no longer suitable to extract the sensitivity maps associated with the multiple receivers. Also, a second bottleneck exposed in Section . . is due to the dependence of the Lipschitz constant to those maps. In contrast, calibration-less methods do not require this prior knowledge to perform image reconstruction in the parallel imaging setting. This makes them more appealing for online reconstruction purposes. Calibration-less methods may operate either in the k-space domain or in the image or any transform (e.g. wavelet) domain. In the context of non-Cartesian acquisition, the existing k-space calibration-less methods in the literature do not perform as well as the self-calibrating ones mainly owing to the gridding step which is usually applied to perform regularization in the k-space. For this reason, we will investigate domain-based calibration-less methods and look at their advantages for imposing some structured sparsity across channels. Then we will evaluate the impact of those regu- From the noisy under-sampled data (y ) 1≤ ≤L , the goal is to retrive L MR images stacked in X = [x 1 , . . . , x L ] ∈ C N ×L such that each x ∈ C N is associated with the th channel of the phased array coil.

X = arg min

X∈C N ×L L =1 1 2σ 2 F M x -y 2 2 + R(ΨX) . ( . ) Hereabove, R ∈ Γ 0 (C N Ψ ×L
) is a regularization function composed with a linear operator

Ψ ∈ C N Ψ ×N
, with the aim to enforce sparsity of the solution within a given multiscale decomposition (e.g., wavelet transform). In the following we suppose that R have a closed form proximity operator. This formulation enables the use of over-complete dictionaries [Elad ]. We will assume that Ψ decomposes the stack of

L images X ∈ C N ×L into a stack of coefficients ΨX ∈ C N Ψ ×L with C scales. Each scale c ∈ {1, . . . , C} is composed of S c sub-bands. Each sub-band s ∈ {1, . . . , S c } has K s(c) coefficients, so that finally N Ψ = C c=1 Sc s=1 K s(c)
. We also introduce Z as the sparsifying decomposition

Z = ΨX = [z 1 . . . z ] ∈ C N Ψ ×L of the muti-channel image X = [x 1 . . . x ] ∈ C N ×L .
The resolution of Problem ( . ) delivers L channel images ( x ) 1≤ ≤L , stacked in X. To retrieve the magnitude and phase information many methods exist as discussed in Section . , in this work we will use the square-root of the sum-of-squares (sSOS) [Roemer ], x sSOS = L =1 x 2 2 , to form a single magnitude image and the virtual coil method [Parker ] for phase image.

. . Optimization algorithm

In order to solve Problem ( . ) we adapted a primal-dual algorithm proposed in [Condat , Vũ ] and recalled Section . , leading to the Algorithm

Algorithm : Condat-Vú algorithm Set τ > 0, κ > 0, X 0 ∈ C N ×L , Z 0 ∈ C N Ψ ×L ;
for k = 0, . . . , T do

X k+1 := X k -τ (∇f (X k ) + Ψ * Z k ); W k+1 := Z k + κ Ψ (2X k+1 -X k ); Z k+1 := W k+1 -κ prox R/κ W k+1 κ ; end
According to [Condat , Theorem . ], the sequence (X k ) k∈N generated by Algorithm weakly converges to a solution of Eq. ( . ) as soon as

1 τ -κ|||Ψ||| 2 ≥ β 2 , with β = L l=1 (σ 2 ) -1 |||F M ||| 2
the Lipschitz constant of the gradient of the data consistency term. In practice, the hyper-parameters of this algorithm are set as follows: τ := 1 β ,

Calibrationless multi-channel reconstruction

κ := β 2|||Ψ||| . Note that when Ψ defines a basis, we get |||Ψ||| = 1. The main advantage of Algorithm is that it does not involve the computation of prox R•Ψ . The latter does not usually have closed form, in particular when Ψ is overcomplete (e.g. undecimated wavelet transform), and would require the use of an inner iterative solver [Combettes ].

. From group-sparsity . .

Least Absolute Shrinkage and Selection Operator

The 1 -norm has been successfully used for promoting sparsity of the solution in many different fields such as generalized linear models, Bayesian statistics and dictionary learning.

The success of such norm in the regression/regularization domain is mainly due to its good with respect to the sparsifying dictionary (as it is the case of Sparkling trajectories), the 1 -norm can be perfectly applied to multi-channel image reconstruction. Hence, in the calibration-less setting the regularization term involvoing the Wavelet coefficients Z can be summarized as follows: ) with λ being a positive hyper-parameter that needs to be set. We recall that the proximity operator of the LASSO corresponds to a soft-thresholding step [Bertsekas ]:

∀Z ∈ C N Ψ×L , λ ∈ R + , R LASSO (Z) = λ Z 1 = λ N Ψ n=1 L =1 |z n, | ( .
prox R LASSO (Z) n, =      z n, 1 -λ |z n, | , if |z n, | ≥ λ 0 , otherwise ( . )
Although the good performances of the LASSO have been demonstrated in many domain, such as dictionnary learning [Mairal ], this regularization also has some well known drawbacks such as its poor behavior in the case of highly correlated variables [Zhao , In the context of multi-channel acquisition, the k-space seen by the different channels is actually highly correlated due to the coil configuration. Therefore in the following we will describe alternatives to LASSO regularization and their application to the multi-channel reconstruction.

. . From group-sparsity

. . Elastic-Net

The first alternative that has been proposed is the Elastic-net regularization [Zou ].

[Zou ] presented Elastic-net as the right counterpart to the LASSO penalty when variables are highly correlated. Elastic-net regularization performs a linear interpolation between 1 (LASSO) and 2 (Ridge regression [Hoerl ] or Tikhonov regularization [Tikhonov ]) norms and is defined as follows:

∀Z ∈ R N Ψ ×L , (λ 1 , λ 2 ) ∈ R 2 + , R Elastic-net (Z) = λ 1 Z 1 + λ 2 2 Z 2 2 ( . )
The proximity operator of such norm is easy to compute, although the authors [Zou ] noticed that the naive version of Elastic net causes a double amount of shrinkage which led to poor predictions. Hence they proposed to rescale the thresholded coefficients by a factor of (1 + λ 2 /2) and thus the proximity operator is reduced to:

prox R Elastic-net (Z) n,l =      z n, 1 - λ 1 (1+λ 2 )|z n, | × 1 + λ 2 2 , if |z n, | ≥ λ 1 1+λ 2 0 , otherwise ( . )
Although the LASSO and Elastic-net regularization terms do not explicitly exploit redundant information provided by the different coils, in the following we will explore various extensions that handle the structural information and redundancy across channels.

To begin with, we study the group-LASSO penalty.

. .

Group-LASSO

Definition and application to image reconstruction

In order to overcome the above mentioned drawbacks associated with LASSO regularization and select a group of variables, the group-LASSO [Yuan ] has been proposed.

Given a set of groups G the group-lasso Ω G group is basically defined as a 2,1 mixed-norm:

∀w ∈ C p , Ω G group (w g ) = g ∈G w g 2 ( . )
with G being the non-overlapping groups that need to be specified a priori. In the context of MR image reconstruction, a group gathers the same wavelet coefficients over all channels [Majumdar a], the regularization can thus be defined as:

∀ Z ∈ C N Ψ ×L , λ ∈ R + R gLASSO (Z) = λ Z 2,1 = λ N Ψ n=1 L =1 |Z n, | 2 ( . )
The proximity operator of the 2,1 mixed-norm is explicit and equal to:

prox R gLASSO (Z) n, =          z n,   1 - λ L =1 |z n, | 2   , if L =1 |z n, | 2 ≥ λ 0 , otherwise ( . )
Mixed norm in other fields. The group penalty has arisen a strong interest in the inverse problem community as it has been used in many applications such as spatiotemporal MEG/EEG source reconstruction [Strohmeier , Bekhti ], where the authors proposed an iterative scheme to approach the 2,0.5 -pseudo norm.

Overlapping Group-LASSO. An overlapping version of the group-LASSO also has been introduced by [Jacob ]. The regularization term is formulated as follow:

w ∈ C p , Ω G overlapp (w) =    g∈G w g 2 : supp(w g ) ⊆ g, w = g∈G w g    ( . )
with G the set of all groups. The major bottleneck of the use of such norm is related to the definition of the groups as the latter need to be defined a priori.

Expected improvement for multi-channel image reconstruction

The theory behind multi-channel sparse recovery problems has been extensively studied by the CS community as multiple-measurement vector (MMV) [Chen , Eldar a, Eldar b, Duarte , Cotter ] and uses the multi-channel basis pursuit:

∀Y ∈ C p×L , min Z∈C N ×L Z 2,1 s.t. Y = AZ ( . )
Although in the worst case scenario [Chen ] the use of the group-sparsity norm may not be beneficial in terms of recovery guarantees, [Eldar b] has proposed an average case study and shown the potential benefit of adding the group structure with a probability of recovery failure decaying exponentially with the number of channels.

In the MRI context, [Chun ] has developed similar results. The authors have shown that for the case of disjoint sensitivity matrices the recovery guarantees of the multi-channel basis pursuit can be divided by a factor of L, with L the total number of channels. However those results were given for the case of disjoint sensitivity matrices, i.e.

when the sensitivity profiles do not overlap and thus bring complementary information, which is not true in the real life Section . . Therefore the benefits of such norms must be mitigated for real case applications.

.

To structured sparsity

In order to exploit the redundant information provided by each coil, the group-LASSO regularization requires a prior definition of the group-structure. While a straight-forward prior for calibrationless reconstruction consist of defining a group by the same coefficients for all the channels, a more appropriate approach would be to consider different groupstructure that leverages the group definition to maximize the final image quality. The following section will be dedicated to more sophisticated regularization terms that infer the structure of the groups. In particular we focus on two competitors, the first one presented in Section . . is based on a Clustering algorithm, and finds the structure promoting a pairwise ∞ -norm, while the second one in Section . . based on the 2norm defines the sets of all possible overlapping groups with a cardinality at most equal to k.

. . Octagonal Shrinkage and Clustering Algorithm for Regression

In the context of highly correlated variables, the first method one could think of consist In what follows, we propose four choices for function R in Problem ( . ) relying on OSCAR norm, with the aim to perform efficient calibration-less MR image reconstruction.

The main difference between these formulations lies in the way the OSCAR norm is applied to the sparsifying decomposition Z of the muti-channel image X.

General definition

Let z ∈ C p with p ≥ 1. We introduce the magnitude sorting operator S p : C p → C p such that vector S p (z) = (S p (z) j ) 1≤j≤p contains the p entries of z sorted in decreasing order in magnitude, i.e. such that

|S p (z) 1 | ≥ |S p (z) 2 | ≥ • • • ≥ |S p (z) p |. ( . )
Then, OSCAR norm is defined as follows:

Ω λ,γ (z) = λ z 1 + γ 1≤j<k≤p max(|z j |, |z k |). (

Proximity operator

The proximity operator of z ∈ C p can be efficiently computed thanks to the following algorithm as shown in [Zeng a, Sec.III A]:

Algorithm : Proximity operator of the OWL norm.

Input: z ∈ C p /{0}, w ∈ R p ; n = |z|/z; Let P ∈ R p×p s.t. S p (n) = P n; Return: prox Θw (z) = n P PAV(S p (n) -w);
Hereabove, PAV refers to the Pool Adjacent Violator Algorithm [Mair ]. For z is equal to zero, then the proximity operator of the OWL norm at z is also equal to zero.

The proximity operator of OSCAR can thus be easily deduced by setting the appropriate value for w mentioned above, so that OSCAR and OWL match together.

In the following we apply different versions of OSCAR regularization to the multi-

channel wavelet coefficients Z ∈ C N Ψ ×L .

Global-OSCAR regularization

The most straightforward way to implement OSCAR-based regularization consists of flattening all wavelet coefficients and thus discarding the multiscale and multi-channel structure (e.g. in scales, sub-bands, coefficients and channels) of Z. For that reason, we call this version global OSCAR (g-OSCAR) regularization. The wavelet coefficients are stacked together, leading to a single but large vector with entries (z j ) 1≤j≤LN Ψ , where we remind that

N Ψ = S C c=1 K c and N Ψ = N when Ψ is orthogonal. The g-OSCAR regularizer then reads: R g-OSCAR (Z) = Ω λ,γ (Z) = LN Ψ j=1 λ|z j | + γ 1≤j<k≤LN Ψ max{|z j |, |z k |} = LN Ψ j=1 (λ + γ(LN Ψ -j)) |S LN Ψ (z) j | . ( . )
with S LN Ψ , the magnitude decreasing sorting operator of dimension LN Ψ .

Scale-wise OSCAR regularization

We now propose a scalewise formulation, where OSCAR norm is applied to each specific scale c of the wavelet decomposition, hence to each vector z c,: separately where z c,: gathers the LSK c wavelet coefficients across all channels in a specific scale c ∈ {1, . . . , C}. This leads to the so-called s-OSCAR regularizer:

R s-OSCAR (Z) = C c=1 Ω λ,γ (z c,: ) = C c=1 LSKc j=1 (λ + γ(LSK c -j)) |S LSKc (z c,: ) j | , ( . )
with S LSKc , the LSK c -dimensional magnitude decreasing sorting operator. In that way, the wavelet coefficients can be clustered together regardless the subband they belong to, their position and their channel dependence. Thus, C sorting operations are required, each of them involving LSK c parameters. As the s-OSCAR regularization is separable by scales, the computation of its proximity operator can be performed efficiently using parallelization over scales.

Band-wise OSCAR regularization

The present formulation applies OSCAR regularization to each specific subband of the wavelet decomposition, hence to each vector z cs,: separately where z cs,: gathers the K c L wavelet coefficients across all channels in a given subband s of scale c. The band-wise (b-OSCAR) regularization thus reads:

R b-OSCAR (Z) = C c=1 S s=1 Ω λ,γ (z cs,: ) ( . ) = C c=1 S s=1 KcL j=1 (λ+γ(K c L -j)) |S LKc (z cs,: ) j | ,
with S LKc , the LK c -dimensional magnitude decreasing sorting operator. Here again, the separability of the regularizer can be exploited for an efficient implementation of the proximity operator. as it operates separately on each pixel in the transform domain. However, instead of implicitly assuming constant noise level over all receivers by taking the 2 -norm, the c-OSCAR regularization allows one to weight these receivers differently thanks to the sorting step and to adapt regularization to stay robust to a space varying noise.

Calibrationless multi-channel reconstruction

Coefficient-wise OSCAR regularization

. . k-support norm

The 1 norm has been intensively and successfully used for solving sparse estimation problems. Its usage has often been justified as being the tightest convex envelope of the 0 pseudo-norm. However this holds only in the case of vectors with upper bounded entries by as it is easy to show that ∀z ∈

C p , z 1 ≤ z ∞ z 0 . [Argyriou ]
claims that exploiting this relation between the 1 -norm and the ∞ norm might be suboptimal (the relation with the ∞ norm has been studied in Section . . with OSCAR).

Instead of exploiting the ∞ -norm they define a new surrogate based on the 2 -norm such as: z 1 ≤ z 2 z 0 and studied the Gauge function associated with the following convex hull:

C k := conv(S k ) = conv(z, z 0 ≤ k, z 2 ≤ 1). ( . )

General definition

Although the formulation . could be difficult to solve [Argyriou ] proposed a variational formulation defined as follows:

∀z ∈ C p , z sp k := min{ I∈G k v I 2 : supp(v I ) ⊆ I, I∈G k v I = z} ( . )
where G k denotes the set of all subsets of {1, • • • , p} of cardinality at most k . To compute the k-support norm the authors proposed an algorithm with a complexity of O(p log p).

We recall the proposition hereafter:

The Gauge function γC k associated with a convex set C k is defined as

γC k (x) = inf{λ ∈ R+ : x ∈ λ C k }
The k-support norm is therefore defined as the norm whose unit-ball is the convex-hull of the set of vectors of cardinality at most equal to k and with an 2-norm not greater than . Proposition . . . For z ∈ C p , and k ∈ N, k ≤ p,

z sp k =    k-q-1 j=1 (S p (|z|) j ) 2 + 1 q + 1   p j=k-q S p (|z|) j   2    1 2 , ( . )
with: S p the magnitude sorting operator, S p : C p → C p such that vector S p (z) = (S p (z) j ) 1≤j≤p

contains the p entries of z sorted in decreasing order in magnitude and q being the unique integer in {0, • • • , k -1}, and z 0 = ∞, that satisfies:

S p (|z|) k-q-1 ≥ 1 q + 1 p j=k-q S p (|z|) j ≥ S p (|z|) k-q ( . )
when k = 1 then q = 0 and the remaining term is the second one wich correspond to the 1 -norm, in a similar way when k = p then q should be equal to 0 to verify the following inequality: S p (|z|) p-1 ≥ S p (|z|) p ≥ S p (|z|) p leading to the definition of the 2 norm. Relation to the overlapping group-LASSO. Interestingly the variational definition of the k-support norm in Eq. ( .) shows similarities to the overlapping group-LASSO in Eq. ( .) when the sets of groups are chosen to be p k . Since the potential number of groups is exponential this makes the overlapping group-LASSO untractable in the case of all the possible subsets of cardinality at least k. On the other hand, the k-support norm can be efficiently computed with a complexity of O(p log(p)).

Dual norm.

Although the variational formulation in Eq. ( . ) seems complex, the dual of the k-support norm is surprisingly easy to compute and corresponds to the 2 -norm of the k largest entries:

u sp k := max{ z, u : z sp k ≤ 1} = max{( i∈I |u| 2 i ) 1 2 : I ∈ G k } = k i=1 S p (|u|) 2 i 1 2
( . ) The dual of the the k-support norm could be seen as an interpolation between the ∞ and 2 since for k = p we get the dual of the 2 norm and for k = 1 the dual takes the largest entries of the u vector which corresponds to the ∞ norm (i.e. the dual of the 1 norm).

Proximity Operator [Argyriou , McDonald ] proposed to use the k-support norm as a sparse regularization for a learning task , in particular they consider the square of the k-support norm, which led to the following regularizer:

∀z ∈ C p , R λ 2 • sp k 2 (z) = λ 2 min{ I∈G k v I 2 2 : supp(v I ) ⊆ I, I∈G k v I = z} ( . )
where G k denotes the set of all subsets of {1, • • • , p} of cardinality at most k. According to [McDonald , Corollary ] the proximity operator of the k-support norm is defined

prox λ 2 • sp k 2 (z) i = θ i z i θ i + λ ( .
)

θ i =            1, if α|z i | > λ + 1 α|z i | -λ, if λ + 1 ≥ α|z i | ≥ λ 0, if λ > α|z i |
The k-support norm encourages z to be a sum of limited number of vectors with small support.

where α is chosen such that p i=1 θ i (α) = k. Even though the computation of this proximity operator is complicated it can be done efficiently using Algorithm .

Algorithm : Proximity operator of

λ 2 • sp k 2 . Input: z ∈ C p /{0}, λ ∈ R ; Output: x = prox λ 2 • sp k 2 (z) ; Compute α 2i , α 2i+1 p i=1 = λ |z i | , λ+1 |z i | p i=1 ; Sort α i such that α i ≤ α i+1 ;
Find i such that:

p j=1 θ j (α i ) ≤ k and p j=1 θ j (α i+1 ) ≥ k ; Interpolate α * between α i and α i+1 such that p j=1 θ j = k; θ i (α * ) =            1, if α * |z i | > λ + 1 α * |z i | -λ, if λ + 1 ≥ α * |z i | ≥ λ 0, if λ > α * |z i | x i = prox λ 2 • sp k 2 (z) i = θ i z i θ i +λ ;

Application to MR image reconstruction

We implemented the straightforward version of the k-support norm regularization, that consists of flattening all wavelet coefficients and thus discarding the multi-scale and multichannel structure (e.g. in scales, sub-bands, coefficients and channels) of Z. Following the wording used for OSCAR-norm regularization, this implementation corresponds to the global version of k-support norm and it models the signal as the g-OSCAR regularization (cf Section . . ). The wavelet coefficients are stacked together, leading to a single but large vector with entries (z j ) 1≤j≤LN Ψ , where we remind that N Ψ = S C c=1 K c and N Ψ = N when Ψ is orthogonal. The k-support regularization then reads:

R λ 2 • sp k 2 (Z) = λ 2 min{ I∈G k v I 2 2 : supp(v I ) ⊆ I, I∈G k v I = Z} ( . )
.

Experimental setting & Results

To Hereafter, we describe the experimental setup used for validation purposes.

. . Experimental setup

All numerical experiments were conducted on D k-space data even though the proposed framework could be extended to D imaging (i.e. D k-space data) directly. Hence, we only report results on slices.

Prospective non-Cartesian acquisition

We consider the reconstruction of an ex vivo human brain with an in plane resolution of 0.39 × 0.39mm 2 and a slice thickness of mm for different acquisition scenarios prospectively accelerated on the same T MR system and with a L = 32-channel coil (Nova Medical). Ex vivo imaging offers the best imaging conditions to remain insensitive to motion when validation is performed on truly accelerated acquisitions. In the literature, for the sake of simplicity, most of reconstruction algorithms are usually validated on retro-

Reconstruction parameters

Algorithm was run until the maximum number of iterations (T = 150) or a tolerance threshold on the image solution was reached, which appears sufficient to converge.

Moreover, we used for Ψ a Daubechies orthogonal wavelet transform (OWT) with C = 4 decomposition scales (i.e., N Ψ = N ). Note that MR image quality could probably be improved using redundant multiscale transforms such as undecimated bi-orthogonal wavelet transforms or curvelets as shown in [Cherkaoui , Ma ] but those transforms may increase the time and memory requirements of the overall algorithm. The hyperparameters were set using a grid-search procedure such that the SSIM [Wang ] score was maximized. In the next section we also report the peak Signal to Noise Ratio (p-SNR)

and the normalized root mean square error (NRMSE) scores.

. . Results

Before presenting in more details the results corresponding to the different regularization schemes, we first computed the zero-order solution (iFT). Since we used non-Cartesian trajectories, the forward model relies on the non-Uniform or Non-equispaced Fast Fourier

Transform (NUFFT or NFFT) [Fessler , Keiner ], which is not invertible, a conjugate gradient descent was performed on the data-fidelity term Eq. ( .). In what follows, we call the corresponding estimator the least squares solution. Although not perfect, those results presented in Table . will be taken as a baseline for comparison purposes. The similarity scores (SSIM and p-SNR) decreases as the under-sampling factor increases for both Sparkling and Spiral reconstructions, which is not surprising as the problem gets harder with fewer measurement. We illustrate the results provided by the least square solution on the -fold accelerated scans for both Sparkling and spiral acquisition on

Comparison of coefficient-based regularization schemes

Additionally, we compare the group-LASSO and the c-OSCAR regularization as both impose a group-structure on the same wavelet coefficient across channels. Table . summarizes the results of the magnitude image obtained using the square root of the sum of square. Once again the hyper-parameters were optimized such that the SSIM index was set to its maximal value, however we also reported the p-SNR and NRMSE scores. Although for all tested sampling patterns the best results were obtained using the c-OSCAR regularization, the difference is not significant enough to notice a notable impact on the final images as presented in Figure . -.

Comparison of OSCAR-norm regularization schemes

As regards the different versions of OSCAR-based regularization, Table . summarizes their respective performances.

For the five Sparkling data sets, the best scores were achieved with the b-and c-OSCAR versions, where each subband and each wavelet coefficient are treated separately, respectively. This is not really surprising as the other versions mix up different orientation or scale details together. The fact that the coefficient-wise version performs slightly better than the subband-wise one suggests that localized regularization in space preserves much better details as the acceleration factor increases. Similar findings were replicated on an ex-vivo baboon brain (results not shown). For the five spiral data sets, the best SSIM scores were yielded in the vast majority of cases by the b-and c-OSCAR regularization schemes illustrating the consistency between the two readout scenarios.

Calibrationless multi-channel reconstruction

Comparison with other reconstruction methods Tab. . summarizes the quantitative performances (SSIM, p-SNR and NRMSE) of competing reconstruction approaches, computed over the magnitude images. These competitors embed i) zero-order inverse NFFT (no regularization, i.e. LS solution), ii) CALM solution which implements a group-LASSO regularization and self-calibrating techniques, namely 1 -ESPIRiT as well as k-space based calibration-less method AC-LORAKS. As in sparkling and spiral imaging, the k-space center is heavily sampled, 1 -ESPIRiT performs well for extracting the sensitivity maps (i.e. low resolution information) which are then used for solving the 1 -norm regularized CS-SENSE reconstruction problem. The bottleneck in LORAKS method lies in the gridding step required to project the non-Cartesian sampling pattern onto a binary Cartesian mask. This step degrades the image quality in the end and demonstrates that AC-LORAKS may not be really compliant with non-Cartesian sampling.

Table . : Quantitative comparison of state-of-the-art methods for non-Cartesian PI-CS reconstruction using either a calibration-less or self-calibrating approach. The sparkling and spiral solutions on the red frame are shown in Figure . -and Figure . -, respectively.

In terms of SSIM

. Conclusion

Non-Cartesian trajectories have been used for many different applications and have recently met a renew of interest to accelerate scan time. However, in the multi-receiver coil context most of reconstruction methods are more suited for Cartesian acquisitions than for non-Cartesian ones. This is especially the case of k-space based calibration-less techniques.

Indeed, we formulated CS-PI image reconstruction as an inverse ill-posed problem and made use of NFFT for dealing with non-Cartesian trajectories and structured sparse induced norm for regularization purposes. We instantiated different versions of structured sparse norm regularization to exploit the redundant information provided by each coil. Then, we relied on state-of-the-art convex nonsmooth optimization tools, namely the Condat-Vù algorithm with sound convergence properties to compute the global minimizer of the derived cost function. The main advantages of this algorithm are first that it can efficiently deal with analysis-based regularization which is known to provide better results than synthesis-based priors. Second, this algorithm is highly flexible, as it can be used the same way for a large range of penalty terms and the parameters that control its Bottom: Wrapped phase images extracted using the virtual coil approach [Parker ] .

Calibrationless multi-channel reconstruction

. Discussion

In the future, we will explore alternative automatic hyper-parameters settings either based on (generalized) cross-validation [Hastie ] across slices for instance as a larger data set is necessary to implement this strategy, or on statistical inference (e.g. SURE estimator [Stein ]), which requires the knowledge of noise statistics. The overall SSIM scores may also be improved using redundant wavelet transforms either considering undecimated multiscale decompositions or more sophisticated transforms such as curvelet or shearlet transforms [Cherkaoui , Ma ]. Although the aforementioned approaches might improve MR image quality for each reconstruction, it could also increase the computation time and the memory loading that can be critical especially for D high-resolution multi-channel reconstruction.

This work did not address the case of accelerated Cartesian acquisition and reconstruction. In this setup, the proposed approach can be useful as far as low acceleration factors are considered. However when the acceleration is pushed further and the sampling is kept over the Cartesian grid, the proposed calibration-less methodology -including CALM and CLEAR methods -may perform worse than state-of-the-art self-calibrating methods ( 1 -ESPIRiT) or domain based calibration-less (P-LORAKS). Nevertheless one could think of a self-calibrating approach with a more sophisticated regularization (e.g. OSCAR norm or group-LASSO) as done in [Chun ].

In order to speed-up reconstruction, many different method have been proposed by the Deep-Learning community, in particular unrolled optimization [Kamilov , Yang , Adler ] methods seems the most promising one. With this perspective, [Wen ] generalizes the structured sparsity norm to use it on deep neural network. However deeplearning methods requires a large database of raw data (i.e. high resolution multi-channel complex k-spaces unavailable until recently [Zbontar ]), whereas traditional methods such as the one proposed are not. Importantly the proposed method has theoretical guarantees of convergence and is made of fully interpretive steps which make it more reliable and robust to small perturbations than recent deep learning approaches (such as [Zhu , Mardani ]) proposed for medical image reconstruction as it has been

shown by [Antun ].

In terms of applications, the proposed method can be easily extended to D imaging and used as such in isotropic high-resolution susceptibility weighted imaging (SWI). In this context, the scan time is lengthy and in the clinical realm SWI acquisitions will benefit from highly accelerated non-Cartesian encoding schemes to reach µm in a scan time of -min. Noticeably, in SWI imaging, post-processing is applied to phase information in order to reveal potential alteration of the microvascular brain network. The fact that our calibration-less regularized reconstruction better preserves phase information is thus an asset for its utilization in SWI application for diagnostic purposes.

. Introduction

Instead of performing offline image reconstruction by minimizing a sparsity promoting regularized objective function as most of existing works [Lustig , Guerquin-Kern , Chaâri ] do, we introduce an online image reconstruction approach so that acquisition and reconstruction processes become interleaved. This method is well-suited for high-resolution MR acquisition. Indeed, the short lifespan of the MR signal makes the spatial encoding in large k-spaces (i.e. high-resolution setting) necessarily segmented. The ultimate goal is actually to reduce the cumulative time of acquisition and reconstruction.

To this end, we adopt a mini-batch formalism, that consist of stacking together collected k-space samples to form a so-called mini-batch. Once a mini-batch is available, the reconstruction algorithm is run from an incomplete k-space data set for a few iterations.

Hereafter, the current image solution is then used as a warm restart for the next mini-batch processing. This new framework allows us to reach a decent image quality by the end of acquisition and can be implemented throughout the Gadgetron [Hansen ] project. It thus enables us to display the resulting images directly to the MR system console, making online CS implementations available beyond any vendor solution.

The rest of the Chapter is organized as follows. We first recall the link between previous work related to online image reconstruction in the context of dynamic MRI.

Then we derive the problem for anatomical imaging in a simplified setup, when then data is collected over a single channel receiver coil. The proposed method is evaluated on a retrospective Cartesian acquisition as well as on prospective non-Cartesian k-space data, we scanned on a Magnetom 7 T scanner (Siemens Healthineers , Erlangen, Germany) an ex-vivo baboon brain using a D high-resolution ( µm and mm slice thickness) T * 2weighted imaging with single channel birdcage coil. Although of pedagogical interest, the practical application of this preliminary study remains limited for high-resolution imaging as we need multiple receivers coil to boost the input SNR in this context. For that reason, the online formalism is then extended to deal with multiple-channel or multireceiver coil. To this end, we actually rely on the calibrationless image reconstruction approach presented in Chapter . The proposed online reconstruction method is evaluated on prospective non-Cartesian k-space data collected with acceleration factor up to in time, where an ex-vivo human brain was scanned using Sparkling trajectories at T for D high-resolution ( µm in plane and . mm slice thickness) T * 2 -weighted imaging using a multi-receiver phased array. Using the same ideas [Dietz ] has proposed an online approach combining CS and Principal Component Analysis, in particular they have demonstrated the benefit of this method for tracking lung cancer.

Related works in dynamic MRI

However this work does not rely on the same assumption as for dynamic MR reconstruction. In the case of dynamic MRI, there is a trade-off between spatial and temporal resolution, whereas in high-resolution anatomical imaging the only one that matters is spatial resolution. Moreover due to the resolution and in order to improve the signal to noise ratio acquisition is necessary segmented which makes the acquisition even lengthier.

The purpose of this chapter is to take advantage of this time to start the reconstruction and to deliver a decent image as fast as possible while maintaining the same image quality compared to the offline reconstruction.

. Extension to anatomical imaging: general online problem statement

To take advantage of the sequential aspect of segmented acquisition, the reconstruction process starts during acquisitions' dead times associated with the time intervals separating consecutive shots (i.e. TR when the magnetization grows back). The acquired shots are stacked to form a mini-batch and the reconstruction starts from this incomplete data set. Hence, the two processes get interleaved and partial feedback may be delivered to the MR technician along the exam with the aim to provide a nearly optimal image by the end of acquisition . Two acquisition setups will be studied hereafter. First, we will consider the case of a single-channel receiver coil. Though simpler, this case is rarely met in the high resolution context since multi-channel receiver coils are traditionally used to boost the input SNR. The second setup corresponds to online reconstruction from data collected over a multi-channel receiver coil. While state-of-the art methods often rely on the coil sensitivity information that requires an estimation step for each scan as described in Chapter , we will rely on a calibration-less method proposed in Chapter that tackles the image reconstruction problem with the great advantage of not requiring the prior knowledge of sensitivity matrices.

In the sense that the recovered image at this stage is close to the final one obtained from the whole data set after complete convergence of the reconstruction algorithm.

. . Acquisition considerations

Let us first recall some acquisition considerations presented in Chapter . In segmented acquisition, multiple radio-frequency (RF) pulses are used to tip the global magnetization of a given volume. More precisely, at each repetition time (TR), a new RF pulse is delivered and a "new shot is collected" [Bernstein , Chap. . ].

In the following, S represents the number of shots used to fill the k-space hence, the total scan time is equal to S×TR. For each i ∈ {1, . . . , S}, we refer to the k-space support of the i th -shot as Γ i and y Γ i ∈ C C gathers C samples measured over this shot.

Next, let k ∈ {1, . . . , S}, we define Ω k as the concatenation of the k first collected shots.

Hence . .

Ω k = ∪ k i=1 Γ i ,

Variational formulation of image reconstruction

In standard offline approaches, image reconstruction from the k-space data is performed by minimizing an objective function that sums a data consistency term f Ω S (depending on the trajectory support Ω S ) and a regularization term g, which usually promotes sparsity in a given multiscale decomposition Ψ ∈ C N Ψ ×N such as a wavelet transform [Pustelnik ].

The general offline image reconstruction problem aims at finding:

x ∈ arg min x∈C N f Ω S (x) + g(Ψx), ( . )
where the following standard assumptions are made: (i) g ∈ Γ 0 (C N Ψ ) with a closed form proximity operator, and (ii) f Ω S is convex, differentiable on C N and its gradient ∇f Ω S is β S -Lipschitz i.e.:

(∀(x, x ) ∈ C N ) ∇f Ω S (x) -∇f Ω S (x ) ≤ β S x -x . ( . )
Problem ( . ) can be efficiently solved using iterative solvers as presented in Section . . Note that for over-complete dictionaries [Elad ] the use of primal-dual approaches is more efficient in terms of computation time since these methods do not require the use of an inner iterative solver [Combettes ] to compute the proximity operator of g • Ψ.

For online reconstruction, the data consistency term is progressively filled during the scan and reconstruction starts with incomplete data. Concretely, incomplete versions of . . Extension to anatomical imaging: general online problem statement Problem ( . ) are solved in order to compute:

x k ∈ arg min x∈C N f Ω k (x) + g(Ψx), ( . )
with k and x k being respectively the number of available shots so far and x k the resulting solution. For the processing of the next shot (i.e. k + 1), x k will be used as initialization.

In an ideal case, the shots are processed one by one (i.e., k is incremented by one at each step), however, such strategy is only feasible in practice if TR is longer than the time T it needed to solve SubProblem ( . ), which is rarely met in the non-Cartesian multi-channel acquisition. Therefore to adapt online reconstruction to this specific case, we adopt a mini-batch formulation presented in Algorithm where b s consecutive spokes are stacked together to form a mini-batch. Then, a warm restart strategy is implemented once an complete mini-batch is available, i.e. SubProblem ( . ) is not solved for any integer k but only for multiples of b s . In this setting, to comply with the online timing requirements, the maximum number of allowed iterations n b per subproblem in Algorithm is given by:

n b × T it ≈ b s × TR. ( . )
In Algorithm , the result of n b iterations of the chosen proximal algorithm for solving SubProblem ( . ), initialized with z 0 , is denoted by A k,n b (z 0 ). The variable z k usually includes an approximation x k to x k , but it may also embody additional information e. initialize k ≤ S, (x k,0 , v k,0 );

κ k = β k 2|||Ψ||| 2 ; τ k = 1 β k ; for t = 1, 2, . . . , n b do x k,t = x k,t-1 -τ k (∇f Ω k (x k,t-1 ) + Ψ * v k,t-1 ); w k,t = v k,t-1 + κ k Ψ (2x k,t -x k,t-1 ); v k,t = w k,t -κ k prox g/κ k w k,t κ k ; end return (x k,n b , v k,n b );
In Sections . . and . . we will discuss the choice of the mini-batch size b s and the setting of hyper-parameters for respectively single and multi-channel receiver coils.

Hereafter, we first derive the reconstruction problem for the single channel coil acquisition before analyzing the more challenging one, namely the multi-channel acquisition setup.

.

Single-channel receiver coil

. .

Problem statement

In the single-channel context, Algorithm is applied to the following mini-batch formulation of SubProblems ( . ):

(∀k ∈ {1, . . . , S}), xk = arg min

x∈C N S 2k F Ω k x -y Ω k 2 F + λ Ψx 1 , ( .
)
where λ is the positive hyper-parameter that controls the sparsity level and can be efficiently set in an offline manner. In the case of Cartesian acquisition, the following relation describes the operator

F Ω k = Ω k F
where Ω k is the under-sampling binary mask and F the fast Fourier transform (FFT). In the case of non-Cartesian acquisition, F Ω k refers to non-equispaced or non-uniform FFT [Keiner , Fessler ]. While the Lipschitz constant is β k = S/k for Cartesian acquisition, its value has to be estimated for each value of k when dealing with the non-Cartesian case. This can be performed offlinei.e. prior to the acquisition -using the power iterative method [Press ]. Importantly, all these computations can be done offline and prior to running the acquisition on the scanner.

Note that when Ψ defines an orthonormal basis, we get |||Ψ||| = 1. Nevertheless, SubProblem ( . ) makes use of an analysis-based prior which enables the use of overcomplete dictionaries and guarantees enhanced image quality at reconstruction at the cost of longer reconstruction times [Florescu , Cherkaoui ].

. . Parameters setting

Acquisition parameters

An ex-vivo baboon brain was scanned on a Tesla MR system (Siemens Healthineers, Erlangen, Germany) using a birdcage Tx/ Rx coil. All animal studies were conducted in accordance with the European convention for animal care and the NIHs Guide for the

. . Results

Retrospective Cartesian sampling

Prospective Non-Cartesian acquisition

The non-Cartesian case was also tested using -fold accelerated Sparkling trajectories as shown in

Y Ω k = [y Ω k ,1 , . . . , y Ω k ,L ] ∈ C kC×L .
As described in Chapter , the multi-channel reconstruction methods can be split in two classes. The first one models the coil sensitivity maps (S ) 1≤ ≤L where S ∈ C N ×N represents the th channel sensitivity profile. In these techniques, a single image is recovered from the combination of all channels [Chaâri , Florescu ]. The reconstruction problem is often formulated as follows:

x = arg min

x∈C N 1 2 L =1 F Ω S S x -y Ω S , 2 2 + λ Ψx 1 , λ > 0. ( . )
However, the coil sensitivities depend on the scanned subject, therefore they have to be calibrated for each subject. The calibration step can be performed beforehand [Uecker ,

El Gueddari b] and then some estimates ( S ) 1≤ ≤L can be injected in Eq. ( . ). Alternatively, the reconstruction can be viewed as a blind bilinear inverse problem [She ] where the optimization alternates between sensitivity profile estimation and image recovery steps. While the second approach is more computationally demanding, the first one does not really fit the constraints of online reconstruction. Indeed, as any shot is able to entirely cover the center of k-space, low frequency information associated with spatially smooth sensitivity profiles cannot be extracted in a straight manner. Moreover, a second weakness of this formalism lies in the fact that the gradient Lipschitz constant β k depends on the coil sensitivities and cannot be computed in advance if we assume that no calibration scan has been run beforehand. For the aforementioned reasons, we will use a domain-based calibrationless formulation that reads as follows: Here, we can introduce the straightforward extension of our online formulation to the matrix case. The online calibrationless reconstruction Subproblem ( . ) thus reads:

X = arg min X∈C N ×L 1 2 F Ω S X -Y Ω S 2 2 + R(ΨX). (
(∀k ∈ {1, . . . , S}), Xk = arg min

X∈C N ×L S 2k F Ω k X -Y Ω k 2 2 + R b-OSCAR (ΨX). ( . )
On the one hand, the gradient Lipschitz constant β k only depends on the sampling scheme Ω k , hence it can be computed offline and loaded at the beginning of each scan. On the other hand, the proximity operator of R b-OSCAR is explicit and can be computed efficiently [Zeng a, Eq. ].

The performance of this new multi-channel image reconstruction formulation and its ability to comply with online processing constraints will be discussed hereafter.

. . Parameters setting

Acquisition parameters

An ex vivo human brain was scanned using a Tx/ Rx (Nova Medical Inc., Wilmington, MA, USA), i.e. L = 32 with µm 2 in-plane resolution and a slice thickness of . mm.

Reconstruction parameters.

The decimated Symmlet Wavelet Transform with decomposition scales was used as sparsifying transform. As mentioned previously, the hyper-parameters (λ, γ) in Eq. ( .)

have been set retrospectively so as to maximize the SSIM score [Wang ] to the reference image once all data are available (i.e. offline scenario). Note that in a more clinically to guarantee that convergence was reached.

Computing parameters setting

All experiments were run on a machine with GB of RAM and an -core ( . GHz)

Intel Xeon E -v Processor. All the codes have been developed in Python using the PySAP package . The values of computing times have been obtained using epochs of Condat-Vú iterations each, the mean and standard deviation being summarized in Table . . Using multiple epochs permits to account for potential variability in computing times due to concomitant processes running on the workstation.

Table . : Computing time for one iteration estimated on epochs of iterations each using timeit for Sparkling with, first row: S = 34 (presented Fig. . -(b)) and second row S = 64 (presented Fig. . -(c)).

. . Results

As multi-channel online reconstruction is more demanding from a computational viewpoint, the way constraint ( . ) can be satisfied requires a larger number of shots. Because Algorithm assumes uniform batch sizes over iterations, we decided to explore a drastically less accelerated (i.

X∈C N ×L S 2k F Ω k X -Y Ω k 2 2 . ( . )
In this formulation, the computation time is dedicated to the evaluation of the gradient step as presented Table . . This permits the batch size deflation down to spokes for the acquisition schemes under consideration.

. Conclusion

Acceleration of MRI acquisition has been a major area of research for the last decade, especially in the high-resolution context where motion becomes critical. In this context, k-space data is collected in a segmented manner through multiple consecutive shots separated by times of repetition. Among the different acceleration strategies proposed in the literature (e.g. partial Fourier, parallel imaging, etc.), Compressed Sensing has been the most appealing as it enables larger acceleration factors without degrading the image quality at the reconstruction step. However, the price to pay lies in a long iterative process for image reconstruction.

Our online reconstruction approach relies on a mini-batch formulation which consists of aggregating the acquired shots in multiple subsets over time on mini-batches processing.

We take advantage of the period of time between the mini-batches acquisitions to launch channels using OSCAR-norm regularization. This calibrationless formulation presents the advantage of being well-suited to non-Cartesian sampling and does not require prior knowledge on the sensitivity profiles associated with the multiple channels. Importantly, the proposed online approach relies on convex optimization and its convergence to the global solution is guaranteed once all the spokes have been acquired.

The proposed framework was validated on single and multi-channel acquisitions at 7 T. We considered ex vivo imaging to avoid any motion-related concern. In the singlechannel acquisition setup, a retrospective Cartesian sampling was implemented using variable density over phase encoding lines with an acceleration factor of . . When the shots are split in small batches, the online reconstruction showed a significant gain in terms of speedup as the SSIM score of the image solution computed by the end of acquisition was very close to that of the final image delivered by offline reconstruction. This finding also holds for non-Cartesian acquisitions performed with a single-channel coil.

In the multi-channel acquisition setup, online reconstruction gets more challenging as both the computing time per iteration and memory usage increase. Hence, it gets harder to comply with the timing constraints using small batch sizes without significantly decreasing the number of iterations per mini-batch processing. We still observed some advantages for MR online reconstruction in -fold accelerated Sparkling acquisitions using a -receiver coil. In contrast, in highly accelerated acquisitions (i.e. -fold Sparkling) we showed that the current approach was not able to recover artifact-free images by the end of acquisition as only large batch sizes can be managed in this acquisition scenario. This . . Outlook intermediate conclusion calls for exploring new directions to overcome this difficulty. We already investigated an imperfect solution that is based on solving the data-fidelity term as the non-Cartesian Fourier transform involved in the gradient step is poorly conditioned and hence requires several iterations to converge to a descent image quality. Then the regularization term is added to the reconstruction as the problem is ill-posed thus needs regularization.

. Outlook

A first attempt for improving numerical efficiency of online reconstruction would consist of improving the initialization x 0 using a density compensated regridded solution [Fessler ].

A second direction that could be relevant to explore is related to coil compression. With the current trend of enlarging the number of channels in reception coils in order to continuously improve the input SNR (e.g. channel Nova coil now available on Tesla systems), coil compression could drastically limit memory usage and computing time [Buehrer ].

However, the question of optimal combination of coil compression with calibrationless reconstruction is still an open issue and there is no contribution in the literature that addresses this point. A third direction would be to consider stochastic implementation [Combettes ], where the global gradient is approximated by the sum of a subset of partial gradients.

Overall, the proposed framework enables to reduce the reconstruction times in CS MRI applications. Importantly, for practical implementation, it is compatible with the Gadgetron system [Hansen ] that provides the right interface to make those reconstruction techniques available to the scanner. A description of the Gadgetron system is provided Section C. . This may be really helpful to guide the physician for selecting further pulse MR sequences in order to sharpen his/her medical diagnosis.

General Conclusions and Perspectives

At the end of these three years, we hope we gave a boost to the clinical application of Compressed Sensing for high-resolution MRI. The aim of this thesis was not merely to propose an efficient reconstruction method for highly accelerated non-Cartesian acquisition particularly suited for high-resolution imaging, but also to speed up the reconstruction in order to bring those trajectories into clinical practice.

Context

Improving the spatial resolution in MRI while limiting the acquisition time has been a crucial challenge for the MRI community over the last two decades. Therefore, a major research field has emerged with multiple contributions dedicated to scan time accelerations.

In the early s parallel imaging has been proposed as a key step and first results were obtained using dedicated phased array receiver coils and deterministic undersampling,

Contributions & Limitations

High-resolution imaging requires a certain level of input signal to noise ratio (SNR), hence the acquisition is often performed with multi-channel receiver coils, to ensure the signal quality. Before the beginning of this work none of the reconstruction algorithm in the laboratory was considered for the multi-channel case, and the SPARKLING acquisitions Calibrationless methods have been preferred since they do not require any prior knowledge on the coil sensitivity profiles and do not involve any gridding step -which are known to deteriorate the final image quality. Consequently, the quest was to design a new regularization scheme that achieves at least the same results as the self-calibrating methods. The outcomes were beyond our expectations since calibrationless reconstruction, with a dedicated structured sparsity norms, outperforms the self-calibration techniques, both with respect to the image magnitude and more importantly regarding the recovery of phase information. The proposed calibrationless reconstruction method constitutes the second contribution of this thesis. Finally, the ultimate goal was to extend this calibrationless method to the online framework. However due to the computational load, the first online results were disappointing, consequently we investigated a faster solution for which only the data-fidelity term was considered and optimized over the scan time. Using this new strategy, we improved the image quality of the estimate obtained by the end of acquisition.

Online reconstruction methods definitely constitutes the last contribution of this thesis.

Besides those scientific contributions, the goal was also to propose a new package for image reconstruction. Hence with the help of CosmoStat laboratory, we have developed a new python package dedicated to the reconstruction named pySAP. The pySAP package was also adapted to be included within the Gadgetron project and thus enable the reconstruction to be performed in the scanner room and the results sent back to the scanner console.

Although the results of our calibrationless method seem promising, the set-up of hyper-

parameters remains an open question, constituting one of the main limitations of these methods. On top of these considerations, the major limitation remains the heavy computation load associated with the non-Cartesian Fourier operator. Although we used GPU accelerations for Non-Cartesian Fourier operations in d imaging we did not fully implement the whole algorithm in GPU. This might be helpful in the future, especially for d imaging combined with multiple receiver phased arrays, for which the increasing number

General Conclusions and Perspectives

of channels directly affects the computation cost.

Perspectives

Many more developments could be thought of using the calibrationless approaches proposed in this PhD thesis. For instance, structured sparsity norms could be adapted to deal with multi-contrast or dynamic MRI. Moreover, susceptibility weighting imaging (SWI) could also benefit from the good performances of calibrationless reconstruction, especially owing to the well preservation of phase information involved in the post-processing of SWI images. This work is currently addressed in the MANIAC project where patients and healthy volunteers are scanned on a 3 T MR system using SPARKLING acquisitions and the proposed calibrationless reconstruction method.

Regarding the computational aspects, the non-Cartesian Fourier transform could benefit from a faster implementation on GPUs. Also to avoid the time loss during the communication between the GPU and CPU, our calibrationless method could be made scalable by implementing it in a pure GPU framework such as the CuPy module.

Last, the reconstruction methods proposed in this thesis are inherently slow due to In general the Gadgetron reconstruction process (Gadget stream) consists of three main components: "Readers", "Writers", and "Gadgets". These components are assembled and controlled by a GadgetStreamController, which is also responsible for maintaining the connection with the client. The role of the Reader is to receive and deserialize data from the client, and the Writer serves the (opposite) role of serializing and returning data to the client. Once the data has been received, it is passed down the chain of Gadgets. Each Gadget can modify or transform the data before passing it on to the subsequent Gadget.

Adding Python modules

To make the transition from prototyping to deployment easier, the Gadgetron is designed so that the functionality of individual Gadgets can be implemented using Python as a scripting language. The Python scripting functionality can be accessed by the user through a PythonGadget, and the framework provides certain hooks so that the user can implement callback functions for in a Python module. The user specifies which Python module

General Architecture

The PySAP project is organized around a core architecture and plugins for end-user ap- operator, but also multi-channel data collected over phased array coils in parallel imaging.

In this context, both two-step self-calibrating and one-step calibrationless reconstruction methods have been developed in this plugin. For the sake of conciceness, we remind here that the former perform sensitivity matrix extraction prior to the reconstruction of a single full FOV image while the latter reconstructs as many images as receiver channels. Among all available features, a key aspect is to guarantee the reliability of the proposed classes and methods. Hence unit-tested have been written for each and every proposed operator as well as for the overall reconstruction methods. This makes continuous integration with Travis and CircleCI smoother. For instance, the non-uniform Fourier operator and its adjoint have been tested so as to check that the adjoint relation Eq. (C. ) is fulfilled.

Also, in the simpler case when the collected k-space samples lie on the Cartesian grid, we guarantee that the output of NFFT and FFT operators provide the same result.

Reconstruction en ligne

I

  would also like to express my gratitude to Dr. Alexandre Vignaud for his kind and regular advice. I am also really grateful to Dr. Nicolas Boulant, Pr. Jean-Christophe Pesquet and Dr. Emilie Chouzenoux for your advice during our insightful meetings. Thank you for your time, for sharing your knowledge with me. I deeply thank my COSMIC colleagues especially Dr. Jean-Luc Starck, Dr. Samuel Farrens, Zac Ramzi and last but not least Dr. Antoine Grigis for sharing their knowledge and their vision on reproducible science. I am thankful to Dr. Franck Mauconduit for his help on the Gadgetron project. Thank you Aurelien Massire, Chai GR, Mathilde Rippart and Carole Lazarus for your help on the MANIAC project. In addition, I also wanted to thank the jury for the interest they have shown in my work, and the time they have spent reviewing and examining the manuscript.

.c

  Correspondence between SSIM score ans Mean Opinion Scores (MOS) between an image and its reference . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Reported computation time for D-Cartesian reconstruction for different methods including: Sparse-BLIP, J-SENSE, Sparse-SENSE and IRGN-TV. . . . . . Ex-vivo baboon and in-vivo human brain data collected using various slice thickness and acceleration factors A . . . . . . . . . . . . . . . . . . . . . . . . SSIM values within the ( , ) range computed over three data sets. The larger the SSIM score the better the image quality. Bold font indicates the best score . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Acquisition parameters used for prospective CS ex vivo T * 2 acquisition on the same human brain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Performances of the least squares (LS) solution (zero-order reconstruction) on the dataset Table . . The results in red box are displayed in Figure . -. Reconstruction of -fold accelerated scans acquired with sparkling and spiral trajectories. Comparison of global regularization functions (LASSO, Elasticnet, g-OSCAR and k-support norm) with hyper-parameters chosen to as to maximize the SSIM score. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Quantitative comparison of the c-OSCAR and group-LASSO . . . . . . . . . . Quantitative comparison of the four implementations of OSCAR norm regularization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Quantitative comparison of the proposed calibration-less results to the stateof-the-art methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Computing time for one iteration estimated on epochs of iterations each using timeit. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Computing time for the different operator involved in a single iteration of the multi-channel reconstruction optimizer . . . . . . . . . . . . . . . . . . . . . . Speed of light in a vacuum 299 792 458 m • s -1 (exact) Planck constant (reduced) 1.055 × 10 -34 J • s k B Boltzmann constant 1.381 × 10 -23 J • K -1 γ H 1 H gyromagnetic ratio 267.5 × 10 6 rad • s -1 • T

agnetic

  Resonance Imaging (MRI) is the key imaging modality for exploring the most fascinating organ in the human body, namely the brain. With a huge diversity of contrasts, MRI can leverage multiple information to ease medical diagnosis. Besides diagnosis, longitudinal studies take advantage of MR non-invasiveness to understand the brain or predict the outcome of diseases using appropriate biomarkers. In order to improve accuracy in the diagnosis and help uncover structural abnormalities at an early stage, high resolution MRI could be beneficial at the expense of scan time. Long scans (> min) are more sensitive to subject's motion which impedes the resolution and increases the number of discarded exams. Although dedicated acquisition strategies combined with motion correction algorithms may reduce artifacts, the scan time remains too long to fit within the min dedicated to the whole clinical MRI exam, during which, multiple acquisitions are performed to collect complementary information on tissues. Consequently, specific acceleration and encoding strategies (e.g. parallel imaging techniques, navigator echoes, propeller scheme, Partial Fourier, etc.) have emerged over the last two decades to limit scan duration and gain in robustness to motion. Yet clinical imaging remains limited to millimetric resolution. During the last decade, the Compressed Sensing theory has made a major breakthrough to accelerate high-resolution scans. This theory combines efficient acquisition schemes and reconstruction algorithm. Recently, the MR community has seen the emergence of new efficient sampling patterns, first using Cartesian straight lines drawn pseudo-randomly, then with simple non-Cartesian sampling such as spiral or radial spokes that allow to reach higher acceleration factors compared to Cartesian acquisition. Even more recently, a new mathematically driven and hardware-adaptive sampling framework, called SPARK-LING, has been proposed. SPARKLING stands for Spreading Projection Algorithm for Rapid K-space samplING, it is the outcome of years of collaborative research between mathematicians (P. Weiss and J. Kahn, CNRS/ITAV Toulouse), physicists (A. Vignaud and F. Mauconduit, CEA/NeuroSpin), a researcher in signal processing for neuroimaging (P. Ciuciu, CEA/NeuroSpin) and former PhD students (C. Lazarus and N. Chauffert, CEA/NeuroSpin). The collaboration led to an efficient optimization algorithm which takes into account MR hardware constraints specific to the magnetic field gradients and a target sampling distribution. The outcome of this algorithm is a set of sampling trajectories whose temporal derivatives, namely the gradient profiles, can then be injected into a MR pulse sequence. SPARKLING has been proved to drastically accelerate scan time at Tesla by a factor of in D and for D acquisitions with a sub-millimetric isotropic resolution. Despite the tremendous acceleration provided by the latter trajectory, SPARK-LING acquisition are unsuitable for clinical routines. The very reason is that several issues related to image reconstruction are still pending in order to deliver reliable MR images shortly after the end of acquisition.

  Outline . Nuclear Magnetic Resonance phenomenon . . . . . . . . . . . . . . . Magnetization. . . . . . . . . . . . . . . . . . . . . . . . Excitation and Relaxation . . . . . . . . . . . . . . . . . . Image formation. . . . . . . . . . . . . . . . . . . . . . . . . . Spatial Encoding . . . . . . . . . . . . . . . . . . . . . . Reception: Detection, Conversion & Demodulation . . . . . . . . . Reconstruction . . . . . . . . . . . . . . . . . . . . . . Multi-channel receiver coil and Parallel Imaging . . . . . . . . . . . . . Signal-to-Noise Ratio increase . . . . . . . . . . . . . . . . . Parallel MRI: first MR accelerations . . . . . . . . . . . . .

  and . . . It is only three decades later that the works of Paul Lauterbur [Lauterbur ] and Sir Peter Mansfield [Mansfield a] allowed to encode the Nuclear Magnetic Resonance (NMR) signal spatially, thus enabling NMR imaging, or MRI. For this, they shared the Nobel Prize in Medicine and Physiology in . Spatial encoding is the subject of Section . .

Figure

  Figure . -: Effect on macroscopic magnetization of a B + 1 field rotating at the Larmor frequency ω 0 viewed (a) from laboratory frame (XYZ) where B 0 is aligned with the Z axis. (b) It is generally convenient to represent it in a rotating frame of reference (xyz) of frequency ω rf , where B + 1 lies along with the x axis. The z axis is aligned with Z. From [Tomi-Tricot ].

Figure

  Figure . -: Free induction decay signal representation, showing the signal that oscillates at the carrier frequency (128 MHz at 3 T), weighted by a decreasing exponential caused by T * 2 decay. No values are shown on the time axis, as in this simple representation, for the sake of readability, the frequency is orders of magnitude lower than in reality.

Figure

  Figure . -: Diagrams of the B z component of the D magnetic field variation of a) G x -, b) G y -, and c) G z -gradients. The field is aligned with the direction of the arrows. The scale of the arrows indicates the magnitude of the magnetic field strength at that point. From [Poole ].

  signals in MR imaging arise from multiple voxels recorded all at once. These voxels contain different materials with different spin densities and relaxation times. In order to discriminate spatial location of the different voxels Lauterbur and Mansfield [Mansfield b] have proposed to change spatially the magnetic field, making the Larmor frequency a function of space. The local changes of the magnetic field is due to the use of three gradients (G x , , G y , G z , one for each spatial dimension as shown in Figure . -).

Figure

  Figure . -: Encoding the k-space line by line. From [Le Ster].

G

  Phase . Hence , to acquire a full d k-space, we have to repeat the G Part → G Phase → G Read pattern N Part ×N Phase times, with N Part the number of partitions to be encoded and N Phase the number of lines. As shown later on, the imaged volume will be recovered by applying an inverse d Fourier transform. One drawback of d imaging is its sensitivity to patient motion, as the whole d kspace must be filled -which takes time -before reconstruction. Any spatial information corrupted by movement occurring during acquisition affects the quality of all slices within the volume. In comparison, in d imaging, each slice is reconstructed separately: motion during the acquisition of one slice has no effect on the rest of the FOV.

  Figure . -illustrates the reception chain. After the low-noise amplification, the signal is split into two channels. In the first one, the signal is multiplied by an signal in-phase (cosinusoid) and in the second one by a signal in quadrature (sinusoid). The in-phase component leads to the real part of the MR signal while the in-quadrature component is considered as its imaginary part. Before being digitally converted the signal is low-pass filtered to get rid of any aliasing artifacts due to digital conversion. Then each channel is converted separately and recombined as real and imaginary components of k-space samples.

Figure

  Figure.-: Block diagram of the reception chain. The signal from the coil is amplified by a low-noise amplifier before it gets split into two channels. Each channel will be multiplied by an in-phase and in quadrature signal generated by a frequency synthesizer, leading to a complex multiplication of the signal. The low-pass filter corresponds to an anti-aliasing step which prevents aliasing before the Analog-to-Digital conversion with a frequency being at least twice the maximal sampling frequency provided by the ADC. Adapted from[Décorps ].

Figure

  Figure.-: Grid-driven interpolation, the value at any Cartesian grid point (+) is computed as the interpolation (e.g. bilinear) of the four surrounding non-Cartesian points (o). The interpolation kernel is a user-defined parameter of grid-driven methods, tat can be either linear, cubic, or based on the nearest neighbor algorithm.

  Eqs. ( . ) and ( . ). This library uses only O |I N | log |I N | + | log ε| d M instead of O (M |I N |), floating point operations, with ε being the desired computation accuracy of the approximation. The key idea of the library is to use standard FFTs (more precisely Background in Magnetic Resonance Imaging the FFTW [Frigo ]) in combination with an approximation scheme that is based on a windowing function ϕ. The latter must be mutually well localized in spatial and frequency domains. Two parameters control the accuracy of the NFFT: the oversampling factor σ and a truncation parameter m.

  Figure . -: (a): Single channel birdcage coil used on the 7 T scanner at Neurospin. (b)channel array coil in transparency, the different channels are color-coded (from [Paolini ]).

  Figure . -: Comparison of , , and -channel phased array coils at 3 T. (a) SNR map for acquisitions with the different receiver coils. Large number of channels with small diameter are more sensitive to the region next to the coils. (b) SNR profile along the central line in the brain for receiver coils with different number of channels, while the SNR gain is less significant in the middle of the brain. From [Wiggins ].

  Figure . -: (a) Four coil phased array showing overlapping coils, (b) their associated sensitivity maps. From [Gruber ].

  Figure . -: (a) Fully sampled Cartesian k-space. (b) Acquired k-space for parallel imaging with an acceleration factor R = 2. Dashed lines represent the missing samples while solid lines correspond to the acquired k-space values. Horizontal axis represents the frequency encoding direction and vertical axis is the phase-encoding one.

Figure

  Figure . -: Comparison of the image quality of SENSE and GRAPPA reconstructions with an acceleration factor R = 3. Left: SENSE reconstruction with accurate coil sensitivity extraction. Right: GRAPPA reconstruction. Visually, no difference can be seen between both methods. From [Blaimer ].

Figure

  Figure . -: Comparison of artifacts in SENSE and GRAPPA reconstructions at high acceleration factors R (chosen close to the number of coils). In this example, R = 4 leads to a four-fold acceleration of the scan time for an equal number of coils (L = 4). Left: The SENSE image shows a local noise enhancement due to non-ideal conditioning for the reconstruction. Right: The noise enhancement in GRAPPA is distributed more evenly over the FOV. Additionally, aliasing artifacts can be seen due to inaccurate calculation of missing k-space lines. From [Blaimer ].

  strated yet . Due to the mentioned impediments of the 0 -pseudo norm, problem ( . ) is often relaxed and its tightest convex envelop [Boyd ] is considered insteadi.e. the tightest convex hull of the 0 -pseudo-norm is the 1 -norm x 1 = n j=1 |x i | -leading to the following optimization problem: z = arg min z∈C p z 1 subject to y = Az. ( . ) Based on the restricted isometry property (RIP), which characterizes matrices nearly orthonormal, the first results of the Compressed Sensing theory [Candès b, Donoho , Candès ] guarantee exact recovery of s-sparse signals. The RIP conditions state that

  with Ψ is the sparse representation such as wavelet transform (synthesis operator) and Ψ H its adjoint (conjugate-transpose as defined Page xv) operator that performs image reconstruction from a set of wavelet coefficients. This is often called the "coherence barrier"[Adcock ].

  (n)). Hence the required number of measurements to get an exact reconstruction in MRI gets closer to O(s log(n) 2 ). Although Theorem . . allows to partially break the coherence barrier, if we look closer at the numerical values for an image size of 356 × 256 we get log(n) 2 123, which makes the above result less attractive. If we assume on top that C ≥ 1 and s 0.2 • n, then m gets larger than n. However, in practice we may observe at least on simulations that a number of measurements m ≈ 2•s provides excellent reconstruction results [Adcock ]. The reason is related to the underlying structure of the signal, more precisely its structured sparsity.

Figure . -

 . Figure.-illustrates the structured sparsity pattern of brain images decomposed within the Symmlet basis with J = 3 decomposition scales, each scale being composed of four subbands. Normally, the scale is inversely proportional to the spatial frequency.So, j = 1 refers to the high-frequency details and corresponds to the subbands depicted in the bottom row and top right corner of Figure.-, while j = 3 is associated with low-resolution information and contains both the low-frequency approximation as well as the three detail sub-bands located both in the top left corner of Figure.-.Unfortunately, for many applications -including MRI -the sampling matrices cannot be constructed at random and the collected measurements are dependent one another.Therefore, concrete applications of CS are often based on sampling strategies that may strongly deviate from theory. Despite having no solid theoretical foundation, these heuristic strategies work very well in practice. Recently[Boyer ], derived CS theoretical results for structured acquisitions and structured signals (i.e. structured sparsity), with a recovery probability that explicitly depends on their support. The main contribution can be formulated informally as follows:

Figure

  Figure . -: Wavelet decomposition of a brain phantom on scales (J = 3). The lowfrequency or higher scale of decomposition are less sparse than the lower or high-frequency scales. This illustrates the structured sparsity pattern of brain images. From [Lazarus a].

  ments, called shots. The latter sequentially fill the considered k-space, either on the grid for Cartesian acquisitions (spin-warp) or out of the grid (radial, spiral) for non-Cartesian sampling. Therefore the global acquisition time (TA) in segmented acquisitions is equal to the product of the number of shots (n s ) and the value of TR, the time of repetition: TA = n s • TR. At each and every TR, a radio-frequency (RF) pulse is used to tip the global magnetization for a given slice ( D imaging) or volume ( D imaging) and then a new shot is collected during T obs [Bernstein , Chap. . ].

  points in D [Lustig ]. Hence, this strategy can only be implemented on D acquisitions as illustrated in Figure . -. Nonetheless this paradigm has been intensively used in retrospective d CS validations.Its advantage is threefold: i) Poisson-disk actually achieves locally uniform coverage as the k-space samples stay far apart from each other, the distance between them being tuned by the disk radius; ii) it actually fully samples the central region of k-space [Levine , Vasanawala ];

Figure

  Figure . -: D acquisition with d Poisson disk sampling. From [Chauffert ]. The readout is performed along the horizontal lines wih black arrows, namely orthogonally to the d cross-section along which isolated points have been drawn.

  Figure . -. However, both D VDS on Cartesian lines and D Poisson disk sampling are hampered by the same drawback: they do not exploit the full dimensions of the ambient space. So their performances are limited in terms of acceleration factor. To go one step further, we summarize in the next part the non-Cartesian sampling strategies that fully implement VDS strategies in all dimensions of the native space, D or D VDS in respectively D and D imaging. The major difficulty lies in the ability to be compliant with the gradient constraints.

Figure

  Figure . -: An example of reconstructions of a × MR image from block of measurements. (a) Sampling pattern with horizontal lines -i.e. variable density in the phase encoding direction -with an under-sampling factor m/n ≈ . (b) Corresponding reconstruction using the 1 -minimization. (c) A zoom on a part of the reconstructed image. (d) Image obtained by using the pseudo-inverse transform. (e) A zoom o a part of the image. From [Boyer ].

  Figure . -: Examples of some non-Cartesian trajectories used within CS theory to accelerate the scan time(a) Originally proposed by [Lauterbur ], Radial trajectory inherently follows a variable density. (b) Spiral trajectories are one of the most efficient curve as it can be implemented so as to take fully advantage of the encoding gradients (i.e. saturation of S max in the k-space center and G max in the periphery). (c) TWIRL: the first combination of radial and spiral trajectories for D acquisitions [Jackson ]. (d) TWIST: the D extension of the TWIRL method [Boada ]. (e) FLORET: recent advances on the design of D non-Cartesian sampling patterns based on the combination of radial and spiral.

Figure

  Figure . -: Generated SPARKLING trajectories for different initializations and target densities. Top tow: The SPARKLING method applied to Cartesian lines with uniform density. Mid row: The SPARKLING method applied to radial spokes with variable radial density. Bottom row: The SPARKLING method applied to the centered-out Archimedean spiral initialization with variable radial density. For each initialization the Point Spread Function (PSF) gives a rough idea on the artifact structure for the under-sampled factor. From [Lazarus a].

Figure

  Figure . -presents the outputs of SPARKLING algorithm for different initializations, in particular the Point Spread Function (PSF) gives an idea on the performances of the different sampling patterns. While for the Cartesian sampling presented in Figure .the spikes in the PSF predict the aliasing artifacts associated with under-sampling, the SPARKLING output is not affected, i.e. its PSF presents a pure central peak without any small replica.

Figure

  Figure . -: Comparison of different acquisition strategies on an T * 2 ex-vivo human brain acquired at 7 T for various acceleration factor AF = 10, 15 and 20. From [Lazarus a].

A

  basic known result states that if Ψ is an orthonormal transformi.e. Ψ H Ψ = ΨΨ H = I N and Ψ a Ψ s = I N ψ with necessary N ψ = N -then both analysis and synthesis formulations Eqs. ( . ) and ( . ) are equivalent. The same property holds more generally for invertible transforms Ψ. However, in other situations, the two formulations provide different solutions [Elad ]. Noticeably, for over-complete transformsi.e. a transform that maps the vectorized image x ∈ C N to a coefficient vector z ∈ C N Ψ with N < N Ψ -the analysis formulation is preferred as it has been shown in [Majumdar b, Cherkaoui , Selesnick ]. The choice of an appropriate fixed decomposition has been intensively studied in the literature [Selesnick , Majumdar b, Cherkaoui ]. Importantly, it has been proved that over-complete dictionaries such as curvelets [Starck ],

Figure

  Figure . -: Set of atoms learned from the CamCan database [Taylor b] where only magnitude information has been learned. From our work published in [Carrié ].

  )); end One of the most popular techniques associated with Nesterov acceleration scheme is the Fast Iterative Soft-Thresholding Algorithm (FISTA) [Beck ] which has a convergence rate of O 1/k 2 . The convergence proof of the iterates of FISTA haven't been demonstrated yet, however [Chambolle ] proposed a variant of FISTA with a convergence proof of the iterates. On drawback of the Nesterov acceleration schemes is the presence of ripples in the trace of the objective function which slows down the convergence as illustrated in Figure . -. In [O'donoghue ] the authors exhibit two regimes depending on the value of the momentum and noticed that the ripples arise when the momentum term was too high. Therefore, by simply restarting the momentum termi.e. setiting it to zero -they could accelerate the different algorithms. The methods proposed by [O'donoghue ] is a very general and could potentially be applied to different optimization algorithms that present the same behavior.

  Figure . -: Comparison of the convergence speed of different proximal algorithms to solve Eq. ( . ). The Nesterov ripples can be seen from the FISTA-CD proposed by [Chambolle ]. From [Ramzi ].

Figure

  Figure . -: A view of the maximum under-sampling factor allowing SSIM scores above . as a function of image size (or resolution as the FOV is fixed) and input SNR. Orange line: Available input-SNR on the 7 T scanner (Siemens-Healthineers, Erlangen, Germany). [Lazarus c].

First

  let us consider methods that extract the sensitivity information beforehand and use it as a prior in the optimization of Eq. ( . ). As originally proposed by[McKenzie ],extra k-space central lines can be added to the accelerated scan and used to extract lowresolution coil images which are then used as sensitivity profiles. While those lines can easily be added to Cartesian acquisitions, it is less obvious for non-Cartesian sampling. However because most of non-Cartesian schemes (such as radial [presented in Figure.-) are heavily sampled in the k-space center region,

  Figure . -: Examples of non-Cartesian trajectories: (a) radial sampling, (b) spiral sampling, (c) PROPELLER sampling. Oversampling the center of k-space brings a gain in robostness to the patient's motion.

  of a calibration matrix. First let us introduce how the calibration matrix A is built up. Considering a calibration region AC with non-missing k-space values y AC Ω , the method proposed by [Uecker ] constructs a block-Hankel matrix A (also called calibration matrix). Block-Hankel matrices are known to possess well-defined subspaces [Heinig ] which in this case can be used to extract sensitivity information. The construction of matrix A is based on data concatenation along a D sliding window that moves over the AC region and across all channels as illustrated in Figure . -.

Figure

  Figure.-: Left: Set of k-spaces obtained within the multi-coil setting. Right: the calibration matrix A is obtained by sliding a window (small green frame) on the calibration data y AC Ω (blue frame).

  Consequently, they have proposed to simplify the eigenvalue decomposition of the operator W by directly applying it in the image domain. While the ESPIRiT method represents an elegant way of extracting the sensitivity matrices, this method remains suited for Cartesian sampling due to the calibration matrix A. For non-Cartesian sampling, the construction of this matrix requires a gridding step which might degrade the extraction and thus the final image quality. While fast estimation of coil sensitivity profiles could be thought of in the context of online reconstruction, the Lipschitz constant of the forward model depends on the estimated sensitivity profiles. Hence, it cannot be computed beforehand. This constitutes a major bottleneck for online MR image reconstruction.

  the original paper [Uecker ], the authors have proposed to decrease the value of α and β parameters over iterations, leading at the end to the minimization to a nonregularized, ill-conditioned problem. To overcome this issue, [Knoll ] have recommended to fix it during the minimization. Based on a previous work [Knoll ], they have also proposed to replace Tikhonnov regularization R im in the image domain by total variation (TV) [Rudin ], which promotes piece-wise constant images (cf Eq. ( . a)) or its generalization, namely the total generalized variation (TGV), which assumes non piece-wise constant images (cf Eq. ( . b)). This eventually leads to the following penalty terms: TV: R im (x) = α ∇x 1 ( . a) TGV: R im (x) = α min v ∇x -v 1 + ∇v + ∇v H 1 ( . b)

  Sobolev regularization on the coil sensitivity profiles, [She ] has proposed the Sparse BLind Iterative Parallel imaging reconstruction (Sparse-BLIP) method based on TV regularization of the sensitivity profiles as defined in Eq. ( . b). On top of this new regularization, following the seminal paper in sparse MRI reconstruction [Lustig ], they have also regularized the image x by combining its 1 -norm in the wavelet domain (i.e. Ψx) with the TV penalty in the image domain as shown in Eq. ( . a):

.

  . Coil-by-coil MR image reconstruction gence of algorithms make them less appealing for online image reconstruction purposes..Coil-by-coil MR image reconstructionWhile multi-channel MR image reconstruction methods based on coil sensitivity profiles are well suited for non-Cartesian acquisition, the extraction of sensitivity profiles is a critical step with a severe impact on the final reconstructed image. In some cases, such as the Cartesian scan with calibration lines, the estimation problem is easier. However, it gets harder when no-calibration data is available or for non-Cartesian acquisitions. In the last case, blind reconstruction methods have been proposed to estimate both the final image and the coil sensitivity profiles, using priors such as sparsity for the combined image and smoothness for the sensitivity profiles. However, due to the non-convexity of the underlying objective function, blind bi-linear methods might have unexpected artifacts depending for instance on the initial estimate of the coil sensitivity profiles. In what follows, we introduce coil-by-coil MR image reconstruction methods. The latter can be split in two different categories depending on the considered optimization domain either the optimization lies in the k-space or in the image or sparse transform domain. Since those methods reconstruct one image per coil, for comparison purposes we will use the square root of the sum-of squares to combine all coil-specific images. . . K-space driven methods promoting data low-rankness K-space methods were originally motivated by the good performances of GRAPPA [Griswold ] reconstruction (presented in Section . . ). While early methods are based on the robust extraction of an interpolation kernel, k-space driven reconstructions have met a shift in their principle. In fact, instead of finding an explicit interpolation kernel to generate the missing samples, state-of-the-art methods enforce the low-rankness of k-space data. This was made possible due to the following assumptions: i) All the sensitivity maps are band-limited hence they can be represented by a small number of Fourier coefficients. Because of the forward model (element wise multiplication in the image domain) this corresponds to a convolution in the k-space with a small kernel.

  a regularization step that enforces the low-rankness of Id C . While the previously mentioned method requires a calibration region to extract the kernel C, recent reconstruction algorithms such as SAKE [Shin ] or p-LORAKS [Haldar ] , have relaxed this assumption and proposed to use a low-rank constraint over a constructed block-Hankel matrix. For the sake of simplicity, we will rely on the SAKE formalism. The construction of this matrix is based on the vectorization of a sliding window that spans the whole k-space. This matrix is actually constructed like the ESPIRiT calibration matrix (presented in Figure . -) however the sliding window is not restricted to the calibration data (blue frame in Figure . -) but covers the entire k-space.If we denote by H the linear operator that generates the block-Hankel matrix A from the k-space Y , and H H its adjoint, then the reconstruction problem solved by SAKE (or p-LORAKS) is defined as:

  R the regularization term that promotes the joint-sparsity. Dictionary Ψ can be either the identity matrix, a fixed sparsifying decomposition (e.g. wavelet transform) or any data-driven dictionary. In their work, [Trzasko ] have proposed to use a locally low-rank formulation for image reconstruction called Calibration-free locally Low-rank EncourAging Reconstruction (CLEAR). They have introduced the use of a patch-based nuclear norm in the image domain to solve the multi-channel reconstruction problem. The nuclear norm A * of a matrix A corresponds to the 1 -norm of the singular value of A. The CLEAR method actually defines Ψ as the image patch selection operator P s,p requiring two other hyper-parameters to be set: the stride s and the patch size p, hence the regularization reads as follows: R (ΨX) = λ P s,p X * . The operator P s,p is illustrated in Figure . -. In the original work, [Trzasko ] have proposed to use a stride equal to the patch size (i.e.

Figure

  Figure . -: Illustration of the P s,p X operation performed by the CLEAR method, implemented on our ex vivo human brain data set.

  Table.: SSIM values within the ( , ) range computed over three data sets. The larger the SSIM score the better the image quality. Bold font indicates the best score.The reconstruction proposed by the Sparse-SENSE is very similar to the Cartesian reference scan, while the 1 -ESPIRiT method presents a smooth background with vertical lines, particularly visible on the zoom depicted in Figure.-(h). Regarding the coilby-coil p-LORAKS reconstruction method the results seems affected by high frequency noise, although no blurring effect is visible. The reconstruction method is hampered by the gridding step required to build the calibration matrix A. This can be seen when comparing the output result of p-LORAKS with the SOS of the gridded k-space displayed in Figure.-. Both figures present the same noise pattern and p-LORAKS does not improve the reconstruction. During this experiment two major impediments were reported regarding p-LORAKS: the first one reflects its poor performance on non-Cartesian acquisitions while the second one concerns the construction of matrix A for high-resolution data collected over a large number of coils (e.g. L = 32) which requires a huge memory load. Hence given our experimental set-up and the code provided by the authors [Kim b],

  Figure . -: Comparison of different parallel MRI reconstruction methods on a prospectively -fold accelerated Sparkling acquisition on the ev-vivo baboon brain ( mm slice thickness). (a) Cartesian reference, (b) Sparse-SENSE reconstruction with the proposed sensitivity maps extraction method, (c) 1 -ESPIRiT reconstruction (d) p-LORAKS solution after the gridding step and (e) CLEAR reconstruction and (f-j) their respective zoom over the red frames.

  Figure . -: (a) Square root of the Sum-Of-Squares (sSOS) of the gridded k-space reconstruction, (corresponding to the result before p-LORAKS reconstruction), (b) zoom in the red frame.

  Figure . -: SSIM score of the sSOS results for the CLEAR method (a) for of = 2 (corresponding to a stride s of 8 × 8 and a patch size p of 16 × 16) and (b) for of = 4 (corresponding to a stride s of 4 × 4 and a patch size p of 16 × 16).

  ded to various other contrast such as T 1 -weighting [Lazarus a, Chapter . ] and also to D imaging with a direct application to susceptibility weighted imaging (SWI) [Laz-

  larization terms on the final image quality and propose an extended comparison to other domain-based reconstruction methods, mainly, Let us first define the general calibration-less multi-channel coil image reconstruction problem. We recall n being the resolution of the sought image per channel and N = n d the image size where d is the dimension , L the number of channels used to acquire the NMR In D imaging d = 2, whereas in D imaging d = 3. . . Calibration-less reconstruction: general problem statement signal and M the number of k-space measurements per channel, with M N . For the sake of compactness, we denote the complete data set Y = [y 1 , . . . , y L ] ∈ C M ×L which stacks all vectors y ∈ C M , the latter gathering the k-space samples collected in the th channel.

  properties as variable selection operator and regularization term as it has been shown [Tibshirani ] via the Least Absolute Shrinkage and Selection Operator (LASSO). In the context of Compressed Sensing the usage of such norm has often been justified as being the tightest convex envelope of the cardinality function (i.e. the 0 pseudo-norm) [Lustig , Donoho , Candès ]. When the sampling pattern presents large incoherence

  this setup, the regularization tends to select one of those variables and discard the others, whereas one can expect to treat all the correlated variables the same way. To overcome this issue many regularization terms have been proposed, such as Elastic-Net [Zou ], group-LASSO [Yuan ], or based on clustering methods such as octagonal shrinkage and clustering algorithm for regression (OSCAR) [Bondell ].

  of clustering those variables before the selection/shrinkage operation. This idea has been exploited in a new regularization called Octagonal Shrinkage and Clustering Algorithm for Regression (OSCAR) [Bondell ]. In [Argyriou ] the authors proved a link between the OSCAR norm and the Ordered Weighted 1 (OWL) norm [Zeng b, Bogdan ]. Appendix B recalls the derivation of this link and the computation of its proximity operator.

  . ) with λ, γ two positive hyper-parameters. The 1 -norm term in Ω λ,γ promotes the sparsity of z while the second term, corresponding to a pairwise ∞ -norm, encourages the equality of each pair of entries in z. As pointed out in [Zeng a, Sec. II. A.], OSCAR norm has a closed relation with the OWL norm defined below: Θ w (z) = p j=1 w j |S p (z) j |. ( . )with w ∈ R p + a vector of hyper-parameters such that w 1 ≥ • • • ≥ w p ≥ 0. More precisely, OWL and OSCAR become equivalent if one sets the OWL weights as follows: w j = λ + γ(p -j) for j = 1, . . . , p.

Figure

  Figure . -: 1 , Elastic net, Ridge and OSCAR ball in R 2 .

Figure

  Figure . -: OSCAR ball in R 2 for various pairs of (λ, γ).

Figure

  Figure . -compares the unit ball of the 1 , Elastic-net and OSCAR norms. When γ is small, the OSCAR norm corresponds to the 1 norm. In contrast when λ → 0 the OSCAR norm tends to the pair-wise ∞ norm, Figure . -depicts the unit-balls for different pairs of (λ, γ) hyper-parameters. Note: in the D case the pair-wise ∞ -norm correspond to the global ∞ norm.

Finally

  we propose to apply OSCAR norm to each wavelet coefficient separately and thus to each vector z csk,: , where z csk,: gathers the L wavelet coefficients across all channels for coefficient k in subband s of scale c. The coefficient-wise (c-OSCAR) regularization thus reads: γ(L -)) |S L (z csk,: ) | , with S L , the L-dimensional magnitude decreasing sorting operator. This formulation is the closest to the usual application of the group-LASSO penalty [Chun , Majumdar a]

  [McDonald ] generalized the k-support norm using the box norm, and derived the same computation for the norm. Using the variational formulation . . , one can notice that the k-support norm for k = 1 is equivalent to the 1 -norm (unit ball with p = 3 is presented Figure . -a), and for k = p it is equivalent to the 2 -norm (unit ball with p = 3 is presented Figure . -c). When k / ∈ {1, p}, its unit ball differs from a linear combination between the 1 and 2 norm, i.e. elastic-net regularization. Taking for instance p = 3, Figure . -(a)-(c) present the unit balls for k = 1, 2, 3 Figure . -(d)-(e) depicts elastic-net and OSCAR norms. A comparison with Elastic-net unit ball shows that the 2-support norm has more rounded corners. Relation to elastic-net. [Argyriou ] shows that the k-support norm provides tighter relaxation bound compared to elastic-net by exactly a factor of √ 2. For prediction purposes, [Gkirtzou ] demonstrated good performances of the k-support norm compared to the elastic-net regularization in functional MRI data analysis (supervised learning task).

  Figure . -: Top row: Unit balls of the k-support norm in R 3 for different k values. (a) presentq the unit ball of the -support norm in R 3 , which corresponds to the LASSO regularization, (b) represents the unit ball of the -support norm, (c) shows the -support unit ball which correspond to the ridge regularization. Bottom row: Elastic-net regularization and unit ball of OSCAR norm in R 3 .

  assess the performances of the different regularization penalties in terms on MR image quality, we used prospectively accelerated non-Cartesian real acquisitions and proceed in four steps. First, we compared the global regularization regardless of the wavelet decomposition structure, in particular we compared the LASSO Section . . , Elasticnet Section . . , the global version of OSCAR (g-OSCAR) and the k-support norm regularization. A particular attention was paid to the performances of the ∞ vs 2 norm assumptions behind respectively OSCAR and k-support norm. Second, we made some comparisons between the coefficient formulation of OSCAR (c-OSCAR) and group-LASSO (i.e. CALM method [Majumdar a]), since both rely on the same assumptions in terms of structured sparsity. Third, we compared the four OSCAR-based formulations, described in Section . . . Finally, an overall comparison was conducted between selfcalibrating and calibration-less state-of-the-art methods, such as 1 -

  spectively under-sampled k-space in vivo data. However in this context their robustness to actual acquisition set ups is overlooked. Spiral trajectories[Lee ] were generated for different under-sampling factors as detailed in Tab. . . More advanced sampling schemes were also designed using Sparkling (Spreading Projection Algorithm for Rapid K-space samplING) [Lazarus a, Boyer ]. Sparkling method generates physically plausible trajectories with improved sampling efficiency and robustness to gradient imperfections and lower sensitivity to off-resonance artifacts as demonstrated in [Lazarus a]. The resulting sampling schemes are known to reach higher image quality for a given scan time, compared to state-of-the art trajectories (e.g. spiral or radial in D imaging). We explain hereafter while in Tab. . we specifically achieve in the non-Cartesian setting an higher acceleration in time as compared to the under-sampling factor expressed in the number of measurements over the number of image pixels. We recall the difference between acceleration factor (AF) in time and under-sampling factor (UF) in the number of measurements: AF = n/n c where n c is the number of shots. For instance, we used n = 512 and varied n c between down to corresponding to acceleration in time from to , as reported in Tab. . . In the same time, an oversampling (with a factor of ) along each shot was performed (n s = 3, 072 samples per shot) as usually done in non-radial non-Cartesian trajectories, hence the under-sampling factor defined as UF = N/M with M = n c × n s varied in a more limited range from . to . . The other acquisition parameters were set as follows: FOV = 200 × 200 mm 2 , TR = 550 ms (for slices), TE = 30 ms, BW = 100 kHz, T obs = 30.72ms (long readout, i.e. each shot is traversed in T obs ) and FA= • . In Tab. . ), we summarized the different data sets we considered by varying the under-sampling or acceleration factor for both Sparkling Table . : Acquisition parameters used for prospective CS ex vivo T * 2 (trajectories are displayed Figure . -). A fully sampled Cartesian dataset was acquired and the corresponding image was reconstructed using inverse FFT for each receiver coil and the square root of the sum of squares was eventually computed over all channels to form a single image Figure . -(a) which will serve as ground truth using the same sequence parameters (matrix size: N = 512 × 512 or n = 512).

  Figure . -: Non-Cartesian in-out variable density sampling schemes used in prospective CS at Tesla. (a): -fold accelerated-in-time Sparkling trajectories ( shots). (b): fold accelerated-in-time spiral trajectories with the same number of shots and measurements as Sparkling. For details about their generation, see [Lazarus a] for Sparkling and [Lee ] for spiral imaging. The trace of a single in-out shot is highlighted.

  Figure . -. Despite the high image quality metrics for Spirals compared to Sparkling (for same acceleration factor), Spiral MR reconstructions are contaminated by severe off-resonance artifacts (cortical ribbon highlighted with red arrow in Figure . -(c)) whereas Sparkling ones are not. These results illustrate the much better point spread function (PSF) associated with Sparkling sampling schemes as shown in [Lazarus a] and highlighted Figure . -.

  Figure . -: Top: (a) Cartesian reference image. Zero-order Fourier reconstructions based on the minimization of the least square data consistency term for -fold accelerated-in-time ( shots) (b) Sparkling ( shots) and (c) spiral trajectories. Bottom: Respective zoom over the red frame. Red arrows show the presence of gradient-related artifacts.

Figure

  Figure . -: Evolution of the SSIM score computed between the reconstructed image fromfold accelerated Sparkling trajectory using Elastic-net regularization and the reference image. The black cross displays the hyper-parameters that maximize the SSIM score.

  Figure . -: Reconstructed MR images (square-root of the sum-of-squares) from -fold accelerated Sparkling (first two rows) and spiral acquisitions (last two rows) using different regularization. (a) Cartesian reference. (b) inverse Fourier Transform. Reconstruction based on (c) group-LASSO regularization R group-LASSO . (d) coefficient version of OSCAR R c-OSCAR . Second and Forth rows: Sparkling and Spiral respective zooms in the red frame .

  and p-SNR scores, CALM regularization outperforms all other techniques. However, CALM performances remain below those of OSCAR, whether it is the subband-or coefficient-wise version. This confirms that extending the group-LASSO penalization using a pairwise ∞ -norm between channels instead of the global the 2norm brings more flexibility to account for varying SNR across channels. Last but not least, in Figs. . --. -we show that calibration-less reconstruction techniques such as OSCAR and CALM better preserve phase information in comparison with 1 -ESPIRiT, AC-LORAKS. In this comparison we also provide the phase information related to each reconstructions. As there is neither sensitivity map extraction nor gridding step in domainbased calibration-less approaches, any phase-related error (e.g. phase shift or wrapping) in those steps does not propagate in the image reconstruction itself. Such artifacts are particularly visible in 1 -ESPIRiT for spiral imaging (see red arrows in Fig. . -). Phase images based on AC-LORAKS are also corrupted by artifacts that may be due to the . . Conclusion gridding operation, especially for Sparkling imaging (see red arrows in Fig. . -).

  Figure . -: Top: Reconstructed MR images (magnitude) from -fold accelerated Sparkling acquisitions using different methods. (a) Cartesian reference. (b) Zero-order reconstruction obtained by solving the least square criterion (data consistency term) involving the NFFT operator. (c) Calibration-less regularized reconstruction based on the subband-wise OSCAR formulation. (d) Calibration-less regularized reconstruction based on k-support penalization. (e) Self-calibrating 1 -ESPIRiT reconstruction. (f) Auto-calibrated (AC) LORAKS reconstruction. Middle: Respective zooms in the red frames. Reconstructed MR images (phase). Bottom: Wrapped phase images extracted using the virtual coil approach.

  Figure . -: Top: Reconstructed MR images (magntiude and phase) from

  Nowadays, real time MR image reconstruction is mainly implemented for dynamic imaging and more precisely for cardiac MRI. Indeed, the idea of real-time MR image reconstruc-. . Extension to anatomical imaging: general online problem statement tion is not new and has raised a huge interest in the dynamic MRI community. Especially [Majumdar c] has proposed a real-time dynamic reconstruction based on the reconstruction of the difference between two consecutive frames. Yet the authors mentioned a reduction of the accuracy compared to the offline methods. Over the last years, dynamic MR image reconstruction has been solved by applying a sparse prior to the images and a low-rank one on the temporal variations from one scan to the next [Otazo ].

  and the k-space measurements associated with Ω k read y Ω k ∈ C kC . Using the aforementioned notation, we recall the definition of the acceleration and undersampling factors (AF and UF, respectively) with respect to the Cartesian reference as respectively the ratios n/S and N/(SC). Note that in the CS context, we get both S < n and SC < N . Although AF and UF usually evolve similarly (i.e. AF = UF), the Sparkling sampling scheme [Lazarus a] introduced in Section . . breaks down this relation as AF > UF.

  g. a dual variable delivered by the optimization algorithm. Note that when the last mini-batch is acquired (k = S), n b must be relaxed and set large enough to ensure convergence. The basic online formulation is recovered for a batch size b s = 1. Algorithm : Online mini-batch reconstruction algorithm for solving Problem . . initialize k = b s , z 0 ; while k ≤ S do z k = A k,n b (z k-bs ); k ← k + b s ; end It is worth noting that, since SubProblem ( . ) is convex, any good optimization algorithm is guaranteed to converge to a global solution whatever its initialization. However, the time to converge to the solution may be very sensitive to the initial guess. Therefore the warm restart procedure is beneficial to reduce the overall reconstruction time. To solve Subproblem ( . ), we propose to make use of the primal-dual Condat-Vù [Condat , Vũ , Combettes ] approach summarized in Algorithm . This choice was originally motivated by the usage in Eq. ( . ) of an analysis-based regularization term, which allows us to consider overcomplete dictionaries for Ψ such as undecimated wavelet transforms [Cherkaoui ]. If a synthesis-based prior had been chosen instead, one would have probably considered a fast proximal gradient method such as the POGM [Taylor a] or greedy FISTA algorithms [Liang ] as improved convergence rates can be reached by those algorithms. Note that β k denotes here the Lipschitz constant of ∇f Ω k . According to [Condat , Theorem . ], Algorithm converges to a solution of SubProblem ( . ). Algorithm : Condat-Vù algorithm A k,n b (x k,0 , v k,0 ) for solving SubProblem ( . ).

  Care and Use of Laboratory Animals. The acquisition parameters were set as follows: TR= ms (for collecting slices), TE= ms and flip angle FA= • with an in-plane resolution of µm for a field of view (FOV) of . cm which leads to a matrix size of N = 512 × 512, to maintain a high SNR, we considered a slice thickness of mm and excitation have been averaged. A fully sampled Cartesian reference scan composed of slices, each slice consisting of lines with samples each, was collected in an acquisition time of . s . In what follows, we will restrict our numerical experiments to the central slice, i.e. to D imaging. The image was reconstructed using an FFT and used as reference for computing image quality scores (e.g.. SSIM metric). Fig. .shows (a) the Cartesian reference, (b) the retrospective Cartesian variable density sampling mask and (c) the prospective -fold accelerated non-Cartesian Sparkling [Lazarus a]pattern. While the prospective reconstruction is more challenging due to potential B 0 inhomogeneities and real artifacts, the retrospective reconstruction is not affected, hence the retrospective validation is addressed first and the prospective one, which is more challenging, is presented afterwards.

  Figure . -: Single-channel coil acquisition: (a) ex vivo baboon brain Cartesian reference image, (b) retrospective Cartesian under-sampling scheme (AF = UF = 2.9), (c) prospective non-Cartesian Sparkling sampling scheme (implemented with a golden-angle ordering, two consecutive shots are highlighted in black) (AF = , UF = 2.5), the first collected being the horizontal one and the rotation being counterclockwise.

First

  Fig. . -(top-row) depicts partial reconstructions for the tested mini-batch sizes by the end of acquisition (i.e. before taking into account the last mini-batch). While large minibatches show aliasing artifacts (see Fig. . -for b s = 44, and b s = 88), the small batches deliver pretty accurate images, which tends to demonstrate the benefits of using small batch sizes for online reconstruction purposes. The reason for which we observed aliasing artifacts in large batch sizes is the varying amount of available k-space data in the top row of Fig. . -. Indeed, by the end of the acquisition (i.e. before the process of the last mini-batch i.e. at T = TR -ε) only a single shot is missing for b s = 1 while for larger b s (such as for b s = 88) a large number of shots are missing in the data-fidelity term (when b s = 88 half of the spoke are not processed).

  Figure . -: Single-channel retrospective Cartesian reconstruction: (a) evolution over time of cost function in Eq. ( . ) and (b) of the SSIM score for different batch sizes b s . The dark dashed line marks the end of acquisition.

  Figure . -: Online reconstruction of -fold (S = 34 shots) retrospectively accelerated Cartesian scan of ex vivo baboon brain in a single-channel coil acquisition setup. Images reconstructed for increasing values of mini-batch size b s at the end of acquisition (top) and at convergence (bottom).

  Figure . -: Online reconstruction of -fold (S = 34 shots) prospectively accelerated Sparkling scan of ex vivo baboon brain in a single-channel coil acquisition setup. Evolution over time of (a) the cost function in Eq. ( . ), and of the SSIM score in (b) for batch sizes b s corresponding to the primary factors of S. The dark dashed line marks the end of acquisition.

  . ) In Eq. ( . ), R refers to a structured sparsity promoting term. As shown in Section . , this formulation is able to recover images with a quality competitive with state-of-the art methods [She ], especially for non-Cartesian sampling schemes. Subsequently, the Octagonal Shrinkage and Clustering Algorithm for Regression (OSCAR) [El Gueddari b] based calibrationless reconstruction method is implemented, where the wavelet transform Ψ decomposes the stack of images into a stack of coefficients C = ΨX = c 1 , . . . , c L = Ψx 1 , . . . , Ψx L ∈ C N Ψ ×L . Each vector of wavelet coefficients c with ∈ {1, . . . , L} is made up of B bands. The b ∈ {1, . . . , B} band contains P b coefficients (c ,b,j ) 1≤j≤P b . The OSCAR penalty is expressed as follows: ,j | + γ < max{|c ,b,j |, |c ,b,j | , λ > 0, γ > 0. ( . )

  The donor gave his written consent before death to the donation program of University of Tours, France. The reference slice is displayed in Fig.. -(a). Prospectively accelerated Sparkling acquisitions were considered with AF = 8 and AF = 15 corresponding to S = 64 and S = 34 shots, respectively and to acquisition times of . s and . s for the slices as compared to . s for the Cartesian reference. Both are shown in Fig. . -(b)-(c). The -fold accelerated sampling pattern was the same as the one used in single-channel acquisition. We also investigated less accelerated acquisition (AF= ) using radially-initialized Sparkling trajectories to see to what extent this may impact online reconstruction. A larger number of shots actually offers more degrees of freedom to define the mini-batch size b s in online reconstruction.

  Figure . -: Multi-channel coil acquisition: (a) ex vivo human brain Cartesian reference image, (b) prospective Sparkling sampling pattern with S = spokes (i.e. AF = 8) and C = points (i.e. UF = 1.3), (c) prospective Sparkling sampling pattern with S = spokes (i.e. AF = 15), each one being composed of C = points (i.e. UF = 2.5).

  Fig. . -shows that there is only a little advantage for small batch sizes (i.e. b s = 8 and b s = 16) in the context of multi-channel receiver coils. This results from the fact that the overall computing time per iteration and constraint ( . ) allow to perform only n b = 1 (n b = 2, respectively) iteration(s) when s b = 8 (s b = 16, respectively). The zoom in the curve of SSIM score in Fig. . -(b) confirms this tiny advantage for small batch sizes. Stronger evidence in favor of small values of b s is shown in Fig. . -(top) where the ventricles (especially the right one) appear darker in b s = 8 as compared to b s = 32. Images obtained at convergence (Fig. . -(bottom)), i.e. once n b = 200 iterations have been run from k = S mini-batches, confirm the convergence to the Cartesian reference. Overall, in this multi-channel acquisition scenario the main bottleneck remains the computation cost per iteration. With a diminished cost, we could perform a larger number of iterations in particular for small batch sizes (for instance n b = 4 for b s = 8) and get better MR images by the end of acquisition.

  Figure . -: Online reconstruction of -fold (S = 34 shots) prospectively accelerated Sparkling scan of ex vivo human brain in a multi-channel coil acquisition setup using b s = 17. (a) Partial solution obtained by the end of acquisition and (b) at convergence.

  Figure . -: Split of the contributions to the global cost function (a) for offline reconstruction using the -fold accelerated Sparkling acquisition, and (b) evolution of the contributions ratio over iterations.

FiguresFigure

  Figures . -and . -, show the improvement of this new approach on the reconstruction of the -fold accelerated Sparkling trajectory (depicted in Figure . -(a)). Compared to Figure . -, the partial solutions show an improvement in terms of image quality especially for smaller batch sizes (b s = 2, 4). The same results hold of course for -fold accelerated scans, i.e. with Sparkling S = 34 (illustrated in Figure . -(b)). (a) (b)

Figure

  Figure . -: Fast-online reconstruction of -fold (S = 64 shots) prospectively accelerated Sparkling scan of ex vivo human brain in a multi-channel coil acquisition setup. Images reconstructed for increasing values of mini-batch size b s by the end of acquisition (top) and at convergence (bottom). There is a clear advantage for short batch sizes..

  Figure . -: Fast-online reconstruction of -fold (S = 34 shots) prospectively accelerated Sparkling scan of ex vivo human brain in a multi-channel coil acquisition setup. Images reconstructed for increasing values (only possible values here as is only divisible by and ) of mini-batch size b s by the end of acquisition (top) and at convergence (bottom).

Future

  work will be devoted to the adaptation of this framework to other imaging contrasts (e.g., T 2 -weighted imaging) and pulse sequences (e.g., turbo spin echo) for which larger times of repetition (e.g. s) give us the opportunity to successfully apply the proposed approach in the multi-channel acquisition setting. A pending challenge remains how we can extend this approach to D imaging in a time-efficient way. A nice feature we have recently implemented in PySAP to deal with D imaging applies to the situation when the k-space data are collected along parallel and uniformly-spaced planes, as done in z-variable density SPARKLING [Lazarus b]. In this case, we can restrict the usage of the NFFT to D only and make use of the FFT along the third dimension. This typically allow us to decrease the memory usage by a factor of and the computation time by a factor of .

  while the reconstruction relying on the spatial discrimination of the coil sensitivity profiles. However the maximum accelerations factor was fairly low. Over the last years, the development of compressed sensing theory and more specifically, the gain in acquisition times and the recovery guarantees on the reconstruction made a breakthrough in the MRI field. Although the practical implementation slightly departs from theory, the resulting acceleration factors exceeded the ambition, up to in D imaging and in D imaging for the SPARKLING trajectories proposed by [Lazarus a]. Despite this acceleration, the high-resolution imaging context remains unsuited for clinical purposes, the reason for that is rather simple: the long acquisition time was transferred to an expensive reconstruction and thus the gain in acquisition provided by CS theory was not converted into a global acceleration in the whole acquisition/reconstruction pipeline. The purpose of this PhD thesis was thus to reduce this reconstruction time and therefore to bring the high-resolution MRI to the clinical practice. For this purpose, we have developed a new approach that interleaves acquisition and reconstruction called online reconstruction.

  were hampered by the lack of SNR. Therefore the first step was to implement a simple yet efficient method to estimate the coil sensitivity profiles involved in the SENSE-like multichannel reconstruction problem. The proposed method has been used for many contrasts and multiple sampling strategies, additionally it has been extended to the d imaging case and has been used as the reference reconstruction technique involved in the d and d SPARKLING papers [Lazarus a, Lazarus b]. However the sensitivity estimation constitutes the major limitation for the application to online reconstruction. Hence, after a deep review of the literature on image reconstruction for non-Cartesian sampling and multiple trials, calibrationless reconstruction appeared as the most appealing technique.

  their iterative nature, therefore to speed up reconstruction, a promising avenue lies in the usage of Deep Learning (DL) approaches instead. In particular unrolled modelreason for this is rather simple, they are less mysterious than datadriven counterparts.However DL methods may require a large database of raw data (i.e. high resolution multi-channel complex-valued k-space measurements unavailable until recently[Zbontar ]). Although the application of these methods are promising, it is worth mentioning that they are sensitive to small perturbations or adversarial attacks as it has been proved in[Antun ].All in all, we have proposed a new framework to adapt MR image reconstruction algorithms to the context of high-resolution highly accelerated non-Cartesian and parallel imaging, and deploy them in a clinical environment.stream is configured dynamically at runtime using a plugin architecture. More specifically, the Gadgets of a particular reconstruction stream can reside in multiple independently compiled binaries (shared libraries) and are loaded at runtime based on a Gadget stream configuration in an extensible markup language (XML) file. Data are passed from a given client into the Gadgetron streaming architecture through a Transmission ControlProtocol/Internet Protocol (TCP/IP) connection through which the resulting images are eventually returned to the client application.

Figure

  Figure C. -: Gadgetron communicates with the client through a TCP/IP protocol. The client application sends data to the Gadgetron and associated with each data package is a Message ID. Based on the Message ID, control of the socket is handed over to a specific Reader, which is capable of deserializing the incoming data package. The data package is converted to message blocks that are added to the first Gadget's queue. Data are then passed down the Gadget stream where each Gadget can modify and transform the data. Any Gadget can return images (or partially processed data) to the Gadgetron framework. Based on the MessageID of this return data package, the control of the socket and the data are handed to a particular Writer, which is responsible for writing the return message to the client. From [Hansen ].

  Figure C. -: Schematic view of the Gadgetron installation at the T scanner.

  Figure C. -: Schematic view of the PySAP architecture.

  ∀x, y ∈ H, Ax , y = x , A H y (C. )Concerning the optimization part, we have tested the gradient descent of all optimizers for the smooth part in the cost functions. We have also checked that the value of the results should be equal to the inverse Fourier transform for Cartesian sampling regardless of the Fourier operator used. Last, The presence of unit-test permits to increase the code reliability and provide use-case to end-users.Introduction à l'IRMUn scanner IRM est organisé autour d'un puissant aimant caractérisé par la puissance de son champ magnétique. Lorsqu'un volume est placé dans un champ magnétique statique noté B 0 , les spins de protons d'hydrogène composant le volume s'alignent ou s'anti-alignent le long du champ, laissant néanmoins une légère propension pour l'état de plus faible énergie, créant ainsi à l'échelle macroscopique un moment magnétique, appelé M 0 . Afin d'obtenir un signal variant dans le temps, l'aimantation M 0 est basculée grâce à un champ radio-fréquence (RF), noté B + 1 transverse au champ statiques B 0 , c'est le phénomène d'excitation. Une fois l'aimantation basculée, l'impulsion RF le champs B + 1 est arrêtée et l'aimantation revient à son état d'équilibre en tournant autour de B 0 (phénomènes de relaxation et de précession). La fréquence de rotation -aussi appelée fréquence de Larmor-dépend de la composition du tissu mais également de l'intensité du champ magnétique généré. Durant la relaxation un signal appelé FID est mesuré. Celui-ci peut-être caractérisé par deux constantes de temps : -T 1 , correspondant à la constante de temps décrivant le retour à l'état d'équilibre de la composante longitudinale ; -T 2 , correspondant à la constante de temps décrivant la disparition de la composante transverse. Ces deux constantes sont caractéristiques du tissu sondé. Afin de discriminer spatialement le signal acquis, il est important d'encoder spatialement celui-ci. Pour cela, des gradients d'encodage, permettent de faire changer localement la fréquence et phase encodant ainsi le domaine conjugué de l'espace de Fourier aussi appelé espace k. Afin de reconstruire l'image à partir du signal magnétique électrique, la transformée de Fourier y est appliquée. Ainsi pour acquérir une image, il est d'usage de répéter l'excitation pour différentes valeurs de fréquence et de phase. De ce fait l'IRM est une modalité lente car les constantes de temps impliquées sont à minima de l'ordre de la dizaine de millisecondes. Pour accélérer les acquisitions, l'imagerie parallèle a été créée. En effet en sous-échantillonnant l'espace k de manière à acquérir une ligne sur deux, le temps d'acquisition peut se voir en effet réduit. Cependant, lors de la reconstruction le théorème de Shannon-Nyquist n'étant plus vérifié, du repliement est observé sur l'image. Pour s'affranchir de ce repliement, l'idée est de se baser sur l'information spatiale apportée par les canaux de réception de l'antenne en résolvant un système linéaire. Bien que les premières méthodes d'accélération développés fassent maintenant parti des séquences utilisées en clinique, les facteurs d'accélération correspondant restent faibles. Pour réduire d'avantage ces temps, l'échantillonnage compressif permet d'accélérer l'acquisition tout en garantissant des preuves de reconstruction

  Figure D-: Comparaison des différentes méthodes de reconstructions d'images acquises à l'aide d'une trajectoire SPARKLING -dont le facteur d'accélération est de en temps -et d'antenne multicanaux. (a) La référence Cartésienne, (b) reconstruction Sparse-SENSE avec notre méthode d'extraction des matrices de sensibilité, (c) la reconstruction 1 -ESPIRiT, (d) La solution p-LORAKS après la projection sur la grille cartésienne. (e-h) leurs zooms respectifs.

  Figure D-: Première ligne : Boule unitaire de la norme à k-support en D pour les différentes valeurs de k. (a) présente la boule unitaire pour k= -correspondant à la régularization 1 , (b) présente la boule unitaire pour k= , (c) présente la boule unitaire pour k= -correspondant à la régularisation 2 . Seconde ligne : Régularisation Elastic-net -combinaison linéaire entre la norm 1 et 2 -et la boule unitaire d'OSCAR en D.

  b-OSCAR (d) Norme k-support (e) 1 -ESPIRiT (f) AC-LORAKS Figure D-: Première ligne : Magnitude des images acquises à l'aide des trajectoires SPARKLING avec un facteur d'accélération égal à . (a) Référence Cartésienne. (b) La transformée de Fourier inverse obtenue grâce à la méthode des moindres carrés. (c) Reconstruction obtenue par la méthode sans calibration basée sur la norme OSCAR. (d) Reconstruction obtenue par la méthode sans calibration basée sur la norme à k-support. (e) Reconstruction obtenue par la méthode auto-calibrée 1 -ESPIRiT. (f) Reconstruction obtenue par la méthode auto-calibrée (AC) LORAKS basée sur la régularisation de l'espace k. Deuxième ligne : Zooms respectifs de la partie encadré en rouge. Dernière ligne : Images de phase extraites pour chacune des méthodes grâce à l'algorithme du canal virtuelle [Parker ].

  Figure D-: Reconstruction rapide en ligne de données SPARKLING accélérées fois (S= segments) d'un cerveau humain ex vivo acquise par un réseau d'antennes. Première ligne : Les images reconstruites à la fin de l'acquisition, pour de tailles de paquets différentes ( et segments par paquets). Seconde ligne : Les images obtenues à convergence.

  represented, one needs to introduce a density compensation step to correct for potential bias arising from high concentration of samples in those regions. Akin to grid-data methods, the reconstruction accuracy is a trade-off between the interpolator complexity and k-space oversampling. A fair compromise consists in choosing a Kaiser-Bessel kernel with

	a limited (e.g. from . to . -times) oversampled grid, see [Beatty	] for details.
	The non-uniform Fast Fourier Transform.	Others existing data-driven grid-
	ding method have been proposed in the literature [Dutt	, Fessler	, Greengard	,
	Keiner	], hereafter we will only remind the nonequispaced Fast Fourier Transform
	(NFFT) introduced by [Keiner	].
	Let us consider a d-dimensional nonequispaced Discrete Fourier Transform defined by
	a set of arbitrary spatial node χ and a frequency bandwidth vector N . Each nodes x j
	in the nonequispaced sampling set χ

  Cartesian acquisition with an image field of view of FOV x × FOV y and a matrix size of N x × N y the image resolution is defined by ∆ x = FOVx Nx and ∆ y =

				The Nyquist-Shannon the-
	orem states that the spectrum of a discretely sampled signal is replicated in the Fourier
	conjugate domain. It provides a prescription for the nominal sampling rate required to
	avoid aliasing. It may be stated as follows:
	Theorem . . . A band-limited continuous-time signal can be sampled and perfectly re-
	constructed from its discrete samples if the waveform is sampled at least twice as fast as
	its highest frequency component.
	Spatial encoding in MRI might not deviate from this rule, hence to prevent wrap-
	around artifact or aliasing the theorem should be applied to the k-space domain. Let us
	consider a D FOVy Ny .
	To comply with the Nyquist-Shannon criterion and get the desired resolution, the k-space
	sampling should respect the following relations:
	∆k x ≤	1 FOV x	and ∆k y ≤	1 FOV y	( . a)

  = |p | e -iϕ ref / k |p k |, I is the reconstructed image of the th coil, ϕ ref is the phase of the weights used to avoid SNR loss, ϑ is the coil phase, θ is the true tissue phase and η is a random component due to the thermal noise in the individual RF coil measurement. While θ and ϑ are supposed to be slowly varying, i.e. smooth in space, the phase noise η is expected to be uncorrelated from one voxel to another. Because the coil phase ϕ ref varies slowly in space and in order to reduce measurement errors the

	|e ı(ϑ +θ+η )	( . )
	where the weights w

Table . :

 . Correspondence between SSIM score ans Mean Opinion Scores (MOS) between an image and its reference.

	SSIM	MOS	Quality	Impairments
	≥ 0.99		Excellent	Imperceptible
	[0.95, 0.99(		Good	Perceptible but not annoying
	[0.88, 0.95(		Fair	Slightly annoying
	[0.5, 0.88(		Poor	Annoying
	≤ 0.5		Bad	Very annoying
	with either Gaussian blurring, wavelet compression, Rician noise [Gudbjartsson	] or
	motion artifacts were presented to two neuroradiologists who carefully read and rank the
	images. The physicians were asked to score overall image quality on a -Likert scale

that was determined by consensus on a training set before scoring the images included in the results. Correlations between each image quality metrics and the physician' scores were measured using the Pearson linear correlation coefficients(PLCC) and the Spearman rank order correlation coefficients (SROCC). According to

[Mason ]

the best quality metric should be linear with respect to the Likert index of the physician. Fortunately, for the SSIM score above . this was the case. Hence in the rest of the thesis, we will stick

  Parts of this chapter (Section . . ) were presented at International conferences as: L. El Gueddari, C. Lazarus, H. Carrié, A. Vignaud and P. Ciuciu. Self-Calibrating Nonlinear Reconstruction Algorithms for Variable Density Sampling and Parallel Reception MRI. In IEEE th Sensor Array and Multichannel Signal Processing Workshop (SAM), pages -,

	July	.
	Chapter L. El Gueddari, C. Lazarus, H. Carrié, A. Vignaud and P. Ciuciu. Self-calibrating nonlinear
	MR image reconstruction algorithms for variable density sampling and parallel imaging. In
	Proceedings of the th Annual Meeting of ISMRM, page	, Paris, France,	.
		Compressed Sensing reconstruction
		from multi-channel acquisition
		8 8 8 Methodological Developments Chapter Outline The method developed in this chapter (Section . . ) was used as reconstruction tech-nique in the following works : C. Lazarus, P. Weiss, N. Chauffert, F. Mauconduit, L. El Gueddari, C. Destrieux, I. Zemmoura, A. Vignaud and P. Ciuciu. SPARKLING: variable-density k-space filling curves for accelerated T * 2 -weighted MRI. Magnetic resonance in medicine, vol. , pages -, . 8 8 Part II C. Lazarus, P. Weiss, L. El Gueddari, F. Mauconduit, C. Destrieux, I. Zemmoura, A. Vig-

. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . Parallel MR image reconstruction based on sensitivity maps . . . . . . . . Coil sensitivity extraction . . . . . . . . . . . . . . . . . . . Blind bi-linear reconstruction . . . . . . . . . . . . . . . . Coil-by-coil MR image reconstruction . . . . . . . . . . . . . . . . . K-space driven methods promoting data low-rankness . . . . . . . . Domain based methods . . . . . . . . . . . . . . . . . . . Experiments & Results . . . . . . . . . . . . . . . . . . . . . . . Experimental setup . . . . . . . . . . . . . . . . . . . . . Results . . . . . . . . . . . . . . . . . . . . . . . . . Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . naud and P. Ciuciu. D SPARKLING trajectories for high-resolution T * 2 -weighted Magnetic Resonance imaging. NMR biomedicine, . Under review.

Table . :

 . Reported computation time for D-Cartesian reconstruction (acquired with anchannel coil) with a target image resolution equal to 256 × 256 and an under-sampling factor

	of 2.78, for different methods including: Sparse-BLIP [She	], J-SENSE [Ying	],
	Sparse-SENSE [Liang	] and IRGN-TV [Knoll	]. Originally presented in [Shin	].
	Sparse-BLIP	J-SENSE	Sparse-SENSE	IRGN-TV
	s	s	s		s

Table . :

 . Ex-vivo baboon and in-vivo human brain data collected using various slice thickness and acceleration factors A.

		Object	Slice thickness	AF
	.		mm	( shots)
	.	Baboon brain (ex-vivo)	mm	( shots)
	.		mm	( shots)
	. .	Human brain (in-vivo)	mm	( shots)

  the context of this thesis, i.e. non-Cartesian high-resolution imaging, SENSE-like me-

				Calibrationless multi-channel reconstruction
	reconstruction for compressed sensing and parallel MR imaging. In Proceedings of the th
	Annual Meeting of ISMRM, page	, Montreal, QC, Canada,	.
	L. El Gueddari, P. Ciuciu, E. Chouzenoux, A. Vignaud and J.-C. Pesquet. Calibrationless
	thods and coil-by-coil image reconstructions are the only potential competitors. However, owing to the presence of coil sensitivity profiles in SENSE-like formulation, this kind of reconstruction methods is hardly online compliant. A potential alternative could be the coil-by-coil image reconstruction promoting joint sparsity across the different channels, such as CLEAR. In this thesis we have not addressed the choice of regularization para-meters which may be a tricky problem in various scientific fields (such as data-science, statistical inference, and image restoration or reconstruction) depending on the overall 8 8 8 Calibrationless multi-channel reconstruction: from group-sparsity to structured sparsity Chapter Outline OSCAR-based image reconstruction in compressed sensing parallel MRI. In IEEE th 8 8 Chapter International Symposium on Biomedical Imaging (ISBI ), Venice, Italy, . IEEE ISBI.
		. .	Octagonal Shrinkage and Clustering Algorithm for Regression . .
		. .
	.	Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . .
	.	Discussion . . . . . . . . . . . . . . . . . . . . . . . . . .

similarity/global cost addressed. Instead, we chose to maximize the SSIM score, however sometimes we found some inconsistencies between visual perception of reconstruction images and the corresponding SSIM scores. Since no artifact related to the sampling scheme was present on the coil-specific images, one can wonder whether sensitivity information really matters for non-Cartesian reconstruction. In the rest of this thesis we review, propose and evaluate new regularization schemes based on structured sparsity for coil-by-coil high-resolution non-Cartesian MR image reconstruction. . Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . Calibration-less reconstruction: general problem statement . . . . . . . . . General problem statement . . . . . . . . . . . . . . . . . . Optimization algorithm . . . . . . . . . . . . . . . . . . . From group-sparsity . . . . . . . . . . . . . . . . . . . . . . . . Least Absolute Shrinkage and Selection Operator . . . . . . . . . Elastic-Net . . . . . . . . . . . . . . . . . . . . . . . . . Group-LASSO . . . . . . . . . . . . . . . . . . . . . . To structured sparsity . . . . . . . . . . . . . . . . . . . . . k-support norm . . . . . . . . . . . . . . . . . . . . . . Experimental setting & Results. . . . . . . . . . . . . . . . . . . . Experimental setup . . . . . . . . . . . . . . . . . . . . . Results . . . . . . . . . . . . . . . . . . . . . . . . The methods and principles contained in this chapter were also presented at international conferences as: L. El Gueddari, P. Ciuciu, E. Chouzenoux, A. Vignaud and J.-C. Pesquet. OSCAR-based

Table . :

 . Performances of the least squares (LS) solution (zero-order reconstruction) on the dataset Table.. The results in red box are displayed in Figure.-.

		Idx	.	.	.	.	.	.	.	.	.	.
		SSIM	.	.	.	.	.	.	.	.	.	.
	iFT	p-SNR	.	.	.	.	.	.	.	.	.	.
		NRMSE	.	.	.	.	.	.	.	.	.	.

Table . :

 . Reconstruction of -fold accelerated scans acquired with sparkling and spiral trajectories. Comparison of global regularization functions (LASSO, Elastic-net, g-OSCAR and k-support norm) with hyper-parameters chosen to as to maximize the SSIM score. Best values appear in bold font. The images associated with the scores highlighted in red boxes are displayed in Figure.-.

	Idx		LASSO			Elastic-Net			g-OSCAR			k-support norm
		SSIM p-SNR NRMSE SSIM p-SNR NRMSE SSIM p-SNR NRMSE SSIM p-SNR NRMSE
	.	.	.	.	.	.	.	.	.	.	.	.	.
	.	.	.	.	.	.	.	.	.	.	.	.	.
	.	.	.	.	.	.	.	.	.	.	.	.	.
	.	.	.	.	.	.	.	.	.	.	.	.	.
	.	.	.	.	.	.	.	.	.	.	.	.	.
	.	.	.	.	.	.	.	.	.	.	.	.	.
	.	.	.	.	.	.	.	.	.	.	.	.	.
	.	.	.	.	.	.	.	.	.	.	.	.	.
	.	.	.	.	.	.	.	.	.	.	.	.	.
	.	.	.	.	.	.	.	.	.	.	.	.	.

Table . :

 . Quantitative comparison of coefficient based regularization schemes (c-OSCAR and group-LASSO) with hyper-parameters set to maximize the SSIM score. Best values appear in bold font.

	Idx		group-LASSO			c-OSCAR
		SSIM p-SNR NRMSE SSIM p-SNR NRMSE
	.	.	.	.	.	.	.
	.	.	.	.	.	.	.
	.	.	.	.	.	.	.
	.	.	.	.	.	.	.
	.	.	.	.	.	.	.
	.	.	.	.	.	.	.
	.	.	.	.	.	.	.
	.	.	.	.	.	.	.
	.	.	.	.	.	.	.
	.	.	.	.	.	.	.

Table . :

 . Quantitative comparison of OSCAR-norm regularization schemes with hyperparameters set to maximize the SSIM score. Best values appear in bold font.

	Idx		g-OSCAR			s-OSCAR			b-OSCAR			c-OSCAR
		SSIM p-SNR NRMSE SSIM p-SNR NRMSE SSIM p-SNR NRMSE SSIM p-SNR NRMSE
	.	.	.	.	.	.	.	.	.	.	.	.	.
	.	.	.	.	.	.	.	.	.	.	.	.	.
	.	.	.	.	.	.	.	.	.	.	.	.	.
	.	.	.	.	.	.	.	.	.	.	.	.	.
	.	.	.	.	.	.	.	.	.	.	.	.	.
	.	.	.	.	.	.	.	.	.	.	.	.	.
	.	.	.	.	.	.	.	.	.	.	.	.	.
	.	.	.	.	.	.	.	.	.	.	.	.	.
	.	.	.	.	.	.	.	.	.	.	.	.	.
	.	.	.	.	.	.	.	.	.	.	.	.	.

  Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . El Gueddari, E. Chouzenoux, A. Vignaud, J.-C. Pesquet and P. Ciuciu. Online MR image reconstruction for compressed sensing acquisition in T * imaging. In Wavelets: Applications in Signal and Image Processing XVIII. International Society for Optics and Photonics, .T his chapter propose a new way for accelerating MR image reconstruction in the context of CS-accelerated acquisitions.

	8 8 8 Online MR image reconstruction Chapter Outline . . Parameters setting . . . . . . . . . . . . . . . . . . . . . . Results . . . . . . . . . . . . . . . . . . . . . . . . Multi-channel receiver coil . . . . . . . . . . . . . . . . . . . . . . Problem statement . . . . . . . . . . . . . . . . . . . . . . Parameters setting . . . . . . . . . . . . . . . . . . . . . . Results . . . . . . . . . . . . . . . . . . . . . . . . . . Online multi-channel reconstruction: improvements . . . . . . . Parts of this chapter were presented at international conference in : L. El Gueddari, P. Ciuciu, E. Chouzenoux, A. Vignaud and J.-C. Pesquet. Online compressed sensing MR image reconstruction for high resolution T * 2 imaging. In Proceedings of the th 8 8 Chapter Annual Meeting of ISMRM, page , Montreal, QC, Canada, .

. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . Extension to anatomical imaging: general online problem statement . . . . . Acquisition considerations . . . . . . . . . . . . . . . . . . . Variational formulation of image reconstruction . . . . . . . . . Single-channel receiver coil . . . . . . . . . . . . . . . . . . . . . Problem statement . . . . . . . . . . . . . . . . . . . . . L

This relation holds only for small objects compared to the wavelength, in the other cases the expression of the signal is still debated.

8 8

http://cosmic.cosmostat.org/ http://www.cosmostat.org/ http://github.com/cea-cosmic/pysap/ https://github.com/CEA-COSMIC/pysap-mri

Acknowledgements

The work presented in this thesis was only possible with the support and the help of many people. I would like to take the opportunity to express my gratitude.

First I would like to express the deepest appreciation to my supervisor Dr. Philippe

Chapter : Online MR image reconstruction

Finally the online reconstruction framework is introduced in this final contributed chapter. A general problem proposes a data splitting into mini-batches scheme, and develops the optimization. Then, we adopt a pedagogical approach to first address the simpler case of single-channel data acquisition. In this case, a single image is reconstructed, either from Cartesian and non-Cartesian data. In particular, we explore the optimal way for tuning the batch size in order to deliver a reliable solution by the end of acquisition. Last, we extend the online approach to multi-channel acquisition by considering the calibrationless method we described in Chapter . For given computational power, we show how the cost of the proximity operator associated with structured sparsity may prevent us from fulfilling the online reconstruction constraints. Therefore we introduced some online trick in the reconstruction process to remain compatible with these online constraints. Finally, we validate the proposed method on the same data set as in Chapter .

Sparkling sampling pattern (Fig. . -(c)) was generated with S = 34 shots from a radial initialization, each shot being composed of C = 3072 samples leading to UF = 2.5 and AF = 15. The total scan time for the slices using Sparkling was reduced to . s. The acquisition was segmented in time according to the golden-angle ordering scheme, as it was proven [Winkelmann ] to increase the maximize the k-space coverage during the scan. Hence, between two consecutive shots a rotation of • was applied as highlighted in Figure . -(c).

Reconstruction parameters

The decimated Symmlet Wavelet Transform with decomposition scales was used as sparsifying transform. Regarding hyper-parameter λ in Eq. ( .) it has been set retrospectively so as to maximize the structural similarity score (SSIM) [Wang ] to the reference image once all data are available (i.e. offline scenario). Note that in a more clinically plausible setting, the value of λ needs either to be estimated using the noise statistical properties [Stein ] based on a pre-calibration scan or set using cross-validation across slices in D imaging [Hastie ]. For non-Cartesian acquisition, the GPU [Lin ] implementation of the NUFFT [Fessler ] was used. In Section . . , we compare the performances of several batch sizes with a given number of iterations defined by Eq. ( . ).

Once all shots are considered in Algorithm , i.e. k = S, then the number of iterations is set to n b = 200 to guarantee that convergence was reached.

Computing parameters setting

All experiments were run on a machine with GB of RAM and an -core ( . GHz)

Intel Xeon E -v Processor. All the codes have been developed in Python using the PySAP package . The values of computing times have been obtained using epochs of Condat-Vú iterations each, the mean and standard deviation being summarized in Table . . Using multiple epochs permits to account for potential variability in computing times due to concomitant processes running on the machine. Although the NUFFT is usually slower than the FFT, here we observed the converse as we used the GPU implementation of the NUFFT and the CPU implementation of the FFT. A more fair comparison might be achieved using the GPU version of FFT but this was not mandatory for the present online study as the FFT was not the computational bottleneck.

Table . : Computing time for one iteration estimated on epochs of iterations each using timeit.. Let z ∈ C p with p ≥ 1. We introduce the magnitude sorting operator S p : C p → C p such that vector S p (z) = (S p (z) j ) 1≤j≤p contains the p entries of z sorted in decreasing order in magnitude, i.e. such that

For the sake of compactness, we renamed the sorting operator S p as the permutation matrix P ↓ , by definition, matrix P ↓ is orthogonal. We recall the definition of the OSCAR norm: with w ∈ R p + a vector of hyper-parameters such that w 1 ≥ • • • ≥ w p ≥ 0. More precisely, OWL and OSCAR become equivalent if one sets the OWL weights as follows: w j = λ + γ(p -j) for j = 1, . . . , p. Proposition B. . . The OSCAR norm defined in Eq. (B. ) is equivalent to the OWL norm in Eq. (B. ) with weights w j = λ + γ(p -j) for j = 1, . . . , p.

Proof. Set z ∈ C p . First, let us remark that, for any permutation matrix P ∈ R p×p , Ω λ,γ (z) = Ω λ,γ (P z). Let us consider now the particular permutation matrix P ↓ such that Proximity Operator details of the Ordered Weighted 1 norm P ↓ z = z ↓ . Then, Ω λ,γ (P ↓ z) = Ω λ,γ (z ↓ ) and we get:

with the above mentioned definition of (w j ) 1≤j≤p .

Proximity operator

By definition, the proximity operator of the OSCAR norm at z reads:

Since the P ↓ operator is orthogonal and using the property of composition with semiorthogonal linear transform as defined in [Combettes ,

Proof. By definition, matrix P is orthogonal, i.e. P P = P P = I p , the p-dimensional identity matrix. Moreover Ω λ,γ (z) = Ω λ,γ (P z). Using the property of composition with semi-orthogonal linear transform as defined in [Combettes , Table . ] we get prox Ω λ,γ (z) = z + P (prox Ω λ,γ (P z) -P z) and finally the result as P is orthogonal.

prox

with w defined as Proposition B. . . The proximity operator of OWL norm can be solved efficiently using the following proposition.

Proposition B. . . Let z in C p and the OWL norm defined in Eq. (B. ). Then, the following equality holds:

with being the Hadamard product ( i.e. element-wise), sign(v) = v |v| , ∀v ∈ C * and sign(0) = 0.

Proximity Operator details of the Ordered Weighted 1 norm

Proof. Let w ∈ R p + being a vector of coefficients sorted in non-increasing order, prox Θw (z ↓ ) = arg min • Facilitate cross platform Development: The platform is operating-system independent such that the users can then deploy it on their platform of choice, which may be dictated by the medical imaging device or their local computing environment.

• Flexibility: The project has been designed to be compatible over several vendors, and is not only restricted to one modality. It can be used either for CT, PET or MRI modality.

• Modular: Different modules are common to all the modalities and can then be shared across vendors.

• Facilitate deployment: Once the communication between the Gadgetron host and the scanner is set the deployment of the new algorithm can be easily performed and since the framework is used directly in the scanner loop, it allows online reconstruction but also online post-processing to be completely transparent to the user view point.

• Various programming languages: One of the key strengths of such method is its independence with respect to the development language hence C++, Matlab or even Python code can be run during the same pipeline.

• Free: The platform is freely available to all researchers.

Software architecture

In order to optimize the development of new reconstruction algorithms the idea of the Gadgetron project is to pool common functions. Each of these functions or modules (called Gadget) has an input fed with data and an output. Multiple Gadgets are wired together (connected) to assemble a reconstruction process (a Gadget stream). The Gadget

PySAP-gadgetron

The PySAP-gadgetron plugin proposes an interface between the MR scanner and the The generated trajectories led to acceleration factor up to 20 in 2D and 60 for 3D-acquisitions on highly resolved T * 2 weighted images acquired at 7 Tesla. Those accelerations were only accessible thanks to the high input Signal-to-Noise Ratio delivered by the usage of multi-channel reception coils. However, such high under-sampling factors were coming at a price of long and complex reconstruction. In this thesis the objective is to propose an online approach for non-Cartesian multi-channel MR image reconstruction. To achieve this goal we rely on a progressive online approach where the reconstruction starts from incomplete data. Hence acquisition and reconstruction are interleaved, and a partial feedback is given during the scan. After exposing the CS theory, we present state-of the art methods dedicated to multi-channel coil reconstruction. In particular we first focus on self-calibrating methods that present the advantage to be adapted to non-Cartesian sampling and we propose a simple yet efficient method to estimate the coil sensitivity profiles based on the data itself, i.e. without any precalibration scan. However, owing to its dependence to userdefined parameters, this two-step approach (extraction of sensitivity maps and then image reconstruction) is not compatible with the timing constraints associated with online reconstruction. Then we study the case of calibration-free reconstruction methods and splits them into two categories, the k-space based and the domain-based. While the k-space calibration-free methods are sub-optimal for non-Cartesian reconstruction, due to the gridding procedure, we retain the domain-based reconstruction and prove their interests and performances for online purposes. Hence in the second part, we first prove the advantage of mixed norm to improve the recovery guarantee in the parallel MRI setting. Then we investigate the impact of structured sparse induced norm on the reconstruction multichannel purposes, where we adapt different penalty to handle those highly correlated images. Finally, the retained method is applied to online purposes. The entire pipeline, is compatible with an implementation through the Gadgetron pipeline to deliver the reconstruction at the scanner console.

Motivations et contextes