INTRODUCTION

The introduction of magnetic resonance imaging (MRI) in the clinical and research community started a new era of medical imaging. Compared with X-ray, nuclear medicine, ultrasound-based imaging technologies, MRI technology, which is based on nuclear magnetic resonance (NMR), can provide images with a high soft-tissue contrast. Another technology based on NMR, which is called magnetic resonance spectroscopy (MRS), can be used to provide metabolic information. Both in MRI and MRS, there is an increasing interest in quantitative approaches.

An important feature of MRI is that the image contrast is fundamentally multiparametric, primarily based on T1, T2 relaxation time and proton density. Furthermore, it was soon recognized that the diffusion and flow of water molecules in tissues can play an important role in the contrast of MRI images. One of the MRI methods exploiting the water diffusion in tissues is the intravoxel incoherent motion (IVIM) imaging, that can be used to assess simultaneously diffusion and perfusion. IVIM imaging has been proved of great value in diagnosis, staging and prognosis, and its clinical application in the human body such as head and neck, prostate, liver, kidney and other parts is increasing rapidly. However, the quantification of IVIM perfusionrelated parameters still suffers from large variability and low repeatability. This is particularly true in organs such as liver, which contains a rich vascular network. The numerous blood vessels of different sizes in liver represent a confounding factor in the IVIM quantification of tissue perfusion.

Between MRI and MRS, the MR spectroscopic imaging (MRSI) is an imaging technology that combines to some extent the metabolic information of MRS with the spatial coverage of MRI. The MRSI is also known as the chemical shift imaging (CSI). Compared with IVIM, which is used to quantify the parameters describing physical phenomena such as diffusion and perfusion in the tissue, CSI is mainly used to quantify the metabolite content in different organs such as brain, prostate, muscle, etc.. In the CSI technology, there are two mainstream solutions for providing reference signals for quantitative quantification: in vitro physical phantom and in vivo reference signal. However, in vivo reference signals can be easily affected by pathology, and the use of physical phantom outside the body has a number of drawbacks, such as for instance the loss of CSI spatial resolution.

I

This dissertation focused on these quantitative techniques in IVIM and CSI, specifically:

1. Sparsity-based All-Voxel Tri-Exponential IVIM (SAVTE-IVIM) algorithm

In view of the current clinical interest in the quantification of the parameters of the IVIM diffusion-weighted images of liver, and the vascular structure with its confounding effect in the liver, we proposed a novel method, a sparsity-constrained all-voxel tri-exponential IVIM (SAVTE-IVIM) algorithm, that can automatically identify the existence of potential blood vessels in the target region-of-interest (ROI). In addition to the main goal of identifying potential blood vessels in a given ROI, the algorithm can simultaneously quantify the IVIM parameters of all voxels in the ROI to evaluate the diffusion, perfusion, and blood vessel confounding effect in each voxel. Specifically, we proposed a tri-exponential model based on sparse constraint to describe all voxels at the same time. Typically, the IVIM parameters are evaluated voxel-by-voxel (voxel-wise), while SAVTE-IVIM can quantify all the voxels in the ROI simultaneously. In addition, in order to solve the proposed new model, an optimization algorithm, based on the idea of Alternating Direction Multiplier (ADMM) together with the use of Levenberg Marquardt algorithm to deal with nonlinear problems, was proposed. Two strategies for the inherent non-negative constraints were also introduced.

Implementation and Comparison of Five Fitting Algorithms for IVIM Quantification on Vertebral Bone Marrow

Since Marchand et al. successfully applied the IVIM method to the quantification of bone marrow in 2014, there has been an increasing interest in bone marrow IVIM in recent years. However, there are still issues with the current image quality of bone marrow IVIM; furthermore, little attention has been paid to the investigation of optimal algorithms for IVIM quantification of bone marrow. It is worth noting that an optimal algorithm can to some extent compensate for the low image quality. In view of the above considerations, i) we applied a recently proposed protocol that improves the image quality in bone marrow IVIM, using the RESOLVE (readout segmentation of long variable echo train) sequence and ii) we implemented five algorithms for the parameter quantification of vertebral bone marrow IVIM. Four algorithms, One-Step, Two-Step, Three-Step and Fixed-D * algorithms are based on the idea of least squares (LSQ), and the fifth one is a Bayesian-based algorithm. A comparison among these algorithms was conducted. Furthermore, maps of the IVIM parameters were generated and compared.

Virtual phantom chemical shift imaging (ViP CSI)

Based on some shortcomings of the mainstream solutions for providing reference signals in CSI technology (including physical phantom in vitro and internal reference): the demand for an additional MR scan, decrease of the CSI resolution, etc., we proposed to use the virtual phantom technique, which was extended from ERETIC (Electronic REference To access In vivo Concentrations) technology to design a virtual phantom, to provide in CSI the reference signal for quantification. The amplitude and frequency can be custom designed, to simulate the FID signal which is going to be acquired with CSI to provide a reference signal. In the classic CSI technology, with different phase encoding gradients, a repeated NMR signal acquisition process is required. For example, to generate an 8 × 8 matrix, data acquisition needs to be repeated 64 times. However, the proposed scheme can generate the same reference signals in all 64 voxels II by transmitting the ViP signal only once in one of the 64 data acquisitions. In this way, an additional MR scan is no more needed, simplifying the CSI imaging protocol. It can also be further applied to different data acquisition schemes, such as elliptical CSI or weighted CSI. In addition, the ViP reference signal can be customized in its amplitude and frequency and displays an excellent spatial uniformity and time stability for all voxels in the sample.

III IV

Contents

INTRODUCTION I

Table of Contents V

Notations XVII

Notations

x scalar

x vector

x (n) n-th component of the N -dimensional vector x X matrix Magnetic Resonance Imaging and Diffusion MRI

X (n 1 ,n 2 ) (n 1 , n 2 )-th entry of the (N 1 × N 2 ) matrix X -→ M M M magnetization vector - → B B

Magnetic resonance imaging

Magnetic resonance imaging is the result of a series of scientific discoveries made all throughout the 20th century, starting with the early days of quantum physics. In 1900, Joseph Larmor described the precession of the magnetic moment of electrons and protons [Larmor, 1903, Larmor, 1900]. Wolfgang Pauli discovered the spin of electrons in 1924 [Pais, 1991] and the phenomenon of nuclear magnetic resonance is highlighted in the 1930s by Isidor Isaac Rabi with the nuclear spin measurement of Sodium by the method of Molecular Beams [Rabi andCohen, 1934, Rabi et al., 1938]. However, using a molecular beam in a vacuum makes individual nuclei isolated from each another and their environment. Until 1945, independent teams led by Felix Bloch at Stanford and Edward Mills Purcell at the Massachusetts Institute of Technology, simultaneously demonstrated NMR and they jointly received the Nobel Prize for Physics in 1952 [Bloch, 1946, Block et al., 1946, Purcell et al., 1946, Purcell, 1952]. In 1966, Richard R. Ernst (awarded the Nobel Prize in Chemistry in 1991) brought the use of the Fourier Transform to NMR [START_REF] Ernst | Application of fourier transform spectroscopy to magnetic resonance[END_REF]. Then in 1971, Raymond Vahan Damadian highlighted the differences in relaxation time between different biological tissues, and proposed to use NMR relaxometry to detect tumors [Damadian, 1971]. Finally in 1973, the concept of magnetic resonance imaging (MRI) was proposed by Paul Lauterbur [Lauterbur, 1973]. Then In 1977, Paul Lauterbur and Peter Mansfield presented the first image in the human body obtained by MRI [Geva, 2006]. They received the Nobel Prize in Medicine for their work in 2003. This section is going to review the principles of magnetic resonance imaging and the formation of an MR image. Magnetic resonance imaging is an imaging modality based on the phenomenon of nuclear magnetic resonance (NMR).

Nuclear Magnetic Resonance (NMR)

1.1.1.1 Spin and free induction decay (FID) NMR technology originated from observation of atomic nuclei with non-zero spin quantum numbers. As early as the 1940s, the technology related to accurate measurement of NMR proposed by Felix Bloch of Stanford University aroused people's interest in NMR technology. At the same time, Felix Bloch proposed the theoretical model of NMR: the spin-charged particles like hydrogen nuclei will generate a magnetic field, and the distribution of this magnetic field is very similar to that of a strip magnet [START_REF] Hashemi | MRI: the basics: The Basics[END_REF]. The spin motion of the nucleus can be regarded as a positive charge rotating at high speed and generating a rotating ring current. At the same time, the rotating ring current also generates the corresponding magnetic field, so that the atomic nucleus has a magnetic dipole moment µ, as shown in Figure 1.1(a). In NMR, we take the magnetic dipole moment µ originating from the spin of the nucleus as a theoretical model. In nature, the nucleus of some specific elements can generate NMR signals and can be used as the research object of NMR, but we usually choose hydrogen nucleus 1 H as the research object in NMR technology, for two reasons: 1) the spin quantum number S of 1 H is 1/2, and it has uniform charge distribution leading to a relatively narrow NMR spectrum, which is easier to detect ; and 2) among several nuclei with S equal to 1/2, the hydrogen nucleus 1 H has the highest gyromagnetic ratio γ, leading to a stronger magnetic resonance signal in comparison with other nuclei.

We know that for a spinning gyro, when its rotation axis deviates from the direction of gravity, the gyro will spin along its own axis while rotating around the direction of gravity. When a nucleus with a magnetic dipole moment µ is placed in an external magnetic field B 0 , a similar phenomenon will occur due to the magnetic moment: the nucleus spins along its own axis, also rotates around the direction of the external magnetic field B 0 , as shown in Figure 1.1(b).

A central concept in NMR is the precession of the spin magnetization around the magnetic field at the nucleus, with the angular frequency ω = γB

(1.1)

where ω = 2πν relates to the oscillation frequency ν and B is the magnitude of the field. Thus, the precession of non-equilibrium magnetization in the applied magnetic field B 0 occurs with the frequency ω L = 2πν L = γB 0 , called Larmor frequency.

There is a large amount of water (H 2 O) and fat (-CH 2 -) in the human body, both of which contain a large number of hydrogen nuclei. Before the external magnetic field B 0 is applied to the human body, the magnetic dipole moments of hydrogen nucleus protons in the body are randomly arranged, and the magnetization vector (the sum of the magnetic dipole moment vectors in a unit volume) is zero, externally not magnetic, see Figure 1.2 (left). The nuclei are arranged in disorder without an external magnetic (left), and rearranged in order after an external magnetic B 0 applied (right).

When a uniform external magnetic field B 0 is applied, the original energy level is divided into 2S + 1 energy levels (S is the spin quantum number), and the corresponding nuclei are rearranged along the direction of the external magnetic field. For hydrogen nuclei, there will be two energy levels (2S + 1 = 2), as shown in Figure 1.2 (right). The spin nucleus in the higher energy state has the opposite orientation of magnetic dipole moment µ as B 0 , and the lower energy state has the same orientation as B 0 . In the thermal equilibrium state, the number of particles in the higher energy state and the lower energy state is not equal and it is given by the Boltzmann distribution. Thus, the proton system in a uniform external magnetic field will form a magnetization vector M 0 along the direction of the external magnetic field.

After applying an external magnetic field B 0 to a proton system such as the human body, the system will be in an equilibrium state with a longitudinal magnetization vector M 0 that does not change with time. After a radiofrequency (RF) pulse is applied, the equilibrium state of the spin-mass subsystem is perturbed, and the nuclei will absorb a certain amount of energy and be excited to the higher energy level. At this time, the system is in an unbalanced state.

Chapter 1

Bloch equations

To calculate the variation of the nuclear magnetization vector -→ M M M = (M x , M y , M z ) over time when relaxation times T1 and T2 are present. Here the relaxation times T1 and T2 will be introduced later in this subsection. The Bloch equations were introduced by Felix Bloch in 1946 [Bloch, 1946].

There are two factors that change the magnetization vector -→ M M M : 1) the moment -→ M M M × -→ B B B and 2) the spontaneous relaxation process inside the spin system. This phenomena can be modeled by the following equation:

d -→ M M M dt = γ( -→ M M M × - → B B B ) + (relaxation term) (1.2)
represented in the form of matrices:

d dt   M x M y M z   = γ - → i - → j - → k M x M y M z B x B y B z -   1 T2 0 0 0 1 T2 0 0 0 1 T1     M x M y M z -M 0   (1.3)
The rotating coordinate system used in NMR is a coordinate system that rotates around the z-axis (the direction of B 0 ) at the Larmor frequency ω of the RF field. Thus, the Bloch equations in the rotating coordinate system can be written as:

d dt   M r x M r y M r z   = - → i - → j - → k M x M y M z γB r x γB r y γB 0 -ω -   1 T2 0 0 0 1 T2 0 0 0 1 T1     M x M y M z -M 0   (1.4)
When the RF pulse disappears, -→ B B B = B 0 -→ k , reaching resonance, γB 0 -ω = 0, we can get:

d dt   M x M y M z   = -    Mx T2 My T2 Mz-M 0 T1    (1.5)
for xy plane, M xy = M 2 x + M 2 y , and the solution is:

M z (t) = M 0 (1 -e -t T1 ) + M (0)e -t T1 M xy (t) = M xy (0)e -t T2 (1.6)
This simple solution of the Bloch equations is based on the assumptions: 1) the nuclear magnetization is exposed to a constant external magnetic field B 0 in the z direction and 2) the RF pulse is applied at the Larmor frequency ω 0 = γB 0 . It is worth noting that these are often necessary conditions in practical applications of NMR.

We call the magnetization vector in the direction of the external magnetic field (z direction) as the longitudinal magnetization vector M z , and the magnetization vector perpendicular to the When an RF pulse is applied, the longitudinal magnetization vector M z decreases from M 0 to M z , and the transverse magnetization vector M xy changes from zero to M xy . Thus, under the action of the RF pulse, the magnetization vector M 0 is flipped by a certain angle θ in the xy plane, as illustrated in Figure 1.3. We call this kind of RF pulse, that makes the magnetization vector M 0 flip a θ angle, a θ angle pulse. In NMR technology, the 90°pulse is often used and the corresponding longitudinal magnetization vector M 0 will be totally flipped to the xy plane. When the RF pulse disappears, due to the energy exchange between particles or between the particles and the environment, the nucleus will recover from the higher energy level to the lower energy level, and the system will return to the original equilibrium state. The above phenomenon is called a relaxation process, and is divided into a T1 or longitudinal relaxation process and a T2 or transverse relaxation process. Taking the 90°pulse as an example, the macroscopic behavior of the corresponding longitudinal relaxation process is that the longitudinal magnetization vector M z is restored from zero to the maximum value M 0 , and the time of the longitudinal relaxation process is usually expressed by T1. With respect to the transverse relaxation process, its behavior is that the transverse magnetization vector M xy gradually decreases from M x,y to zero, and the time of the transverse relaxation process is usually expressed by T2.

As shown in Figure 1.4, given the T1 and T2 500 ms and 50 respectively, the longitudinal magnetization (red) and transverse magnetization (blue) increases and increases during the relaxation process.

In Table 1.1 examples of T1 and T2 values of different tissues at 1.5 T are shown [START_REF] Plewes | Physics of mri: a primer[END_REF]. Still taking the 90°pulse as an example, considering the example from the 3-D perspective, the magnetization vector M 0 is first flipped to the xy plane, and the transverse magnetization vector M xy = M 0 at this time. Then M xy rotates around the direction of the external magnetic field, the rotation frequency is the Larmor precession frequency (also the resonant frequency), and the M xy starts to decay from M 0 according to an exponential law, as shown in Figure 1.5. The rotating transverse magnetization vector M xy will generate an induced current in the receiving coils, and this induced current is the NMR signal of the corresponding sample tissue.

Formation of a magnetic resonance (MR) image

The NMR signal is acquired in the frequency domain, which is so-called K-space [START_REF] Schick | Proton relaxation times in human red bone marrow by volume selective magnetic resonance spectroscopy[END_REF], through the application of magnetic fields whose amplitude varies linearly in a direction of space and superimposed on the main magnetic field. When applied, these magnetic field gradients will cause a local variation both in the Larmor frequency and in the phase of the spins according to their location in space. The information of the NMR signal acquired in space is spatially encoded both in frequency and in phase. The image is then reconstructed via an Inverse Fourier transform. The reconstructed image is an assembly of rectangular parallelepipeds, which are called voxels.

Spatial information encoding

In a uniform magnetic field the NMR signal does not contain spatial information. The encoding of spatial information is done by superimposing three magnetic field gradients, which are orthogonal, on the main magnetic field B 0 . These gradients are used to select the slice and its thickness, and then to encode the spatial information of the voxels of the sample via frequency and phase encoding.

Slice selection

Taking an example of an axial slice, to encode spatial information in the z-direction, a gradient is applied in the z-direction, resulting in a linear change in Larmor frequencies as described above. The first gradient, is so-called slice selection gradient. A scheme showing how the slice selection gradient works is shown in Figure 1.6. The slice selection gradient is applied simultaneously with the excitation RF pulse, which selectively excites the spins of a certain slice. The position of the excited slice depends on the center frequency of the RF pulse, and the thickness of the excited slice depends on the bandwidth of the RF pulse and the amplitude of the slice selection gradient. Taking Figure 1.6 as an example, the center frequency of the RF pulse is γ (B 1 -B 2 ) 2 and the bandwidth of the RF pulse is γ(B 1 -B 2 ).

Frequency encoding and phase encoding

The basic idea of frequency encoding is to apply a second gradient in the x-direction during the signal acquisition process. As shown in Figure 1.7(a), since the main magnetic field is changing linearly in the x-direction because of the gradient G x , each vertical line of the slice along the x-axis will have its own Larmor frequency (f L = γB) and the amplitude of the signal should be proportional to the combination of the number of protons in the vertical strip, and related to the T1 and T2 relaxation [START_REF] Mcrobbie | Seeing is believing: introduction to image contrast[END_REF]. Phase encoding uses the same idea as frequency encoding in that a phase encoding gradient G y is applied in the y-direction, resulting in a linear change in the Larmor frequency in the y-direction (see Figure 1.7(b)). However, in phase encoding the gradient G y is turned on for a period of time, and then turned off. The result is that regions of higher Larmor frequencies will have net magnetizations rotating over a larger angle than those regions of lower Larmor Chapter 1 1.1. Magnetic resonance imaging frequencies. These angles, corresponding to the phases, which are dependent on the gradient G y and thus space, can be recorded. If we then apply the phase-encoding gradient for a longer amount of time we will have a longer average phase (see Figure 1.8). In this way we can obtain different average phases to encode the spatial information in the y-direction.

Figure 1.8: Effect of phase encoding gradient duration and Larmor frequency on phase Phase encoding along the y-axis using the gradient G y , with different activating time to record the corresponding phase information thus the spatial information.

K-Space and inverse Fourier transform

The K-space is essentially a plot of raw MRI data. The elements of K-space are stored as complex numbers. K-space holds raw data before reconstruction. To reconstruct an MR image, the key is to sample K-space. For example, in a Cartesian sampled K-space, a K-space line is acquired at one time using the phase encoding gradient with the same activation time of the gradient G y . The acquisition is repeated as many times as the lines to sweep the K-space. The filled K-space data is then used to reconstruct the MR image through an inverse Fourier transform.

Pulse sequences

A sequence defines the time dependent combination of the applied RF pulses and gradient pulses. The classification of MRI sequences is complex and continues to grow and expand with the development of technology. In this section, we are going to review two main types of sequences: the Spin Echo (SE) Sequence and the Gradient Recall Echo (GRE) Sequence.

Spin Echo Sequence

An example of a Spin Echo (SE) pulse sequence diagram is depicted in Figure 1.9. SE pulse sequences use a 90 • excitation pulse followed by a 180 • refocusing pulse to recover the T2 decay relaxation [START_REF] Bitar | Mr pulse sequences: what every radiologist wants to know but is afraid to ask[END_REF]. First a 90 • excitation pulse tips and phases all the precessing protons into the xy plane. After a period of time τ the precessing protons will dephase because , leading the slower protons ahead of the faster protons along the precessional trajectory. After another period of time τ , all the protons will rephase and recover some fraction of the original transverse magnetization. This is known as the spin echo and at this point the time TE is defined as shown in Figure 1.9, the time between the excitation pulse and the echo. The TR is defined as the repetition time (or time between excitation pulses) as illustrated in Figure 1.9. The G z and G x refer to the slice selecting gradient and frequency encoding gradient (or reading gradient) respectively. The phase encoding gradient G y is usually displayed as a step function where every TR it moves up by a step. The signal readout (or analog/digital (A/D) readout) occurs during the echo with the reading gradient G x . The acquired signal amplitude A echo can be described by [START_REF] Liang | Principles of magnetic resonance imaging: a signal processing perspective[END_REF]]

A echo = k • ρ 0 (1 -e -TR T1 )e -TE T2 (1.7)
where ρ 0 is the proton density (PD), k is a constant related to the experimental settings. The spin echo pulse sequences can be manipulated by varying TE and TR: e.g. 1) with short TE and long TR, the two exponential components (1 -e -TR T1 and e -TE T2 ) are similar to one in equation (1.7), the signal is PD-weighted; 2) with short TE and short TR the signal is T1-weighted; and 3) with long TE and long TR the signal is T2-weighted.

For in vivo applications, the SE sequence is often used together with accelerated acquisition strategies such as RARE, Turbo-RARE, SE-EPI sequences [START_REF] Bernstein | Handbook of MRI pulse sequences[END_REF].

Gradient Recall Echo Sequence

The gradient recall echo (GRE) sequence was developed for rapid imaging. In GRE sequences, the transverse magnetization M xy is tilted by an angle α (0 • < α < 90 • ) in the xy plane. Applying a flip angle less than 90 • allows for the reduction of the repetition time (TR) and thus the acquisition time, compared to a SE sequence. Instead of applying a 180 • refocusing pulse, GRE uses an additional reading gradient (x-direction) in opposite direction to generate the echo. An example of the GRE sequence diagram is depicted in Figure 1.10. GRE uses an α (0 • < α < 90 • ) RF pulse, followed by a bipolar gradient pulse (G x in Figure 1.10). This bipolar gradient pulse consists of a gradient with negative polarity followed by a rapid switch to a gradient with positive polarity [Blink, 2017]. In this case, the reversal of the gradient generates the rephasing of the protons spins. The GRE sequence is much faster than the SE sequence. Thus, it can be used for certain clinical scans that require a much quicker scanning time, e.g. perfusion MR imaging [Blink, 2017]. In contrast to the SE sequence, the GRE sequence rephases only a part of transverse magnetization, and the measured signal is weighted by a transverse relaxation time called T2 * (see Figure 1.10). The T2 * dephasing is faster than T2 dephasing. As a matter of fact, the T2 * dephasing is affected by the static B 0 inhomogeneities which can be compensated by the 180 • pulse in the spin echo sequence.

The classification of MRI sequences is complex and continues to grow with the development of technology. Nevertheless, the SE and GRE are still the two primary sequences in the MRI pulse sequences family because many other sequences are generally the variants or combinations of the two basic sequences.

Diffusion MRI

Diffusion-weighted magnetic resonance imaging (DWI or DW-MRI) is an MRI technology that uses specific MRI sequences, which utilize the diffusion of water molecules to construct signal intensities sensitive to the water diffusion [Le Bihan and Breton, 1985, Merboldt et al., 1985, Taylor and Bushell, 1985]. This MRI technology allows for non-invasive mapping of the diffusion process of water molecules in the biological tissues. Molecular diffusion in tissues is not free, but reflects interactions with many obstacles, such as macromolecules, fibers, and membranes. Water molecule diffusion patterns can therefore reveal microscopic details about tissue architecture, either normal or in a diseased state. DW-MRI was developed in the 1980s, to generate the NMR signal from molecular diffusion through the diffusion sensitizing gradients in the sequence. DW-MRI provides the possibility to quantify the molecular diffusion via the calculation of the apparent diffusion coefficient (ADC). The first in vivo application of diffusion imaging was in 1986 by Le Bihan [Le Bihan et al., 1986]. Following the demonstration of its clinical potential for the early detection of ischemia in the early 1990s [START_REF] Moseley | Diffusion-weighted mr imaging of acute stroke: correlation with t2-weighted and magnetic susceptibility-enhanced mr imaging in cats[END_REF], Chien et al., 1992], diffusion imaging became a routine clinical practice very quickly. However, the use of diffusion imaging was limited to brain imaging and was less favorable on other organs for a long time. Throughout the years, technical improvements of MRI sequences and gradients allowed diffusion imaging applied on other anatomical regions, e.g. liver, kidney, spleen, pancreatic, gallbladder, prostate and etc..

Brownian motion and diffusion in vivo

Molecular diffusion corresponds to the random movement of molecules, caused in nature by their molecule thermal motion. To better understand the diffusion, we start from the Brownian motion. Motion trajectories (blue lines) of one single particle diffusing by 2-D and 3-D Brownian motion models.

Brownian motion

In 1827, while looking through the microscope at pollen immersed in water the botanist R. Brown discovers the phenomenon that the movement of small pollen particles is disordered and erratic. The Brownian motion denotes the random or disordered motion of particles suspended in a liquid originating from their collision with the fast-moving molecules in the liquid [Einstein, 1905, Einstein, 1906].

According to Einstein's theory [Einstein, 1905, Einstein, 1906], the mathematical formulations of the Brownian motion, can be written as:

   X 2 = 2Dt for 1-D motion model X 2 = 4Dt for 2-D motion model X 2 = 6Dt for 3-D motion model (1.8)
where X denotes the motion displacement, D is the diffusion coefficient representing the diffusivity, and t is the time. The corresponding cases of the motion models are shown in Figure 1.11 and Figure 1.12. Each model records 1000 motions of one object.

Diffusion in vivo and ADC

Considering the water molecules diffusing in a tissue, rather than in a free medium, they would also collide with the tissue elements (membranes, organelles, macromolecules and etc.). Thus the distance traveled during the diffusion in tissue would be reduced compared to in a free medium. The diffusion coefficient in vivo is so-called the apparent diffusion coefficient (ADC). The ADC depends on the organization of the tissue of interest. In tissues with a high cellularity, e.g. tumor tissue, diffusion of water molecules is restricted. On the other hand, in tissues where cellularity is lower, e.g. the necrotic tissues, the diffusion of water molecules is much freer. Measurement of the ADC therefore can provide information not only of the tissue structures, but also information of the tissue pathophysiological conditions.

Taking the liver cirrhosis as an example (Figure 1.13), the ADC decreases in the liver cirrhosis compared to the normal liver which means the measurements of ADC can be useful for diagnosing liver cirrhosis [START_REF] Gourtsoyianni | Respiratory gated diffusion-weighted imaging of the liver: value of apparent diffusion coefficient measurements in the differentiation between most commonly encountered benign and malignant focal liver lesions[END_REF]. In [START_REF] Jafar | Diffusion-weighted magnetic resonance imaging in cancer: Reported apparent diffusion coefficients, in-vitro and in-vivo reproducibility[END_REF], a review of 115 selected articles in which healthy extra-cranial body DW-MRI was presented, reported the ADC values in different organs and discussed their variation across different studies. The ADC values are listed in Table 1.2 with different ranges for each tissue involved in the 115 studies [START_REF] Jafar | Diffusion-weighted magnetic resonance imaging in cancer: Reported apparent diffusion coefficients, in-vitro and in-vivo reproducibility[END_REF]. Note that the range only covers the mean ADC reported by these studies. The detailed ADC values and corresponding study can be found in [START_REF] Jafar | Diffusion-weighted magnetic resonance imaging in cancer: Reported apparent diffusion coefficients, in-vitro and in-vivo reproducibility[END_REF] 1 Uterus contains endometrium, myometrium and cervix 2 FBT: fibroglandular breast tissue

Diffusion weighted imaging

Diffusion imaging is an MRI method that produces in vivo MR images of biological tissues sensitized with the local characteristics of molecular diffusion (generally water) [START_REF] Posse | Human brain: proton diffusion mr spectroscopy[END_REF].

Sensitize MRI images to diffusion

A notion called "field gradient pulse" was initially proposed by [START_REF] Stejskal | Spin diffusion measurements: spin echoes in the presence of a time-dependent field gradient[END_REF] to sensitize MR signal to diffusion. This field gradient pulse is activated to provide a magnetic field strength varying linearly. Thus, with the gradient pulse, the protons begin to precess at different rates, resulting in dispersion of the phase and signal loss. A second gradient pulse with the same magnitude but opposite direction is applied to refocus the protons spins. In this case, the signal acquired for the protons would be reduced related to the amount of diffusion:

S S 0 = e -γ 2 G 2 δ 2 (∆-δ 3 )•D (1.9)
where S 0 is the signal intensity without the diffusion weighting, S is the signal with the gradient activated, γ is the gyromagnetic ratio, G is the strength of the gradient pulse, δ is the duration of the pulse, ∆ is the time interval between the two gradient pulses, and D is the diffusion coefficient. For simplification, Le Bihan suggested a "b" value to represent all the gradient terms as following [Le Bihan and Breton, 1985]:

b = γ 2 G 2 δ 2 ∆ - δ 3 (1.10)
so that the signal attenuation becomes:

S S 0 = e -b•ADC (1.11)
Note that the b value depends only on the acquisition parameters. The diffusion coefficient D is replaced by an apparent diffusion coefficient, ADC, indicating that the diffusion is not free in tissues. In this way, the MR images are "weighted" by the diffusion. In the DW images, the signal is more attenuated the faster the diffusion is and the larger the b value is. However, those diffusion-weighted images are still also sensitive to T1 and T2 relaxation contrast, but it is beyond the scope of this thesis.

The basic pulse sequence of DWI

Edward Stejskal and John Tanner developed the pulsed gradient spin echo (PGSE) technique for diffusion weighted NMR signals, in the mid-1960s [START_REF] Stejskal | Spin diffusion measurements: spin echoes in the presence of a time-dependent field gradient[END_REF]. A simple diagram of the PGSE sequence, together with the phase changing of the spins (stationary spins and diffusing spins are included) are depicted in Figure 1.14. In the slice selection direction, double diffusion-sensitizing gradient (DG) are applied on either side of the 180 • inversion pulse, symmetrically. The aim of the double DG is to affect the phase of diffusing spins with the stationary spins unaffected.

Immediately following the second DG, an acquisition module is played out. In modern DW sequences, it is typically an echo-planar imaging (EPI) sequence. EPI is designed for rapid imaging, which uses rapidly oscillating phase and frequency gradients to generate multiple gradient echoes. The modern implementations of diffusion weighted imaging are still based on the PGSE technique, and with certain modifications. In a normal spin echo sequence, two diffusion-sensitizing gradients (DG) are applied on either side of the 180 • inversion pulse, symmetrically. In this case, the stationary spins are unaffected by the paired DG while the diffusing spins are dephased, leading to a signal loss. It should be noted that, the 180 • inversion pulse lasts for a very short time and the phase is reversed immediately. For a better understand of the phase changing, this period of time is enlarged.

Actually based on the pioneering work of [START_REF] Stejskal | Spin diffusion measurements: spin echoes in the presence of a time-dependent field gradient[END_REF], for a long time it has been a challenging task to integrate the diffusion encoding gradients to the conventional sequences. After the availability of EPI in the 1990s, the diffusion weighted imaging could be implemented in the field of clinical imaging [START_REF] Turner | Echo-planar imaging of intravoxel incoherent motion[END_REF], Le Bihan, 2014]. The diffusion sequences based on EPI were fast and solved the problems of motion artifacts. Typically, the single shot EPI sequences have been used over the years and in recent years a new sequence, the RESOLVE (readout segmentation of long variable echo train) sequence [START_REF] Porter | High resolution diffusion-weighted imaging using readout-segmented echo-planar imaging, parallel imaging and a two-dimensional navigator-based reacquisition[END_REF], was introduced. The segmented readout of the RESOLVE allows for an improved spatially resolution and the shortening of the EPI echo train. However, the disadvantage of the RESOLVE is the increase of the acquisition time in comparison to the single shot approach.

In Figure 1.15, an example of the K-space sampling schemes using a single-shot EPI sequence and a multi-shot (three shots) EPI sequence is shown. The RESOLVE allows for multi-shot readout segmented EPI acquisitions of the K-space (see Figure 1.15(b) as an example for three shots). Comparing to single-shot EPI, RESOLVE provides MR images with higher spatial resolution and reduced distortions while the cost is a longer acquisition time. IntraVoxel Incoherent Motion

IVIM: model and quantitative analysis

In this section, we are going to review the intravoxel incoherent motion (IVIM) imaging technique, including the basics of diffusion and perfusion in Section 2.1.1, the IVIM bi-exponential model in Section 2.1.2 and the quantitative analysis of IVIM in Section 2.1.3.

Diffusion and perfusion

IntraVoxel Incoherent Motion (IVIM) imaging was initially introduced and developed by Le Bihan [Le Bihan et al., 1986, Le Bihan et al., 1988, Le Bihan, 2008] for quantitative assessment of the microscopic translational motions of water molecules that could contribute to the signal acquired with DW-MRI. In Figure 2.1 an model of a biological tissue voxel is shown. In this model, the biological tissue contains two types of environments: molecular diffusion of water in the tissue (so-called 'true diffusion', or apparent diffusion), and microcirculation of blood in the capillary network (perfusion, or pseudo-diffusion). First, the diffusion in the tissue, characterized with apparent diffusion coefficient (ADC), is generally denoted by D. The true diffusion D can provide a mean to probe the tissue microstructure [Le Bihan and Breton, 1985]. Second, the perfusion in capillaries (obviously, perfusion is much faster than tissue diffusion), characterized with a pseudo-diffusion coefficient, is generally denoted by D * . Thus, D * > D.

Besides the IVIM technique, there are different techniques there are different techniques for measuring perfusion with MRI, such as the dynamic contrast-enhanced (DCE) MRI, dynamic susceptibility contrast imaing (DSC) and arterial spin labelling (ASL) [START_REF] Jahng | Perfusion magnetic resonance imaging: a comprehensive update on principles and techniques[END_REF]. Compared to the DCE-MRI and DSC, the IVIM and ASL techniques have the advantage of not relying on an injected contrast agent. The ASL technique was developed in 1992 [Koretsky, 2012] to quantify blood perfusion in brain. In other organs, it is rather challenging to apply the ASL technique. The advantage of IVIM is that it can be applied in a straightforward way in different organs of the body.

The true diffusion ADC can be used to probe the tissue microstructure and represents a Chapter 2 diagnostic biomarker for a number of diseases. The perfusion has become another important physiological parameter which is sensitive to tissue disease status. As shown in Figure 1.13, still taking the liver as an example, the liver cirrhosis decreases not only the true diffusion D but also the pseudo-diffusion (perfusion) D * . Actually a large number of IVIM-MRI studies have been performed recently on liver to investigate pathologies such as nonalcoholic fatty liver disease (NAFLD) and steatohepatitis (NASH), which are nowadays reaching epidemic levels worldwide [START_REF] Loomba | The global NAFLD epidemic[END_REF], Andreou et al., 2013, Dyvorne et al., 2013a, Hayashi et al., 2013, Pasquinelli et al., 2011, Guiu et al., 2012, Leporq et al., 2015].

Figure 2.1: The random motion of water in a voxel model Schematic representation of random water motion in a simple biological tissue voxel, where free diffusion component (in green, described by the apparent diffusion coefficient) is complemented by blood flowing in capillaries (in red, described by the pseudo-diffusion coefficient).

IVIM bi-exponential model

Using diffusion weighted MRI, with different diffusion sensitizing gradient parameters, denoted by b values (introduced above in Section 1.2.2), a set of images can be obtained, as illustrated in Figure 2.2(a). For one voxel of this region, the observed intensities for different b values can be modeled as a bi-exponential model, as shown in Figure 2.2(b), in two components: a slow decay component (in blue) and a fast decay component (in red). They are related to the diffusion and perfusion respectively. This model allows for the measurements of these parameters: the ADC D and the pseudo-diffusion coefficient D * and their weights (Amplitudes, A 1 and A 2 ). In mathematical formulation, the IVIM bi-exponential model is:

S(b) = A 1 e -bD + A 2 e -bD * (2.1)
where S(b) denotes the signal intensity for the chosen voxel with a certain b value, and the A 1 , A 2 are contributions of the diffusion and perfusion respectively.

Using the bi-exponential IVIM model in equation (2.1), clinical diagnosis is based on the quantification of three IVIM parameters. In addition to D and D * , the perfusion fraction PF (denoted by f , which is defined in equation (2.2)) can also be a biomarker of a number of tissue pathologies. The perfusion fraction characterizes the contribution of the vascularization to the signal intensity of a given voxel.

f = A 2 A 1 + A 2 × 100% (2.2)
Once the IVIM parameters are estimated, the maps for each parameter can be created, the parametric maps of the IVIM model are of great interest from diagnosis point of view. Bihan et al., 1986, Le Bihan et al., 1988], the microscopic perfusion effects have been largely ignored for over twenty years. However it still costs over 30 years to apply the IVIM imaging to the MRI clinical community, because of the lack of high image quality acquisition. Resulting from the wide application of IVIM imaging to clinical and research community, work on IVIM quantification algorithms gets more interest recently [START_REF] Suo | Intravoxel incoherent motion diffusion-weighted mr imaging of breast cancer at 3.0 tesla: Comparison of different curve-fitting methods[END_REF], Barbieri et al., 2016, Meeus et al., 2017, Gurney-Champion et al., 2018]. The IVIM model was originally applied in the brain [Le Bihan et al., 1986], and nowadays it has been more commonly used in the abdomen, where relatively high-perfused organs, such as the liver and kidneys, have highlighted its potential [START_REF] Luciani | Liver cirrhosis: intravoxel incoherent motion mr imaging-pilot study[END_REF], Yamada et al., 1999, Chandarana et al., 2012, Gaing et al., 2015].

Quantitative analysis of IVIM

Fitting algorithms for IVIM quantification

Algorithms which are commonly employed to quantify the IVIM parameters include leastsquares (LSQ) methods and "step-wise" techniques considering the diffusion signal at high and low b-values step by step.

Least-squares based algorithms

An important source of least squares problems is data fitting. More precisely, the goal is to find the optimal parameters of a predefined data model such that a lowest model error is Chapter 2 obtained. In the case of the IVIM model, least squares problem can be generally formulated as follows:

arg min

Θ b (S b -fb (Θ)) 2 (2.3)
where fb (Θ) denotes the signal intensities defined in equation (2.1) with certain set of b values, fb (Θ) is the predefined model of the signal intensity with the estimates parameters Θ (in this case, Θ = {D, f, D * }). According to equation (2.1), the above optimization problem is clearly non-linear as the objective function to be minimized is non-linear in any of the model parameters (i.e. D, f and D * ). To solve such non-linear optimization problem, adequate optimization algorithms such as Gauss-Newton, Levenberg-Marquardt (LM) [Levenberg, 1944, Marquardt, 1963] or Trust-Region (TR1 ) algorithms [START_REF] Byrd | Approximate solution of the trust region problem by minimization over two-dimensional subspaces[END_REF], Branch et al., 1999, Nocedal and Wright, 2006] can be used.

The LM algorithm has been employed in the IVIM quantification by several researchers [START_REF] Wittsack | Statistical evaluation of diffusion-weighted imaging of the human kidney[END_REF], Thoeny and De Keyzer, 2011, Mazaheri et al., 2012]. The LM algorithm can simultaneously estimate the three IVIM parameters f , D and D * for each voxel by a LSQ fit. This algorithm does not allow for boundary constraints of the fitted parameter values.

Similarly to LM, the TR algorithm is an LSQ algorithm and is used to simultaneously determine the three IVIM parameters. However, the search space is restricted to a subset of the domain of the cost function, which is also called trust region. Thus, with TR algorithm, boundary constraints of the parameter values can be incorporated. Actually in physical phenomena, the LM and TR are commonly used for without and with constraints situations for the exponential analysis [START_REF] Istratov | Exponential analysis in physical phenomena[END_REF].

As shown in Figure 2.3, an example of the fitting results with LM (A) and with TR (B) is shown. Generally in IVIM applications, the performances of LM and TR are comparative.

Step-wise fitting strategy

However, the most frequently used algorithm in IVIM quantification is a step-wise based method. Based on the fact that the pseudo-diffusion coefficient D * is much larger than the ADC D, the impact of the perfusion fraction on the DW MR signal at high b values (b ≥ b th , the threshold b value b th depends on the specific application and employed sequences) is considered negligible. In the case of Figure 2.4, 200 s/mm 2 was chosen to be the threshold that b th = 200 s/mm 2 . Thus, the algorithm starts from fitting the model

S(b ≥ b th ) = S intercept • e -bD
(2.4) to the DW MR signals at b values which are larger than b th . Here the S intercept is the intercept in the vertical axis, as shown in Figure 2.4, the green curve denotes the pure diffusion contribution, and hence the f can be calculated directly by:

f = S 0 -S intercept S 0 (2.5)
Figure 2.3: Comparative fitting performances of LM and TR An example of the fitting results using LM algorithm in (A) and using TR algorithm in (B). The input data is marked in blue squares and the fitting curve is marked in red. For this certain example, the ground truth for the parameters are f = 0.2, D = 0.001, D * = 0.05, SNR = 50.

And finally, a TR algorithm is employed to determine the pseudo-diffusion coefficient D * with the values of f and D fixed. The final fitting curve in Figure 2.4 is marked in red. The input data is marked in blue squares. The final fitting curve is marked in red and the curve in green denotes the contribution from pure diffusion. For this example, the ground truth for the parameters are f = 20%, D = 0.001 mm 2 /s, D * = 0.05 mm 2 /s.

Besides the introduced algorithms based on LSQ, the approaches based on Bayesian inference also got attention in recent years. The Bayesian-based fitting algorithm will be introduced in Section 4.

IVIM: model and quantitative analysis

Chapter 2

Parametric maps

Once the IVIM parameters are estimated, the maps for each parameter can be created. In Figure 2.5 an example of the parametric maps for human liver is shown [START_REF] Gambarota | Eliminating the blood-flow confounding effect in intravoxel incoherent motion (IVIM) using the non-negative least square analysis in liver[END_REF].

In Table 1.2, the estimates of ADC values reported in healthy volunteers were listed for several different organs. Note that the estimation values in Table 1.2 were based on the monoexponential model, where the perfusion was not taken into account. Now with the IVIM biexponential model, examples of previously reported IVIM parameters (D, f and D * ) in some healthy tissues and in some pathological tissues, are given in Table 2.1 [Le Bihan et al., 2018]. It should be noted that, the accurate and precise quantification of the perfusion coefficient D * is still a pending problem in the IVIM framework. Thus, the estimates of D * reported by previous studies are not listed in Table 2.1.

Figure 2.5: An example of parametric maps

These results were reported for one healthy volunteer in the liver with the IVIM model 1.2 [START_REF] Gambarota | Eliminating the blood-flow confounding effect in intravoxel incoherent motion (IVIM) using the non-negative least square analysis in liver[END_REF].

Chapter 2 2.1. IVIM: model and quantitative analysis

IVIM in liver

Liver disease is a major cause of illness and death worldwide [START_REF] Blachier | The burden of liver disease in europe: a review of available epidemiological data[END_REF], Byass, 2014, Nader et al., 2014, Wang et al., 2014, Pimpin et al., 2018, Asrani et al., 2019]. For example in China, liver diseases, primarily viral hepatitis (predominantly hepatitis B virus (HBV)), nonalcoholic fatty liver disease, and alcoholic liver disease, affect approximately 300 million people [START_REF] Wang | The global burden of liver disease: the major impact of china[END_REF]. As reported in [START_REF] Bray | Global cancer statistics 2018: Globocan estimates of incidence and mortality worldwide for 36 cancers in 185 countries[END_REF], among 36 cancers in 185 countries, the liver cancer has been the fourth leading cause of cancer death.

In this section, we are going to introduce the liver fibrosis (Section 2.2.1) and the assessment of liver tumors with IVIM imaging (Section 2.2.2).

Liver fibrosis

Fibrosis is a pathological wound healing in which connective tissue replaces normal parenchymal tissue to the extent that it goes unchecked, leading to considerable tissue remodelling and the formation of permanent scar tissue [Wynn, 2004, Birbrair et al., 2014]. Fibrosis can occur in many tissues in human body, e.g. lungs, liver, brain and etc..

Stages of liver fibrosis

Supposing that the liver parenchyma is damaged (caused by viral hepatitis, alcoholic or nonalcoholic steatohepatitis, or an autoimmune-related disorder), the fibroblasts will proliferate and produce collagen (inhibiting the molecular exchange between the sinusoids and hepatocytes) around the sinusoids. The changes in molecules diffusivity can be observed with DW-MRI or IVIM imaging, for the diagnosis of pathological status. Histopathologic assessment of liver fibrosis is usually divided into stages with the F score (with the METAVIR scoring system), ranging frome F0 to F4 [Group andBedossa, 1994, Bedossa andPoynard, 1996]:

• F0, no fibrosis.

• F1, portal fibrosis without septa;

• F2, portal fibrosis with few septa (significant fibrosis);

• F3, numerous septa without cirrhosis (severe fibrosis);

• F4, cirrhosis.

A model of the stages of liver fibrosis is depicted in Figure 2.6. The hepatocirrhosis causes liver failure which is typically irreversible, whereas liver fibrosis in early stages can be treatable [START_REF] Marcellin | Regression of cirrhosis during treatment with tenofovir disoproxil fumarate for chronic hepatitis b: a 5-year open-label follow-up study[END_REF]. Thus, it is crucial to detect the liver fibrosis in early stages for planning the therapy. Staging of hepatic fibrosis is important in the management of patients with chronic liver disease [START_REF] Tsukuma | Risk factors for hepatocellular carcinoma among patients with chronic liver disease[END_REF], Nakagawa et al., 2009]. The gold standard for staging liver fibrosis is liver biopsy followed by histopathology. However, biopsy may cause severe complications, e.g. hemorrhage and infection. There are also some inherent drawbacks in biopsy such as sampling error because of the small specimen sample, interobserver disagreement, and low patient acceptance [START_REF] Bravo | Liver biopsy[END_REF], Regev et al., 2002, Standish et al., 2006]. Therefore, noninvasive methods (e.g. imaging-based methods) have been developed for assessing hepatic fibrosis.

IVIM for staging liver fibrosis

Generally, the ADC decreases with the progression of liver fibrosis. Several studies [START_REF] Bedossa | An algorithm for the grading of activity in chronic hepatitis c[END_REF], Koinuma et al., 2005, Lewin et al., 2007, Taouli et al., 2007, Sandrasegaran et al., 2009, Do et al., 2010] have reported the reduced ADC values in liver cirrhosis compared to in normal liver. However, these reports were based on mono-exponential model fitting, in which the perfusion effects were not taken into account. It should be noted that the liver is a highly vascularized organ and hence the perfusion should not be neglected when evaluating the tissue diffusivity with diffusion imaging. As introduced previously, the IVIM imaging is bi-exponential model which is used for the estimation of both tissue diffusion and perfusion. The histopathological structures of the liver can be altered by fibrosis leading to heavy effects of the diffusivity measurements with IVIM. Thus, compared to the mono-exponential based DWI method, the IVIM DWI is thought to be better for precisely estimating the three parameters D, D * and PF (f ) to assess the liver fibrosis.

Assessment of liver lesions with IVIM

Besides evaluating the stages of liver fibrosis, IVIM is also suitable for assessment of focal liver lesions (FLL). FLL are generally divided into two types: benign or malignant lesions. The most commonly benign lesions in liver are hemangiomas, cysts, focal nodular hyperplasia (FNH) and hepatocellular adenoma (HCA) [START_REF] Kaltenbach | Prevalence of benign focal liver lesions: ultrasound investigation of 45,319 hospital patients[END_REF]. On the other hand, the most encountered malignant liver lesions are primary liver cancers, such as hepatocellular carcinoma (HCC), intrahepatic cholangiocarcinoma (ICC) and combined "biphenotypic" hepatocellular-cholangiocarcinoma (cHCC-ICC) [START_REF] Massarweh | Epidemiology of hepatocellular carcinoma and intrahepatic cholangiocarcinoma[END_REF]. To assess FLL, the current clinical standard method is multiphasic contrast-enhanced MRI (CEMRI) [START_REF] Coenegrachts | Improved focal liver lesion detection: comparison of single-shot diffusion-weighted echoplanar and single-shot t 2 weighted turbo spin echo techniques[END_REF]. However, there is an increasing interest in noninvasive and noncontrast methods, e.g. DW-MRI. With the diffusion imaging, IVIM has attracted more and more attention for detection and characterization of liver lesions. Based on technical considerations for acquisition of IVIM imaging, a summary of the recommendations suggested Chapter 2 by [START_REF] Hectors | Assessment of liver tumors with ivim diffusion-weighted imaging[END_REF] Ranges of the IVIM parameters estimates (mean values), reported by previous studies for both benign and malignant liver lesions, see [START_REF] Doblas | Determination of malignancy and characterization of hepatic tumor type with diffusion-weighted magnetic resonance imaging: comparison of apparent diffusion coefficient and intravoxel incoherent motion-derived measurements[END_REF], Yoon et al., 2014, Watanabe et al., 2014, Zhu et al., 2015, Wang et al., 2016, Klauss et al., 2016, Luo et al., 2017b, Luo et al., 2017a, Jerjir et al., 2017] 1 FNH: Focal Nodular Hyperplasia 2 HCA: HepatoCellular Adenoma 3 HCC: HepatoCellular Carcinoma 4 ICC: Intrahepatic CholangioCarcinoma

D (×10 -3 m 2 /s) D * (×10 -3 m 2 /s) PF (%) Benign (All) 1.

Limitations of IVIM

In state-of-the-art IVIM techniques, the estimation of the perfusion coefficient D * (sometimes also called D p , subscript p denoting "perfusion" or "pseudo"-diffusion) suffers from a weak repeatability and reproducibility, appearing especially prone to error [START_REF] Kakite | Hepatocellular carcinoma: short-term reproducibility of apparent diffusion coefficient and intravoxel incoherent motion parameters at 3.0 t[END_REF], Klauss et al., 2016, Joo et al., 2016, Hectors et al., 2016, Luo et al., 2017b, Luo et al., 2017a]. Based on considerations of the more complex model and more prone to errors in calculations, some studies ( [START_REF] Zhu | A comparative study of apparent diffusion coefficient and intravoxel incoherent motion-derived parameters for the characterization of common solid hepatic tumors[END_REF], Doblas et al., 2013, Klauss et al., 2016]) indicated that the IVIM is inferior to ADC in liver differentiation. Besides the weakness of precisely quantification of D * , the major limitation of IVIM is a combination of physiologic and hemodynamic variability and inherent technical limitations, such as the nonstandardization of scanning program, the acquisition technique, the b values' selection and also a variety of fitting algorithms. With respect to the variability of physiologic and hemodynamic, measurements on the left lobe of liver could be affected by cardiac or respiratory motions and bowel peristalsis [START_REF] Kwee | Influence of cardiac motion on diffusion-weighted magnetic resonance imaging of the liver[END_REF], Schmid-Tannwald et al., 2013, Lee et al., 2015].

Despite the flaws, as an imaging technique that can simultaneously probe the tissue diffusion and perfusion, IVIM still has great potential for a wide range of applications such as lesion characterization, assessment of tumor grade, prognostic evaluation and etc.. Furthermore, with the development of denoising algorithms and the progress of optimized and standardized protocols, the accuracy of the IVIM parameters estimation is steadily improving.

Blood vessel confounding effect

Gambarota [START_REF] Gambarota | Eliminating the blood-flow confounding effect in intravoxel incoherent motion (IVIM) using the non-negative least square analysis in liver[END_REF] and Cercueil [START_REF] Cercueil | Intravoxel incoherent motion diffusion-weighted imaging in the liver: comparison of mono-, bi-and triexponential modelling at 3.0-T[END_REF] investigated the presence of large blood vessels in liver and the confounding effect originating from them for generating the signal for IVIM images. Both [START_REF] Gambarota | Eliminating the blood-flow confounding effect in intravoxel incoherent motion (IVIM) using the non-negative least square analysis in liver[END_REF] and [START_REF] Cercueil | Intravoxel incoherent motion diffusion-weighted imaging in the liver: comparison of mono-, bi-and triexponential modelling at 3.0-T[END_REF] suggested to consider a tri-exponential IVIM-MRI model for post processing of the IVIM data. A tri-exponential IVIM model is introduced in Chapter 3 taking all voxels together into account. The SAVTE-IVIM algorithm is proposed for the post data processing to quantitatively evaluate the IVIM parameters. Besides the estimation of the D, D * and PF (f ), the proposed SAVTE-IVIM is also capable of identifying and assessing the blood vessels in a chosen ROI. Experiments and simulations were designed to evaluate the performances of the proposed SAVTE-IVIM algorithm.

IVIM in bone marrow

Bone marrow is the primary organ of new blood cell production (or haematopoiesis) and an important immune organ [START_REF] Arikan | Haematology of amphibians and reptiles: A review[END_REF]. Bone marrow is also one of the largest organs in human bodies, accounting for approximately 4% to 6% (for healthy grown-ups) of the total body weight [Vogler 3rd andMurphy, 1988, Hindorf et al., 2010]. In adult human bodies, bone marrow is primarily located in the ribs, vertebrae, sternum, and bones of the pelvis [Katherine, 2013]. As a semi-solid tissue, it can be found within the spongy or cancellous portions of bones [Farhi, 2009].

Besides the supportive vascular sinuses, collagen fibers, reticular cells, macrophages, and lymphocytes, bone marrow is primarily composed of hematopoietic tissue and adipocytes. An overview of the bone anatomy is depicted in Figure 2.7 [START_REF] Board | Chronic lymphocytic leukemia treatment (pdq®)[END_REF]. There are two types of bone marrow: red marrow and yellow marrow. It can be distinguished macroscopically on the basis of their colors. The red marrow is mainly responsible for hematopoiesis with a higher proportion of hematopoietic cells (≈ 60%, and adipocytes ≈ 40%) [START_REF] Ricci | Normal age-related patterns of cellular and fatty bone marrow distribution in the axial skeleton: Mr imaging study[END_REF]. The main component of yellow marrow is fat (adipocytes ≈ 80% in yellow marrow) [START_REF] Ricci | Normal age-related patterns of cellular and fatty bone marrow distribution in the axial skeleton: Mr imaging study[END_REF].

Figure 2.7: Anatomy of the bone [START_REF] Board | Chronic lymphocytic leukemia treatment (pdq®)[END_REF] The bone is made up of compact bone, spongy bone, and bone marrow. Compact bone makes up the outer layer of the bone. Spongy bone is found mostly at the ends of bones and contains red marrow. Bone marrow is found in the center of most bones and has many blood vessels.

In this section, we review the up-to-date development of IVIM in human bone marrow.

Evaluation of bone marrow

Both aging and the hematopoietic needs can influence the relative balance of red and yellow marrow across lifetime [START_REF] Custer | [END_REF]Ahlfeldt, 1932, Dunnill et al., 1967]. In addition, aplastic anemia, malignancies and infections can damage or displace the normal bone marrow architecture. Furthermore, bone marrow can also be affected by various form of leukemia [START_REF] Bonnet | Human acute myeloid leukemia is organized as a hierarchy that originates from a primitive hematopoietic cell[END_REF].

For the evaluation of bone marrow, much information can be gleaned directly by testing the blood itself (drawn from a vein through phlebotomy). However, with respect to a number of conditions that involves bone marrow, such as leukemia, multiple myeloma, lymphoma, anemia, and pancytopenia, a bone marrow examination (the pathologic analysis of samples of bone marrow obtained via biopsy and bone marrow aspiration) is needed for the diagnosis.

Considering the increasing interest in the non-invasive diagnostic technologies, imaging plays a more and more important role in the bone marrow's evaluation.

Compared with CT and X-rays, MRI is more sensitive and specific for assessing bone composition and MRI has become the main imaging modality for the assessment of bone marrow. For example, in routine clinical practices, the conventional anatomical T1-weighted, T2weighted and contrast-enhanced imaging are generally employed for the diagnosis of bone marrow lesions and for the monitoring of patients treatment response [Luciani andRahmouni, 2017, Morone et al., 2017]. DW-MRI can also be used to acquire additional information at the tissular and cellular level within the bone marrow, for both detection and characterization of focal lesions. However, as shown in Figure 2.7, the capillary network in bone marrow can not be neglected [Board, 2020, Nombela-Arrieta andManz, 2017]. As a result, the fast movement of the water molecules within the bone marrow, induced by its vascularized stroma, will directly impact the signal observed in DWI, especially for low b values [START_REF] Yeung | Bone marrow diffusion in osteoporosis: evaluation with quantitative mr diffusion imaging[END_REF]. Therefore it is mandatory to apply the IVIM model to the evaluation of bone marrow and the first successful IVIM measurements in bone marrow was conducted in 2014 by [START_REF] Marchand | Mri quantification of diffusion and perfusion in bone marrow by intravoxel incoherent motion (ivim) and non-negative least square (nnls) analysis[END_REF].

Thanks to its ability to evaluate both diffusion and perfusion contributions at the cellular level, the IVIM DW-MRI, using a bi-exponential analysis and with multiple b values, has been recently applied for bone marrow analysis, with promising clinical applications [START_REF] Karampinos | Quantitative mri and spectroscopy of bone marrow[END_REF].

It should be noted that, certain pathologies can affect perfusion in bone marrow, e.g. the proliferation of tumors is generally accompanied by the angiogenesis (regeneration of blood vessels). A recent study on IVIM DWI [START_REF] Bourillon | Intravoxel incoherent motion diffusion-weighted imaging of multiple myeloma lesions: correlation with whole-body dynamic contrast agent-enhanced mr imaging[END_REF] has shown that microvascular density increases in patients with myeloma and decreases in remission, which can be a result of angiogenesis. Some dynamic MRI studies [START_REF] Griffith | Vertebral bone mineral density, marrow perfusion, and fat content in healthy men and men with osteoporosis: dynamic contrastenhanced mr imaging and mr spectroscopy[END_REF], Griffith et al., 2006] on the bone marrow also reported that in people with osteoporosis, the loss of bone mineral density (BMD) is related to the decrease of perfusion.

Thus, medullary perfusion can be a biomarker of the pathophysiological state of the bone marrow and is of interest for prognosis and patient follow-up.

Clinical applications of IVIM in bone marrow

Bone marrow IVIM in healthy volunteers

To date, there are few literature data on normal quantitative IVIM parameters observed in bone marrow:

The first study that takes the perfusion into account was carried out by Yeung et al. [START_REF] Yeung | Bone marrow diffusion in osteoporosis: evaluation with quantitative mr diffusion imaging[END_REF] in 2004. This study was conducted on 20 healthy subjects (with a mean age of 28 years old). Single-shot spin-echo echo-planar imaging (SSSE-EPI) sequence was employed and 10 b values (0,20,40,60,80,100,200,300,400,500) s/mm 2 were used. In this study, the perfusion-related parameters (D * and perfusion fraction PF) were not quantified.

Ten years later in 2014, the second study using IVIM on normal bone marrow was reported by [START_REF] Marchand | Mri quantification of diffusion and perfusion in bone marrow by intravoxel incoherent motion (ivim) and non-negative least square (nnls) analysis[END_REF]. In this study, all IVIM parameters were successfully quantified for the first time. Fifteen healthy female volunteers with the age ranging from 18 to 30 years old were observed. The protocol comprised the SSSE-EPI sequence with five b values (0, 50, 100, 200, 600) s/mm 2 . Data analysis was carried out with the non-negative least-squares (NNLS) method.

In 2015, Ohno et al. studied the IVIM quantification on 11 healthy volunteers (7 males and 4 females, with a mean age of 25.8 years old) [START_REF] Ohno | Evaluation of perfusion-related and true diffusion in vertebral bone marrow: a preliminary study[END_REF]. The sequence SSSE-EPI was used again, with 16 b values (0,10,20,30,50,80,120,180,250,350,500,800,1200,1800,2200,2600) s/mm 2 were used this time. In this study, Ohno et al. reported a positive correlation between the BMD and pseudo-diffusion D * and a negative correlation between fat fraction and diffusion coefficient D or perfusion fraction PF.

More recently, Lasbleiz et al. did IVIM quantification on six healthy volunteers (24.2 ± 4.3 years old), using the multishot, readout-segmented (RESOLVE) EPI sequence, with seven b values (0, 50, 100, 150, 400, 800, 1000) s/mm 2 [START_REF] Lasbleiz | Measurements of diffusion and perfusion in vertebral bone marrow using intravoxel incoherent motion (ivim) with multishot, readout-segmented (resolve) echo-planar imaging[END_REF]. The effects of different fat suppression techniques on the IVIM parameters measurements were also assessed.

All the four studies were carried out with an 1.5 T MR scanner. The quantified values for the IVIM parameters are illustrated below in Table 2.4.

For diseased bone marrow

Compared to the normal bone marrow in health, more (actually, only a little bit more) studies were carried out on abnormal bone marrow in disease.

For example, for the osteoporosis, in [START_REF] Yeung | Bone marrow diffusion in osteoporosis: evaluation with quantitative mr diffusion imaging[END_REF], besides the 20 healthy subjects, Yeung et al. also reported a population of 44 old women (with a mean age of 70 years old) with IVIM. The diffusion coefficient (ADC with mono-exponential or D with bi-exponential) values tended to decrease with the decay of bone marrow density. The authors suggested that the progressive fatty replacement of bone marrow (commonly observed in elderly patients) is generally accompanied by a decrease in the molecular diffusion.

Besides, IVIM are also applied on some other diseases such as diffuse bone marrow a In [START_REF] Yeung | Bone marrow diffusion in osteoporosis: evaluation with quantitative mr diffusion imaging[END_REF], an ADC value of 0.43 ± 0.08 × 10 -3 mm 2 /s within the bone marrow was calculated with mono-exponential model, the value of D was quantified with bi-exponential, but the values of D * and f were not quantified. b In [START_REF] Marchand | Mri quantification of diffusion and perfusion in bone marrow by intravoxel incoherent motion (ivim) and non-negative least square (nnls) analysis[END_REF], values listed in this table were calculated by NNLS method; In this study, STIR fat suppression was used. c In [START_REF] Ohno | Evaluation of perfusion-related and true diffusion in vertebral bone marrow: a preliminary study[END_REF], the detailed values were not provided, so a range (can be found in the figures depicted in this study) of the quantified parameters are listed. d In [START_REF] Lasbleiz | Measurements of diffusion and perfusion in vertebral bone marrow using intravoxel incoherent motion (ivim) with multishot, readout-segmented (resolve) echo-planar imaging[END_REF], SPAIR fat suppression was used.

infiltration [START_REF] Shah | Intravoxel incoherent motion imaging for assessment of bone marrow infiltration of monoclonal plasma cell diseases[END_REF], Bourillon et al., 2015, Niu et al., 2017] and focal bone lesions [START_REF] Bourillon | Intravoxel incoherent motion diffusion-weighted imaging of multiple myeloma lesions: correlation with whole-body dynamic contrast agent-enhanced mr imaging[END_REF], Baik et al., 2017, Park et al., 2017, Lim et al., 2018].

Furthermore, IVIM can also play a role for therapeutic response evaluation, which has been explored in several studies [START_REF] Bourillon | Intravoxel incoherent motion diffusion-weighted imaging of multiple myeloma lesions: correlation with whole-body dynamic contrast agent-enhanced mr imaging[END_REF], Reischauer et al., 2017, Lee et al., 2017a].

Technical challenges

Single-shot spin-echo EPI (SSSE-EPI) sequence is designed for optimal IVIM imaging thanking to its inherent ability to reduce motion artifacts. Three out the four IVIM studies on normal bone marrow have employed SSSE-EPI sequence as introduced previously. SSSE-EPI has been the most commonly used sequence in IVIM imaging. But for the application of bone marrow, technical challenges originating from the anatomical specificities of bone marrow need to be overcome.

• The off-resonance effects

It is well known that, bone marrow is one of the largest organs in human bodies, but another fact is that bone marrow is contained within a relatively small bone medullar cavity, surrounded by cortical bone and spread out over the entire skeleton. As a result, one needs a high spatial resolution to differentiate the bone marrow signal from adjacent cortical bones signal. In this way, the off-resonance effects, related to the magnetic susceptibility difference between the bone marrow and the cortical bones, are going to increase, leading to spatial distortion or blurring.

• Signal-to-Noise Ratio (SNR) Also, a higher spatial resolution means a lower SNR in diffusion MR images.

To reduce the artifacts resulting from off-resonance effects, parallel imaging can be used or try to reduce the field of view. Both of them can increase the bandwidth to improve the image quality [Bammer, 2003]. The latter one is mostly used in imaging spinal cord.

To achieve a sufficient SNR for DW images, the diffusion gradient settings are greatly important. For example, the number of excitations for higher b values can be increased to compensate for the loss of signal intensities.

Future outlook

For the imaging of bone marrow, IVIM is a promising technique. IVIM imaging can better reflect the changes observed within the bone marrow at the tissular and cellular level.

Technical challenges still have to be overcome in order to reduce off-resonance effects and increase the SNR of bone marrow images. To date, however, IVIM robustness and scan times are still insufficient for it to be fully integrated in daily routine practice.

An example of the IVIM DW-MR images is depicted in Figure 2.8.

Figure 2.8: An example of bone marrow images with IVIM DW-MRI

The diffusion images are acquired on the spine of a healthy volunteer using the RESOLVE sequence with the b values equal to 0, 400 and 1000 s/mm 2 (from left to right).

To investigate the feasibility of different fitting algorithms on the IVIM parameters quantification in bone marrow and to provide nice parametric maps of IVIM in bone marrow, we compared five methods on both simulated IVIM data and real IVIM data within bone marrow. The details are introduced in Chapter 4.

Chapter 3

Sparsity-based All-Voxel Tri-Exponential IVIM (SAVTE-IVIM) algorithm

Introduction

Given the great clinical interest in the IVIM-MRI method for quantifying tissue diffusion and perfusion, especially for ROIs where the identification of blood vessels is not evident due to insufficient contrast level, a new algorithm for automatic identification of the potential presence of blood vessels in a given ROI is proposed in this chapter. In addition to this identification, which is the main goal of this algorithm, contribution originating from the tissue diffusion, perfusion and from the blood vessels to the acquired signal intensity, is quantified for all voxels in the ROI, in one single step. This is accomplished with the all-voxel IVIM-MRI model introduced in Section 3.2. The proposed algorithm does not require any manual operations, such as image contrast enhancement commonly used in such context, which are operator-dependent and prone to errors. The proposed algorithm takes into account (i) the assumption of a sparse representation of the spatial distribution of large blood vessels in the considered ROI leading to the Sparsity-based All-Voxel Tri-Exponential IVIM (SAVTE-IVIM) algorithm; and (ii) the non-negativity constraint that is inherent to the IVIM-MRI model. To deal with this non-negativity constraints, two possible strategies are introduced: (i) a rough strategy leading to the SAVTE-IVIM R method presented in Section 3.3.1; and (ii) an embedded strategy leading to the SAVTE-IVIM E one presented in Section 3.3.2.

The performance of the proposed SAVTE-IVIM algorithm (in its two variants SAVTE-IVIM R and SAVTE-IVIM E ) is evaluated and compared to the one of the NNLS-based method [START_REF] Gambarota | Eliminating the blood-flow confounding effect in intravoxel incoherent motion (IVIM) using the non-negative least square analysis in liver[END_REF] in terms of identifying blood vessels and quantifying tissue diffusion and perfusion. A statistical analysis of obtained results in both identifying blood vessels and quantifying tissue diffusion and perfusion, is also performed. Both realistic simulated and real DW-MRI images of liver acquired from six healthy volunteers are used for the aforementioned numerical analysis.

The all-voxel tri-exponential IVIM model

The conventional bi-exponential model

Assume, for a given volunteer, that a set of M DW-MR images acquired using M different diffusion-sensitizing gradient strengths (denoted here by b m , m ∈ {1, • • • , M }) is available. Then, for a given voxel in a region-of-interest (ROI) of size (N 1 ×N 2 ) taken in the m-th DW-MR image, the conventional IVIM-MRI model is given by [Le Bihan et al., 1986]:

s (n 1 ,n 2 ) (b m ) = L=2 =1 a (n 1 ,n 2 ) e -bmd (n 1 ,n 2 ) + (n 1 ,n 2 ) (b m ) (3.1)
where s (n1 ,n 2 ) denotes the acquired signal intensity of the

(n 1 , n 2 )-th voxel, 1 ≤ n 1 ≤ N 1 , 1 ≤ n 2 ≤ N 2 , and (n 1 ,n 2 ) is an additive Rician noise. Coefficients d (n 1 ,n 2 ) 1 and d (n 1 ,n 2 ) 2 with d (n 1 ,n 2 ) 1 < d (n 1 ,n 2 ) 2
stand respectively for the ADC 1 (characterizing the slow exponential decay, around 10 -3 mm 2 /s) and the D * (the pseudo-diffusion coefficient, characterizing the fast exponential decay, around 10 -2 to 10 -1 mm 2 /s, originating from the blood perfusion in capillaries); a

(n 1 ,n 2 ) 1
, and a

(n 1 ,n 2 ) 2
are their corresponding amplitudes. According to equation (3.1), the IVIM-MRI model provides a means to simultaneously assess the ADC and the perfusion fraction (PF), denoted here by f , which is defined as a normalized amplitude of the fast diffusion component [START_REF] Gambarota | Eliminating the blood-flow confounding effect in intravoxel incoherent motion (IVIM) using the non-negative least square analysis in liver[END_REF] such that for the (n 1 , n 2 )-th voxel we have:

f (n 1 ,n 2 ) = a (n 1 ,n 2 ) 2 L =1 a (n 1 ,n 2 ) .
It should be noted that, the modern IVIM analysis of the tissue diffusion and perfusion on the DW-MR images is based on this kind of voxel-wise bi-exponential model fitting. However, this conventional bi-exponential IVIM-MRI model does not take into account the potential contribution of a large blood vessel in the considered

(n 1 , n 2 )-th voxel (1 ≤ n 1 ≤ N 1 , 1 ≤ n 2 ≤ N 2 ) to the acquired signal.
Indeed, as shown in Figure 3.1, where an example of retrohepatic vena cava anatomy on liver [START_REF] Araújo | Retrohepatic vena cava lesion: which we cannot forget?[END_REF] is depicted, fast blood flow in hepatic veins and arteries can contribute, in the voxel they affect, as a confounding effect factor to the acquired voxel intensity.

The issue of this confounding effect was deeply investigated in [START_REF] Gambarota | Eliminating the blood-flow confounding effect in intravoxel incoherent motion (IVIM) using the non-negative least square analysis in liver[END_REF] and [START_REF] Cercueil | Intravoxel incoherent motion diffusion-weighted imaging in the liver: comparison of mono-, bi-and triexponential modelling at 3.0-T[END_REF]. Both studies have considered the contribution of potential blood vessels. In [START_REF] Gambarota | Eliminating the blood-flow confounding effect in intravoxel incoherent motion (IVIM) using the non-negative least square analysis in liver[END_REF], the authors proposed to use the Non-Negative Least Squares (NNLS) algorithm to detect for a single voxel if there is a confounding effect from the blood vessels within a framework of tri-exponential model. From now on, the NNLS-based terminology will be used to refer to the approach proposed in [START_REF] Gambarota | Eliminating the blood-flow confounding effect in intravoxel incoherent motion (IVIM) using the non-negative least square analysis in liver[END_REF]. Diaphragmatic view of the liver and its three-dimensional anatomical relationships involving the vena cava, hepatic pedicle and its branches, biliary tree, and suprahepatic veins [START_REF] Araújo | Retrohepatic vena cava lesion: which we cannot forget?[END_REF]. The inferior vena cava (IVC) and three hepatic veins (HV) are the largest veins in liver.

3.2.2

The NNLS-based approach [START_REF] Gambarota | Eliminating the blood-flow confounding effect in intravoxel incoherent motion (IVIM) using the non-negative least square analysis in liver[END_REF] In recent years, the NNLS algorithm [START_REF] Lawson | Solving least square problems[END_REF] has been widely used to fit the bi-exponential IVIM-MRI model [START_REF] Marchand | Mri quantification of diffusion and perfusion in bone marrow by intravoxel incoherent motion (ivim) and non-negative least square (nnls) analysis[END_REF], Keil et al., 2017]. Recently, Gambarota et al. [START_REF] Gambarota | Eliminating the blood-flow confounding effect in intravoxel incoherent motion (IVIM) using the non-negative least square analysis in liver[END_REF] have employed the NNLS algorithm to detect the presence of potential confounding effect factor in the IVIM measurements of liver. This is possible by minimizing, based on equation (3.1) for an arbitrary L, the following non-negative least squares optimization problem:

P N : arg min a (n 1 ,n 2 ) ||s (n 1 ,n 2 ) -E (n 1 ,n 2 ) a (n 1 ,n 2 ) || 2 2 s.t. a (n 1 ,n 2 ) > 0 1 ≤ n 1 ≤ N 1 , 1 ≤ n 2 ≤ N 2 (3.2)
where the L-dimensional vector a

(n 1 ,n 2 ) = [a (n 1 ,n 2 ) 1 , a (n 1 ,n 2 ) 2 , • • • ., a (n 1 ,n 2 ) L ] T , E (n 1 ,n 2 ) (M × L) is a matrix whose (m, )-th entry, E (n 1 ,n 2 ) m
, is equal to e -bmd (n 1 ,n 2 ) , and

s (n 1 ,n 2 ) = [s (n 1 ,n 2 ) (b 1 ), s (n 1 ,n 2 ) (b 2 ), • • • , s (n 1 ,n 2 ) (b M )] T .
A brief description of the NNLS algorithm [START_REF] Lawson | Solving least square problems[END_REF] is given in Appendix A.

In practice, for the above optimization problem to be solved, the the number L of exponential terms and hence E (n 1 ,n 2 ) in the IVIM model should be fixed. Therefore, a set of diffusion coefficients d

(n 1 ,n 2 ) , 1 ≤ ≤ L, whose cardinal is equal to an overestimated value of L, denoted here by L init , is defined (here we choose L init = 250). Now with the L init candidates of diffusion coefficients being evaluated using the NNLS algorithm, an "NNLS spectrum" is generated as shown in Figure 3.2 [START_REF] Gambarota | Eliminating the blood-flow confounding effect in intravoxel incoherent motion (IVIM) using the non-negative least square analysis in liver[END_REF].

Figure 3.2: NNLS spectrum for detection of a third component Signal decays of three voxels, (a) a voxel close to a blood vessel, (b) and (c) two voxels in two regions where no macroscopic blood vessels appeared to present, are displayed on the left side; and the corresponding NNLS spectrum are displayed on the right side. With the NNLS spectrum, first voxel was determined affected with a confounding effect while the second and third ones were determined negative. It is normal to set a threshold of a peak's amplitude to decide if certain component exists (e.g., 10% of the signal intensity with b = 0 s/mm 2 for the third peak) [START_REF] Gambarota | Eliminating the blood-flow confounding effect in intravoxel incoherent motion (IVIM) using the non-negative least square analysis in liver[END_REF].

Each observed peak in this spectrum corresponds to a diffusion component in the IVIM model with a peak area standing for its related amplitude. Consequently, for each voxel being analyzed with the NNLS algorithm, the number of observed peaks in the NNLS spectrum reveals the number L of exponential decays in the IVIM-MRI model. According to [START_REF] Gambarota | Eliminating the blood-flow confounding effect in intravoxel incoherent motion (IVIM) using the non-negative least square analysis in liver[END_REF], the NNLS spectrum of a voxel affected by the presence of a blood vessel (Figure 3.2(a)) comprises, compared to the one of a non-affected voxel (Figure 3.2(b) and (c)), three peaks contributing to the voxel intensity. This means that a tri-exponential (e.g. L = 3 in equation 3.1) rather than a bi-exponential IVIM-MRI model is to be considered for each voxel affected by the presence of blood vessel. Now for each voxel labeled as affected by a blood vessel (a third exponential decay being detected using the NNLS algorithm), Gambarota et al. [START_REF] Gambarota | Eliminating the blood-flow confounding effect in intravoxel incoherent motion (IVIM) using the non-negative least square analysis in liver[END_REF] proposed a heuristic strategy to deal with the detected confounding effect. Indeed, the data point (s (n 1 ,n 2 ) (b 0 )) related to the contribution of blood vessels to the signal intensity (i.e. the component with a very fast exponential decay) is discarded from subsequent analysis. Then, a bi-exponential fitting using the Trust-Region (TR) approach is performed to estimate both the tissue diffusion and perfusion.

Algorithm 1 below summarizes, for the (n 1 , n 2 )-th voxel, 1 ≤ n 1 ≤ N 1 , 1 ≤ n 2 ≤ N 2 , the main steps of the NNLS-based approach proposed in [START_REF] Gambarota | Eliminating the blood-flow confounding effect in intravoxel incoherent motion (IVIM) using the non-negative least square analysis in liver[END_REF]:

Algorithm 1 Main steps of the NNLS-based approach [START_REF] Gambarota | Eliminating the blood-flow confounding effect in intravoxel incoherent motion (IVIM) using the non-negative least square analysis in liver[END_REF] INPUT:

s (n 1 ,n 2 ) , L init , a set of diffusion coefficient condidates d (n 1 ,n 2 ) , 1 ≤ ≤ L init , COUNT = 0. STEP 1: Create the matrix E (n 1 ,n 2 ) .
STEP 2: Solve the P N problem (3.2) using the NNLS algorithm (see Algorithm 4 in Appendix A).

STEP 3: Generate the NNLS spectrum.

STEP 4: For each observed peak in the NNLS spectrum, DO:

If the peak amplitude > threshold COUNT = COUNT + 1 If COUNT ≥ 3, GO TO STEP 5 END if END if STEP 5: Remove s (n 1 ,n 2 ) (b 0
) and GO TO STEP 6.

END For.

STEP 6: Do a bi-exponential fitting based on the Trust-Region (TR) technique.

The all-voxel tri-exponential model

As discussed in the above subsection, the contribution of a blood vessel to the measured signal intensity of a voxel can be modeled as a very fast exponential decay term that can be added to the conventional bi-exponential IVIM model leading to a tri-exponential IVIM one:

s (n 1 ,n 2 ) (b m ) = L=3 =1 a (n 1 ,n 2 ) e -bmd (n 1 ,n 2 ) + (n 1 ,n 2 ) (b m ) (3.3)
Thus, in addition to the slow (i.e., d

(n 1 ,n 2 ) 1
for ADC) and the fast (i.e., d

(n 1 ,n 2 ) 2
for D * ) exponential decays, the tri-exponential IVIM-MRI model (equation (3.1)) comprises a third very fast exponential decay (i.e., d

(n 1 ,n 2 ) 3 > 0.2 mm 2 /s) with related amplitude a (n 1 ,n 2 ) 3
. Despite the efficiency, to some extent, of methods proposed in [START_REF] Gambarota | Eliminating the blood-flow confounding effect in intravoxel incoherent motion (IVIM) using the non-negative least square analysis in liver[END_REF], Cercueil et al., 2015] in dealing with this confounding effect, none of them are able to consider possible prior information regarding the spatial distribution of the model parameters which might considerably improve the identification of blood vessels in the considered ROI. This is mainly due to the fact that the employed tri-exponential IVIM-MRI model (e.g., L = 3 in equation (3.1)) is, conventionally, a single voxel-wise model. To cope with this limitation, the latter model is extended, as shown hereafter, to an all-voxel tri-exponential IVIM-MRI one. With this extended version being considered, not only spatial prior regarding the model parameters can be employed, but also a simultaneous processing of all voxels in the considered ROI is henceforth possible. For convenience purposes, a vectorized version of any defined ROI in the DW-MR image will be adopted in the sequel. Such vectorization is performed using the vec operator which maps a matrix T (I × J) to an IJ-th dimensional vector whose i + (j -1)I-th component stands for the (i, j)-th entry of T . Consequently, its inverse function, denoted by unvec, is defined such that unvec(vec(T )) = T . Accordingly, based on equation (3.3) and for a given ROI in the m-th (m ∈ {1, • • • , M }) DW-MRI image, the all-voxel tri-exponential IVIM-MRI model can be written as:

s(b m ) = L=3 =1 a e -bmd + (b m ) (3.4)
where a = vec(A ) and d = vec(D ), ∀ ∈ {1, 2, 3} with matrices A l , D l of size (N 1 × N 2 ) denoting respectively the spatial distribution of a and d in the considered ROI, (b m ) denotes the N 1 N 2 -th dimensional vector of MRI noise related to the m-th image. Since spatial distributions A , D , ∈ {1, 2, 3} are invariant through the M DW-MR images [Le Bihan et al., 1986], the all-voxel tri-exponential IVIM-MRI model defined over the M DW-MR images can then be written as:

s = L=3 =1 (1 M ⊗ a ) e -b⊗d + (3.5)
where s is a N M -dimensional vector gathering the acquired signals of the N = N 1 N 2 voxels of the considered ROI over the M DW-MR images. The noise vector , is defined such that

= [ (b 1 ) T , (b 2 ) T , • • • , (b M ) T ] T , b = [b 1 , b 2 , • • • , b M ] T .

The SAVTE-IVIM algorithm

A reliable clinical diagnosis based on a given ROI in the liver, when the presence of blood vessels is questionable, requires an efficient way of detecting blood vessels. As clinical experts are generally interested in investigating the tissue status, delineating an ROI that is, to a large extent, free from the confounding blood vessel effect is mandatory. In this sense, a chosen ROI is considered as informative in terms of both ADC and PF when the spatial distribution of blood vessels, in this ROI, is sparse. More precisely, the spatial distribution, A 3 , of a 3 is a sparse matrix. As the spatial distribution of blood vessels in the ROI reflects its vascularization, the sparsity pattern of A 3 is invariant through the M DW-MR images. Therefore, identifying blood vessels in the considered ROI can be performed by solving the following model identification problem:

P1 : arg min a 1 ,a 2 ,a 3 ,d 1 ,d 2 ,d 3 ψ(a 1 , a 2 , a 3 , d 1 , d 2 , d 3 ) 2 2 +λ a 3 1 s.t. a > 0 and d > 0, ∀ ∈ {1, 2, 3} (3.6) where ψ(a 1 , a 2 , a 3 , d 1 , d 2 , d 3 ) = s -L=3 =1 (1 M ⊗ a )
e -b⊗d denotes the model fidelity, λ stands for a penalty parameter. For the sake of readability, arguments in the function ψ(a 1 , a 2 , a 3 , d 1 , d 2 , d 3 ) will be omitted in the sequel. Note that a LASSO (Least Absolute Shrinkage and Selection Operator) penalty term (i.e., 1 -norm) is used in P1 (3.6) to describe the sparsity of the spatial distribution A 3 of a 3 .

It is noteworthy that in clinical settings MR images are typically characterized by relatively high SNR values. In this case, least square like methods are still applicable in the case of a Rician noise, since the Gaussian distribution approximates well the Rician distribution for relatively high SNR values. However, in the case of relatively low-SNR MR images, treating the Rician distribution as a Gaussian distribution is no more valid and a specific treatment such as the one based on the Majorization-Minimization framework is necessary [Varadarajan andHaldar, 2015, El-Hajj et al., 2020].

The above constrained minimization problem, P1 (3.6), can be solved following the spirit of the well-known ADMM method which is based on the augmented Lagrangian technique [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF] (A brief description of this technique is given in Appendix B). Regarding the non-negativity constraints, rough and embedded strategies are proposed hearafter leading, as mentioned previously, to the SAVTE-IVIM R and the SAVTE-IVIM E algorithms, respectively.

The rough strategy

The P1 problem (3.6) is solved here using the ADMM method [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF] which introduces the latent variable z and replaces the P1 problem (3.6) by:

P2 : arg min a 1 ,a 2 ,a 3 ,d 1 ,d 2 ,d 3 ,z ψ 2 2 + λ z 1 s.t. z = a 3 , a > 0, d > 0, ∀ ∈ {1, 2, 3}
(3.7)

A typical scheme to solve P2 (3.7) using the ADMM method is based on the minimization of its associated augmented Lagrangian function, denoted here by L. This leads to rewrite P2 (3.7) as follows:

P3 : arg min a 1 ,a 2 ,a 3 ,d 1 ,d 2 ,d 3 ,z,y L = ψ 2 2 + λ z 1 + y T (a 3 -z) + ρ 2 a 3 -z 2 2 s.t.a > 0, d > 0, ∀ ∈ {1, 2, 3} (3.8)
where the N -dimensional vector y stands for the Lagrangian multiplier and ρ > 0 denotes a penalty parameter. The estimation of a , ∀ ∈ {1, 2, 3}, can be then performed in a least squares sense by computing the stationary points of L in a , ∀ ∈ {1, 2, 3}. Regarding the estimation of d 1 , d 2 and d 3 , it is performed using the Levenberg-Marquardt (LM) algorithm [START_REF] Nocedal | Numerical Optimization[END_REF]. Note that instead of using the LM algorithm, other alternative non-linear methods such as the Gauss-Newton (GN) [START_REF] Boyd | Convex optimization[END_REF]] can be employed instead, as shown in [Liu et al., 2017]. However, the LM algorithm yielded a higher performance in terms of the estimation quality compared to the GN algorithm [Liu et al., 2017]. Therefore, the SAVTE-IVIM algorithm using the LM will be discussed in this chapter.

As far as the estimation of the latent variable z and the Lagrange multiplier y is concerned, a proximal operator [START_REF] Komodakis | Playing with Duality: An overview of recent primal-dual approaches for solving large-scale optimization problems[END_REF] is used to estimate the former while the dual ascent scheme [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF] is used to estimate the latter. Mathematical derivations of all update rules employed to solve the P3 problem (3.8), in the framework of the ADMM method, are given in Appendix C. At each iteration of the SAVTE-IVIM R algorithm, variables are updated in an alternative way. Indeed, each variable is estimated while keeping the other variables fixed to their last estimates. The algorithm stops when the model error exhibits, between two successive iterations, a value that is smaller than or equal to a predefined threshold or when a maximum number of iterations is reached. As for the non-negativity constraints in the optimization problem P2/P3, a rough strategy is used to deal with it. More precisely, potential negative components in a 1 , a 2 and a 3 encountered during iterations are set to a value in the neighborhood of zero (e.g., 10 -5 ). The choice of the latter threshold value is justified since no physiological prior information regarding the value of a , ∀ ∈ {1, 2, 3} is to be considered. Regarding possible negative components of d 1 and d 2 , they are projected back to their respective physiological ranges d II. While the stop criterion is not fulfilled or the maximum number of iterations is not reached, do:

(1) Update J a = ∂ψ ∂a T , ∀ ∈ {1, 2, 3} using equation (C.5) (2) Update a 1 using equation (C.2)

• if a (n) 1 < 0 then a (n) 1 ← 10 -5 , ∀n ∈ {1, • • • , N } (2) Update a 2 using equation (C.3) • if a (n) 2 < 0 then a (n) 2 ← 10 -5 , ∀n ∈ {1, • • • , N } (3) Update a 3 using equation (C.4) • if a (n) 3 < 0 then a (n) 3 ← 10 -5 , ∀n ∈ {1, • • • , N } (4) Update d using equation (C.6) • if d (n) 1 < 0 then d (n) 1 ← 10 -5 else if d (n) 1 > 0.01 then d (n) 1 ← 0.01, ∀n ∈ {1, • • • , N } • if d (n) 2 < 0.01 then d (n) 2 ← 0.01 else if d (n) 2 > 0.2 then d (n) 2 ← 0.2, ∀n ∈ {1, • • • , N } • if d (n) 3 < 0.2 then d (n) 3 ← 0.2, ∀n ∈ {1, • • • , N }
(5) Update z using equation (C.7) (6) Update y using equation (C.8)

The embedded strategy

The SAVTE-IVIM E algorithm deals with non-negativity constraints in P1 (3.6), by resorting to a change of variable into square such that: [START_REF] Coloigner | Line search and trust region strategies for canonical decomposition of semi-nonnegative semi-symmetric 3rd order tensors[END_REF]. As the sparsity patterns are invariant with respect to the Hadamard product (e.g a 3 and ã3 have the same sparsity pattern), P1 (3.6) can then be reformulated as follows:

a = ã ã = ã 2 , d = d d = d 2 with ã , d ∈ R * , ∀ ∈ {1, 2, 3}
P4 : arg min ã1 , ã2 , ã3 , d1 , d2 , d3 , z ψ 2 2 + λ z 1 s.t. z = ã3 (3.9) where ψ = s -L=3 =1 (1 M ⊗ ã 2 ) e -b⊗ d 2 .
Similarly to the SAVTE-IVIM R approach, the ADMM method is employed to solve the above optimization problem. To do so, P4 (3.9) is reformulated as follows:

P5 : arg min ã1 , ã2 , ã3 , d1 , d2 , d3 , z, ỹ L 1 = ψ 2 2 + λ z 1 + ỹT ( ã3 -z) + ρ 2 ã3 -z 2 2 (3.10)
where L 1 is its associated Lagrangian function. Parameters a , ∀ ∈ {1, 2, 3}, are estimated in a least squares sense by computing the stationary points of L 1 in a , ∀ ∈ {1, 2, 3}. Regarding the estimation of d , ∀ ∈ {1, 2, 3}, the LM method is also employed. Like the SAVTE-IVIM R , parameter estimation is also performed in an alternative way. Indeed, at each iteration, each parameter is updated while keeping the other ones fixed to their last estimates. The algorithm stops when the model error exhibits, between two successive iterations, a value that is smaller than or equal to a predefined threshold or when a maximum number of iterations is reached.

Algorithm 3 below summarizes the different steps in the proposed SAVTE-IVIM E algorithm. Mathematical derivations of all update rules figured in are given in Appendix D.

Numerical complexity

The numerical complexity of the proposed SAVTE-IVIM algorithm in its two variants, the SAVTE-IVIM R and SAVTE-IVIM E , is expressed in numerical flop. A flop is defined as a multiplication followed by addition. However, since in practice more multiplications than additions are encountered, only multiplications are taken into account here. Thus, the numerical complexity per iteration of the SAVTE-IVIM R algorithm is equal to (11M + 30)N 3 + (6M + 9)N 2 + (39M + 5)N + 5 flops while the one of the SAVTE-IVIM E is equal to (22M + 29)N 3 + (6M + 9)N 2 + (39M + 12)N + 5 flops, where as mentioned previously, N denotes the number of voxels in the considered ROI, and M stands for the number of b-values.

Simulations

Simulations on both fully simulated data and realistic simulated data generated from in vivo data are conducted to evaluate the behavior of the proposed SAVTE-IVIM algorithm in its two variants, the SAVTE-IVIM R and the SAVTE-IVIM E , compared to the NNLS-based method [START_REF] Gambarota | Eliminating the blood-flow confounding effect in intravoxel incoherent motion (IVIM) using the non-negative least square analysis in liver[END_REF] recently proposed to deal with the blood vessel confounding effect in the IVIM-MRI model.

Selection of penalty parameters

The penalty parameters ρ and λ in equations (3.8) and (3.10), are chosen in a trial-and-error manner such that a smallest model error is obtained. Indeed, the model error is evaluated over a grid of (ρ, λ)-values with (ρ, λ) ∈ {0.1, 0.2, 0.5, 1, 10} × {10 -8 , 10 -6 , 10 -4 , 10 -2 , 1, 10}, for both noise-free and noisy data and for different ROIs taken from the six volunteers in our data set. Regarding the noisy data configuration, three different values of the Signal-to-Noise Ratio (SNR) are considered (i.e., 50, 100 and 150). Note that, from an MR point of view, the SNR is defined as the ratio of the mean intensity of the chosen ROI to the standard deviation of the intensity outside the tissue [Lipton, 2008]. According to our numerical experiments, the couple (ρ, λ) that fulfils the smallest model error condition is (0.2, 10 -6 ) and hence it is retained for all subsequent experiments. It is worth noting that the quality of the identification results is highly dependent on the choice of the couple (ρ, λ).

To illustrate this fact, Figure 3.3 shows the model error per voxel, e v = ||ψ|| 2 /N as a function of (ρ, λ) values and for the aforementioned different situations of noise level (noise-free and noisy data with SNR values equal to 50,100 and 150). More precisely, the model error is evaluated for different realistic ROIs of size (16 × 16) taken from real MR images with a bifurcated blood vessel being simulated as shown in subsection 3.4.3. The model error is then averaged over the considered ROIs. Besides, for noisy data, results are also averaged over 25 Monte-Carlo (MC) trials. According to this figure, there is a set of (ρ, λ) values for which the model error per voxel, e v , exhibits relatively small values. Among the latter set and whatever the noise level being studied here is, the smallest averaged model error per voxel is obtained for (ρ, λ) = (0.2, 10 -6 ), as illustrated in dark square.

Fully simulated data

In this study, sixteen ROIs of size (16 × 16) voxels were designed. For each voxel in a designed ROI, a tri-exponential IVIM model was simulated, with simulation parameters chosen in the normal physiological ranges (i.e., d 1 < 0.01, 0.01 < d 2 < 0.2 and d 3 > 0.2), in order to obtain parametric maps of ADC (d 1 ), D * (d 2 ) and PF (f ) as depicted in Figure 3 , 10, 20, 40, 60, 80, 100, 200, 300, 400, 600, 800} s/mm 2 is used).

Regarding the contribution of the third exponential decay to the IVIM model, its amplitude a 3 was chosen randomly with values between 10% to 30% of the maximum voxel intensity in such a way a bifurcated blood vessel is modeled as shown in Figure 3.5. Rician noise was then added to each voxel [START_REF] Gudbjartsson | The Rician Distribution of Noisy MRI data[END_REF]. More precisely, the noise variance was adjusted to obtain the desired SNR values of 50, 100 and 150. Experiments were conducted over 15 Monte Carlo (MC) trials and the results of these trials were averaged to obtain the estimates of the simulated parametric maps and bifurcated blood vessels.

The above procedure (parametric map generation and bifurcated blood vessel simulation) was repeated 16 times to simulate 16 different noisy ROIs (e.g., 16 "simulated" volunteers). It should be noted that, the designed ground truth including both the parameters distribution and amplitude is changing for the 16 simulated ROIs.

Realistic simulated data

Besides the fully simulation, a realistic simulation based on in vivo data was conducted to evaluate the behavior of the proposed SAVTE-IVIM R and SAVTE-IVIM E approaches. For each volunteer, M realistic noisy DW-MR images are generated from the M available in vivo DW-MR images by performing a conventional (voxel-by-voxel) bi-exponential fitting. Then, a third exponential decay effect with a very fast diffusion coefficient (i.e. d The IVIM DW-MR images were acquired from six healthy volunteers (aging from 23 to 28, five males and one female). Table 3.1 lists the parameter settings for the data acquisitions which were conducted on a 3.0 T MR scanner (GE Discovery MR 750, GEHC, Milwaukee, WI).

Evaluation criteria

Voxel localization error

As the main objective of the proposed SAVTE-IVIM algorithm is to identify blood vessels in the considered ROI, the identification quality is evaluated here in terms of the Voxel Localization Error (VLE) criterion, which provides a measure of similarity between the original and the estimated configuration. Inspired from [START_REF] Becker | SISSY: an efficient and automatic algorithm for the analysis of EEG sources based on structured sparsity[END_REF], the VLE is defined here as:

V LE= 1 2Q t∈I min w∈ Î r t -r w 2 + 1 2 Q w∈ Î min t∈I r t -r w 2 (3.11)
where I and Î denote, respectively, the ground truth and the estimated sets of indices of voxels affetcted by the third diffusion decaying component; Q and Q are the cardinals of I and Î, respectively; and r , ∈ {t, w} denotes the position of the -th voxel.

Normalized mean square error

In addition to blood vessel identification, the SAVTE-IVIM allows not only for a quantification of the blood vessel contribution to the IVIM model but also for a quantification of the tissue diffusion and perfusion. Hence, the quantification quality of tissue diffusion and perfusion is evaluated here using the Normalized Mean Square Error (NMSE) defined by:

NMSE = H h=1 p -ph 2 2 H p 2 2 (3.12)
where p and ph denote, respectively, the ground truth of the target parameter (ADC (d 1 ), D * (d 2 ) or PF (f )) and its estimate computed at the h-th (h ∈ {1, ..., H}) MC trial.

Results

The performance of the proposed SAVTE-IVIM algorithm (in its two variants SAVTE-IVIM R and SAVTE-IVIM E ) is evaluated and compared to the one of the NNLS-based method [START_REF] Gambarota | Eliminating the blood-flow confounding effect in intravoxel incoherent motion (IVIM) using the non-negative least square analysis in liver[END_REF] in terms of identifying blood vessels and quantifying tissue diffusion and perfusion. In this section, the results of the simulations designed in Section 3.4 are described in terms of the identification of blood vessels (Section 3.4.5.1) and quantification of IVIM parameters (Section 3.4.5.2). Furthermore, a statistical analysis of obtained results in both identifying blood vessels and quantifying tissue diffusion and perfusion, is also performed. All the results with respect to the identification, quantification and statistical analysis, are illustrated for fully simulated data and for realistic simulated data separately.

Chapter 3

Identification of blood vessels

Since the NNLS-based algorithm [START_REF] Gambarota | Eliminating the blood-flow confounding effect in intravoxel incoherent motion (IVIM) using the non-negative least square analysis in liver[END_REF] provides only an identification of the blood vessel confounding effect, the spatial distribution of a 3 is binary. Thus, for a fair comparison of the latter method with the proposed SAVTE-IVIM R and SAVTE-IVIM E ones, which provide an identification and a quantification of the blood vessel confounding effect, a binary spatial distribution of a 3 obtained using the SAVTE-IVIM R and the SAVTE-IVIM E algorithms is created. To this end, the (n 1 , n 2 )-th voxel is labelled as affected when the estimated a 3.4. Simulations SNR values. Indeed, for a relatively low SNR value (SNR= 50), the NNLS-based method [START_REF] Gambarota | Eliminating the blood-flow confounding effect in intravoxel incoherent motion (IVIM) using the non-negative least square analysis in liver[END_REF] shows a high false negative rate in terms of identifying affected voxels while the proposed methods succeed in providing a localization map that is consistent with the ground truth. For higher SNR values (i.e., 100 and 150), better identification of blood vessels is to be noticed for the three considered methods, but with higher performance of the proposed approaches over the NNLS-based one [START_REF] Gambarota | Eliminating the blood-flow confounding effect in intravoxel incoherent motion (IVIM) using the non-negative least square analysis in liver[END_REF]. These results can also be confirmed using the VLE criterion (3.11) as depicted in Figures 3.8. As can be observed, compared to the NNLS-based method [START_REF] Gambarota | Eliminating the blood-flow confounding effect in intravoxel incoherent motion (IVIM) using the non-negative least square analysis in liver[END_REF], smaller VLE values are generally obtained using the proposed methods for all SNR values, for all 16 ROIs. As shown in Figures 3.9, the two proposed SAVTE-IVIM R and SAVTE-IVIM E approaches show better blood vessels identification results compared to the NNLS-based one for all SNR values. Indeed, for a relatively low SNR value (SNR= 50), the NNLS-based method [START_REF] Gambarota | Eliminating the blood-flow confounding effect in intravoxel incoherent motion (IVIM) using the non-negative least square analysis in liver[END_REF] shows a high false negative rate in terms of identifying affected voxels while the proposed methods succeed in providing a localization map that is consistent with the ground truth.

(n 1 ,n 2 ) 3 , 1 ≤ n 1 ≤ N 1 , 1 ≤ n 2 ≤ N 2 is

Quantification of the IVIM model parameters

In addition to the identification of blood vessels, the proposed SAVTE-IVIM approach allows to quantify the contribution of (i) detected blood vessels and (ii) the tissue perfusion and diffusion to the IVIM model.

The comparison between the proposed SAVTE-IVIM method and NNLS-based one, with respect to the quantification of tissue diffusion and perfusion (specifically the three IVIM parameters ADC, D * and PF), is conducted in the form of parametric maps and the NMSE defined in equation (3.12).

Case of fully simulated data

In terms of fully simulated data, Figure 3.11 shows the estimated parametric maps of ADC (d 1 ), D * (d 2 ) and PF (f ) obtained using the SAVTE-IVIM R , SAVTE-IVIM E and the NNLSbased algorithms, for the lowest (i.e., 50) and the highest (i.e., 150) SNR value. According to Figure 3.11, the three considered methods show generally comparable quantification results of ADC (d 1 ) and D * (d 2 ) which are, to some extent, consistent with the ground truth (the fourth column). Regarding the quantification of PF, the NNLS-based method [START_REF] Gambarota | Eliminating the blood-flow confounding effect in intravoxel incoherent motion (IVIM) using the non-negative least square analysis in liver[END_REF] shows, contrary to the proposed algorithms, a high number of spurious voxels with high PF values. This is mainly due to the fact that when the NNLS-based method [START_REF] Gambarota | Eliminating the blood-flow confounding effect in intravoxel incoherent motion (IVIM) using the non-negative least square analysis in liver[END_REF] fails in detecting the presence of blood vessels, the PF is then calculated on the basis of the bi-exponential IVIM-MRI model. In this case, the blood vessel effect will contribute as a confounding factor to the estimation of D * (d 2 ). As a result, higher values of a 2 are to be expected leading to high PF values in those voxels. It is worth noting that a similar behavior of the three considered methods was also observed for the quantification of ADC (d 1 ), D * (d 2 ) and PF (f ) in the other ROIs.

As mentioned previously, the above quantification results can be evaluated also in terms of the NMSE defined in equation (3.12). Figure 3.12 shows the NMSE of the estimated ADC (d 1 ), D * (d 2 ) and PF (f ) in the case of fully simulated data. According to this figure, the proposed methods show better quantification results for the D * and the PF parameters compared to the NNLS-based one especially for relatively low SNR value (e.g., SNR = 50). In addition, comparable results of the three methods are to be noticed for the quantification of ADC (d 1 ). The IVIM parameters contains ADC (d 1 , mm 2 /s) (top row), the pseudo-diffusion coefficient D * (d 2 , mm 2 /s) (middle row) and the perfusion fraction PF (f , bottom row). Parametric maps were obtained using the SAVTE-IVIM R , the SAVTE-IVIM E and the NNLS-based methods for SNR=50 of (a) volunteer 1 and of (b) volunteer 2.

According to Figure 3.13, the three considered methods show generally comparable quantification results of ADC (d 1 ) and D * (d 2 ) which are, to some extent, consistent with the ground truth. Regarding the quantification of PF (f ), the NNLS-based method [START_REF] Gambarota | Eliminating the blood-flow confounding effect in intravoxel incoherent motion (IVIM) using the non-negative least square analysis in liver[END_REF] shows, contrary to the proposed algorithms, a high number of spurious voxels with high PF values. A similar behavior of the three considered methods was also observed for the quantification of ADC (d 1 ), D * (d 2 ) and PF (f ) for all six volunteers for realistic simulation. This is consistent with the results of the fully simulation.

Similarly, the NMSE is also used to assess the quantification of the IVIM parameters, as depicted in Figure 3.14. According to this figure, the NMSE calculated for the three IVIM parameters ADC (d 1 ), D * (d 2 ) and PF (f ) are represented for several SNR values (50, 100 and 150). As can be observed in Figure 3.14, for all considered parameters, better NMSE values are obtained when the SNR increases from 50 to 150. Regarding the quantification of ADC (d 1 ), the proposed methods shows higher quantification quality compared to the NNLSbased method for relatively low SNR value (SNR= 50) while comparable results are observed for the three considered methods for higher SNR values and for all volunteers in our data set except for volunteer 5 where the NNLS-based method shows better quantification of the ADC (d 1 ). Regarding the quantification of D * (d 2 ), we note generally comparable results for the three methods, with a slight superiority of the SAVTE-IVIM R over the SAVTE-IVIM E and the NNLS-based one. As far as the quantification of PF (f ) is concerned, the SAVTE-IVIM R and the SAVTE-IVIM E generally outperform the NNLS-based one especially for relatively low SNR. This behavior is confirmed for all volunteers in our data set except for volunteer 5 where the NNLS-based method shows relatively a higher quantification quality of the PF parameter.

Statistical analysis

In order to evaluate the statistical significance of the obtained results in the context of blood vessel identification and the quantification of tissue diffusion and perfusion, the Wilcoxon signed-rank test is employed. This test is a non-parametric test, and can be used to determine whether two dependent samples were selected form populations having the same distribution [McDonald, 2009].

Indeed, the Wilcoxon signed-rank test is conduct on each paired of correlation series calculated from the results provided by the proposed SAVTE-IVIM R and SAVTE-IVIM E algorithms and NNLS-based [START_REF] Gambarota | Eliminating the blood-flow confounding effect in intravoxel incoherent motion (IVIM) using the non-negative least square analysis in liver[END_REF] one. In addition, the statistical analysis is done for the case of blood vessel identification and the case of IVIM parameters quantification separately.

Case of blood vessel identification

As described above, in order to evaluate the statistical significance of the blood vessel identification results using the Wilcoxon signed-rank test, a correlation series per method (e.g., the SAVTE-IVIM R , the SAVTE-IVIM E and the NNLS-based [START_REF] Gambarota | Eliminating the blood-flow confounding effect in intravoxel incoherent motion (IVIM) using the non-negative least square analysis in liver[END_REF]) is first constructed leading to three correlation series denoted here by γ R,TM , γ According to Figure 3.15, results obtained using both the SAVTE-IVIM R and the SAVTE-IVIM E algorithms are more correlated to the ground truth than the NNLS-based approaches for both fully simulated and realistic simulated data configurations. Obtained p-values (2.92 × 10 -4 and 0.0017 for the case of fully simulated data, 2.3 × 10 -3 and 1.9 × 10 -3 for realistic simulated data) confirm this fact. Furthermore, no difference between the results of the SAVTE-IVIM R and those of the SAVTE-IVIM E is to be noticed (p-value of 0.63 for fully simulated data and 0.91 for realistic simulated data). However compared to the NNLS-based approach [START_REF] Gambarota | Eliminating the blood-flow confounding effect in intravoxel incoherent motion (IVIM) using the non-negative least square analysis in liver[END_REF], blood vessel identification results obtained using both the proposed SAVTE-IVIM R and the SAVTE-IVIM E algorithms are of high statistical significance (p-values of 2.92 × 10 -4 and 1.7 × 10 -3 in the case of fully simulated data and 2.3 × 10 -3 and 1.9 × 10 -3 in the case of realistic simulated data, respectively).

Case of IVIM parameters quantification

Regarding the quantification of tissue diffusion and perfusion, specifically the parameter ADC (d 1 ), D * (d 2 ) and PF (f ) a correlation series per parameter is built leading to γ θ R,TM , γ θ E,TM , γ θ NNLS,TM with θ ∈ {ADC, D * , PF}. For example, each element of γ ADC R,TM denotes the correlation coefficient computed, for one ROI, between the true parametric map (the GT or the TM) related to the ADC and its estimation using the SAVTE-IVIM R method. Then, a Wilcoxon signed-rank test is applied to the three pairs of correlation series (γ According to Figure 3.16, the proposed SAVTE-IVIM R and the SAVTE-IVIM E algorithms show comparative or better performances in terms of the quantification of IVIM model parameters, compared to the NNLS-based algorithm [START_REF] Gambarota | Eliminating the blood-flow confounding effect in intravoxel incoherent motion (IVIM) using the non-negative least square analysis in liver[END_REF], both for the case of fully simulated data and realistic simulated data:

θ R,TM , γ θ E,TM ), (γ θ R,TM , γ θ NNLS,TM ) and (γ θ E,TM , γ θ NNLS,TM ), θ ∈ {ADC, D * , PF}. Box plots of the differences γ θ 1 = γ θ R,TM -γ θ E,TM , γ θ 2 = γ θ R,TM -γ θ NNLS,
• Case of fully simulated data (Figure 3.16(a))

-ADC: the SAVTE-IVIM R shows higher performance compared to the SAVTE-IVIM E and to the NNLS-based approaches. This fact is confirmed by the obtained p-value (4.3 × 10 -4 and 1.7 × 10 -3 , respectively). In addition, no statistical significance in the obtained results using the SAVTE-IVIM E and the NNLS-based approach is to be noted (p-value of 0.2); -PF: the proposed SAVTE-IVIM R and SAVTE-IVIM E methods show higher performance compared to the NNLS-based one as confirmed by the obtained p-values (1.0 × 10 -3 and 2.7 × 10 -3 , respectively). Furthermore, higher quantification of the PF is obtained using the SAVTE-IVIM R compared to the SAVTE-IVIM E (p-value of 8.4 × 10 -2 ); -D * : no statistical significance is to be noted in the obtained results using the proposed SAVTE-IVIM R method compared to the SAVTE-IVIM E one (p-value of 1.0) and also to the NNLS-based algorithm (p-value of 0.4). However, higher quantification of the D * ( is to be noted for the SAVTE-IVIM E method compared to the NNLS-based one. This fact is confirmed using the obtained p-value (2.8 × 10 -2 ).

• Case of realistic simulated data (Figure 3.16(b)) -ADC: all methods show comparable results. This fact is also confirmed by the obtained p-values (0.94, 0.96 and 0.96); -PF: all methods show again comparable results as also confirmed by the obtained p-values (0.99, 0.99 and 0.96); -D * : compared to the NNLS-based approach, the SAVTE-IVIM E algorithm shows higher statistical significance (p-value of 1.5 × 10 -4 ) while no statistical significance is to be noted from the SAVTE-IVIM R (p-value of 0.55).

Application to in vivo data

To explore the feasibility of the considered algorithms on in vivo data, two comparative studies are considered hereafter. First, the behavior of these algorithms is investigated, for a given volunteer (i.e., volunteer 1) in our data set, as a function of the chosen ROI. Second, these methods are evaluated using the DW-MR images of the other volunteers (e.g., volunteer 2, 3, 4, 5 and 6 ) where only one informative ROI per volunteer is considered. The ROI selection is performed in the following way: (i) regions where macroscopic blood vessels are present, as assessed by direct visual inspection of the images, are excluded. Indeed, these areas where the blood vessel confounding effect is evident are typically excluded by clinicians when evaluating tissue diffusion and perfusion; (ii) the image contrast of the DW-MR image at b = 0 s/mm 2 is manually enhanced, to reach nearly saturation levels. This step is performed in order to identify areas characterized by a limited number of voxels with higher signal intensities compared to the neighboring voxels; this is an indicative of partial volume between blood vessels and liver parenchyma [START_REF] Gambarota | Eliminating the blood-flow confounding effect in intravoxel incoherent motion (IVIM) using the non-negative least square analysis in liver[END_REF]. Regarding the first study, four ROIs of size (16 × 16) are chosen from volunteer 1 using manually enhanced image contrast. Figure 3.17(a) shows the real DW-MR image (top left) with a standard contrast of a human liver. A zoom-in of these four chosen ROIs (white squares) is shown in Figure 3.17 It should be noted that the behavior of the manual enhancement of image contrast is operatordependent. To illustrate this fact, the ROI 4 is shown in Figure 3.18 with three different contrast levels leading to three different visualization aspects. Thus, the contrast-enhancing approach can not be employed as an unbiased mean to identify blood vessels. The spatial distribution of blood vessels quantified using the SAVTE-IVIM R , the SAVTE-IVIM E and the NNLS-based methods for four different ROIs.

We stress on the fact that contrary to the NNLS-based approach wherein the provided spatial distribution of this third exponential component is a binary map, the proposed approaches result in simultaneous identification and quantification of this component. Variations of the contrast level in the spatial distribution of a 3 obtained using the SAVTE-IVIM R and the SAVTE-IVIM E algorithms confirm this fact.

Regarding the quantification of ADC (d 1 ), D * (d 2 ) and PF (f ) shown in Table 3.2, the three methods show competitive results which are consistent with the values found in previous studies [START_REF] Gambarota | Eliminating the blood-flow confounding effect in intravoxel incoherent motion (IVIM) using the non-negative least square analysis in liver[END_REF], Barbieri et al., 2016, Leporq et al., 2015]. Four different ROIs were chosen from volunteer 1 as shown in Figure 3.17

Note that since no constraints regarding the physiological values of D * are to be considered neither in the SAVTE-IVIM E algorithm nor in the NNLS-based one, some estimates of D * , for certain voxels, are expected to be higher than its largest physiological value (e.g., 0.2). In such situation, corresponding voxels are excluded from the subsequent within-ROI averaging operations, as suggested in [START_REF] Barbieri | Impact of the calculation algorithm on biexponential fitting of diffusion-weighted mri in upper abdominal organs[END_REF].

As far as the second study is concerned, for each of the five volunteers (i.e. volunteer 2, 3, 4, 5 and 6), one informative ROI is selected (Figure 3.20, first column). The indicative distribution of blood vessels in each chosen ROI is obtained by a manual contrast enhancement (Figure 3.20, second column).

The spatial distribution of blood vessels (e.g., a 3 ) in each ROI using the SAVTE-IVIM R , the SAVTE-IVIM E and the NNLS-based methods is depicted in Figure 3.20. We note that while the NNLS-based approach provides only a binary map reflecting only an identification of blood vessels, the proposed algorithms provide simultaneous identification and quantification of blood vessels. Furthermore, we note from Figure 3.20 that all methods generally succeed in identifying voxels which seem to be affected by blood vessels. However, the NNLS-based approach provides generally a spatial distribution of blood vessels that is, to some extent, less consistent with the available indicative distribution of blood vessels. For instance, compared to the latter, in the case of volunteer 5, the NNLS-based method seems to provide a higher false negative identification rate than the proposed algorithms. Regarding the quantification of ADC (d 1 ), D * (d 2 ) and PF (f ), Table 3.3 shows that the three considered methods provide generally, for the five volunteers, quantification results that are consistent with those found in previous studies on human liver [START_REF] Gambarota | Eliminating the blood-flow confounding effect in intravoxel incoherent motion (IVIM) using the non-negative least square analysis in liver[END_REF], Barbieri et al., 2016, Leporq et al., 2015].

However, Table 3.3 shows some relatively higher PF values using the NNLS-based method compared to the ones obtained using the SAVTE-IVIM E and the SAVTE-IVIM R algorithms. This is probably due, as discussed previously, to the blood vessel confounding effect that systematically appears when a bi-exponential IVIM model fitting is to be considered for each voxel where a false negative identification of blood vessel has occurred. Five ROIs were chosen from volunteer 2 -6 as shown in Figure 3.20(a)

Discussion and Conclusion

In this Chapter, two algorithms, the SAVTE-IVIM R and the SAVTE-IVIM E , were proposed to deal with the problem of the confounding blood vessel effect in the IVIM-MRI model in liver. These algorithms were evaluated and compared to the NNLS-based one recently proposed to deal with such problem. Our comparative study was conducted using simulated and real DW-MR images of liver acquired from six healthy volunteers. The results on simulated data showed the superiority of the proposed methods over the NNLS-based one in terms of blood vessel identification. Indeed, the NNLS-based method showed higher false negative rate in terms of identifying affected voxels, especially for low SNR, as confirmed by the VLE measure and by the statistical analysis. In terms of quantifying ADC (d 1 ) and D * (d 2 ), all methods showed competitive results. However, the proposed algorithms showed higher performance in terms of quantifying PF (f ). In fact, contrary to the NNLS-based method, the SAVTE-IVIM R and the SAVTE-IVIM E methods provided PF maps that are more consistent with the ground truth. Furthermore, a high rate of spurious voxels with high PF values were observed especially for NNLS-based method. This is mainly due to the false negative identification of affected voxels. More precisely, for a false negative detection of affected voxels, a bi-exponential fitting instead of a tri-exponential one is to be systematically considered in the NNLS-based method. Indeed, with the bi-exponential fitting being employed, the confounding blood vessel effect will highly contribute to the estimation of the tissue perfusion. Consequently, high PF values are to be expected. Regarding the experiments on real DW-MR images, all methods gave, to some extent, comparable spatial distribution of blood vessels in the ROI. However the NNLS-based one showed higher PF quantification values compared to the ones obtained using the proposed methods. This phenomena is, as discussed above, mainly due to the problem of false negative identification of affected voxels.

It is noteworthy that as the spatial distribution of blood vessel in liver is not uniform, there exists some ROIs for which the sparsity assumption of blood vessels on which the proposed approach relies, is fulfilled. Besides, it is worth noting that the sparsity rate is highly dependent on the size of the chosen ROI. For example, the sparsity assumption can easily be violated for small ROIs ( i.e., 2 × 2). But, as the typical size of an ROI delineated by clinical expert is, in general, not less than (10 × 10), one can always find an ROI for which the sparsity assumption of the spatial distribution of blood vessels holds true. In such case the defined ROI is well-called informative.

In conclusion, a new approach SAVTE-IVIM to cope with the confounding blood vessel effect in the IVIM-MRI model in liver was proposed. This is accomplished by resorting to an automatic identification of the potential presence of large blood vessels in the ROI, especially in situations where visual identification of such blood vessels is not evident. In addition, two different strategies have been proposed to take into account the non-negativity constraints in the SAVTE-IVIM model: (i) a rough strategy, where potential negative values of parameter estimates were dealt with by taking into account prior information regarding the physiological ranges of the model parameters; and (ii) an embedded strategy where a change of variable into square was employed during the optimization process. In addition to the automatic blood-vessel identification, the SAVTE-IVIM approach, in its two proposed variants, provides a quantification of (i) the confounding blood vessel effect of each voxel; (ii) the apparent diffusion coefficient and (iii) the tissue perfusion. Our results on both realistic and real DW-MR images of six healthy volunteers have shown the efficiency of the proposed algorithm. 65 After 2014, there has been an increasing interest in the bone marrow IVIM:

In 2015, Ohno et al. followed up immediately, and performed IVIM measurement on 11 patients [START_REF] Ohno | Evaluation of perfusion-related and true diffusion in vertebral bone marrow: a preliminary study[END_REF].

Between 2015 and 2017, five new publications have applied the IVIM technique to do the bone marrow assessment, see [START_REF] Park | Intravoxel incoherent motion diffusion-weighted magnetic resonance imaging of focal vertebral bone marrow lesions: initial experience of the differentiation of nodular hyperplastic hematopoietic bone marrow from malignant lesions[END_REF], Lee et al., 2017a, Baik et al., 2017, Niu et al., 2017, Yoon et al., 2017].

Between 2017 and 2020, twelve new publications appeared [START_REF] Karampinos | Quantitative mri and spectroscopy of bone marrow[END_REF], Lasbleiz et al., 2019, Wu et al., 2019, Switlyk, 2019, Tan et al., 2019, Li et al., 2020a, Chen et al., 2019[START_REF] Li | Correlation of Intravoxel Incoherent Motion Parameters and Histological Characteristics From Infiltrated Marrow in Patients With Acute Leukemia[END_REF], Li and Schwartz, 2020, Fathi Kazerooni et al., 2020, Fan et al., 2020, Minutoli et al., 2020].

To date, there have been approximately 20 publications that have conducted relevant research on bone marrow IVIM. In general, the results of these studies indicate that IVIM is a promising technique to observe changes at tissular and cellular within the bone marrow.

As introduced in section 2.3, bone marrow comprises two marrow components: red marrow, which is composed of hematopoietic cells (approximate to 60%) and adipocytes (fatty tissue, approximate to 40%), and yellow marrow, which consists mainly of fatty tissue (approximate to 80%) [START_REF] Ricci | Normal age-related patterns of cellular and fatty bone marrow distribution in the axial skeleton: Mr imaging study[END_REF]. The large amount of lipids1 represents a confounding factor in the measurement of water diffusion [Mulkern andSchwartz, 2003, Biffar et al., 2010] and tissue perfusion in bone marrow [START_REF] Biffar | Quantitative analysis of vertebral bone marrow perfusion using dynamic contrast-enhanced mri: Initial results in osteoporotic patients with acute vertebral fracture[END_REF].

In a recent study, Lasbleiz et al. have assessed the effects of different fat suppression techniques on the measurement of the diffusion and perfusion IVIM parameters [START_REF] Lasbleiz | Measurements of diffusion and perfusion in vertebral bone marrow using intravoxel incoherent motion (ivim) with multishot, readout-segmented (resolve) echo-planar imaging[END_REF]. In this work, the authors have also developed an optimized IVIM protocol in bone marrow using the RESOLVE (readout segmentation of long variable echo train) EPI sequence, which was introduced in 2009 [START_REF] Porter | High resolution diffusion-weighted imaging using readout-segmented echo-planar imaging, parallel imaging and a two-dimensional navigator-based reacquisition[END_REF]. Compared with singleshot EPI, which was employed in all previous studies on bone marrow, the RESOLVE provides MR images with higher spatial resolution and with reduced distortions at the cost, however, of a longer acquisition time [START_REF] Lasbleiz | Measurements of diffusion and perfusion in vertebral bone marrow using intravoxel incoherent motion (ivim) with multishot, readout-segmented (resolve) echo-planar imaging[END_REF]. The results of the study conducted by Lasbleiz et al. indicate that the RESOLVE-EPI sequence, with the SPAIR fat suppression technique, is the current optimal choice for bone marrow IVIM.

It should be pointed out that the confounding factor of the lipid signal could partly explain the relatively broad range of values for the IVIM parameters reported in previous studies [START_REF] Marchand | Mri quantification of diffusion and perfusion in bone marrow by intravoxel incoherent motion (ivim) and non-negative least square (nnls) analysis[END_REF], Park et al., 2017, Lee et al., 2017a, Baik et al., 2017, Niu et al., 2017, Yoon et al., 2017, Lasbleiz et al., 2019]. Other factors responsible for this broad range of values include: different bone marrow locations (spine, pelvis, etc.), different imaging protocols, and different health conditions of the subjects investigated.

Despite the progress achieved over the past years, to date IVIM in bone marrow is still far from being a "plug-and-play" method. In section 2.1.3.2 the IVIM parametric maps were introduced to play a vital role in the medical evaluation and diagnosis process. However, out of the 20 bone-marrow IVIM publications so far, only 6 of them have provided the parametric maps; moreover the quality of these maps is debatable. As a matter of fact, little attention has been paid to the data post-processing of bone-marrow IVIM, i.e. the quantification of the 4.1. Introduction IVIM diffusion and perfusion parameters. As introduced previously in section 2.1.3.1, the most commonly used algorithms are generally based on the idea of least squares (LSQ).

In the current study, we implemented five different fitting algorithms for the IVIM quantification on bone marrow. Four out of the five are based on the LSQ, called respectively One-Step, Two-Step, Three-Step, and Fixed-D * algorithms. The fifth one is based on Bayesian inference [START_REF] Neil | On the use of bayesian probability theory for analysis of exponential decay date: An example taken from intravoxel incoherent motion experiments[END_REF], Bretthorst et al., 2005, Dyvorne et al., 2013b], named Bayesian-based algorithm.

Overall, the aim of this study was to evaluate and compare five different algorithms for IVIM quantification of human bone marrow, with the goal of generating parametric maps of diffusion and perfusion.

Fitting algorithms

In this section, the five fitting algorithms for the IVIM parameters quantification are introduced. The methods belong to two classes: deterministic and probabilistic. For the deterministic approaches, such as the One-Step, Two-Step, Three-Step and the Fixed-D * algorithms, they solve the parameters quantification problem in a least squares (LSQ) sense. And for the probabilistic one, it is based on Bayesian inference [START_REF] Neil | On the use of bayesian probability theory for analysis of exponential decay date: An example taken from intravoxel incoherent motion experiments[END_REF] to do the parameters estimation.

To better understand the algorithms, let us first review the bi-exponential IVIM model for a single voxel. For the target voxel, the recorded intensity signal for different b values S b , can be expressed following the bi-exponential model:

S b = S 0 f e -bD * + (1 -f )e -bD (4.1)
where S 0 stands for the signal intensity taken at b = 0 s/mm2 , f denotes perfusion fraction PF and D, D * 2 stand, respectively, for the diffusion and perfusion coefficients.

Algorithms based on LSQ

An important source of least squares problems is data fitting. More precisely, the goal is to find the optimal parameters of a predefined data model such that a lowest model error is obtained.

In the case of the IVIM identification, least squares problem can be generally formulated as follows:

arg min

Θ b S b -fb (Θ) 2 (4.2)
where fb (Θ) is the predefined model of the signal intensity with its set of parameters Θ (i.e.. Θ = {D, f, D * }). According to equation (4.1), the above optimization problem (equation (4.2)) is clearly non-linear as the objective function to be minimized is non-linear in any of the model parameters (i.e., D, f and D * ). To solve such non-linear optimization problem, adequate optimization algorithms such as Gauss-Newton, Levenberg-Marquardt (LM) or Trust-Region (TR) algorithms [START_REF] Nocedal | Numerical Optimization[END_REF]] can be used.

In our study, the TR algorithm is adopted.

Step-wise strategies

For IVIM parameters (including D, f and D * ), a simultaneous nonlinear fitting can easily produce some uncertainties, especially in the case of low perfusion fraction, limited quality of images [START_REF] Istratov | Exponential analysis in physical phenomena[END_REF]. Thus, in the IVIM practices, a more commonly used fitting approach is in a "segmented" strategy in which the parameters were analyzed separately [START_REF] Patel | Diagnosis of cirrhosis with intravoxel incoherent motion diffusion mri and dynamic contrast-enhanced mri alone and in combination: preliminary experience[END_REF], evaluating the diffusion coefficient D individually using the data of high b values.

Here, we summarize the two common strategies of the "segmented" algorithms of data fitting, together with the algorithm of simultaneous parameters quantification, as "step-wise" algorithms, which are One-Step, Two-Step and Three-Step algorithms.

One-Step strategy

This algorithm solves the IVIM identification problem in a direct manner [Thoeny and De Keyzer, 2011, Wittsack et al., 2010, Mazaheri et al., 2012] by simply fitting the bi-exponential model (equation (4.1)) to the acquired signal intensities of each voxel in the Region-Of-Interest (ROI).

It is worth noting that, the quantification of IVIM parameters is generally [START_REF] Patel | Diagnosis of cirrhosis with intravoxel incoherent motion diffusion mri and dynamic contrast-enhanced mri alone and in combination: preliminary experience[END_REF] on the average value of the signal in the target ROI, which to some extent improves the SNR.

The optimization problem defined in equation (4.2) with three variables (D, f and D * ) can be written for as following:

arg min f,D,D * b S b -S 0 (f e -bD * + (1 -f )e -bD ) 2 (4.3)

Two-Step strategy

With respect to the Two-Step algorithm, the IVIM quantification problem is solved in two successive steps.

In the first step, the value of the apparent diffusion coefficient D is first estimated using a mono-exponential fitting of data points corresponding to b values which are greater than 200 s/mm 2 and the estimate of D is denoted as D. This is done by solving the following non-linear optimization problem:

arg min D b≥200
(S b -S 0 e -bD ) 2 (4.4)

In the second step, with the D parameter being fixed to its estimated value D, parameters f and D * are then estimated by solving the following non-linear optimization problem:

arg min f,D * b S b -S 0 (f e -bD * + (1 -f )e -b D) 2 (4.5)
where solving this double-variable (f and D * ) optimization problem is performed over all the b values used to acquire the IVIM data.

An example showing the scheme of "segmented" analysis is depicted in Figure 4.1.

Three-Step strategy

For Three-Step algorithm, three steps are used to solve the IVIM parameters quantification problem:

Step 1. Similarly to the Two-Step algorithm, an estimate of the D parameter, denoted as D, is obtained using a mono-exponential fitting (see equation (4.4)) of data points corresponding to b values greater than 200 s/mm 2 .

Step 2. With the D value being fixed to its estimate D in Step 1, the signal interception, denoted here by S intercept , is computed with a mono-exponential model (i.e., equation (4.4)) using the entire set of b values (see Figure 4.1 the green curve). Then, an estimate of the perfusion fraction, denoted as f , is calculated directly as follow [START_REF] Pekar | On the precision of diffusion/perfusion imaging by gradient sensitization[END_REF]:

f = 1 - S intercept S 0 (4.6)
Step The first step of the Two-Step or Three-Step algorithm described above relies essentially on the fact that the perfusion coefficient (D * ) is much larger than the diffusion one (D). Thus, the contribution of the exponential term related to D * in the bi-exponential IVIM model (equation (4.1)) can be omitted for high b values (i.e., b ≥ 200 , s/mm 2 ) compared to the contribution of the one related to the D coefficient. Therefore, a mono-exponential fitting of those data points corresponding to high b-vales can result in a reliable estimation of the D parameter [START_REF] Pekar | On the precision of diffusion/perfusion imaging by gradient sensitization[END_REF], Chandarana et al., 2011]. Once this parameter is fixed, estimating the other parameters (f and D * ) can be done easily using a non-linear optimization scheme.

The Fixed-D * algorithm

Instead of fixing the value of the D parameter with a mono-exponential fitting, this algorithm relies on specifying the parameter D * [START_REF] Wittsack | Statistical evaluation of diffusion-weighted imaging of the human kidney[END_REF], Lemke et al., 2010, Heusch et al., 2013]. In this study, the D * parameter was fixed to (15 × 10 -3 mm 2 /s) reliably obtained from previous studies [START_REF] Lasbleiz | Measurements of diffusion and perfusion in vertebral bone marrow using intravoxel incoherent motion (ivim) with multishot, readout-segmented (resolve) echo-planar imaging[END_REF], Bourillon et al., 2015]. Then, estimating the other parameters in the bi-exponential IVIM model is performed by solving the following double-variable (D and f ) non-linear optimization problem:

arg min f,D b S b -S 0 (f e -bD * F + (1 -f )e -bD ) 2 (4.8)

Bayesian-based algorithm

Besides the very commonly used non-linear least squares algorithms, there has been an increasing interest in Bayesian approaches for the IVIM identification problems [START_REF] Schmid | A bayesian hierarchical model for the analysis of a longitudinal dynamic contrastenhanced mri oncology study[END_REF], Orton et al., 2014, While, 2017, Gustafsson et al., 2018]. The key point of Bayesian framework is that each parameter of interest is described through its probability density function. Parameters estimation in Bayesian framework is solved by maximizing the joint posterior probability of each parameter of interest. Regarding the case of the bi-exponential IVIM, the parameter estimation problem, for a single voxel, can be written as follows: The uncertainty (the model error) of the IVIM model can be assumed to follow a Gaussian distribution with a standard deviation characterized using Jeffreys' prior [Jeffreys, 1946, Jeffreys, 1998]. Thus we can get

P ( S|θ, I) ∝ (Q/2) -M/2 (4.12)
where the model error is denoted as

Q = b (S b -S 0 (f e -bD * + (1 -f )e -bD )
) 2 and M is the number of b values.

In addition, the prior probability P (ν|I), l ν ≤ ν ∈ θ ≤ h ν , follows a Gaussian distribution N (µ ν , σ ν ) of a mean, µ ν , and a standard deviation, σ ν where l ν and h ν are denoting the lower and higher bounds of the parameter ν ∈ θ.

A Monte Carlo Markov Chain (MCMC) approach [Neil andBretthorst, 1993, Neal, 1993] is then employed to solve the maximization problem defined in equation (4.11) where the Metropolis-Hastings (MH) sampling method [Yildirim, 2012] is used to iteratively generate the Markov Chain. At each iteration of this procedure, the mean µ ν , ∀ν ∈ θ, is updated while the standard deviation, σ ν , is kept fixed to a value equal to (h ν -l ν )/3. An initial guess of µ ν , is set according to table 4.1 obtained from previous studies [START_REF] Lasbleiz | Measurements of diffusion and perfusion in vertebral bone marrow using intravoxel incoherent motion (ivim) with multishot, readout-segmented (resolve) echo-planar imaging[END_REF], Bourillon et al., 2015].

For each implementation process, the maximum iteration number was set to 1000, and Monte Carlo repetition was set to 50.

ν f (%) D(×10 -3 mm 2 /s) D * (×10 -3 mm 2 /s)
Initial guess of µ ν 10 0.5 15

Lower bound (l ν ) 1 0.001 5

Higher bound (h ν ) 40 2 50

The set of parameters characterizing the prior distribution P (ν|I) = N (µν , σν ), ν ∈ θ{f, D, D * } with σν = (hν -lν )/3

DW-MRI and simulation settings

All experiments were conducted in accordance with procedures approved by the local Institutional Review Board. Written informed consent was signed by each volunteer before the measurements. In the current study, all MRI data were measured using a 1.5T MAGNETOM Aera system (Siemens Healthcare, Erlangen, Germany). The IVIM measurements were carried out on the lumbar spine of six healthy volunteers (age range 18-29 years, mean age 26.1 ± 4.0 years, three women and three men).

DW-MRI: Sequence parameters

In this study, the IVIM measurement was performed in sagittal orientation using the RESOLVE sequence. During the data acquisition, the spectral attenuated inversion recovery (SPAIR) technique provided by the system vendor was used to suppress the lipid signal.

In Table 4.2, the parameters of the sequence for IVIM DW-MRI are listed. According to the fitting curves in Figure 4.3 and the N res (defined as the norm of the residual between the fitted and input signal intensities) given in Table 4.3, in terms of the averaged data of an ROI (here is the ROI of L1, which contains approximately 100 voxels), the five algorithms which are One-Step, Two-Step, Three-Step, Fixed-D * and Bayesian-based approaches have comparable performances on the IVIM parameters quantification. Part of the reason is that the average operation on an ROI can increase the SNR to some extent (theoretically times √ N with N denoting the number of voxels in the ROI). In other words, an ROI can be regarded as a magnified version of one voxel.

Simulations and Evaluation Criteria

Two sets of experiments were designed. One set is dedicated for the numerical analysis, to evaluate and compare different algorithms in terms of the IVIM parameters quantification. The other set is designed to assess the IVIM parametric maps provided by different algorithms.

-Numerical analysis

According to the most recent bone marrow IVIM study on healthy subjects [START_REF] Lasbleiz | Measurements of diffusion and perfusion in vertebral bone marrow using intravoxel incoherent motion (ivim) with multishot, readout-segmented (resolve) echo-planar imaging[END_REF], a set of IVIM parameters is chosen to be the ground truth, which is 48 × 10 -5 mm 2 /s for D, 13% for PF and 18 × 10 -3 mm 2 /s for D * . Together with the IVIM parameters, a set of b values (0, 50, 100, 150, 400, 800, 1000 s/mm 2 ), which is the same as the set of b values for in vivo data acquisition, is used to generate the noise-free DW-MR signals. Then the Rician noise was add to the data to generate noisy data of different SNR (10, 20, 50 and 100) through the adjustment of the noise variance. This kind of designed DW-MR data was generated for each SNR value (out of 10, 20, 50 and 100), for 20 times, to simulate 20 Monte Carlo (MC) trials.

-Parametric maps

The quantification of perfusion-related parameters, especially the pseudo-diffusion D * , is well known of weak reproducibility/repeatability and of large variance. To investigate the sensitivity of the IVIM parametric maps provided by different algorithms to the perfusion-related parameters, a set of simulated data with changing D * and changing PF was designed. This data set simulates 16 × 36 = 576 voxels, distributed in 4 × 4 = 16 square blocks, to simulated 16 square ROIs. Each block consists of 6 × 6 = 36 voxels with the same ground truth of the IVIM parameters, in which way 36 MC trials were simulated.

Relative Error (RE)

To numerically assess the performances of different algorithms in terms of IVIM parameters quantification, the quantification results given by each algorithm can be used to calculate a relative error (RE) compared with the ground truth. The RE is defined as follows:

RE = estimation -groundtruth groundtruth (4.13)
For each algorithm, each SNR value and each MC trial, a RE value was calculated. The results are displayed in the form of box plots in Figure 4.5.

Results

Results

Numerical analysis

For the first group of simulation, the ground truth for the IVIM parameters are 48 × 10 -5 mm 2 /s for D, 13% for PF and 18 × 10 -3 mm 2 /s for D * . In The advantage of using box plots is that they can more intuitively describe the accuracy and 79

Parametric maps

In Figure 4.7, the parametric maps originating from the second group of simulation, provided by different algorithms, are shown for the parameters (a) D, (b) PF and (c) D * .

With respect to the ground truth values settings for all the 576 simulated voxels, the parameter D is keep fixed to one value 0.48 × 10 -3 mm 2 /s, PF has four ground truth values 5%, 10%, 15%, and 20%, while D * has also four ground truth values, {6, 12, 18, 24} × 10 -3 mm 2 /s. For each block, the ground truth is the same for all voxel in it. For this data set, the SNR is set to 20. It is worth noting that, the acquired in vivo data for spine DW-MRI, the SNR is close to 20.

According to Figure 4.7, the five algorithms (including the four LSQ based algorithms and Bayesian-based algorithm) provide the comparable parametric maps for the IVIM parameter D, and also comparable with the ground truth. With respect to the perfusion-related parameters PF and D * , non-negligible bias of the parametric maps of all five methods can be observed compared to the ground truth maps. But in comparison with the four LSQ based algorithms, the Bayesian-based approach provide parametric maps relatively closer to the ground truth, both for PF and D * . Furthermore, the results for the Bayesian-based algorithm are of smaller variance which can be derived from its smoother maps.

In addition, the five algorithms were also conducted on in vivo DW-MR data acquired from six healthy volunteers. Take a 24 years old male volunteer as an example, the parametric maps are given in Figure 4.8. For better readability, parts of the images that do not need to be focused on were eliminated, and the left part include the five lumbar vertebrae from L1 to L5. The ROI for each lumbar vertebra from L1 to L5 are marked in color red in

Results

Chapter 4

The numerical values of the IVIM parameters corresponding to the parametric maps in Figure 4.8 are given following in Table 4.4. With respect to the perfusion-related parameters PF and D * , the estimates using Bayesian-based algorithm have smaller variance compared to the estimates using other LSQ based algorithms. In general, the quantification of the IVIM parameters are acceptable for all algorithms. The mean values and standard deviations of D, PF, and D * were calculated over the ROI in L1 for one volunteer.

The full results for the six volunteers and for the five lumbar vertebrae L1 to L5 can be found in Appendix E, Table E.1 and Table E.2.

Discussion and conclusion

In 2014, Marchand et al. started the era of bone marrow IVIM by realizing the quantification of the IVIM parameters, including diffusion and perfusion, with a bi-exponential model for the first time [START_REF] Marchand | Mri quantification of diffusion and perfusion in bone marrow by intravoxel incoherent motion (ivim) and non-negative least square (nnls) analysis[END_REF]. Since then, there has been an increasing interest in the bone marrow IVIM technology. Recently, more and more studies applying the IVIM techniques to the evaluation on bone marrow. This is partly due to the improvements of the sequence and hardware. However, the results of the published IVIM studies on bone marrow covers a relatively broad range of values of the IVIM parameters and few work was dedicated to the fitting algorithms used for the IVIM parameters quantification.

Thus, we focused on the quantification algorithms for the post data processing of the bone marrow IVIM. In this study, we implemented five different fitting algorithms, four out of which are based on the LSQ -One-Step, Two-Step, Three-Step, and Fixed-D * algorithms and fifth of which is based on Bayesian inference [START_REF] Neil | On the use of bayesian probability theory for analysis of exponential decay date: An example taken from intravoxel incoherent motion experiments[END_REF], Bretthorst et al., 2005, Dyvorne et al., 2013b] -Bayesian-based algorithm. These algorithms were conducted on both simulated data and in vivo data, based on voxel level. All algorithms show acceptable results and the Bayesian-based algorithm gives the best quantification results and parametric maps. For the simulations, Bayesian-based algorithm shows higher accuracy and precision. For the in vivo data, it gives the quantification of perfusion-related parameters of smaller variance.

It would be better if our work can inspire and bring some help to subsequent researchers. Based on the above work, here are some suggestions:

From the accuracy and precision point of view, the Bayesian-based algorithm is suggested; from the perspective of real-time calculation, the Two-Step and Three-Step algorithms are both good choices. Of course, if there is reliable prior information about the value of D * , the Fixed-D * algorithm is preferable.

Difference between MRS and MRI

MRI is primarily concerned with the generation of anatomic images. MRI scans usually contain a large number of voxels whose signals reflect the bulk magnetic properties (e.g. T1, T2, susceptibility, flow) of the tissues they contain. The signals of interest recorded by MRI come predominantly from protons in water and fat.

The goal of MRS is to obtain metabolic information from molecules thanks to the chemical shift of the resonance frequencies of the nuclear spins in the sample. By excluding the overwhelming signals from water and fat, MRS can provide detailed and quantitative information of molecules, e.g. it can detect small metabolites existing in millimolar (mM) concentrations. These metabolites can be differentiated because they resonate at slightly different frequencies based on their local chemical environments. The degree of frequency separation between two molecular species is characterized by their chemical shift (which will be introduced subsequent in section 5.1.2), a small number displayed on the horizontal axis below their spectra. Since the relative areas under each peak are proportional to the number of spins involved, one can determine composition quantitatively with peak integrals.

In 1980, Ackerman et al. obtained surface-resolved NMR spectra for the first time using surface coils [START_REF] Ackerman | Mapping of metabolites in whole animals by 31 p nmr using surface coils[END_REF]. In 1982, Brown et al. proposed a set of three-dimensional NMR chemical shift imaging (CSI) theory [START_REF] Brown | Nmr chemical shift imaging in three dimensions[END_REF], which provided the theoretical basis for further study on CSI. Since then, the development of related theories and experimental instruments has promoted the further development of body core magnetic resonance spectroscopy (MRS) imaging technology with good positioning accuracy.

The NMR spectroscopy techniques include single-voxel spectroscopy (SVS) and multi-voxel spectroscopy (MVS). SVS uses frequency-selective pulses and it can only choose to excite slices with a certain thickness. The MVS usually uses magnetic field gradients to phase-encode the signals within a delay of free precession, thus the obtained NMR signals contain the spatial information.

Chemical shift imaging (CSI) is such a kind of technology that provides spatial information for MRS imaging.

Chemical shift imaging (CSI)

To understand what is CSI, we have to introduce the chemical shift first. Here let's revise the notion of Larmor frequency:

ω = γB ef f (5.1)
where the B ef f denotes the effective magnetic field strength. Because in the region of interest (ROI), the magnetic field generated by neighboring nuclei and electrons will affect the effective magnetic field strength around the nuclei of interest, resulting in a local magnetic field different from the external main magnetic field. The motion electrons around the nucleus are affected by the external magnetic field B o , will generate an induced current and an induced magnetic field on the nucleus. The superposition of the induced magnetic field and the external magnetic field causes the effective magnetic field strength B ef f around the nucleus to change. The shielding factor σ is used to measure the shielding effect between the nucleus and the external magnetic field:

B ef f = B 0 -σB 0 = (1 -σ)B 0 (5.2)
Thus, placing a specific nucleus with a gyromagnetic ratio γ in a magnetic field of known strength, different resonance frequency signals on the NMR spectrum will be obtained. This phenomenon is called chemical shift because the frequency signal on the spectrum looks like to shift from where it was expected to occur. With this property, different compounds can be distinguished from the NMR spectrum. Obviously, this property is also the basis of CSI technology [Brateman, 1986].

Pulse sequences for CSI

The scheme of 1-D CSI is shown in Figure 5.1. A spin-echo sequence, with interpulse delays T E/2 and T E/2, was employed. T E is the echo time, G slice is the slice gradient and G y is the phase encoding gradient. In this particular sequence, eight phase encoding steps are indicated. In this scheme, multiple 180°pulses can be employed (e.g. two 180°pulses to make it double spin-echo sequence). At the end of the sequence, the CSI signal (FID signal) is acquired. The CSI signals are of different frequencies, corresponding to the chemical shifts. Thus, after the Fourier transform of the CSI signal, the spectrum can show the information of the chemical structure of the sample. the metabolite concentration estimates. The conversion of a peak intensity into a meaningful quantitative metabolite concentration is a critical step that requires the use of a reference peak: a spectral peak from a compound whose concentration is known. A metabolite concentration estimate, then, is obtained by comparing the intensity of the desired metabolite peak with the intensity of a reference peak. That's why we need to assume the concentration of compound A to be known. By comparing the intensity of peak B with the intensity of peak A, the concentration of B can be obtained.

To do the quantification of the peak intensity, a number of methods can be employed. Apparently, peak integration can be used to directly determining the peak intensities. Another very common method is peak fitting.

Peak fitting of MRS

Peak fitting needs to choose a model function that well describes the shape of the peak of interest. An example of the peak fiting method is shown in Figure 5.5. The Lorentzian lineshape is used here to be the model function. After curve fitting, the residual which is defined as the difference between the data and the fitted curve, is employed to assess the goodness of the fit. For the choice of the model functions, the common alternatives are Lorentzian, Gaussian, and Voigt lineshapes. Among the three lineshapes, the Lorentzian one, which is also employed in the work of section 3, is the most commonly used fitting model function in MR spectral quantification. An MRS FID signal can be approximated as a decaying mono-exponential signal in the timedomain. A Lorentzian lineshape is simply a representation of such a decaying in frequency domain. Thus, if we are going to measure the MR signal with a perfect mono-exponential T 2 decay, the Fourier transform of this signal would be a Lorentzian lineshape. For this reason, the Lorentzian function is a natural choice to fit the MRS peaks. In the MRS literature, the term "Lorentzian fitting" generally refers to the process of fitting the data to a Lorentzian curve in the frequency domain, or to a decaying exponential in the time-domain. Mathematically, the Lorentzian function, represented as a decaying exponential in the time-domain can be written as:

s(t) = A • e -i(ω 0 (t+t 0 )+φ) • e -(t+t 0 )/T 2 (5.4)
where A is the signal amplitude, T 2 is the transverse relaxation time, ω 0 is the frequency offset in radians, t 0 is the temporal shift between t = 0 and the top of the echo, and φ is the global phase offset in radians.

Fitting the model to the data can be achieved by adjusting the parameters (A, T 2 , ω 0 , t 0 , and phi) until the fitted curve matches the data curve as closely as possible, or more accurately, until the residual is minimized as follows:

minimize|S -s| 2

(5.5)

where S is the data and s is the model function.

Peak fitting method can also be improved by imposing some restrictions, known as prior knowledge, on some of the model parameters. In the 1990s, the accurate spectral curve fitting was significantly improved using advanced prior knowledge and in 1997, an algorithm called the Advanced Method for Accurate, Robust and Efficient Spectral fitting (AMARES) was proposed [START_REF] Vanhamme | Improved method for accurate and efficient quantification of mrs data with use of prior knowledge[END_REF]. The AMARES method can incorporate a great range of prior knowledge constraints, e.g. choosing a Lorentzian as well as a Gaussian lineshape, imposing upper and lower bounds on the parameters [START_REF] Vanhamme | Improved method for accurate and efficient quantification of mrs data with use of prior knowledge[END_REF]. There are a number of user-friendly software packages available providing state-of-the-art spectral fitting algorithms in an easy-to-use graphical user interface. jMRUI is one of them and it is freely available. jMRUI includes a batch of tools for both spectral peak fitting and spectral simulation, and it supports for both single-voxel MRS and MRS imaging data analysis [START_REF] Naressi | Java-based graphical user interface for the mrui quantitation package[END_REF], Naressi et al., 2001a, Stefan et al., 2009]. Peak fitting in jMRUI is performed using AMARES [START_REF] Vanhamme | Improved method for accurate and efficient quantification of mrs data with use of prior knowledge[END_REF] method. The application of jMRUI with the AMARES approach is highly cited in literature [START_REF] Torriani | Intramyocellular lipid quantification: repeatability with 1h mr spectroscopy[END_REF], Stagg et al., 2009, Near et al., 2011], and the jMRUI is one of the most commonly used peak fitting software packages for in-vivo MRS. Figure 5.6 shows an example of peak fitting in jMRUI, using the AMARES algorithm.

Virtual Phantom (ViP)

The optimization of the CSI methodology, in terms of both data acquisition and data analysis, has been and still is a very active field for research [Lee et al., 2017b, Mlynárik et al., 2008]. With respect to the data analysis, for the quantification of the spectroscopy signal that the CSI provides, a reference signal is necessary. Two approaches are typically employed: the use of an internal reference signal (in brain, for instance, water or creatine) or an external reference signal, which originates from a physical phantom with a well known metabolite concentration.

An alternative approach for providing a reference signal in MR spectroscopy was proposed, for applications to high resolution NMR, by [START_REF] Barantin | A new method for absolute quantitation mrs metabolites[END_REF]. This approach, also referred to as the "Electronic REference To access In vivo Concentrations" (ERETIC), consists in generating a reference signal using radiofrequency (RF) electronics and transmitting it to the receiver coil of the NMR spectrometer during the data acquisition. Many studies have employed the ERETIC technique, with applications to in vitro NMR, for instance 2D-NMR, HR-MAS NMR and solid state NMR [START_REF] Martínez-Bisbal | Determination of metabolite concentrations in human brain tumour biopsy samples using hr-mas and ERETIC measurements[END_REF], Michel and Akoka, 2004, Ziarelli et al., 2007]. However, there have been only a very limited number of studies where the ERETIC technique was implemented on MR systems dedicated to in vivo applications. In other words, only few applications of ERETIC have been employed in techniques that are used in vivo [START_REF] Albers | Evaluation of the ERETIC method as an improved quantitative reference for 1h hr-mas spectroscopy of prostate tissue[END_REF], Franconi et al., 2002, Heinzer-Schweizer et al., 2010, Lee et al., 2010], such as single voxel spectroscopy, CSI and MRI. One technique that was developed to extend the ERETIC concept to in vivo application is the Virtual Phantom (ViP) for MRI [START_REF] Rondeau-Mouro | µ-vip: Customized virtual phantom for quantitative magnetic resonance micro-imaging at high magnetic field[END_REF], Saint-Jalmes et al., 2014, Salvati et al., 2016, Saint-Jalmes et al., 2017].

In Figure 5.7, the a scheme of how to apply the ViP technique, including the design and transmission of the ViP signal together with the data acquisition in an MR experiment is shown [START_REF] Saint-Jalmes | ViP mri: virtual phantom magnetic resonance imaging[END_REF].

It could be advantageous to further extend the ViP technique to CSI. As a matter of fact, in CSI : 1) internal reference signals might not be valid, as they might be affected by the pathology and 2) the acquisition of the reference signal from an external phantom taped to the body might not be practical and would require a change of the FOV parameters, with a consequent loss of the CSI spatial resolution. Thus, it is of interest to consider the generation of a virtual phantom for signal quantification in CSI.

To apply the ViP technique to CSI, we first provide a theoretical analysis of the 2D CSI signal; then, with the designed ViP signal, we show the feasibility of adjusting and fine-tuning the amplitude and frequency of the ViP reference peak. Finally, the evaluation of the spatial uniformity and time stability of the ViP signal is performed.

Overall, the purpose of the current study was to explore the potential of ViP to generate reference signals in 2D CSI and to implement ViP 2-D CSI on an MR system dedicated to preclinical research. The steps of ViP MRI: First, the phantom shape of the virtual phantom is designed and the K-space representation is also generated based on it. Second, the simulated K-space lines is converted into an RF signal using a waveform generator. Third, the signal representing K-space line is transmitted by a dedicated RF coil positioned in the scanner bore. In the end, the ViP RF signal is synchronized with the MR scanner data acquisition by the RF unblank signal from the MR console [START_REF] Saint-Jalmes | ViP mri: virtual phantom magnetic resonance imaging[END_REF].

Theory

The MR signal R MR (t) originating from a physical phantom can be written as: R MR (t) =

x,y S MR (x, y, t)dxdy =

x,y ν ŜMR (x, y, ν)e i2πνt dν dxdy (5.6)

where S MR (x, y, t) is the MR signal of the element volume and ŜMR (x, y, ν) = S MR (x, y, t)e -i2πνt dt is the corresponding MR spectrum. In a 2-D CSI experiment, it is necessary to localize the signal from each individual voxel. The spatial encoding is then achieved by the use of phase encoding gradients, with magnetic strength G x , G y and time duration T. Therefore, the equation (5.6) becomes:

R MR (k x , k y , t) = ν x,y
ŜMR (x, y, ν)e i2π(xkx+yky) dxdy e i2πνt dν (5.7)

where k x and k y are defined as γG x T and γG y T , respectively. Thus, the MR spectrum of the voxel at the spatial position (x, y) can be obtained by a 2D-FT as follows:

ŜMR (x, y, ν) = kx,ky RMR (k x , k y , ν)e -i2π(xkx+yky) dk x dk y (5.8)

where

RMR (k x , k y , ν) = R MR (k x , k y , t)e -i2πνt
dt denotes the FT of the received (time-domain) signal. When considering the signal of the ViP, the same mathematical formulation can be used.Thus, similarly to equation ( 2), the signal of the ViP, R ViP (k x , k y , t), can be written as:

R ViP (k x , k y , t) = ν x,y
ŜViP (x, y, ν)e i2π(xkx+yky) dxdy e i2πνt dν (5.9)

where ŜViP (x, y, ν) denotes the expected ViP reference peak. Now, the question arises: what type of external (ViP) signal needs to be transmitted to the scanner bore in order to obtain a given ViP reference peak in the CSI voxels? Let's consider the case of a reference peak with the same amplitude and frequency in all voxels. In this case, ŜViP (x, y, ν) is independent of the spatial coordinates (x, y), and thus can be written as ŜViP (ν). As a consequence, the integral in equation (5.9) yields a 2D Dirac distribution function in the K-space, δ(k x , k y ), as follows:

R ViP (k x , k y , t) = ν x,y
ŜViP (ν)e i2π(xkx+yky) dxdy e i2πνt dν = S ViP (t)

x,y e i2π(xkx+yky) dxdy = S ViP (t)δ(k x , k y )

(5.10)

where S ViP (t) = ν ŜViP (ν)e i2πνt dν is the designed signal of ViP.

As a result, the ViP peak could be generated by transmitting the ViP signal only once in any of the repetition times of the 2-D CSI sequence. In the current study, we sought to obtain a Lorentzian reference peak ( ŜViP (ν)) and thus a exponential time decay signal (S ViP (t)) was employed. 99 5.3. Materials and methods Chapter 5

Materials and methods

Experiment settings

CSI experiments were performed on a 4.7 Tesla MR scanner (47/40 Biospec, Bruker, Wissembourg, France) using a volume coil as a transmitter/receiver. This coil, of 72-mm innerdiameter, is referred here as to the "scanner coil".

The hardware configuration of the ViP apparatus was as follows: a homemade surface coil (15-mm inner-diameter, the "ViP coil") was used for transmitting the ViP signal; this coil was fixed in the scanner bore 40 cm away from the scanner coil. A low quality factor of the ViP coil ensured no interference with the scanner coil.

Furthermore, a waveform generator (Redstone, Tecmag Inc., Houston, TX, USA) was positioned outside the scanner room to generate the ViP signal. The ViP coil was connected to the output of the waveform generator. An external trigger signal from the MR scanner was used to ensure the synchronization between the ViP signal and the ADC (analog-to-digital converter) of the MR scanner. The time-domain signal representing the ViP reference peak in the CSI spectra was designed using Mathematica (Wolfram Research, Inc.); the numerical values (as magnitude and phase) of this signal served as the input of the waveform generator. Coarse values of the current intensity to be injected in the ViP coil were obtained by electromagnetic modeling, using the Principle of Reciprocity [Hoult, 2011]. These values allowed us to generate a ViP peak with an amplitude comparable to that of a physical phantom. A double spin-echo sequence, with interpulse delays TE/4, TE/2 and TE/4, was employed. TE is the echo time, G slice is the slice gradient and G x and G y are the phase encoding gradients. The signal of the ViP was transmitted in synchronization with the data acquisition CSI experiments were performed on 15 ml cylindrical tube phantoms consisting of pure water or oil-in-water emulsions [Salvati and Gambarota, 2016]. We refer to these phantoms as to the physical phantoms. MR images were acquired with a FLASH sequence for localization purposes and to position the CSI grid on the phantoms. The CSI measurements were performed using a Chapter 5 5.3. Materials and methods 2D-spin-echo CSI sequence. The sequence diagram is shown in Figure 5.8. A double spin-echo module (equivalent to the PRESS scheme [Bottomley, 1987]) was used for signal localization.

The acquisition parameters for 2-D CSI are listed in Tabel 5.1. It should be noted that, the ViP signal, sent by the ViP coil in synchronization with the MR acquisition, was acquired simultaneously with the signal from the physical phantom, as shown in Figure 5.8. 1. Can it be adjusted and fine-tuned for the peak amplitude and peak frequency of the ViP reference?

2. The spatial uniformity and time stability of the ViP reference peaks in different voxels.

Data Analysis

Data analysis was performed using the jMRUI software [START_REF] Stefan | Quantitation of magnetic resonance spectroscopy signals: the jmrui software package[END_REF]. The AMARES algorithm was employed to fit the CSI MR spectra, in order to determine the amplitude and frequency of the resonance peaks. The peak amplitude is here defined as the area under the peak. Since the time-domain signal of ViP was designed as a free induction decay (i.e. exponential decay), the ViP resonance peak was fitted to a Lorentzian lineshape.

After the first successful experiment for a ViP 2-D CSI, we performed additional experiments to explore ViP performances. In the first set of experiments, we investigated the feasibility of adjusting and fine-tuning the amplitude and frequency of the ViP reference peak. With respect to the amplitude of the ViP signal, five measurements were performed with different attenuation factors of the waveform generator while the frequency was kept constant. A linear regression calculation was done, and the coefficient of determination (R 2 ) was obtained in order to assess the goodness of the fit. With respect to the frequency, three ViP 2-D CSI data sets were acquired with an increment of 200 Hz, with a constant attenuation factor.

A second set of experiments was dedicated to assess the spatial un iformity and time stability of the ViP reference peak. The spatial uniformity was evaluated by fitting the ViP peak in each CSI voxel and by calculating the coefficient of variation (CV, defined as the mean divided by standard deviation) of the fitted amplitudes and frequencies. The ViP time stability was determined with 5 repeated measurements. For each measurement, the normalized mean amplitude and mean frequency were obtained and their CV was calculated.

Results

ViP 2-D CSI

Results

The MR spectra of the same ViP 2-D CSI measurement is illustrated in Figure 5.10. The water peak originating from the physical phantom can be observed at the center frequency and the ViP peak is at the left edge of the spectrum, as designed. The amplitude of the ViP peak appears to be uniform in all voxels, as expected. The ViP peak was designed to be toward the edge of the spectrum, approximately at 1 kHz from the center frequency (the spectral bandwidth was 2.5 kHz). The ViP peak appears to be uniform over all voxels. The water resonance peak, originating from the physical phantom, can be observed at the center frequency.

In Figure 5.11, a screen-shot of the CSI visualization/processing tool of the MR Bruker Biospec system is shown. It can be observed that the 64 NMR spectra were recorded in a 2-D image in the form of a 8 × 8 matrix. Each voxel (or pixel) in this matrix displays a NMR spectroscopy of one voxel of the sample.

Discussion and Conclusion

In the current study, we have extended the ViP approach to 2-D CSI methodology. The results of the ViP 2-D CSI experiments indicate that the virtual phantom can provide a reference signal uniform in space and stable in time.

The classical 2-D CSI experiment consists of repeated acquisitions of the NMR signal, for different values of the phase encoding gradients. This allows for the spatial encoding of the signal: for instance, to generate an 8 × 8 matrix, total 64 acquisitions are needed. Based on the theoretical analysis of the 2-D CSI introduced in section 5.2, we developed a strategy to generate a reference peak in all voxels, by sending the ViP signal just in one of the 64 acquisitions. This approach greatly simplifies the experimental protocol; furthermore, it is applicable to different acquisition schemes, such as elliptical and weighted CSI. Since it is sufficient to send the ViPequivalent free-induction-decay only once for generating ViP reference peaks in all voxles, the knowledge of the entire sampling scheme of K-space is not required.

In 2-D CSI experiments, to obtain a reference signal using physical phantoms it is necessary to perform an additional measurement. This also applies to the case of the use of an internal reference, such as water for instance; as a matter of fact, typically 2-D CSI is performed with water-suppression modules, to eliminate the deleterious effects of the water resonance, which is 10 4 times greater than the metabolite resonances. In addition, the internal reference signal might be affected by the pathology, thus it could yield a bias in the quantification.

However, the proposed ViP CSI method eliminates the need of a second measurement; furthermore, in contrast to the physical phantom, it provides a reference signal in all voxels, acquired simultaneously with the signals of interest. There are a number of useful additional features of the ViP approach: e.g. it is possible to adjust and fine-tune the amplitude and frequency of the ViP reference signal (in other words, it can be customized). In the current study, we have also shown the uniformity of the ViP signals over the CSI FOV.

It is important to point out a number of differences between the virtual phantom and the physical phantom. First, they are not interchangeable in all situations; in other words, the virtual phantom cannot always replace the physical phantom. For instance, when developing new pulse sequences, it is of interest to test different RF pulses, gradient schemes and K-space sampling trajectories; to do so, it is necessary to acquire the MR signal originated from a physical phantom, since the ViP signal is transparent to the RF and gradient pulses.

Similarly the physical phantom cannot always replace the virtual phantom. A classical MR measurement consists of two separate processes: the spin excitation, achieved by the use of the MR-system transmission chain, and the signal reception, carried out by MR-system receiver chain. The signal of a physical phantom depends on both MR-system transmission and receiver chain; on the other hand, the signal of ViP, its generation and transmission, is independent from the MR-system transmission chain. Thus, with the ViP signal it is possible to decouple the two processes; as a result, the ViP method can be used to selective probe the MR-system receiver chain.

Conclusion and Perspectives

In this dissertation, we focused on the quantitative techniques in intravoxel incoherent motion (IVIM) imaging and chemical shift imaging (CSI) and have done some novel works. Specifically, 1) we proposed a Sparsity-based All-Voxel Tri-Exponential IVIM (SAVTE-IVIM) algorithm that can automatically identify the existence of potential blood vessels in the target region-ofinterest (ROI) and can simultaneously quantify the IVIM parameters of all voxels in the ROI to evaluate the diffusion, perfusion, and blood vessel confounding effect;

2) we implemented and compared five algorithms for the parameter quantification of vertebral bone marrow IVIM;

3) we proposed to use the virtual phantom (ViP) technique, which was extended from ERETIC (Electronic REference To access In vivo Concentrations) technology to design a virtual phantom, to provide in CSI the reference signal for quantification of the magnetic resonance spectroscopy (MRS).

However, there are still many issues worthy for further discussion:

1. For the ViP CSI, it is worth pointing out that in the current study a single receiving coil was employed in all experiments. In the case of multiple receiving coils, the uniformity of the ViP signal needs to be validated, since it may be affected by the specific vendor-provided reconstruction approach. Indeed, a classical MR measurement consists of two separate processes: the spin excitation, achieved by the use of the MR-system transmission chain, and the signal reception, carried out by the MR-system receiver chain. The signal of a physical phantom depends on both the MR-system transmission and receiver chains; on the other hand, the signal of ViP, its generation and transmission, is independent from the MR-system transmission chain but depending on the MR-syetem receiver chain. Thus, with the ViP signal it is possible to decouple the two processes; as a result, the ViP method can be used to selective probe the MRsystem receiver chain. In this way, the ViP signal behaviour could be used to gain insight into the vendor reconstruction approach.

2. For the bone marrow IVIM, we implemented five algorithms to do the IVIM parameters quantification and to generate the parametric maps. Our work is just a beginning and it would be better if our work can inspire and bring some help to subsequent researchers. Furthermore, more fitting algorithms can be exploited on bone marrow IVIM, e.g. the non-negative least squares (NNLS) algorithm. Especially for the fitting algorithm based on Bayesian inference, in terms of the generation of its Markov chain, the design of its sampling scheme is a direction that could be further studied. Besides, different assumptions about its prior information are also worth exploring.

3. For the SAVTE-IVIM algorithm proposed to cope with the confounding blood vessel effect in the IVIM-MRI model in liver IVIM, it is still a hard task to evaluate real diffusion-weighted magnetic resonance (DW-MR) images since the ground truth regarding the distribution of blood vessels in the ROI is generally unavailable. To cope with this issue, additional MR images such as the MR angiogram could be employed to provide a ground truth in vivo. However it should be noted that: (i) in the current work, we investigate the case of a small (10 -20%) partial volume effect between blood vessels and liver parenchyma; (ii) the IVIM DW-MRI protocol consists typically of data acquisition in free-breathing over a few minutes; (iii) DW-MR images need to be co-registered with those of the MR angiogram. Due to this intrinsic limitation of the co-registration process between free-breathing images acquired on a moving organ in two separate measurements, the approach of an MR angiogram as a reference ground truth for small (10 -20%) partial volume effects was discarded in the current work. Further studies would be necessary to investigate in detail this challenging strategy. In addition, it could be of interest to apply such an approach (the proposed SAVTE-IVIM algorithm) to other organs such as brain, kidney, muscle, etc.. As a matter of fact, there have been a number of studies that have observed tri-exponential features in the IVIM decay measured in these organs. The SAVTE-IVIM model has the potential to reliably quantify the third component observed in these organs.

B The alternating direction method of multipliers (ADMM)

A brief description of the ADMM method (see [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF] and the references therein) used to solve the optimization problems P1 (3.6) and P4 (3.9) is given here. The ADMM algorithm is based on the combination of the dual decomposition and the method of multipliers approaches. It considers a problem of the form:

P A : arg min x,z f (x) + g(z) s.t. Ax + Bz = c (B.1)
where functions f (•) and g

(•) are both convex, x ∈ R Nx , z ∈ R Nz , A ∈ R Nc×Nx , B ∈ R Nc×Nz and c ∈ R Nc .
The solution of the problem P A (B.1) relies on the minimization of the augmented Lagrangian function associated with P A (B.1):

P B : arg min x,z,y L(x, z, y) = f (x) + g(z) + y T (Ax + Bz -c) + ρ 2 ||Ax + Bz -c|| 2 2 (B.2)
where y is a Lagrangian multiplier vector and ρ > 0 is a penalty parameter.

At each iteration, variables x, z, y are updated in an alternative way. Indeed, at the (k + 1)th iteration, each of the latter variables is updated while fixing the other variables to their last estimates:

x k+1 = arg min x L(x, z k , y k ) (B.3) z k+1 = arg min z L(x k+1 , z, y k ) (B.4) y k+1 = y k + ρ(Ax k+1 + Bz k+1 -c) (B.5)
C Technical materials on the solutions of P3 (3.8)

As mentioned in section 3.3.1, the ADMM method [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF] is used to solve the P2 problem (3.7). To do so, the latter is reformulated in terms of minimizing its associated augmented Lagrangian function , L, leading to the P3 problem (3.8) recalled hereafter:

P3 : min a 1 ,a 2 ,a 3 ,d 1 ,d 2 ,d 3 ,z,y L = ψ 2 2 + λ z 1 + y T (a 3 -z) + ρ 2 a 3 -z 2 2 s.t.a > 0, d > 0, ∀ ∈ {1, 2, 3} (C.1)
where z is an N -dimensional latent variable, y ∈ R N is the Lagrangian multiplier and ρ > 0 denotes a penalty parameter. Estimating a , ∀ {1, 2, 3} can be easily computed in a least squares sense as the solution of ∂L ∂a = 0, ∀ ∈ {1, 2, 3}. This leads to:

a 1 = -(J T a 1 J a 1 ) -1 J T a 1 (s + J a 3 a 3 + J a 2 a 2 ) (C.2) a 2 = -(J T a 2 J a 2 ) -1 J T a 2 (s + J a 3 a 3 + J a 1 a 1 ) (C.3) a 3 = ( ρ 2 I N + J T a 3 J a 3 ) -1 ( ρ 2 z - 1 2 y -J T a 3 (s + J a 2 a 2 + J a 1 a 1 )) (C.4)
where I N is the identity matrix of size (N × N ) and

J a = ∂ψ ∂a T = -diag(e -b⊗d )(1 M ⊗ I N ), ∀ ∈ {1, 2, 3} (
d = [J d 1 , J d 2 , J d 3 ] with J d = ∂ψ ∂d T .
Regarding the damping coefficient µ, it is updated following the scheme proposed in [START_REF] Madsen | Methods for nonlinear least square problems[END_REF] with an initial guess being taken as the maximum element in the matrix J T d J d . As far as the latent variable z is concerned, it is computed as follows: [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF]:

z = prox φ, λ ρ (a 3 +y) (C.7)
where the prox function stands for the proximity operator dealing with the non-smooth function φ (here φ(•) = • 1 ) initially proposed in [Ding, 2009]. Regarding the Lagrangian variable y, the dual ascent method is used such that [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF]:

∆y = a 3 -z (C.8)
C.1 Derivation of a 1 (C.2), a 2 (C.3) and a 3 (C.4)

Based on Kronecker product's properties (see [START_REF] Coloigner | Line search and trust region strategies for canonical decomposition of semi-nonnegative semi-symmetric 3rd order tensors[END_REF]) and equation (C.5), we have:

(1 M ⊗ a ) e -b⊗d =diag(e -b⊗d )(1 M ⊗ I N )a =-J a a (C.9) then the all-voxel tri-exponential model (3.5) becomes:

s = L=3 =1 (1 M ⊗ a ) e -b⊗d + = - L=3 =1 J a a + (C.10)
Consequently, we can write ψ = s + J a 1 a 1 + J a 2 a 2 + J a 3 a 3 . Parameters a , ∀ ∈ {1, 2, 3} are then computed as the solution of ∂L ∂a T = 0 T N , ∀ ∈ {1, 2, 3} with:

L = ψ 2 2 + λ z 1 + y T (a 3 -z) + ρ 2 a 3 -z 2 2 (C.11) = ψ T ψ + λ z 1 + y T (a 3 -z) + ρ 2 (a 3 -z) T (a 3 -z) (C.12)
where 0 N is an N -dimensional vector of zeros. Then we have:

∂L ∂a T = ∂L ∂ψ T × ∂ψ ∂a T = 2ψ T ∂ψ ∂a T , ∀ ∈ {1, 2} = 2(s + J a 1 a 1 + J a 2 a 2 + J a 3 a 3 ) T J a (C.13) ∂L ∂a 3 T = 2ψ T ∂ψ ∂a 3 T + y T + ρ(a 3 -z) T (C.14)
which immediately leads to:

a 1 = -(J T a 1 J a 1 ) -1 J T a 1 (s + J a 3 a 3 + J a 2 a 2 ) (C.15) a 2 = -(J T a 2 J a 2 ) -1 J T a 2 (s + J a 3 a 3 + J a 1 a 1 ) (C.16) a 3 = ( ρ 2 I N + J T a 3 J a 3 ) -1 [ ρ 2 z - 1 2 y -J T a 3 (s + J a 1 a 1 + J a 2 a 2 )] (C.17)
D Technical materials on the solutions of P4 (3.9)

As also discussed previously, the ADMM method [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF] is used to solve the P4 optimization problem (3.9). Thus, P4 (3.9) is reformulated as a minimization of its associated augmented Lagrangian function L 1 leading to the P5 problem (3.10) recalled hereafter:

P 5 : L 1 = ψ 2 2 + λ z 1 + ỹT ( ã3 -z) + ρ 2 ã3 -z 2 2 (D.1)
Then, estimating parameters ã1 , ã2 and ã3 is performed by solving ∂L 1 ∂ ã1

T = 0 T N , ∂L 1 ∂ ã2 T = 0 T N and ∂L 1 ∂ ã3
T = 0 T N , respectively. Then the update rules of ã1 and ã2 are given by: ã1 = -2(J T ã1 J ã1 ) -1 J T ã1 (s + J a 3 a 3 + J a 2 a 2 ) (D.2) ã2 = -2(J T ã2 J ã2 ) -1 J T ã2 (s + J a 3 a 3 + J a 1 a 1 ) (D.3) while the one for ã3 is computed by rooting the following equation:

4β ( ã3 3 ) + (ρ • 1 N + 4J T a 3 v) ã3 + (ỹ -ρz) = 0 N (D.4)
where v = s-(1 M ⊗a 1 )e -b⊗d 1 -(1 M ⊗a 2 )e -b⊗d 2 , β = diag -1 (J T a 3 J a 3 ) and where the Jacobian matrix J ã is computed as follows [Coloigner et al., 2014, Lemme 2]: where d = [ d1 T , d2 T , d3 T ] T and the Jacobian matrix J d is given by:

J ã = ∂ψ ∂ ã T = ∂ψ ∂a T × ∂a ∂ ã T = J a × 2 diag( ã ) (D.
J d = ∂ψ ∂ dT = ∂ψ ∂d T × ∂d ∂ dT = J d × 2 diag( d) (D.7)
As mentioned in Appendix C, the damping factor µ is computed as suggested in [START_REF] Madsen | Methods for nonlinear least square problems[END_REF] but with an initial guess corresponding to the maximum entry of J T Following the same gradient computation scheme adopted in Appendix C, we have: then immediately we obtain: ã1 = -2(J T ã1 J ã1 ) -1 J T ã1 (s + J a 3 a 3 + J a 2 a 2 ) (D.12) ã2 = -2(J T ã2 J ã2 ) -1 J T ã2 (s + J a 3 a 3 + J a 1 a 1 ) (D.13) Regarding ã3 we can write J ã3 = ∂ψ ∂ ã3

∂L 1 ∂ ã T = 2ψ T ∂ψ ∂ ã T = 2(s + 3 k=1 J a k a k ) T J ã , ∀ ∈{1, 2} (
= J a 3 • 2 diag( ã3 ). Also, we have a 3 = ã3 2 , and ψ = s + L=3 =1 J a a . Then we can write:

∂L 1 ∂ ã3 T = 2(s + L=3 =1
J a a ) T J ã3 + ỹT + ρ( ã3 -z) T (D.14) = 4a T 3 J T a 3 J a 3 diag( ã3 ) + (ỹ -ρz) T + ρ ã3 T + 4v T J a 3 diag( ã3 ) where v = s + J a 1 a 1 + J a 2 a 2 . Based on properties of the Hadamard product, we can write 4v T J a 3 diag( ã3 ) = 4v T J a 3 ã3 T . Also, according to equation (C.5), we have: which means that J T a 3 J a 3 is a (N × N ) diagonal matrix. In addition, based again on the Hadamard product's properties, the first term on the right-hand side of equation (D.14) can be written as:

4a T 3 J T a 3 J a 3 diag( ã3 ) = 4a T 3 • diag(β) • diag( ã3 ) (D.16) = 4 ã3 T ã3 T β T ã3 T = 4β T ( ã3 3 ) T Also, we have ρ ã3 T = ρ • 1 T N ã3
T . Then, by substituting the latter expression together with (D.16) and the expression of 4v T J a 3 diag( ã3 ) in equation (D.15), we obtain:

∂L 1 ∂ ã3 T = 4β ( ã3 3 ) + (ρ • 1 N + 4J T a 3 v) ã3 + (ỹ -ρz) (D.17)
which means that ã3 can be found by rooting the above equation.

Identification du modèle IRM-IVIM multi-exponentiel: application à la quantification de diffusion et de perfusion de tissue

Contexte et importance de la recherche

Au cours des trente dernières années, en tant que technologie d'application médicale, la technologie de résonance magnétique nucléaire (RMN) a fait de grands progrès. Ces technologies peuvent être grossièrement divisées en deux catégories: l'imagerie par résonance magnétique (IRM) et la spectroscopie par résonance magnétique (MRS). L'introduction de l'IRM dans la communauté clinique et de recherche a ouvert une nouvelle ère de l'imagerie médicale. Par rapport aux rayons X, à la médecine nucléaire et aux technologies d'imagerie par ultrasons, la technologie IRM peut fournir des images avec un contraste élevé des tissus mous. Une autre technologie basée sur la RMN, le MRS, peut être utilisée pour fournir des informations métaboliques. Tant en IRM qu'en MRS, il y a un intérêt croissant pour les approches quantitatives.

Une caractéristique importante de l'IRM est que le contraste de l'image est fondamentalement multiparamétrique, principalement basé sur le temps de relaxation T1, T2 et la densité de protons. En outre, il a été rapidement reconnu que la diffusion et l'écoulement des molécules d'eau dans les tissus peuvent également jouer un rôle important dans le contraste des images IRM. L'une des méthodes d'IRM exploitant la diffusion de l'eau dans les tissus est l'imagerie à mouvement incohérent intravoxel (IVIM) [Le Bihan et al., 1986, Le Bihan et al., 1988, Le Bihan, 2008], qui peut être utilisée pour évaluer simultanément la diffusion et la perfusion. L'imagerie IVIM s'est avérée d'une grande valeur pour le diagnostic, la stadification et le pronostic, et son application clinique dans le corps humain, comme la tête et le cou, la prostate, le foie, les reins et d'autres parties, augmente rapidement [Federau et al., 2014, Kim et al., 2017, Hecht et al., 2017, Meeus et al., 2017, Keil et al., 2017, Le Bihan, 2019, Maximov and Vellmer, 2019]. Cependant, la quantification des paramètres liés à la perfusion IVIM souffre toujours d'une grande variabilité et d'une faible répétabilité. Cela est particulièrement vrai dans les organes tels que le foie, qui contient un riche réseau vasculaire. Les nombreux vaisseaux sanguins de différentes tailles dans le foie représentent un facteur de confusion dans la quantification IVIM de la perfusion tissulaire. Entre l'IRM et le MRS, l'imagerie spectroscopique RM (MRSI) est une technologie d'imagerie qui combine dans une certaine mesure les informations métaboliques du SRM avec la couverture spatiale de l'IRM. Le MRSI est également connu sous le nom d'imagerie par déplacement chimique (CSI). Par rapport à l'IVIM, qui permet de quantifier les paramètres décrivant des phénomènes physiques tels que la diffusion et la perfusion dans les tissus, le CSI est principalement utilisé pour quantifier la teneur en métabolites dans différents organes tels que le cerveau, la prostate, les muscles, etc. Dans la technologie CSI , il existe deux solutions courantes pour fournir des signaux de référence pour la quantification quantitative: le fantôme physique in vitro et le signal de référence in vivo. Cependant, les signaux de référence in vivo peuvent être facilement affectés par la pathologie, et l'utilisation de fantôme physique à l'extérieur du corps présente un certain nombre d'inconvénients, tels que par exemple la perte de résolution spatiale CSI.

Par conséquent, les technologies basées sur la RMN étant largement utilisées dans la recherche clinique et scientifique aujourd'hui, cette thèse a choisi des applications spécifiques dans les domaines de l'IRM et du SRM, que sont l'imagerie IVIM et le CSI. Les techniques quantitatives connexes ont été étudiées et certains des problèmes mentionnés ci-dessus ont été améliorés dans une certaine mesure, ce qui a une signification et une valeur importantes à la fois pour la recherche universitaire et pour l'application clinique.

Le travail principal de cette thèse

Cette thèse portait sur les techniques quantitatives en IVIM et CSI, en particulier:

1. Parcimonie basé tout Voxel Tri-Exponentielle IVIM (SAVTE-IVIM) algorithme Compte tenu de l'intérêt clinique actuel pour la quantification des paramètres des images du foie pondérées en diffusion IVIM, et de la structure vasculaire avec son effet confondant dans le foie, nous avons proposé une nouvelle méthode, un tri-voxel à algorithme exponentiel IVIM (SAVTE-IVIM), qui peut identifier automatiquement l'existence de vaisseaux sanguins potentiels dans la région d'intérêt cible (ROI). En plus de l'objectif principal d'identifier les vaisseaux sanguins potentiels dans un ROI donné, l'algorithme peut quantifier simultanément les paramètres IVIM de tous les voxels dans le ROI pour évaluer la diffusion, la perfusion et l'effet de confusion des vaisseaux sanguins dans chaque voxel. Plus précisément, nous avons proposé un modèle tri-exponentiel basé sur des contraintes clairsemées pour décrire tous les voxels en même temps. En règle générale, les paramètres IVIM sont évalués voxel par voxel (voxel), tandis que SAVTE-IVIM peut quantifier simultanément tous les voxels du ROI. En outre, afin de résoudre le nouveau modèle proposé, un algorithme d'optimisation, basé sur l'idée du multiplicateur de direction alternée (ADMM) avec l'utilisation de l'algorithme de Levenberg Marquardt pour traiter les problèmes non linéaires, a été proposé. Deux stratégies pour les contraintes non négatives inhérentes ont également été introduites.

2. Mise en oeuvre et comparaison de cinq algorithmes d'ajustement pour la quantification IVIM sur la moelle osseuse vertébrale Puisque Marchand et al. appliqué avec succès la méthode IVIM à la quantification de la moelle osseuse en 2014, il y a eu un intérêt croissant pour la moelle osseuse IVIM ces dernières années. Cependant, il existe encore des problèmes avec la qualité d'image actuelle de la moelle osseuse IVIM; en outre, peu d'attention a été accordée à l'étude des algorithmes optimaux pour la quantification IVIM de la moelle osseuse. Il convient de noter qu'un algorithme optimal peut dans une certaine mesure compenser la faible qualité d'image. Au vu des considérations ci-dessus, i) nous avons appliqué un protocole récemment proposé qui améliore la qualité de l'image dans la moelle osseuse IVIM, en utilisant la séquence RESOLVE (lecture de segmentation de long train d'écho variable) et ii) nous avons implémenté cinq algorithmes pour la quantification des paramètres de moelle osseuse vertébrale IVIM. Quatre algorithmes, les algorithmes à une étape, à deux étapes, à trois étapes et à D * fixe, sont basés sur l'idée des moindres carrés (LSQ), et le cinquième est un algorithme bayésien. Une comparaison entre ces algorithmes a été réalisée. De plus, des cartes des paramètres IVIM ont été générées et comparées.

Imagerie par déplacement chimique fantôme virtuel (ViP CSI)

Sur la base de certaines lacunes des solutions grand public pour fournir des signaux de référence en technologie CSI (y compris fantôme physique in vitro et référence interne): la demande d'un scan RM supplémentaire, une diminution de la résolution CSI, etc., nous avons proposé d'utiliser le virtuel technique fantôme (ViP), qui a été étendue à partir de la technologie ERETIC (Electronic REference To Access In vivo Concentrations) pour concevoir un fantôme virtuel, afin de fournir dans CSI le signal de référence pour la quantification MRS. L'amplitude et la fréquence peuvent être conçues sur mesure, pour simuler le signal FID (free induction decay) qui va être acquis avec CSI pour fournir un signal de référence. Dans la technologie CSI classique, avec différents gradients de codage de phase, un processus d'acquisition de signal RMN répété est nécessaire. Par exemple, pour générer une matrice 8 × 8, l'acquisition de données doit être répétée 64 fois. Cependant, le schéma proposé peut générer les mêmes signaux de référence dans les 64 voxels en ne transmettant le signal ViP qu'une seule fois dans l'une des 64 acquisitions de données. De cette manière, un scan RM supplémentaire n'est plus nécessaire, simplifiant le protocole d'imagerie CSI. Il peut également être appliqué à différents schémas d'acquisition de données, tels que le CSI elliptique ou le CSI pondéré. De plus, le signal de référence ViP peut être personnalisé dans son amplitude et sa fréquence et affiche une excellente uniformité spatiale et stabilité temporelle pour tous les voxels de l'échantillon.

La structure de cette thèse

Compte tenu du contenu de recherche principal ci-dessus, cet article est divisé en six chapitres et la structure est la suivante: Indeed, the quantification of the CSI is enhanced in the current work by introducing a virtual phantom as a reference to the acquired image.

Figure 1 . 1 :

 11 Figure 1.1: Spin and precession (a) The spin of the nucleus and the corresponding magnetic dipole moment µ; (b) The precession of the spin magnetization around the external magnetic field B 0 .

Figure 1 . 2 :

 12 Figure 1.2: The arrangement of nuclei: from disorder to orderThe nuclei are arranged in disorder without an external magnetic (left), and rearranged in order after an external magnetic B 0 applied (right).

  resonance imaging direction of the external magnetic field (xy plane) is called the transverse magnetization vector M xy .

Figure 1 . 3 :

 13 Figure 1.3: The flip of magnetization vector with an RF pulse (a) A volume of sample placed in an external magnetic field -→ B 0 B 0 B 0 acquires a longitudinal magnetization -→ M 0 M 0 M 0 parallel to -→ B 0 B 0 B 0 ; (b) an RF pulse flips the magnetization of an angle θ.

Figure 1 . 4 :

 14 Figure 1.4: Longitudinal and transverse relaxation Longitudinal (left) and transverse (right) relaxation of the magnetization for a sample with a T1 relaxation time of 500 ms and a T2 relaxation time of 50 ms.

  In the time-domain, NMR signals are usually the superposition of various sinusoidal oscillations and exponentially decaying signals (Figure 1.5). We call such signals as Free Induction Decay (FID) signals, where the sinusoidal oscillations correspond to the frequency of nuclear precession in the xy plane and the exponential decay corresponds to the transverse relaxation of the nucleus. Sampling and processing of the decaying MR signals leads to different applications, e.g. Magnetic Resonance Spectroscopy (MRS) and Magnetic Resonance Imaging (MRI).

Figure 1 . 5 :

 15 Figure 1.5: Free induction decay Schematic diagram of FID signal generated by transverse magnetization vector.

Figure 1 . 6 :

 16 Figure 1.6: Slice selection gradient How the slice selection gradient G z converts a frequency bandwidth into a physical slice.

Figure 1 . 7 :

 17 Figure 1.7: Frequency encoding and phase encoding Spatial-change in Larmor frequency with the (a) frequency encoding gradient G x and (b) phase encoding gradient G y , given a slice selected by the slice selection gradient.

Figure 1 . 9 :

 19 Figure 1.9: An example of a Spin Echo (SE) pulse sequence

Figure 1 .

 1 Figure 1.10: An example of a Gradient Recall Echo (GRE) pulse sequence

Figure 1 . 11 :

 111 Figure 1.11: Brownian motion model (1-D) 1-D Brownian motion model for one single particle. The position starts from zero at the very first beginning, and 1000 motions in one second are recorded.

Figure 1 .

 1 Figure 1.12: Brownian motion model (2-D and 3-D) Motion trajectories (blue lines) of one single particle diffusing by 2-D and 3-D Brownian motion models.

Figure 1 .

 1 Figure 1.13: 3-D model of normal liver and liver cirrhosis

Figure 1 . 14 :

 114 Figure 1.14: Diagram of the pulsed gradient spin echo (PGSE) techniqueIn a normal spin echo sequence, two diffusion-sensitizing gradients (DG) are applied on either side of the 180 • inversion pulse, symmetrically. In this case, the stationary spins are unaffected by the paired DG while the diffusing spins are dephased, leading to a signal loss. It should be noted that, the 180 • inversion pulse lasts for a very short time and the phase is reversed immediately. For a better understand of the phase changing, this period of time is enlarged.

Figure 1 .

 1 Figure 1.15: K-space sampling with (a) single-shot EPI and (b) multi-shot EPI

Figure 2 . 2 :

 22 Figure 2.2: Diffusion-weighted images (liver) and data analysis (a) An example of the Diffusion Weighted (DW)-MRI in human liver, with a set of M b values. (b) For one voxel, the contributions of apparent diffusion component (in blue) and pseudo-diffusion component (in red) to the total diffusion signal decay. The pseudo-diffusion component is substantially faster than the apparent diffusion component, and so is only observed at low b-values.

Figure 2 . 4 :

 24 Figure 2.4: The step-wise strategy for IVIM data fittingThe input data is marked in blue squares. The final fitting curve is marked in red and the curve in green denotes the contribution from pure diffusion. For this example, the ground truth for the parameters are f = 20%, D = 0.001 mm 2 /s, D * = 0.05 mm 2 /s.

Figure 2 . 6 :

 26 Figure 2.6: Stages of liver fibrosis Liver fibrosis may be assessed by liver biopsy or noninvasive methods, e.g. imaging-based methods. The etiological factors includes alcohol abuse, nonalcoholic fatty liver disease (NAFLD), viral infections and etc.. Regardless of the etiological factors, the liver fibrosis would get through three stages from F0 to F3, and in the end becomes cirrhosis.

Figure 3

 3 Figure 3.1: Retrohepatic vena cava anatomyDiaphragmatic view of the liver and its three-dimensional anatomical relationships involving the vena cava, hepatic pedicle and its branches, biliary tree, and suprahepatic veins[START_REF] Araújo | Retrohepatic vena cava lesion: which we cannot forget?[END_REF]. The inferior vena cava (IVC) and three hepatic veins (HV) are the largest veins in liver.

  0.01, 0.2] and d (n) 3 > 0.2, where d (n) , ∈ {1, 2, 3} stands for the n-th component of the N -dimensional vector d . Indeed, potential negative values of d to 10 -5 , 0.01, respectively. Also, values of d 1 that are greater than 0.01 are set to 0.01 while those of d are greater(smaller) than 0.2 are set to 0.2. Algorithm 2 below summarizes the main steps of the proposed SAVTE-IVIM R algorithm. Mathematical derivations of all update rules figured in Algorithm 2 are given in Appendix C. Algorithm 2 Pseudo code of the SAVTE-IVIM R algorithm I. Initialization of a , d , z and y (∀ ∈ {1, 2, 3}):

Algorithm 3

 3 Pseudo code of the SAVTE-IVIM E I. Initialization of z, ỹ, ã and d ( ∈ {1, 2, 3}): II. While the stop criterion is not fulfilled or the maximum number of iterations is not reached, do:(1) Update J a = ∂ψ ∂a T using equation (C.5) andJ ã = ∂ψ ∂ ãT using equation (D.5), ∀ ∈ {1, 2, 3}(2)Update ã1 using equation (D.2) (3) Update ã2 using equation (D.3) (4) Update ã3 by rooting equation (D.4) (5) Update d using equation (D.6) (6) Update z using equation (D.8) (7) Update ỹ using equation (D.9) (8) a ← ã 2 , d ← d 2

  .4. Simulations were performed with the b-values used for the in vivo data acquisition (see equation (3.5), a set of 12 b-values with b ∈ {0

Figure 3 . 3 :

 33 Figure 3.3: Determination of penalty parametersThe averaged model error per voxel over considered ROIs taken from the available six volunteers in our data set and over four noise levels (noise-free, SNR of 50, 100 and 150), as a function of penalty parameters ρ and λ. For noisy data, results are averaged over 25 MC trials. The smallest model error is marked in dark square.

Figure 3 . 4 :

 34 Figure 3.4: The ground truth of IVIM parametric maps for fully simulated data From left to right, the IVIM parameters are ADC (mm 2 /s), D * (mm 2 /s) and PF (×100%).

Figure 3 . 5 :

 35 Figure 3.5: The simulated confounding effect in blood vessels for fully simulated data The designed ground truth (left) of the spatial distribution of a blood vessel and (right) its binarized representation.

  ≤ n ≤ N ) ranging from 0.2 to 1 mm 2 /s, and amplitudes (i.e., a(n)3 , 1 ≤ n ≤ N ) ranging from 10% to 30% of the voxel intensities, is added to each voxel in the chosen ROI. In this way, the realistic DW-MR images data set, in which the blood vessel confounding effect is taken into account, is generated.

Figure 3 . 6 :

 36 Figure 3.6: An example of the realistic simulated data Left: the realistic DW-MR image (of volunteer 1) at b = 0 s/mm 2 with parameters estimated from a bi-exponential fitting on a real DW-MR image of a liver in axial view. The ROI of size (16 × 16), marked in white square, contains a simulated bifurcating blood vessel. Right: a zoom-in of the considered ROI.

  Figure 3.7 below shows the spatial distribution, A 3 = unvec(a 3 ), of the blood vessel in the considered ROI (out of the 16 simulated ROIs). Similar behavior of the considered algorithms was obtained for the other ROIs. This fact is confirmed in terms of VLE plots as shown in Figure 3.8.

Figure 3 . 7 :

 37 Figure 3.7: Identification of blood vessels (One case of fully simulated data) Binarized estimated spatial distribution of the very fast diffusion component (corresponding to the blood vessle), for one fully simulated ROI, as a function of the SNR using the NNLSbased method and the proposed SAVTE-IVIM R and SAVTE-IVIM E methods.

Figure 3 . 8 :

 38 Figure 3.8: Voxel localization error comparison (fully simulated data) VLE vs. SNR for the NNLS-based method and the proposed SAVTE-IVIM R and SAVTE-IVIM E methods, in the 16 fully simulated ROIs.

Figure 3 . 9 :

 39 Figure 3.9: Identification of blood vessels (two cases of realistic simulated data) The spatial distribution (binary map) of the very fast diffusion component, a 3 (corresponding to the blood vessel effect), for (a) volunteer 1 and (b) volunteer 2, as a function of the SNR and obtained using the NNLS-based method and the proposed SAVTE-IVIM R and SAVTE-IVIM E methods.These results can also be confirmed using the VLE criterion defined in equation (3.11) as depicted in Figures 3.10. One can observe that, compared to the NNLS-based method[START_REF] Gambarota | Eliminating the blood-flow confounding effect in intravoxel incoherent motion (IVIM) using the non-negative least square analysis in liver[END_REF], smaller VLE values are generally obtained using the proposed methods for all SNR values.

Figure 3 . 10 :

 310 Figure 3.10: Voxel localization error comparison (realistic simulated data) VLE vs. SNR for the NNLS-based method and the proposed SAVTE-IVIM R and SAVTE-IVIM E ones for six healthy volunteers. VLE values are averaged over 25 MC trials.

Figure 3 . 11 :

 311 Figure 3.11: Parametric maps (two cases of fully simulated data) The IVIM parameters contains ADC (d 1 , mm 2 /s) (top row), the pseudo-diffusion coefficient D * (d 2 , mm 2 /s) (middle row) and the perfusion fraction PF (f , bottom row). Parametric maps were obtained using the SAVTE-IVIM R , the SAVTE-IVIM E and the NNLS-based methods for (a) SNR = 50, and (b) SNR = 150.

Figure 3 . 12 :

 312 Figure 3.12: Normalized mean square error comparison (fully simulated data) NMSE of the ADC (d 1 ), D * (d 2 ) and PF (f ), using the SAVTE-IVIM R , the SAVTE-IVIM E and the NNLS-based method, for the 16 fully simulated volunteers and the three SNR values (50, 100 and 150).

  Figure3.13 shows, for SNR = 50, the parametric maps (i.e., the spatial distributions) of ADC (d 1 ), D * (d 2 ) and the PF (f ) for two volunteers (similar behavior is obtained for the other

Figure 3 .

 3 Figure 3.13: Parametric maps (two cases of realistic simulated data)The IVIM parameters contains ADC (d 1 , mm 2 /s) (top row), the pseudo-diffusion coefficient D * (d 2 , mm 2 /s) (middle row) and the perfusion fraction PF (f , bottom row). Parametric maps were obtained using the SAVTE-IVIM R , the SAVTE-IVIM E and the NNLS-based methods for SNR=50 of (a) volunteer 1 and of (b) volunteer 2.

Figure 3 .

 3 Figure 3.14: Normalized mean square error comparison (realistic simulated data) NMSE of the ADC (d 1 ), D * (d 2 ) and PF (f ), using the SAVTE-IVIM R , the SAVTE-IVIM E and the NNLS-based method, for different SNR values (150, 100 and 50) and for the six volunteers in our data set. For the sake of readability (i) the logarithm of NMSE values increased by one is considered and (ii) volunteers are irregularly ordered.

  E,TM and γ NNLS,TM . Each component of these series denotes the correlation coefficient computed, for one SNR value and for one ROI (or one volunteer), between a binarized version of the estimated spatial distribution of blood vessels in the considered ROI and a binarized version of the available Ground Truth (GT) or Target Map (TM). Second, the Wilcoxon signed-rank test is applied to each of the pairs (γ R,TM , γ E,TM ), (γ R,TM , γ NNLS,TM ) and (γ E,TM , γ NNLS,TM ). Box plots of the differences γ R,TM -γ E,TM , γ R,TM -γ NNLS,TM and γ E,TM -γ NNLS,TM are shown in Figure 3.15 together with the corresponding p-values (shown in a box).

Figure 3 . 15 :

 315 Figure 3.15: Wilcoxon signed-rank test on blood vessel identification Box plots of the difference between paired correlation series γ R,TM , γ E,TM and γ NNLS,TM used for the Wilcoxon signed-rank test, in terms of the identification results for (a) fully simulated data and (b) realistic simulated data. Obtained p-values are presented in boxes.

  TM and γ θ 3 = γ θ E,TM -γ θ NNLS,TM are shown in Figure 3.16 together with the corresponding p-values (shown in boxes), for both (a) fully simulation and (b) realistic simulation.

Figure 3 . 16 :

 316 Figure 3.16: Wilcoxon signed-rank test on parameters quantification Box plots of the difference between paired correlation series γ θ 1 = γ θ R,TM -γ θ E,TM , γ θ 2 = γ θ R,TM -γ θ NNLS,TM and γ θ 3 = γ θ E,TM -γ θ NNLS,TM for θ ∈ {ADC, D * , PF} used for the Wilcoxon signed-rank test, in terms of the quantification results for (a) fully simulated data and (b) realistic simulated data. Obtained p-values are presented in boxes.

Figure 3 .

 3 Figure 3.17: Four ROIs of volunteer 1 (a) Real DW-MR image of a liver (volunteer 1) in axial view. Four ROIs (white squares) of size (16 × 16) are chosen (from left to right, ROI 1 to ROI 4). (b)-(e) Zoom-in images (ROI 1 to ROI 4) shown with manually enhanced image contrast.

  (b)-(e).

Figure 3 . 18 :

 318 Figure 3.18: Comparison of the visual performances with different contrast levels Three zoom-in images of ROI 4 are shown in three different contrast levels.

Figure 3 .

 3 Figure3.19 shows the spatial distribution (i.e., A 3 = unvec(a 3 )) of the identified third diffusion component (the blood vessel effect) using the SAVTE-IVIM R , the SAVTE-IVIM E and the NNLS-based[START_REF] Gambarota | Eliminating the blood-flow confounding effect in intravoxel incoherent motion (IVIM) using the non-negative least square analysis in liver[END_REF] methods.

Figure 3 . 19 :

 319 Figure 3.19: Identification of blood vessels (volunteer 1)The spatial distribution of blood vessels quantified using the SAVTE-IVIM R , the SAVTE-IVIM E and the NNLS-based methods for four different ROIs.

Figure 3 . 20 :

 320 Figure 3.20: Identification of blood vessels (volunteer 2 -6) (a) Selected ROIs (white squares) from five volunteers (one ROI per volunteer) shown with Standard Contrast (SC), (b) a zoom-in of ROIs shown in manually Enhanced Contrast (EC),(c-e) estimated spatial distribution maps of blood vessels using the SAVTE-IVIM R , the SAVTE-IVIM E and the NNLSbased methods[START_REF] Gambarota | Eliminating the blood-flow confounding effect in intravoxel incoherent motion (IVIM) using the non-negative least square analysis in liver[END_REF] , respectively.

Figure 4 . 1 :

 41 Figure 4.1: An example of the step-wise fitting algorithm In blue squares are the input IVIM data which is the signal intensities to fit. The green curve is reconstructed with a mono-exponential and the diffusion coefficient D which is estimated in the first step using mono-exponential fitting to the high b (≥ 200 s/mm 2 ) values data. The red curve is the final fitted curve reconstructed from all the estimated IVIM parameters with a bi-exponential model. In this example, the ground truth are: f = 20%, D = 0.001 mm 2 /s, D * = 0.05 mm 2 /s and SNR = 50.

S

  3. With parameters D and f being fixed to their estimates D and f , obtained in step 1 and step 2 respectively, the D * parameter is estimated by solving the following single-variable b -S 0 ( f e -bD * + (1 -f )e -b D

  {S b } ∀b is the set of signal intensities, of a single voxel, acquired over all defined b values, θ = {f, D, D * } denotes the set of all model parameters and I stands for the prior information. P (θ| S, I), P ( S|θ, I) and P (θ|I) are, respectively, the parameter posterior probability, the data likelihood and the joint prior probability of the parameters. P ( S|I) is a normalization factor Chapter 4 which is a constant and can be dropped from the calculation. Then equation (4.9) becomes: of statistical independence of the model parameters, f, D and D * , the equation (4.10), can be rewritten as follows: arg max f,D,D * P (f |I)P (D|I)P (D * |I)P ( S|θ, I) (4.11) where P (f |I), P (D|I) and P (D * |I) are the prior probabilities of parameters f , D and D * , respectively.

In

  Figure 4.2, an example of the IVIM DW-MR images on vertebral bone marrow with three different b values (a) 0 s/mm 2 (b) 150 s/mm 2 and (c) 1000 s/mm 2 is shown. In addition, the locations of the five lumbar vertebrae are pointed out with five arrows (in red) from L1 to L5.In addition, we take the L1 as an example, five algorithms were applied to the mean values of the voxels in L1 and the fitting results are shown in Figure4.3. The quantification results of the IVIM parameters, together with the norm of the residual between the fitted and input signal intensities, are listed in Table4.3, in terms of different algorithms.

Figure 4 . 2 :

 42 Figure 4.2: An example of spine IVIM DW-MRI From left to right, the b values (a) b = 0 s/mm 2 , (b) b = 150 s/mm 2 and (c) b = 1000 s/mm 2 .For better visualization, the viewable area here was cut from the whole image (a larger area can be seen in Figure2.8), and the contrast level was manually adjusted.

Figure 4 . 3 :

 43 Figure 4.3: An example of the fitting results with different methods The fitting curves using the five different algorithms, One-Step, Two-Step, Three-Step, Fixed-D * and Bayesian-based methods, over the signal intensities averaged from the ROI in L1.

  Figure 4.4, the box plots illustrate the quantification results, from top to bottom, for different IVIM parameters (a) D, (b) PF and (c) D * .

Figure 4 . 4 :

 44 Figure 4.4: Quantification of IVIM parameters with different methods Estimation results of the parameters (a) D, (b) PF and (c) D * , using five different approaches on the simulated data set. Each box plot draws 20 Monte Carlo trials for one specific SNR out of 10, 20, 50 and 100, and for one method out of the five methods: One-Step, Two-Step, Three-Step, Fixed-D * and Bayesian-based methods. The blue line in each sub-figure indicates the ground truth for each parameter, 48 × 10 -5 mm 2 /s for D, 13% for perfusion fraction (PF) and 18 × 10 -3 mm 2 /s for D * . Different colors were employed to mark different SNR for better visualization.

Figure 4 . 5 :

 45 Figure 4.5: Comparison among the five methods in terms of Relative Error (RE) Comparison among the five methods, One-Step, Two-Step, Three-Step, Fixed-D * and Bayesianbased methods, over the IVIM parameters (a) D, (b) PF and (c) D * . The comparison criterion is the relative error of the estimation results for each parameter and for each SNR. Each boxplot in the figure indicates all trials for one specific SNR out of 10, 20, 50 and 100, and for one parameter. For better visualization, different colors were employed to mark different methods

  Figure 4.8(a). In this figure, maps of the IVIM parameters (b) D, (c) PF and (d) D * gained from five different algorithms (from left to right, they are One-Step, Two-Step, Three-Step, Fixed-D * and Bayesian-based algorithms) are shown.According to Figure4.8, the similar phenomenon can be observed that five algorithms provide comparable maps of diffusion coefficient D and the Bayesian-based algorithm gives smoother maps of the parameters PF and D * compared to the LSQ based algorithms. This can also be derived from the results listed in Table4.4.

Figure 4

 4 Figure 4.7: IVIM MRI parametric maps (simulated data) Estimated parametric maps of the IVIM parameters (a) D, (b) PF

Figure 4

 4 Figure 4.8: IVIM MRI parametric maps (in vivo data) Estimation results of the IVIM parameters (b) D, (c) PF and (d) D * , using the five different approaches on the data in vivo captured from the ROIs in lumbar vertebrae L1 to L5 (a) marked in red. From left to right, the order of the five methods is: One-Step, Two-Step, Three-Step, Fixed-D * and Bayesian-based methods. For better visualization,

Figure 5 . 1 :

 51 Figure 5.1: The pulse sequence for chemical shift imaging (1-D)

Figure 5 . 2 :

 52 Figure 5.2: The pulse sequence for CSI (2-D)

Figure 5 . 3 :

 53 Figure 5.3: The pulse sequence for 2-D CSI: double spin echo

Figure 5 . 5 :

 55 Figure 5.5: NMR spectroscopy quantification: an example of the peak fitting Lorentzian lineshape is used as the model function to fit the data. The goodness of the fit is assessed by the residual.

Figure 5 . 6 :

 56 Figure 5.6: Spectral peak fitting in jMRUI software, using the AMARES algorithm

Figure 5 . 7 :

 57 Figure 5.7: The experimental apparatus for ViP MRIThe steps of ViP MRI: First, the phantom shape of the virtual phantom is designed and the K-space representation is also generated based on it. Second, the simulated K-space lines is converted into an RF signal using a waveform generator. Third, the signal representing K-space line is transmitted by a dedicated RF coil positioned in the scanner bore. In the end, the ViP RF signal is synchronized with the MR scanner data acquisition by the RF unblank signal from the MR console[START_REF] Saint-Jalmes | ViP mri: virtual phantom magnetic resonance imaging[END_REF].

Figure 5 . 8 :

 58 Figure 5.8: The pulse sequence for ViP 2-D CSI A double spin-echo sequence, with interpulse delays TE/4, TE/2 and TE/4, was employed. TE is the echo time, G slice is the slice gradient and G x and G y are the phase encoding gradients. The signal of the ViP was transmitted in synchronization with the data acquisition

  Successful implementations of the VIP 2-D CSI are shown in Figure 5.9, Figure 5.10 and Figure 5.11.

Figure 5 .

 5 Figure 5.9 displays the magnitudes of the 64 time-domain signals of the ViP 2-D CSI experiment.The ViP signal was sent in the first TR only, and for this reason the ViP magnitude needed to be much higher than that of physical phantom, as it can be observed in Figure5.9a. A zoom-in of the first time-domain signal, consisting of both the ViP signal and the signal from the physical phantom, is shown in Figure5.9b for a better visualization of both the magnitude (black line) and the real part (red line). We can observe the contribution of the physical phantom signal, as a form of oscillations, in the initial part of the decay.

Figure 5 . 9 :

 59 Figure 5.9: The time-domain signals acquired in a ViP 2-D CSI measurement With an 8 × 8 matrix size (repetition time = 1200 ms, echo time = 30 ms, acquisition points = 1024, dwell time = 0.4 ms). a. All the 64 signals are displayed in magnitude mode, for easier visualization. The ViP signal was sent only in the first of the 64 acquisitions, with a signal intensity much higher than the signal from the physical phantom. b. A detailed view of the first time-domain signal. Both the magnitude (black curve) and real part (red curve) of the signal are shown

Figure 5 .

 5 Figure 5.10: The MR spectra of the ViP 2-D CSIThe ViP peak was designed to be toward the edge of the spectrum, approximately at 1 kHz from the center frequency (the spectral bandwidth was 2.5 kHz). The ViP peak appears to be uniform over all voxels. The water resonance peak, originating from the physical phantom, can be observed at the center frequency.

Figure 5 .

 5 Figure 5.11: An example of the visualization tool for CSI This is a screen-shot of the CSI visualization/processing tool of the MR Bruker Biospec system. Upper panel: an example of an MR spectrum acquired in the ViP 2-D CSI experiment; the ViP peak as well as the water peak can be observed. The MR spectrum originates from the voxel whose location is illustrated by the hair-cross on the MR image (lower left panel). In the lower right panel, all spectra of the ViP 2-D CSI experiment are shown. It can be noticed that the ViP peak is present in all voxels

Figure 5 .

 5 Figure 5.14: The spatial uniformity and time stability of ViP reference peaks Coefficient of Variation: (a) 0.67%, (b)0.028‰, (c) 0.66%, (d) 0.033‰ and (e) 0.17%. (f) No measurable variation is observed.

  C.5) More details regarding the derivation of equations (C.2), (C.3) and (C.4) are given in the subsequent subsection. The estimation of d , ∀ ∈ {1, 2, 3} is performed using the LM algorithm as follows: ∆d = (J T d J d + µI 3N ) -1 J T d ψ (C.6) where d = [d T 1 , d T 2 , d T 3 ] T , ∆d denotes the difference between two estimates of d taken at two successive iterations, and the matrix J d of size (M N ×3N ) is obtained by concatenating J d 1 , J d 2 and J d 3 as follows: J

  5)A detailed derivation of equations (D.2), (D.3) and (D.4) is given in the subsequent subsection. As far as the estimation of parameters d , ∀ ∈ {1, 2, 3} is concerned, the LM method is employed. Then we have:∆ d = (J T dJ d + µI 3N ) -1 J T dψ (D.6)

  Similarly to equations (C.7) and (C.8), variables z and ỹ, are updated as follows: of ã1 (D.2), ã2 (D.3) and ã3 (D.4)

J T a 3

 3 J a 3 = (1 M ⊗ I N ) T diag(e -2b⊗d 3 )(1 M ⊗ I N )

Chapitre 1 .

 1 Imagerie par résonance magnétique et IRM de diffusion Ce chapitre présente principalement les principes de l'IRM, y compris les phénomènes RMN, la génération d'images RM et les séquences d'imagerie RM. Il introduit également l'imagerie par résonance magnétique pondérée en diffusion (DW-MRI), qui est sensible à la diffusion moléculaire dans les applications d'IRM. C'est le fondement de la technologie IVIM. Chapitre 2. Mouvement incohérent IntraVoxel Dans ce chapitre, nous allons d'abord passer en revue la technique d'imagerie du mouvement incohérent intravoxel (IVIM), y compris les bases de la diffusion et de la perfusion, le modèle bi-exponentiel classique pour IVIM et l'analyse quantitative d'IVIM. De plus, pour deux applications spécifiques, dans le foie humain et dans la moelle osseuse humaine, nous introduisons d'abord la fibrose hépatique et l'évaluation des tumeurs hépatiques par imagerie IVIM. Ensuite, nous passons en revue le développement à jour de l'IVIM dans la moelle osseuse humaine. Chapitre 3. Parcimonie basé tout Voxel Tri-Exponentielle IVIM (SAVTE-IVIM) algorithme Concernant l'application de la technologie IVIM dans le foie humain, nous avons proposé un algorithme IVIM tri-exponentiel tout-voxel basé sur la parcimonie (SAVTE-IVIM). Ce chapitre présente d'abord ce modèle IVIM tri-exponentiel puis présente l'algorithme SAVTE-IVIM utilisé pour résoudre le problème du nouveau modèle, y compris ses deux stratégies différentes pour gérer les contraintes non négatives. Enfin, la faisabilité et les performances de l'algorithme proposé sont vérifiées par une analyse expérimentale de données simulées et de données cliniques, ainsi que par comparaison avec l'algorithme basé sur NNLS. Chapitre 4. Mise en oeuvre et comparaison de cinq algorithmes d'ajustement pour la quantification IVIM sur la moelle osseuse vertébrale Concernant l'application de la technologie IVIM dans la moelle osseuse humaine, nous avons implémenté et comparé cinq algorithmes d'ajustement quantitatif. Ce chapitre présente d'abord les cinq algorithmes, y compris quatre algorithmes basés sur la méthode des moindres carrés, à savoir la méthode en une étape, en deux étapes, en trois étapes et en D fixe, et le cinquième est basé sur l'inférence bayésienne, à savoir bayésienne algorithme basé sur. Par la suite, à travers l'analyse expérimentale de données simulées et de données cliniques, les performances des cinq algorithmes considérés dans la quantification des paramètres IVIM ont été comparées. Chapitre 5. ViP CSI: imagerie par déplacement chimique fantôme virtuel Concernant l'application de la technologie ViP dans CSI, ce chapitre présente le principe de CSI et l'application de la quantification MRS. Dans cette étude, nous avons proposé un schéma ViP-CSI, qui a été vérifié par des expériences de comparaison.Enfin, dans le chapitre Conclusion et perspectives, les travaux de recherche sur le texte intégral sont résumés, quelques problèmes dans la recherche de cet article sont discutés, et le contenu qui pourra être étudié en profondeur dans le futur est prospecté.

Titre:

  Identification du modèle IRM-IVIM multi-exponentiel : application à la quantification de diffusion et de perfusion de tissue Mots clés : IVIM-IRM, modèle multi-exponentiel, parcimonie, optimisation proximale Résumé : Le travail dans cette thèse concerne l'identification du modèle IntraVoxel Incoherent Motion (IVIM) des images IRM de diffusion. Cette identification a pour objective la quantification de diffusion et de perfusion des tissues humaines avec une application aux images IRM du foie et de la moelle osseuse. La première partie de ce travail de thèse porte sur la détection d'une potentielle présence d'un volume partiel dans la région d'intérêt. Ce problème d'identification est formulé sous forme d'un problème d'identification d'un modèle IVIM tri-exponentiel sous contraintes de parcimonie de la distribution spatiale du volume partiel dans la région d'intérêt étudiée et de la non-négativité des paramètres de ce modèle. La prise en compte d'informations spatiales est rendue possible grâce au modèle IVIM multi-voxel introduit dans cette thèse et permettant un traitement simultané de l'ensemble de voxels définis dans la région d'intérêt considérée. L'approche proposée, ne permets pas seulement une bonne identification du volume partiel mais aussi une quantification de la diffusion et de la perfusion du tissue étudié. La deuxième partie de cette thèse porte sur la comparaison de méthodes d'identification du modèle IVIM bi-exponentiel avec une application aux images IRM de diffusion de la moelle osseuse. La dernière partie de cette concerne l'amélioration de la quantification des images métaboliques (Chemical Shift Imaging) par l'introduction d'un fantôme virtuel comme une référence dans l'image acquise. Title: Multi-exponential IVIM MRI model identification: application to the quantification of tissue diffusion and perfusion aKeywords: IVIM-MRI, multi-exponential model, sparsity, proximal optimization Abstract: This PhD work concerns the identification of the IntraVoxel Incoherent Motion (IVIM) diffusion weighted MR images. This is for the quantification of tissue diffusion and perfusion with application to MR images of liver and bone marrow. The first study in this PhD work tackles the detection of potential partial volume in the region of interest (ROI). This identification problem is dealt with as a constrained tri-exponential IVIM model identification problem. Indeed, in addition to the non-negativity constraints of the IVIM model parameters, the spatial distribution of the volume partial is supposed to be sparse over the considered ROI. Incorporating such spatial priors is henceforth possible thanks to an allvoxel IVIM-MRI model introduced in this PhD. The proposed approach allows for an identification of potential partial volume in the ROI together with a quantification of tissue diffusion and perfusion, in only one single step. The second part of this PhD addresses the quantification of the tissue diffusion and perfusion in bone marrow. A full comparative study of two families of approaches: leastsquare based and Bayesian are conducted. The final part of this PhD work is a contribution to the quantification of Chemical Shift Imaging (CSI).

Table 1 . 1 :

 11 Examples of T1 and T2 values at 1.5 T

	Tissue	T1 (ms)	T2 (ms)
	Gray matter	950	100
	White matter	600	80
	Muscle	900	50
	CSF *	4500	2200
	Fat	250	60
	Blood	∼ 1400	∼ [180, 250]

Table 1 . 2 :

 12 Estimates of ADC values reported in healthy human organs

	Organ	No. of Studies	ADC range (×10 -3 m 2 /s)
	Liver	28	0.81 -2.40
	Kidney	15	1.50 -3.76
	Spleen	14	0.63 -2.81
	Pancreatic	13	1.02 -2.61
	Gallbladder	6	2.51 -3.50
	Prostate	13	0.90 -1.99
	Uterus 1	13	1.27 -2.09
	FBT 2	13	1.33 -2.37

Table 2

 2

	.1: Values of IVIM parameters in healthy organs
	Organ	D (×10 -3 m 2 /s)	PF (%)
	Cardiac 1	1.26 -3.33	8 -31
	Pancreatic 1	0.56 -1.65	24 -55
	Kidney 1	1.65 -2.07	8 -37
	Placenta 1	1.40 -1.74	25 -56
	G&W matter 1,2	0.65 -0.84	0.2 -14 3
	Breast 1	1.84 -1.96	5.3 -8.9 3
	1 Estimates of the diffusion coefficient D (ADC) and perfusion fraction in some healthy
	human organs, reported by previous studies with the IVIM bi-exponential model, see
	[Le Bihan et al., 2018], and Chapters 9, 12, 13, 16 and 24 in it
	2 G&W matter:

the Gray matter and White matter 3 The tissues with perfusion fraction less than 10% are often called low perfused tissues

	2.2. IVIM in liver	Chapter 2
	25

Table 2

 2 is listed in Table 2.2. Table 2.2: Recommendations of parameters settings for IVIM imaging on liver lesions .3, ranges of reported IVIM parameters estimates are listed, for both benign and malignant liver lesions. It should be noted that, in early studies, the IVIM estimation and ADC estimation based on mono-exponential are mixed. Thus, in the Table 2.3, the results of ADC (mono-exponential) are removed. Table 2.3: Values of IVIM parameters in liver lesions

	Variable	Recommendation
	b values	More than 10 b values
	SNR	Not less than 40
	Acquisition scheme	Free breathing
	Field strengh	Either 1.5T or 3.0T
	In

Table 2 .

 2 4: Values of IVIM parameters in normal bone marrow D (×10 -3 m 2 /s) D * (×10 -3 m 2 /s)

	PF (%)

Table 3 . 1 :

 31 Parameter settings for DW-MRI

	Pulse sequence	SE-EPI eDWI 1
	Parallel imaing	ASSET 2
	Repetition time (ms)	2000
	Echo time (ms)	54
	Field of view (cm 2 )	40 × 30
	Slice thickness (mm)	8
	Matrix	256 × 256
	Hold breath?	NO

Table 3

 3

	.2: Quantification of IVIM model parameters (volunteer 1)
	ADC (×10 -4 mm 2 /s)	SAVTE-IVIM R	SAVTE-IVIM E	NNLS-based
	ROI 1	9.0 ± 1.3	8.8 ± 1.2	8.6 ± 1.8
	ROI 2	9.7 ± 1.2	9.6 ± 1.2	9.2 ± 1.6
	ROI 3	11.0 ± 2.1	11.0 ± 2.1	10.0 ± 2.5
	ROI 4	8.6 ± 7.5	9.0 ± 8.9	8.6 ± 8.2
	PF (%)	SAVTE-IVIM R	SAVTE-IVIM E	NNLS-based
	ROI 1	19.2 ± 9.1	24.4 ± 7.4	27.4 ± 8.8
	ROI 2	13.7 ± 7.7	15.2 ± 6.2	19.5 ± 9.5
	ROI 3	24.9 ± 13.3	26.7 ± 11.5	31.3 ± 13.2
	ROI 4	49.1 ± 11.8	49.5 ± 10.5	59.0 ± 10.0
	D * (mm 2 /s)	SAVTE-IVIM R	SAVTE-IVIM E	NNLS-based
	ROI 1	0.057 ± 0.022	0.022 ± 0.018	0.024 ± 0.018
	ROI 2	0.047 ± 0.053	0.031 ± 0.034	0.029 ± 0.031
	ROI 3	0.048 ± 0.046	0.038 ± 0.037	0.034 ± 0.035
	ROI 4	0.023 ± 0.036	0.019 ± 0.026	0.021 ± 0.025

Table 3 . 3 :

 33

	ADC (×10 -4 mm 2 /s)	SAVTE-IVIM R	SAVTE-IVIM E	NNLS-based
	Volunteer 2	10.1 ± 1.6	10.4 ± 1.4	4.6 ± 4.7
	Volunteer 3	10.8 ± 1.4	10.8 ± 1.4	10.8 ± 1.6
	Volunteer 4	9.6 ± 3.4	9.8 ± 3.0	7.9 ± 2.8
	Volunteer 5	16.2 ± 2.2	16.5 ± 2.1	15.1 ± 4.6
	Volunteer 6	12.7 ± 2.1	12.8 ± 2.0	12.8 ± 2.3
	PF (%)	SAVTE-IVIM R	SAVTE-IVIM E	NNLS-based
	Volunteer 2	9.8 ± 7.5	12.7 ± 4.2	44.3 ± 24.3
	Volunteer 3	16.1 ± 6.0	16.5 ± 4.1	19.9 ± 15.5
	Volunteer 4	11.7 ± 9.7	14.1 ± 7.3	24.3 ± 16.9
	Volunteer 5	20.3 ± 10.4	19.5 ± 9.7	26.0 ± 17.7
	Volunteer 6	30.8 ± 12.8	29.1 ± 10.7	41.5 ± 14.1
	D * (mm 2 /s)	SAVTE-IVIM R	SAVTE-IVIM E	NNLS-based
	Volunteer 2	0.109 ± 0.073	0.075 ± 0.065	0.093 ± 0.040
	Volunteer 3	0.130 ± 0.031	0.126 ± 0.029	0.125 ± 0.031
	Volunteer 4	0.089 ± 0.067	0.076 ± 0.060	0.055 ± 0.055
	Volunteer 5	0.036 ± 0.022	0.037 ± 0.021	0.038 ± 0.025
	Volunteer 6	0.073 ± 0.056	0.063 ± 0.044	0.068 ± 0.047

Table 4 . 1 :

 41 Parameters settings for M-H sampling method to generate Markov chain

Table 4 .

 4 2: Parameters settings for spine IVIM DW-MRI

	Sequence

Table 4 . 3 :

 43 Quantification of the IVIM parameters of the ROI in L1 for one volunteer Nres denotes the norm of the residual between the fitted and input signal intensities.

	Method	D (×10 -3 mm 2 /s)	f (%)	D * (×10 -3 mm 2 /s)	N res	1
	One-Step	0.39	9.1	11	0.011
	Two-Step	0.41	7.9	13	0.013
	Three-Step	0.41	7.4	13	0.014
	Fixed-D *	0.41	7.8	15	0.014
	Bayesian-based	0.38	9.6	12	0.014
	1

Table 4 .

 4 4: Quantification results of the IVIM parameters for in vivo data

	Method	D (×10 -3 mm 2 /s)	f (%)	D * (×10 -3 mm 2 /s)
	One-Step	0.42 ± 0.14	14.4 ± 11.4	14.0 ± 10.6
	Two-Step	0.49 ± 0.13	11.7 ± 10.1	16.3 ± 10.6
	Three-Step	0.49 ± 0.13	10.1 ± 10.4	16.5 ± 9.1
	Fixed-D *	0.46 ± 0.13	11.6 ± 10.6	15.0
	Bayesian-based	0.43 ± 0.14	14.6 ± 4.3	15.4 ± 2.6

Table 5

 5

	.1: Acquisition parameters for 2-D CSI
	TR	1200 ms
	TE	30 ms
	CSI matrix size	8 × 8
	FOV	4 × 4 cm 2
	Slice thickness	8 mm
	Spatial resolution	5 × 5 mm 2
	Bandwidth	2500 Hz
	No. of data points	1024
	Different sets of experiments were performed to investigate the performances of the ViP
	technique applied to 2-D CSI:

It should be noted that here the "TR" denotes the abbreviation of Trust-Region, and "TR" was also used to represent Repetition Time in Chapter 1. Since it can be easily distinguished by context, it will not be specifically explained again.

For the sake of clarity, variables d1 and d2 are used in this chapter to designate the commonly used nomenclature ADC and D * in the MR community.

eDWI: enhanced diffusion-weighted imaging.

ASSET: array spatial sensitivity encoding technique.

The term "lipid" is sometimes used as a synonym for fats.

It should be noted that here we use back the nomenclature for the diffusion and perfusion parameters, which is typically used by the MR community when the standard bi-exponential approach is applied. More precisely, the apparent diffusion coefficient ADC and pseudo-diffusion coefficient are respectively denoted by D and D * .

Chapter 4

Implementation and Comparison of Five Fitting Algorithms for IVIM Quantification on Vertebral Bone Marrow

Introduction

As pointed out previously in this thesis, it has been more than three decades since Le Bihan introduced the MRI approach, which allows for the measurements of both tissue diffusion and perfusion [Le Bihan et al., 1986]. Recently, the IVIM technique has been "rediscovered" as an attractive method for assessing diffusion and perfusion in different organs [START_REF] Yamada | Diffusion coefficients in abdominal organs and hepatic lesions: evaluation with intravoxel incoherent motion echo-planar mr imaging[END_REF], Le Bihan, 2008, Luciani et al., 2008, Döpfert et al., 2011, Shinmoto et al., 2012, Delattre et al., 2012, Federau et al., 2012, Ichikawa et al., 2013, Liu et al., 2013].

With respect to the application of IVIM in bone marrow, the first measurements were performed by [START_REF] Yeung | Bone marrow diffusion in osteoporosis: evaluation with quantitative mr diffusion imaging[END_REF][START_REF] Yeung | Bone marrow diffusion in osteoporosis: evaluation with quantitative mr diffusion imaging[END_REF]. In this study, Yeung et al. observed two sets of subjects, 20 healthy volunteers and 44 subjects with osteoporosis. The ADC was evaluated with a mono-exponential fitting and the diffusion coefficient D with a bi-exponential IVIM model. However, Yeung et al. did not provide the assessment of perfusion (i.e. the perfusion fraction and the pseudo-diffusion coefficient D * ) because the quality of the acquired data did not allow them to perform a reliable data-fitting of the perfusion parameters.

During a ten years gap (2004 -2014), the evaluation of marrow perfusion was generally assessed using the dynamic contrast-enhanced (DCE) imaging technique [START_REF] Griffith | Vertebral bone mineral density, marrow perfusion, and fat content in healthy men and men with osteoporosis: dynamic contrastenhanced mr imaging and mr spectroscopy[END_REF], Griffith et al., 2006, Patel et al., 2010].

In 2014, Marchand et al. performed an IVIM study on bone marrow of healthy volunteers and the quantification of the IVIM diffusion and perfusion parameters with a bi-exponential model was successfully performed [START_REF] Marchand | Mri quantification of diffusion and perfusion in bone marrow by intravoxel incoherent motion (ivim) and non-negative least square (nnls) analysis[END_REF]. The optimization of data acquisition parameters (in particular, the b-values) and a dedicated data analysis approach, based on the NNLS strategy, contributed to the successful IVIM quantification.

Chapter 4 precision of the results. On the one hand, the proximity of the box's median (the horizontal line inside the box) and the ground truth can represent the accuracy. On the other hand, the height of the box, that is, the degree of dispersion of the result, can characterize the precision. Now we look at Figure 4.4 again, the results are grouped into five parts, from left to right in order corresponding to the five algorithms, which are One-Step, Two-Step, Three-Step, Fixed-D * and Bayesian-based algorithms. In this way, it is more obvious that the Bayesian-based method shows significantly better performance than the least squares based algorithms (One-Step, Two-Step, Three-Step and Fixed-D * algorithms), both in terms of accuracy and precision, especially when the SNR is low (e.g. SNR = 10 and SNR = 20). In addition, for each algorithm, as expected when the SNR increases, the performances of quantification on D and on PF improves. For the quantification of D * , the Bayesian-based algorithm gives the estimates of smaller variance compared to the three step-wise strategies, for all SNR values.

To further compare the Bayesian-based algorithm and the four algorithms based on LSQ, the relative error (RE, the smaller the better) between the result (in Figure 4.4) and the ground truth was calculated and the results of RE comparison is displayed in two figures in different forms, Figure 4.5 and Figure 4.6.

In Figure 4.5, the RE values are grouped according to SNR values (10, 20, 50 and 100). In this way, it is clearer that the performances are better when the SNR is higher. According to this figure, five algorithms give comparable results on quantification of D. With respect to the quantification of perfusion-related parameters PF and D * , for higher SNR (i.e. SNR = 50 and SNR = 100), the results are similar for all methods while for lower SNR (i.e. SNR = 10 and SNR = 20), the Bayesian-based algorithms provides results of much lower relative errors.

In Figure 4.6, the RE values are displayed in scatter plots way. In this way, the Bayesianbased algorithm is compared to each of other four LSQ based algorithms, which are One-Step, Two-Step, Three-Step and Fixed-D * algorithms. For each plot in this figure, the two values of its coordinates, the x-axis value is the RE of Bayesian-based algorithm result for certain SNR (characterized with different colors) and the y-axis value is the RE of one algorithm (out of the four LSQ based methods). The straight line with a slope of one is called equal line. Simply put, if the points are below the equal line, then a higher RE value of Bayesian-based algorithm for the same data set can be observed. This means that Bayesian-based algorithm lose this comparison. According to Figure 4.6, it can be observed that most of the points are above the equal line, which means the Bayesian-based algorithm is generally better than the other four least squares based algorithms.

Chapter 5

ViP CSI: virtual phantom chemical shift imaging

Introduction

In the last three decades, there have been considerable advances in magnetic resonance techniques for medical applications. These techniques can be broadly classified into two categories:

-Magnetic Resonance Imaging (MRI) MRI provides anatomical, structural and functional information on a whole organ, with a high spatial resolution (mm 3 )

-Magnetic Resonance Spectroscopy (MRS)

MRS yields metabolic information on a very limited region of interest within an organ, with a low spatial resolution (cm 3 ).

Between these two categories, we can consider the MR spectroscopic imaging (MRSI), is referred to as chemical shift imaging (CSI) [START_REF] Brown | Nmr chemical shift imaging in three dimensions[END_REF]. With CSI technology, metabolic information can be acquired on a large volume, and with a spatial resolution that is intermediate between MR imaging and MR spectroscopy techniques. CSI applications have nowadays covered the majority of organs and tissues, including brain, prostate, breast, liver, muscle and kidneys [START_REF] Smith | Magnetic resonance spectroscopy in medicine: clinical impact[END_REF].

In this section, we introduce the magnetic resonance spectroscopy in section 5.1.1, and the chemical shift imaging in section 5.1.2.

Quantification of MRS

With the CSI technology, we can obtain for each voxel one NMR spectroscopy. For example, assuming there is a sample that contains two MR-visible compounds, compound A and compound B. With respect to the concentrations of compound A and compound B, suppose that it is known for compound A (e.g. supposing it as 10 mM) and unknown for compund B. Of course the two compounds need to have different chemical shifts and thus can give different peaks in the spectrum. Figure 5.4 shows an example of a processed spectrum that might be obtained from this hypothetical sample. In Figure 5.4, for the horizontal axis, denoting the frequency or the chemical shift, no unit is marked here. For a MR spectroscopy, the spectrum can be illustrated in the scale of Hertz (Hz) directly. However, the absolute frequency shifts in Hz is dependent on the external magnetic field. But if the resonance frequencies ω are assigned relative to a reference frequency ω ref , the derived chemical shift δ can become independent of the applied external magnetic field, which facilitates the comparison of spectra that have been acquired at different scanner B 0 field strengths:

Thus, ppm is more commonly used instead of Hz.

With respect to the vertical axis, the amplitude of the MR spectrum, one can notice that the units are arbitrary units. If denoting the area under a given spectral peak as the intensity of the peak, the peak intensity should be directly proportional to the concentration of nuclei giving rise to that peak. So the spectral quantification generally involves measuring the intensity of a spectral peak or a group of peaks. And then the measured peak intensities can be converted into

Customizable amplitude/frequency of ViP reference peak

The feasibility of adjusting and fine-tuning the amplitude and frequency of the ViP signal are shown in Figure 5.12 and Figure 5.13. Figure 5.12 shows the results of five ViP 2-D CSI measurements. The ViP amplitudeaveraged over all voxels -was obtained for each measurement, and it is plotted in Figure 5.12 as a function of the attenuation factor. An excellent linearity was observed with the coefficient of determination R 2 = 0.9979.

APPENDICES

A The non-negative least squares (NNLS) algorithm

The non-negative least squares (NNLS) algorithm is initially proposed in [START_REF] Lawson | Solving least square problems[END_REF] to solve the following optimization problem:

Algorithm 4 below summarizes the main steps of the NNLS algorithm:

Algorithm 4 Main steps of the NNLS algorithm INPUTS: A, y, .

MAIN LOOP: WHILE R = ∅ and max(w) > DO

• Let j in R be the index of max(w) in w

• Move j from R to P

• Let A P be A restricted to the variables included in P

• Let t be a vector of the same length as x, let t P denote the sub-vector with indexes from P and t R the sub-vector with indexes from R

(3) Move all indices j from P to R such that x j = 0 (4) Set t P = (A P ) T A P -1 (A P ) T y, set t R to zero

L1

One-Step (0.42 0.14) 1 14.4 11.4 14.0 10.6 0.53 0.20 13.5 11.1 12.3 10.2 0.41 0.14 10.4 10.6 12.7 10.3 Two-Step 0.49 0.13 11.7 10.1 16.3 10.6 0.60 0.21 10.0 8.5 15.2 10.6 0.47 0.12 8.6 8.3 14.5 10.8 Three-Step 0.49 0.13 10.1 10.4 16.5 9.1 0.60 0.21 8.5 9.0 15.7 9.4 0.47 0.12 7.3 9.0 15.3 9.2 Fixed-D * 0.46 0.13 11.6 10.6 15.0 0.0 0.59 0.20 9.3 8.8 15.0 0.0 0.44 0.11 8.3 8.9 15.0 0.0 Bayesian-based 0.43 0.14 14.6 4.3 15.4 2.6 0.54 0.19 13.0 4.3 14.6 2.4 0.39 0.12 13.2 3.7 14.9 2.2

L2

One-Step 0.43 0.16 15.5 12.7 13.4 10.4 0.55 0.25 12.8 10.2 12.9 10.0 0.42 0.16 10.4 10.4 13.9 10.8 Two-Step 0.49 0.16 12.5 10.9 15.9 10.6 0.61 0.23 10.0 8.7 15.8 10.2 0.46 0.14 8.9 9.0 14.9 10.9 Three-Step 0.49 0.16 10.7 11.6 15.6 9.0 0.61 0.23 9.0 10.2 15.8 8.8 0.46 0.14 7.7 9.8 15.8 9.3 Fixed-D * 0.47 0.16 11.9 11.4 15.0 0.0 0.59 0.26 9.7 8.9 15.0 0.0 0.44 0.14 8.6 9.5 15.0 0.0 Bayesian-based 0.45 0.16 14.6 4.5 15.4 2.3 0.55 0.23 13.2 4.0 14.9 2.2 0.39 0.14 13.4 3.9 15.1 2.3

L3

One-Step 0.44 0.16 14.3 11.8 13.9 10.9 0.51 0.22 12.6 10.8 12.1 10.3 0.42 0.19 10.5 10.1 12.8 10.3 Two-Step 0.51 0.17 11.3 10.3 16.4 11.0 0.56 0.22 9.8 8.7 13.7 10.1 0.46 0.19 8.9 9.0 14.9 11.0 Three-Step 0.51 0.17 9.4 10.6 16.1 9.5 0.56 0.22 8.7 11.1 14.4 8.7 0.46 0.19 7.6 9.8 14.8 9.3 Fixed-D * 0.48 0.17 10.8 10.7 15.0 0.0 0.55 0.24 9.2 9.5 15.0 0.0 0.45 0.18 8.1 9.2 15.0 0.0 Bayesian-based 0.45 0.16 14.3 4.6 15.3 2.4 0.50 0.22 13.4 4.1 14.8 2.6 0.40 0.19 12.9 3.9 14.9 2.2

L4

One-Step 0.46 0.18 13.6 11.1 13.0 10.2 0.50 0.22 12.2 11.1 13.6 10.5 0.43 0.23 12.6 10.7 16.7 11.2 Two-Step 0.52 0.17 10.7 8.7 15.7 10.5 0.54 0.14 9.5 9.0 15.5 10.3 0.49 0.24 10.3 9.0 19.0 10.6 Three-Step 0.52 0.17 8.9 9.0 16.0 9.2 0.54 0.14 8.7 11.3 15.9 8.9 0.49 0.24 9.3 10.2 18.7 9.4 Fixed-D * 0.50 0.16 10.4 9.3 15.0 0.0 0.54 0.20 9.4 9.5 15.0 0.0 0.45 0.22 11.0 10.1 15.0 0.0 Bayesian-based 0.46 0.16 14.6 4.0 15.2 2.2 0.49 0.18 13.0 4.3 14.7 2.2 0.44 0.22 14.2 4.0 15.9 2.2

L5

One-Step 0.41 0.30 13.5 11.6 15.2 11.1 0.60 0.30 10.9 10.7 13.3 10.2 0.36 0.21 11.7 11.8 13.5 10.9 Two-Step 0.47 0.29 11.3 10.0 17.3 10.8 0.64 0.30 8.9 8.9 14.5 10.1 0.42 0.22 9.6 9.8 15.5 10.9 Three-Step 0.47 0.29 10.0 11.1 17.4 9.4 0.64 0.30 8.0 10.8 15.1 8.6 0.42 0.22 8.6 10.4 15.4 9.0 Fixed-D * 0.43 0.30 11.3 10.9 15.0 0.0 0.63 0.30 8.4 8.9 15.0 0.0 0.38 0.21 9.5 10.0 15.0 0.0 Bayesian-based 0.43 0.26 14.4 4.6 15.8 2.3 0.58 0.28 12.2 4.3 14.7 2.3 0.36 0.19 13.9 4.9 15. Fixed-D * 0.53 0.35 8.5 9.6 15.0 0.0 0.36 0.18 14.3 12.0 15.0 0.0 0.48 0.28 10.6 10.7 15.0 0.0 Bayesian-based 0.50 0.34 12.3 4.5 14.7 2.5 0.36 0.15 16.2 4.6 15.8 2.3 0.45 0.28 13.5 4.5 15.2 2.3 L2 One-Step 0.47 0.24 12.2 10.6 13.7 10.4 0.31 0.26 15.9 13.3 12.6 10.0 0.40 0.24 14.4 12.2 13.3 10.2 Two-Step 0.52 0.24 10.3 9.0 15.7 10.5 0.39 0.28 13.0 11.2 15.1 10.7 0.45 0.25 12.1 11.1 14.8 10.2 Three-Step 0.52 0.24 9.0 10.0 15.1 8.8 0.39 0.28 10.6 11.5 14.9 9.1 0.45 0.25 10.9 12.5 15.4 8.9

Fixed-D * 0.50 0.24 10.0 9.8 15.0 0.0 0.36 0.28 11.8 11.7 15.0 0.0 0.43 0.24 11.8 11.8 15.0 0.0 Bayesian-based 0.47 0.24 13.8 4.4 15.2 2.2 0.37 0.26 15.6 4.3 15.5 2.2 0.41 0.23 14.8 4.8 15.4 2.4

L3

One-Step 0.36 0.19 16.4 12.5 13.1 10.7 0.27 0.22 17.2 13.2 14.0 10.6 0.44 0.23 11.7 11.0 14.3 10.9 Two-Step 0.44 0.20 12.8 10.2 15.4 10.8 0.37 0.25 12.9 10.7 16.8 10.7 0.49 0.23 9.8 9.2 15.6 11.0 Three-Step 0.44 0.20 10.8 10.6 15.8 9.5 0.37 0.25 11.1 11.2 16.9 9.5 0.49 0.23 8.2 9.3 15.7 9.4 Fixed-D * 0.42 0.17 11.9 10.8 15.0 0.0 0.33 0.25 13.0 11.3 15.0 0.0 0.46 0.22 9.5 10.0 15.0 0.0 Bayesian-based 0.40 0.16 15.1 3.9 15.4 2.1 0.37 0.23 15.6 4.7 15.8 2.2 0.42 0.21 13.8 4.2 15.3 2.3 L4 One-Step 0.38 0.17 Three-Step 0.45 0.17 9.6 11.2 15.5 9.4 0.38 0.36 10.3 12.6 16.9 8.9 0.45 0.16 9.2 9.3 15.4 8.9

LIST OF PUBLICATIONS

Liu, J., Gambarota, G., Shu, H., Jiang, L., Leporq, B., Beuf, O., & Karfoul, A. (2020). On the identification of the blood vessel confounding effect in intravoxel incoherent motion (IVIM) Diffusion-Weighted (DW)-MRI in liver: An efficient sparsity based algorithm. Medical Image Analysis, 61, 101637.

Liu, J., Karfoul, A., Jiang, L., Shu, H., Saint-Jalmes, H., & Gambarota, G. (2018). ViP-CSI: Virtual Phantom Chemical Shift Imaging. Applied Magnetic Resonance, 49(4), 369-380.

AKNOWLEDGEMENTS

This dissertation was conducted in the frame of CRIBs (Centre de Recherche en Information Biomédicale sino-français), which is an international associate French-Chinese laboratory (Université de Rennes 1 -France, INSERM -France, Southeast University -China).

First of all, I would like to express my sincere thanks to my supervisors for their invaluable support during these years of the thesis. I thank very deeply Giulio Gambarota for being an excellent mentor to me, for his knowledge and kindness, for his helping on writing scientific articles and doing research. I also thank Ahmad Karfoul for his knowledge of mathematics, for his guidance on writing and researching. I am grateful to them for being with me forever even when I am in China -9000 kilometers away from them.

Besides, I would like to share my appreciation to my Chinese supervisor Huazhong SHU. He provided to me the chance to be a PhD student in LTSI (Laboratoire Traitement du Signalet de l'Image, Université de Rennes 1) to do this thesis in France working with Giulio and Ahmad. When I was in Rennes, he cared about me not only on my work but also on my daily life.

I would like to also thank Longyu JIANG for opening the door of scientific research for me.

I would like to also thank Jean-Marie Bonny and Nathalie Just who accepted to review my thesis, Virginie Callot and Said Moussaoui who accepted to participate in my jury.

In addition, I will not forget all the friendliness and help I got from LTSI and members there, Hervé Saint-Jalmes, Jean-Jacques Bellanger, Louis, Xu HAN, Wentao XIANG. And the "BM-ROOM" will always be with me.

Finally, I would like to share my thank to my parents for always being there standing behind me. The last person to mention is my wife Hao SONG. She motivated and encouraged me much for this thesis. And from now on in all my life, she is the person always being with me.