Changyu Sun

R Gourdon

M Jean

Marc Lancelin

M Gérard Scorletti
email: gerard.scorletti@ec-lyon.fr

H Charles

Gudrun Mme

Bornette

Safia Ait

Chalal Insa

M Lagarde

M Didier

M Johannes Kellendonk

M Laboune

Jean-Yves Buffiere
email: yves.buffiere@insa-lyon.fr

Mega

M Philippe Boisse

Scso Scso

M Obadia

Lionel Sec

Viviane Polsinelli

J Y Toussaint

Dr Jean-Philippe Thiran

Dr François Rousseau

Dr Rachid Deriche

Dr Su Ruan

Dr Nicole Vincent

Patrick

Dr Edward W Hsu

This thesis concerns the reconstruction and description of orientation distribution functions (ODFs) in high angular resolution diffusion imaging (HARDI) such as q-ball imaging (QBI). QBI is used to analyze more accurately fiber structures (crossing, bending, fanning, etc.) in a voxel. In this field, the ODF reconstructed from QBI is widely used for resolving complex intravoxel fiber configuration problem. However, until now, the assessment of the characteristics or quality of ODFs remains mainly visual and qualitative, although the use of a few objective quality metrics is also reported that are directly borrowed from classical signal and image processing theory. At the same time, although some metrics such as generalized anisotropy (GA) and generalized fractional anisotropy (GFA) have been proposed for classifying intravoxel fiber configurations, the classification of the latters is still a problem. On the other hand, QBI often needs an important number of acquisitions (usually more than 60 directions) to compute accurately ODFs. So, reducing the quantity of QBI data (i.e. shortening acquisition time) while maintaining ODF quality is a real challenge.

In this context, we have addressed the problems of how to reconstruct high-quality ODFs and assess their characteristics. We have proposed a new paradigm allowing describing the characteristics of ODFs more quantitatively. It consists of regarding an ODF as a general three-dimensional (3D) point cloud, projecting a 3D point cloud onto an angle-distance map (ADM), constructing an angle-distance matrix (ADMAT), and calculating morphological characteristics of the ODF such as length ratio, separability and uncertainty. In particular, a new metric, called PEAM (PEAnut Metric), which is based on computing the deviation of ODFs from a single fiber ODF represented by a peanut, was proposed and used to classify intravoxel fiber configurations. Several ODF reconstruction methods have also been compared using the proposed metrics. The results showed that the characteristics of 3D point clouds can be well assessed in a relatively complete and quantitative manner. Concerning the reconstruction of high-quality ODFs with reduced data, we have proposed two methods. The first method is based on interpolation by Delaunay triangulation and imposing constraints in both q-space and spatial space. The second method combines random gradient diffusion direction sampling, compressed sensing, resampling density increasing, and missing diffusion signal recovering. The results showed that the proposed missing diffusion signal recovering approaches enable us to obtain accurate ODFs with relatively fewer number of diffusion signals.

Résumé

Changyu SUN II

Reconstruction et description des fonctions de distribution d'orientation en imagerie de diffusion à haute résolution angulaire

Résumé

Ce travail de thèse porte sur la reconstruction et la description des fonctions de distribution d'orientation (ODF) en imagerie de diffusion à haute résolution angulaire (HARDI) telle que l'imagerie par q-ball (QBI). Dans ce domaine, la fonction de distribution d'orientation (ODF) en QBI est largement utilisée pour étudier le problème de configuration complexe des fibres. Toutefois, jusqu'à présent, l'évaluation des caractéristiques ou de la qualité des ODFs reste essentiellement visuelle et qualitative, bien que l'utilisation de quelques mesures objectives de qualité ait également été reportée dans la littérature, qui sont directement empruntées de la théorie classique de traitement du signal et de l'image. En même temps, l'utilisation appropriée de ces mesures pour la classification des configurations des fibres reste toujours un problème. D'autre part, le QBI a souvent besoin d'un nombre important d'acquisitions pour calculer avec précision les ODFs. Ainsi, la réduction du temps d'acquisition des données QBI est un véritable défi.

Dans ce contexte, nous avons abordé les problèmes de comment reconstruire des ODFs de haute qualité et évaluer leurs caractéristiques. Nous avons proposé un nouveau paradigme permettant de décrire les caractéristiques des ODFs de manière plus quantitative. Il consiste à regarder un ODF comme un nuage général de points tridimensionnels (3D), projeter ce nuage de points 3D sur un plan angle-distance (ADM), construire une matrice angle-distance (ADMAT), et calculer des caractéristiques morphologiques de l'ODF telles que le rapport de longueurs, la séparabilité et l'incertitude. En particulier, une nouvelle métrique, appelé PEAM (PEAnut Metric) et qui est basée sur le calcul de l'écart des ODFs par rapport à l'ODF (représenté par une forme arachide) d'une seule fibre, a été proposée et utilisée pour classifier des configurations intravoxel des fibres. Plusieurs méthodes de reconstruction des ODFs ont également été comparées en utilisant les paramètres proposés. Les résultats ont montré que les caractéristiques du nuage de points 3D peuvent être évaluées d'une manière relativement complète et quantitative. En ce qui concerne la reconstruction de l'ODF de haute qualité avec des données réduites, nous avons proposé deux méthodes. La première est basée sur une interpolation par triangulation de Delaunay et sur des contraintes imposées à la fois dans l'espace-q et dans l'espace spatial. La deuxième méthode combine l'échantillonnage aléatoire des directions de gradient de diffusion, le compressed sensing, l'augmentation de la densité de ré-échantillonnage, et la reconstruction des signaux de diffusion manquants. Les résultats ont montré que les approches de reconstruction des signaux de diffusion manquants proposées nous permettent d'obtenir des ODFs précis à partir d'un nombre relativement faible de signaux de diffusion. This motion is known as diffusion. Shown here are sample random paths that can be taken by molecules within a medium that allows isotropic (the same is in every direction) diffusion. As time passes, the net distance traveled by a molecule increases. [START_REF] Huettel | Functional Magnetic Resonance Imaging[END_REF] ................... 16 Fig. 1.10 Diagram shows the diffusion-driven random trajectory (red line) of a single water molecule during diffusion. The dotted white line (vector r ) represents the molecular displacement during the diffusion time interval, between 1 t0  and 2 t . [START_REF] Hagmann | Understanding Diffusion MR Imaging Techniques: From Scalar Diffusion-Weighted Imaging to Diffusion Tensor Imaging and Beyond[END_REF] ............................ 17 [Stejskal, 1965[Stejskal, , 1965] ] diffusion gradient according to their respective position. 180 o pulse and a second diffusion gradient rephrase those spins. Yet moving spins do not have the same position when the second diffusion gradient is applied, there at echo time all spins are not perfectly in phase, resulting in an echo of smaller magnitude. [START_REF] Chabert | Diffusion Signal in Magnetic Resonance Imaging: Origin and Interpretation in Neurosciences[END_REF] are used to represent diffusion displacements. The diffusion is highly anisotropic in fibrous tissues such as white matter and the direction of greatest diffusivity is generally assumed to be parallel to the local direction of white matter. [START_REF] Jing | [END_REF] Left: ODF from q-ball imaging [Descoteaux, 2007] Middle: ODF from q-ball imaging within constant solid angle [Aganj, 2010] Right: FOD from Constrained spherical deconvolution [Tournier, 2007] (b) K=30, the ODF map reconstructed in the noise-free condition. Left: the map reconstructed using the original data (K=30). Right: the map reconstructed from the data recovered using our method. (c), (d) and (e) are the ODFs maps reconstructed from the original diffusion signal (left) and after processing by our proposed method (right), at SNR=30, 20 and 10 respectively. .......... 67 reduced HARDI dataset is used to reconstruct the phantom ("fiber cup") data. ODF on the left: constructed using the original q-ball imaging. ODF on the right: constructed after processing using the proposed approach. 4.6 ODF maps of the circled fiber crossing region of the b=2650s/mm 2 phantom dataset, obtained using our method and analytical q-ball imaging. (a) b0 image of the physical phantom from the b=2650s/mm 2 dataset. The circled region is a fiber crossing region. (b) ODF maps of the circled region obtained using analytical q-ball imaging with complete 64 directions (left), analytical qball imaging with 29 directions (middle), and our method with 29 directions (right). .................. 72 Fig. 4.7 ODF maps of real brain data, obtained using analytical q-ball imaging and our diffusion data recovering combined with analytical q-ball imaging. Left: the analytical q-ball imaging results. Middle: b0 image of the brain with the red box indicating the region where ODFs are reconstructed. Right: our method's results. The ODFs at each voxel is superimposed on a grayscale background modulated by the GFA [Tuch, 2004] at that voxel (black: GFA=0; white: GFA=1). . AQBI, (b) the correct detection rates of a single and two-fiber diffusion pattern by our proposed method in red, the directly using q-ball imaging method in blue and the spherical harmonic based CS method in green for b=3000s/mm 2 at SNR=30, 20 and 10 dB, (c) the average length ratio at the same conditions, (d) the same as above, the average separability, (e) the average uncertainty, The three numbered rows (1, 2, and 3) correspond respectively to Rectangle 1, Rectangle

Content of Figures

Changyu SUN VII

Content of Figures

Introduction Générale

L'imagerie par résonance magnétique de diffusion (IRMd) est un outil puissant pour étudier la microstructure de tissus, par exemple la structure neurologique du cerveau ou l'architecture des fibres myocardiques du coeur ex-vivo ou in-vivo. Parce que la fonction du cerveau ou du coeur est fondamentale liée à son architecture de fibres spécifique, l'étude cette dernière s'avère fondamentale. L'eau étant le constituant majoritaire des tissues biologiques, ce sont les protons des molécules d'eau qui contribuent majoritairement au signal IRM de diffusion. S'il n'y a des obstacles, les molécules d'eau diffusent de manière isotrope. Par contre, dans le corps (cerveau ou coeur), le déplacement des molécules d'eau est restreint par la structure biologique. Les molécules d'eau se déplacent plus librement dans la direction des fibres qu'à travers elles. Les signaux IRM de diffusion permettent donc de mesurer la structure sous-jacente, particulièrement pour la distribution d'orientation des fibres.

Le signal de résonance magnétique peut être rendu sensible à la diffusion par l'utilisation de séquences IRM spécifiques. On obtient des images pondérée en diffusion (DWI) à travers l'atténuation du signal due à la diffusion des molécules d'eau. La diffusivité peut être mesurée au moyen du coefficient de diffusion apparent ADC. Elle peut également être mesurée suivant des directions spécifiées. Dans le cas de la diffusion restreint où la direction de diffusion n'est pas connue a priori, le coefficient scalaire ADC n'est pas suffisant pour décrire la diffusion anisotrope. Pour résoudre ce problème, beaucoup de méthodes sont introduites.

Les modèles qui traitent et extraire les caractéristiques de la diffusion peuvent être classés suivant le critère de la résolution angulaire d'analyse de l'anisotropie de la diffusion.

Parmi les modèles d'analyse à basse résolution angulaire, l'imagerie du tenseur de diffusion (DTI) est la plus connue. Avec ce modèle, il est possible de reconstruire le tenseur de diffusion en chaque point à partir de six images pondérées en diffusion, qui correspondent à six directions de gradient de diffusion. Ceci est dû au fait qu'un tenseur de diffusion étant symétrique, il est suffisant de le représenter par six coefficients. Ainsi, pour calculer le tenseur, un système de six équations linéaires est suffisant. Géométriquement, un tenseur de diffusion est souvent représenté par un ellipsoïde. L'indice d'anisotropie fractionnelle (FA) est une mesure d'anisotropie de la diffusion de l'eau dans un voxel et la mesure moyenne de la diffusion (MA) d'un voxel est aussi proposée.

DTI présente deux avantages principaux: un temps d'acquisition relativement court (donc le patient n'a pas besoin de rester immobile longtemps) et une technique standard de l'IRMd (nombreux résultats de recherche existent). Cependant, à cause des problèmes de croisement des fibres dans le cerveau ou pour vérifier s'il y a croisement des fibres dans le coeur, DTI n'est pas suffisant pour récupérer toute l'information angulaire. L'hypothèse de diffusion gaussienne est une autre limitation majeure.

Afin de pouvoir mesurer l'information angulaire et morphologique de manière plus précise à partir du signal de diffusion IRM, nous avons besoin des modèles d'analyse en imagerie à haute résolution angulaire (HARDI). Le modèle basé sur l'espace q est appelé imagerie spectrale de diffusion (DSI, diffusion spectrum imaging) [START_REF] Wedeen | Mapping Complex Tissue Architecture with Diffusion Spectrum Magnetic Resonance Imaging[END_REF]. Mais DSI nécessite un échantillonnage intensif dans l'espace q et un gradient de diffusion puissant, d'où l'intérêt de recourir à l'imagerie HARDI [DS Tuch et al., 2002] qui utilise une réduction de l'échantillonnage du signal dans l'espace q. Le plus simple modèle de HARDI est un mixte des fonctions gaussiennes [ DS Tuch, 2002] ; il permet de résoudre le problème de croisements des fibres à condition que le nombre des fibres doivent être connus a priori. C'est la limitation principale de cette méthode.

Il y a deux principaux modèles qui n'ont pas besoin de connaitre le nombre de fibres a priori : l'imagerie q-ball (QBI) [ [START_REF] Ds Tuch | [END_REF] et la déconvolution sphérique [Tournier et al., 2004]. De nombreuses études ont été menées sur ces deux types de modèle. [ [START_REF] Ds Tuch | [END_REF][START_REF] Descoteaux | Apparent Diffusion Coefficients from High Angular Resolution Diffusion Imaging: Estimation and Applications[END_REF][START_REF] Descoteaux | Regularized, Fast, and Robust Analytical Q-Ball Imaging[END_REF][START_REF] Hess | Q-Ball Reconstruction of Multimodal Fiber Orientations Using the Spherical Harmonic Basis[END_REF]Tournier et al., 2008;[START_REF] Kezele | Spherical Wavelet Transform for ODF Sharpening[END_REF]Michailovich et al., 2010c]. Les résultats de l''imagerie QBI sont sous la forme de la fonction d'orientation de distribution (ODF). La base utilisée pour représenter le signal de diffusion est une base d'harmoniques sphériques. QBI propose une estimation du profil de la diffusion sans a priori. Normalement QBI a besoin de plus de 60 directions de gradient de diffusion [START_REF] Michailovich | Spatially Regularized Compressed Sensing for High Angular Resolution Diffusion Imaging[END_REF] et des valeurs de b (b-values) supérieures à 3000s/mm 2 , ce qui reste long pour une application en routine clinique. La méthode proposée dans [Michailovich, 2011] a utilisé compressed sensing pour réduire le nombre d'acquisitions. L'indice d'anisotropie fractionnelle généralisée (GFA) est la variance de l'ODF normalisée. [ [START_REF] Ds Tuch | [END_REF] L'estimation par déconvolution sphérique est linéaire et rapide avec une base harmonique sphérique. La déconvolution sphérique peut résoudre le problème de croisements des fibres. Mais elle a besoin de choisir un noyau de déconvolution, et ce choix dépend souvent des données ou se fait via une simulation gaussienne. Les résultats de déconvolution sphérique sont sous forme de l'orientation de distribution de fibre (FOD). En général, dans la même condition d'acquisition, FOD est plus tranchant qu'ODF. Mais le fait d'utiliser les noyaux de déconvolution tous identiques est une hypothèse forte. Dans [Tournier, 2008], la déconvolution sphérique par contrainte permet de résoudre des angles de croisement à 60 o , même avec une b-value de 1000 s/mm 2 . D'autres modèles comme la transformation d'orientation de diffusion (DOT) [Özarslan et al., 2006] ont également été proposés.

L'évaluation des différentes HARDI méthodes est importante. Dans ce domaine, l'erreur angulaire est une mesure souvent utilisée [START_REF] Hess | [END_REF]Descoteaux, 2007;Tournier et al., 2007;Michailovich, 2011]. Mais le problème avec cette mesure est que dans le cas de données réelles, la vérité terrain n'est pas connue. Le taux de détection des croisements est aussi une mesure qui nécessite la connaissance de la vérité terrain pour comparer.

La deuxième sorte de mesure pour comparer les résultats HARDI est la comparaison de la forme de HARDI quantité comme ODF. Dans cette catégorie, l'erreur moyenne quadratique entre la forme d'estimation ODF et la vérité terrain de ODF est une méthode plus souvent utilisé. [START_REF] Hess | [END_REF] Les méthodes sont Kullback-Leibler (sKL) [Zhan et al.] et Dip test [START_REF] Aganj | Reconstruction of the Orientation Distribution Function in Single and Multiple Shell Q-Ball Imaging within Constant Solid Angle[END_REF]. La fonction SIP [START_REF] Jiao | Uncertainty Visualization in HARDI Based on Ensembles of ODFs[END_REF] et Hessian [START_REF] Seunarine | Exploiting Peak Anisotropy for Tracking through Complex Structures[END_REF] sont aussi utilisés. La méthode Boostrap est aussi proposée dans [START_REF] Cohen-Adad | Quality Assessment of High Angular Resolution Diffusion Imaging Data Using Bootstrap on Q-Ball Reconstruction[END_REF] pour évaluer la qualité des données HARDI ainsi que les résultats.

Dans ce contexte, on peut trouver trois problèmes :

Le premier problème : l'évaluation des méthodes HARDI reste encore un problème, malgré l'existence d'un certain nombre de mesures visuelle ou empruntées directement de la théorie classique de traitement d'images, parce qu'elles sont simples et intuitives à calculer. Cependant, ces mesures existantes ne sont pas toujours pertinentes pour caractériser les résultats en HARDI.

Le deuxième problème : La configuration des patterns de diffusion peut refléter la microstructure des tissus. De nombreuses techniques ont été proposées pour analyser de telles configurations. Cependant, la classification des différentes configurations de patterns de diffusion est toujours un problème.

Le troisième problème : l'imagerie QBI est présentée comme une méthode permettant d'analyser la structure des fibres dans un voxel. Bien que cette technique d'imagerie échantillonne sur une sphère et par conséquent réduit déjà de manière significative les signaux à acquérir pour la reconstruction des ODFs, elle nécessite néanmoins encore beaucoup d'acquisitions pour le calcul précis des ODFs (souvent plus de 60 directions de gradients), ce qui apparaît encore loin pour une utilisation en routine clinique.

La présentation de ce manuscrit de thèse est structurée comme suit : Partie I : Etat de l'art Chapitre 1: Une brève présentation des principes de l'IRM. Nous introduisons les notations sur l'espace k, les images pondérées en diffusion, le coefficient de diffusion apparent et le tenseur de diffusion.

Chapitre 2 : Technique de la reconstruction de croisement des fibres Un état de l'art des méthodes de reconstruction des croisements des fibres est présenté dans ce chapitre. Nous avons aussi exposé des mesures d'anisotropie en HARDI et des méthodes d'évaluation des résultats d'HARDI.

Partie II : Contributions Chapitre 3 : l'évaluation des caractéristiques des nuages de points tridimensionnels (3D) en utilisant des mesures morphologiques.

Nous proposons un modèle pour évaluer quantitativement et précisément la fonction d'orientation de distribution (ODF) ou plus généralement d'autres nuages de points 3D via la mesure de leurs caractéristiques morphologiques. L'idée consiste à projeter un nuage de points 3D sur une carte angle-distance (ADM), construire une matrice angle-distance (ADMAT), et calculer de nouveaux paramètres (le rapport de longueur, la séparabilité et l'incertitude). Trois méthodes ont été comparées en utilisant les mesures proposées sur des données simulées et réelles.

Etant donné une fonction de distribution d'orientation (ODF) d'analyse à haute résolution angulaire (HARDI) de diffusion, il est considéré comme un nuage de points 3D qui est ensuite projetée sur une carte angle distance. De celle-ci, une matrice angle-distance est construite et de nouveaux indicateurs, tels que le rapport de longueur, séparabilité et de l'incertitude, sont calculées.

Chapitre 4 : Reconstruction des ODFs de haute qualité par l'interpolation des fonctions de probabilité de distribution Dans ce chapitre, nous proposons une idée pour faire face au problème d'obtenir la fonction d'orientation distribution précise tout en utilisant moins de signaux. La méthode consiste à reconstruire, par interpolation triangulation de Delaunay, des signaux de diffusion dans des directions de gradients de diffusion où l'acquisition de données n'avait pas eu lieu.

Un aperçu de la méthode de reconstruction des signaux de diffusion non acquis par triangulation de Delaunay. Un nombre réduit de données HARDI (des données de fantôme physique "fiber cup") est utilisé. Gauche : images pondérées en diffusion et ODF. Droite : images pondérées en diffusion et ODF après interpolation par triangulation de Delaunay.

Chapitre 5 : Dans ce chapitre, nous proposons une autre idée pour obtenir plus précis ODF en utilisant moins de mesures. Elle a basé sur des approches d'acquisition comprimée en l'espace Q, mais notre objectif est d'obtenir des signaux de diffusion dans plusieurs directions de diffusion par ré-échantillonnage non acquises ou donnée manquant. La deuxième partie : Utilisant notre méthode en un fantôme physique. Gauche : résultats de intravoxel complexité quantification de l'image. Dans le proposé "PEAM" mesure, les pixels plus clairs indiquant une plus grande complexité. Droite : résultats de la classification des voxels en deux classes: les orientations homogenerous et orientations hétérogènes.

        1 1 1 1 1 1 ,, ,, B       K sub P P K P P YY YY           Ec  sub qB         1 1 1 1 1 1 ,, ,, B       K RS RS K RS RS YY YY         Modèle d'
La troisième partie : Utilisant notre méthode en une donnée réelle de cerveau macaque. Gauche : les cartes ODF obtenus par AQBI, la couleur de fond est le résultat de la complexité de introvoxel quantification de ce carte ODF. Droite : les résultats de la classification des pixels en deux classes.

Perspective : L'existence de croisement de fibre dans la structure de fibre myocardique sera vérifiée dans prochain travail. On va utiliser deux HARDI méthodes (CSD et AQBI) dans les données simulation myocardique (images pondérée en diffusion) par l'imagerie en lumière polarisée (PLI) et donnée réel en coeur cochon et coeur humain. Pour montrer la relation entre l'échelle du voxel et le croisement des fibres, les simulations avec différent échelles serons utilisées. Deux modèles des HARDI avec différent noyau mathématique vont être utilisés pour éviter l'influence des modèles des reconstructions. Le voisinage de croisement de fibre doit être vérifié pour comparer directions de fibres adjacentes avec les directions croisements. Les différent tranches des donnée réel doivent être aussi calculé et reconstruit par les deux méthodes pour réduire l'influence de bruit.

Part I Part I Introduction and background

Résumé en français

L'imagerie par résonance magnétique (IRM) est une technique d'image médicale permettant d'obtenir des vues de l'intérieur du corps de façon non invasive. Le principe d'IRM est présenté d'abord. On a décrit brièvement les concepts de précession, de résonance et de relaxation. Afin de bien comprendre comment obtenir une image à partir du signal magnétique, espace réciproque est présenté. L'application successive de gradients de champ magnétique dans chaque direction spatial par un gradient de sélection de coupe, le codage de phase, et le codage par la fréquence. La séquence écho de spin est introduite, parce que la séquence écho de spin est la séquence considérée comme la séquence de base de l'IRM. Nous l'expliquons pour comprendre l'acquisition de signal d'IRM pratique.

L'IRMd a été présenté basée sur l'imagerie par résonance magnétique. Pour bien comprendre l'IRMd, le phénomène de diffusion, en suite le principe de l'IRM est introduit. Deux phénomènes de diffusion (la diffusion libre, diffusion contrainte). Séquence Stejskal et Tanner (PGSE) est présenté pour réaliser de mesurer le déplacement microscopique des spins.

L'imagerie de diffusion est basée de la séquence standard, par exemple écho de spin, et des gradients de diffusion destinés à pondérer le signal par la diffusion, comme PGSE. L'image du coefficient de diffusion apparent (ADC) est introduite pour obtenir les signaux en une direction. Mais l'image d'ADC dépend de la direction du gradient de diffusion.

Pour obtenir information sur directions of microstructure générale, l'imagerie du tenseur de diffusion (DTI) est traité de signal de diffusion en mesurer l'information anisotropies par 3 diffusion directions. Résolution des équations du DTI, diffusivité moyenne et indice d'anisotropie (FA) sont aussi expliqué ici. En revanche, l'hypothèse de tenseur de diffusion de ce modèle, qui suppose la diffusion locale comme étant gaussienne, ne peut pas montrer le plus complexe diffusion configuration comme croisement de fibres.

Magnetic Resonance Imaging (MRI)

Magnetic resonance imaging is a technology that provides high quality images utilizing the magnetic properties. It is a non-invasive method for acquiring high resolution images of the internal structures of the body. This ability makes it as a technique widely use in hospitals for medical diagnosis, staging of disease and for follow-up without exposure to ionizing radiation. [Revett, 2011] In Fig. 1.1, we find that MRI systems can be comprised by magnet for static field, magnets for gradient fields (needed for position encoding), RF transmission coil (for transmission of RF pulses), RF receive coil (for receiving the MR signal) and Reconstructor (for calculating images from the acquired data, it is not shown in Fig. 1.1).

Basis of MRI

1.1.1 Magnetic resonance imaging (MRI) exploits the properties of the body's natural magnetic to make images from the body. For imaging purposes, the hydrogen nucleus (a simple proton) is used because of its abundance in water and fat. In Fig. 1.2, the hydrogen proton is spinning on its axis. They can be seen as a small bar magnet. In normal condition, they spin with their axes randomly.

In Fig. 1.3, when they are placed in a strong magnetic field, for example, an MRI scanner, the protons' axes all line up after a magnetic field applying. This uniform alignment creates a magnetic vector oriented along the axis of the MRI scanner. [Edelman, 1993] Fig. 1.3 With external magnetic field. [Edelman, 1993] In medical MRI, people are interested in looking at the NMR signal from water and fat, the major hydrogen containing components of the human body. The resonance equation is the principle of all magnetic resonance imaging, which shows that the resonance frequency v of a spin is proportional to the magnetic field, 0 B , the frequency v depends on the gyromagnetic ratio of the particle,  . After transmission of the RF pulse cease, the effective magnetic field returns to B 0 . M lies in the xy plane and proceeded about the z-axis (Fig. 1 .5,Left). This precession induces a current in a receiving coil, which the scanner uses to translate an image. As soon as the RF pulse ceases, however, the signal begins to diminish as Mxy (where the signal is at its maximum) returns to Mz. This free induction decay (FID) is shown in Fig. 1.5, Right. K-space is a formalism widely used in magnetic resonance imaging introduced in 1979 by [Twieg, 1983] and [Ljunggren, 1983].

Slice selective technique

For imaging in two-dimension (2D), we choose through the three-dimensional (3D) tissues. Then the signal stems only from this slice. The image is arisen by encoding the two in-plane directions. The resonance frequency of the spins during the z gradient is:

  oz v z = B + G z . (1.2)
After selecting the slice by z direction, strong external linear magnetic field gradients are applied across the imaging region in xy plane. Therefore, different points in space become recognized by different resonance frequencies, which can let the location of spins emitting RF fields to be fixed by their frequency. where  is the gyromagnetic ratio. ()  x, y is the spatial distribution of protons within a sample. The Fourier transform of this equation generates ()  x, y , the MR image [Ljunggren, 1983;[START_REF] Stehling | [END_REF]. [Ljunggren, 1983;Twieg, 1983], as shown in Fig. 1.7, the MR signal can be expressed in terms of spatial frequencies k, the spatial frequency coordinates axes k x and k y (or the read and phase-encode directions.) extend over k-space, the 2D inverse Fourier transform of MR image.

Frequency and phase encoding technique

x y t xx 0 t yy 0 k = G (t)dt k = G (t)dt     , (1.4)

Spin echo 1.1.3

It is known that the scanner is unable to use of the free induction decay (FID), since a majority of signal has been lost as it is being received by the time. Many experiments showed that an echo can make it possible. A spin echo is used to make the out of phases spins go back to phase. Fig. 1.8 illustrates that 90 o pulse is used to make the spins phase. After t=TE/2, a 180 o pulse is applied to produce an echo, which can be repeated several times during the research in the chosen region. After repetition time TR, we reuse this method for reproducing an echo. The signal at the echo is measured for imaging.

However, for generating an image, we should determine the spatial location in the spins providing the signal. In MRI technology, we exploit magnetic gradients to encode the phase and frequency of the spins according to spatial location. Then the Fourier transform is used to convert the frequency and phase signal to a spatial image. Molecular diffusion was first described by Einstein [Einstein, 1956]. In Fig. 1.9, we show the phenomenon of diffusion with time. At all temperatures above absolute zero, thermodynamic effects cause molecules to move randomly. The motion of molecules due to thermodynamic effects is known as diffusion. If we focus on a single water molecule during diffusion, as shown in Fig. 1.10, at the beginning time 1 t0  , there are N water molecules in water. After a given time  (diffusion time interval), for each displacement distance r , The number of the water molecule that displaced to that distance as n can be counted. The diffusion is illustrated by a Gaussian function as shown in Fig. 1.11. For each displacement distance r, there is a corresponding probability n/N, which is the proportion of molecules within a voxel that distance within a time interval  .

Diffusion Weighted Imaging (DWI) Diffusion

1.2.1

According to the type of molecule, the temperature of the medium, and the time for diffusion, the distribution will be narrower or wider. The single parameter for controlling Gaussian distribution is variance ( 2 ). In diffusion here, it can be defined by two parameters: D and  , where D is the diffusion coefficient, it depends on the viscosity of the medium.  controls the diffusion time, diffusion longer, the variance is larger. This relation is described as below:

2 2D
     .

(1.5)

In the human body, the existence of the abundance of water molecules makes it possible to perform diffusion-weighted imaging using MRI.

Principles for DW-MRI 1.2.2

Diffusion weighted imaging is a form of MR imaging based on the diffusion of water molecules in a voxel. It allows the quantification of water translational motion.

As in Fig. 1.9, the random Brownian motion of water molecules has been presented. However, the problem is that the diffusion water inside a voxel of a biological tissue is hindered primarily by cell membrane boundaries. The more the cellularity, the greater the diffusion restriction is. For example, in Fig. 1.12, we show a kind of restricted diffusion that is different from free diffusion. The randomly diffusion can be described by Gaussian distribution; however, the diffusion of water molecules in heterogeneous media is in general not isotropic. The presence of ordered biological membranes can provide obstacles to the diffusing spin, restricting their movement in certain directions. [START_REF] Cleveland | Nuclear Magnetic Resonance Measurement of Skeletal Muscle: Anisotrophy of the Diffusion Coefficient of the Intracellular Water[END_REF] In diffusion imaging studies, the contrast is based on spin dephasing and the rate of diffusion can be measured by the application of diffusion-encoding gradients. [START_REF] Stejskal | Spin Diffusion Measurements: Spin Echoes in the Presence of a Time-Dependent Field Gradient[END_REF]] Therefore, we should give attention to the famous pulsed-gradient spin-echo sequence (PGSE). In Fig. 1.13, the diffusion-weighted pulse sequence is shown. It is also known as the pulsed-gradient spin-echo sequence (PGSE). The amount and orientation of diffusion weighting is controlled by varying the gradient sthength G, duration  , and the time between the onsets of the two gradient pulses  .

Below we will introduce how to meaure diffusion using PGSE in Fig. 1.14. At echo time, normally the 180 o pulse will cancel out the difference of the phase. However, if any spin moves according to diffusion phenomena, this effete of rephasing of this 180 o pulse will not be equal to the first 90 o pulse. This imperfect rephasing will create smaller single magnitude in Fig. 1.14. This single difference depends on the two diffusion gradients G.

According to diffusion MRI, many theories have been proposed for addressing the relation between the diffusion and the shape and size of the restricted volume.

Apparent diffusion coefficient (ADC) 1.2.3

As we have already said, in the case of some fluids such as intracellular water, diffusion will be restricted to a limited volume. If enough time is allowd for diffusion, the expected free diffusion distance may exceed the available range. When this phenomenon occurs, the apprent diffusion coefficient is reduced and becomes a function of the diffusion time and of the geometry of the limiting volume. [START_REF] Bihan | MR Imaging of Intravoxel Incoherent Motions: Application to Diffusion and Perfusion in Neurologic Disorders[END_REF] The apparent diffusion coefficient provides very useful information about the underlying microstructure. When water molecules diffuse without any preferential direction over the space of an image voxel, it is the pattern of isotropic. In this condition, ADC will be all the same for every direction. MR signal attenuation is given by -bD 0 S / S = e . , S is the diffusion weighted signal observed at one of the diffusion direction. 0 S is the signal at the same voxel without the additional diffusion gradients.  is gyromagnetic ratio ,  is the pulse duration, G is the strength of diffusion gradient, and  is the time between the application of the two gradients pulses. However, if restricted diffusion occurs, the ADC will be changed depending on the gradient direction and can be presented by where i S is the diffusion weighted signal at a given voxel with the diffusion-gradients applied along direction i , ADC i is the apparent diffusion coefficient of water along the i th direction of the diffusion gradients. Here the diffusion coefficient D has been replaced by an apparent diffusion coefficient ADC i to indicate that the diffusion process is not free in tissues.

Though ADC present very helpful information about the underlying microstructure, it is measured only in a single direction. For free diffusion, ADC is all the same value in each direction. Nevertheless, the measuring of ADC is depending on the choice of directions, if the scatter patter is anisotropic, for example, the underlying orientation of the white matter or the fiber structure of heart. In this situation, the ADC measure can be replaced by giving rise to the concept of anisotropy of diffusion in three directions.

Diffusion tensor imaging (DTI) 1.2.4

The model of diffusion tensor was proposed in [START_REF] Basser | MR Diffusion Tensor Spectroscopy and Imaging[END_REF]. Diffusion tensor imaging (DTI) is one of the most popular ways to measure the diffusion anisotropy. To achieve a DTI experiment, we must have one image where b=0 s/mm 2 and at least 6 diffusion weighted (DW) images to compute the whole tensor. We will obtain a better result if more than six images should be detected. where the diffusion tensor D is to be estimated.

DTI model

2 2 2 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2

Extraction of diffusion characteristic

Through diagonalizing D from equation (1.9), the eigenvalues 1  , 2  and 3  can be determined. The tensor diagonal can be present by 11 ' 22 33 00 00 00

D                 (1.12)
where 1  , 2  and 3  are the three eigenvectors of D . When the diffusion is isotropic, the three eigenvalues are equal. Conversely, when for anisotropic situation, the eigenvalues will be different and span an ellipsoid, characterizing the major diffusion direction. [Basser, 1994] For example, the major diffusion eigenvector (the direction of greatest diffusivity) is assumed to be parallel to the tract orientation [START_REF] Jing | Enhancement of Fiber Orientation Distribution Reconstruction in Diffusion-Weighted Imaging by Single Channel Blind Source Separation[END_REF] as shown in Fig. 1.15. The fractional anisotropy (FA) is a scalar value between zero and one that denotes the degree of anisotropy of anisotropy of a diffusion process. FA equaling to 0 reflects that diffusion is isotropic and the result of one reflects diffusion occurs only along one axis. FA is described mathematically in [START_REF] Basser | Microstructural and Physiological Features of Tissues Elucidated by Quantitative-Diffusion-Tensor MRI[END_REF]

      22 2 1 2 3 2 2 2 1 2 3 ( MD) MD 2 MD FA              , (1.13)
where

1 2 3 () 3 MD     
is the mean diffusivity, it presents the average diffusivity along three axis x, y and z. Therefore, it is seen as an information overall. However, there is a very clear limitation for DTI: according to ellipsoid (Gaussian assumption) presentation, it can only present a single direction.

Résumé en français

Nous avons déjà présenté dans le chapitre précèdent, l'IRMd (coefficient de diffusion apparent et l'imagerie du tenseur de diffusion) permet d'imager les propriétés de diffusion des molécules d'eau dans les tissus biologiques et donc inférer la microstructure sous-jacente. Mais il est aussi présenté que les multiple fibres ne prouvent pas montrer par les modèles dans première chapitre.

Afin de permettre une mesure angulaire précise de l'intensité du signal de diffusion IRM, l'espace-Q doit être présenté. PGSE étroit avec le temps d'application d'un gradient est très court, donc le temps d'application d'un gradient est négligeable par rapport au temps de diffusion. La théorie de l'espace Q permet de former un nombre des modèles de reconstruction.

L'imagerie du spectre de diffusion (DSI) propose une mesure directe du propagateur de diffusion en se basant de signaux de l'espace Q par Fourier transformation inverse. Bien que la méthode DSI est le plus précise méthode, il requiert un temps d'acquisition très long et il a besoin des gradients de champs magnétique très puissance.

Haute résolution angulaire (HARDI) a été proposé afin de réduire le temps d'acquisition de DSI et permettre de mesurer angulaire précise de signal de diffusion IRM. Haute résolution angulaire réduit l'échantillonnage du signal de diffusion à une unique ou plusieurs sphères dans l'espace Q, mais il peut aussi résoudre le problème de croisement de fibres. Il existe de nombreux types de modèles dans ce groupe, ici on a classé à trois groupes : l'imagerie Q-ball, déconvolution sphérique, et les autres modèles. Pour bien comprendre les modèles HARDI, modèle mixte est présenté avant, il a fait l'hypothèse que le signal de diffusion peut se décomposer en une somme pondérée de modèle de diffusion gaussienne, mais la somme des faisceaux fibres doivent être un nombre précision.

L'imagerie Q-ball a proposé pour extraire la fonction d'orientation de diffusion (ODF). Nombreux travaux sont proposé sur améliorer la qualité de ODF : l'orignal fonctions sphérique utilisé par Tuch est fonction de base radial sphérique, plusieurs travaux ont recommandé l'utilisation d'autre fonctions sphériques ou différent méthodes régularisation. Entre les méthodes proposé, l'imagerie Q-ball utilises de la base des harmonique sphériques est détaillé dans ce chapitre. Aussi la méthode Q-ball angle constant solide est proposé pour améliorer la qualité de ODF. Pour obtenir un plus précise ODF, multiple coquille sphériques méthodes est proposé.

Mais les méthodes proposés ont besoin aussi beaucoup d'acquisition de signaux diffusion, donc les méthodes pour réduire l'acquisition sur un sphérique est suggéré comme l'acquisition comprimé est utilisé dans ce domaine.

L'estimation du signal IRM de diffusion par déconvolution sphérique a été proposée. Il peut utiliser la même manière que par le DTI comme faible valeur de b value par rapport des autres méthodes HARDI. Par contre l'hypothèse que les caractéristiques de diffusion de toutes les faisceaux de fibres présents sont identiques. Déconvolution sphérique contrainte (CSD) ou super-résolue CSD est proposé pour résoudre les angles croisement plus petit.

Les autres méthodes HARDI comme DOT et GQI sont aussi présenté. A partir de l'état de l'art de ces HARDI méthodes de modélisation sur l'IRMd, on peut extraire différent mesures des anisotropies. Comme l'anisotropie fractionnelle (FA) dans DTI, l'indice d'anisotropie fractionnelle généralisée (GFA) est proposé. L'anisotropie généralisée (GA) est défini aussi pour décrire l'anisotropie de résultats des méthodes HARDI.

Pour évaluer les résultats des méthodes HARDI, plusieurs méthodes ont été proposées. Les méthodes peuvent être classées à deux catégories : la première catégorie est l'évaluation des orientations des fibres. Sur cet catégorie, les plus souvent utilisant méthodes sont erreur de angle de croisement et le taux de détection de croisement des fibres. La deuxième catégorie est l'évaluation de formes des résultats HARDI. L'évaluation qualitative comme utilisant donnée synthétique ou fantôme physique avec la vérité terrain est souvent utilisé pour comparer les formes des résultats. Donnée réel avec très bien qualité est aussi utilisé. Les méthodes quantitatives peuvent être divisées à deux groupes. La première est l'évaluation de forme globale comme RMSE ou NMSE (en anglai). Les autres méthodes comme Kullback-Leibler (sKL) et Dip test sont proposé pour évaluer les résultats comme ODF en globale. La deuxième classe est mesure quantitative, comme Jensen-shannon divergence (JSD) ou mesures d'incertitude des ODFs en Hessian ou SIP functions.

Q-space

The process of signal acquisition in the time domain with a gradient is designated k-space mapping. k-space is conjugate to the image space by the Fourier transformation. However, in Fig. 2.1, the narrow-pulse approximation PGSE was used, excepting that the gradient pulses as sufficiently narrow and the  here can be neglected. For distinguishing the magnitude of the gradient pulse of upper case symbol in k-space, g is to represent the magnitude.

      exp 2 S k r i k r dr    , ( 2
In the narrow gradient pulse approximation,

          ' ' ' -, exp - s S g r P r r i g r r dr dr       , (2.2)
where   Sg here is defined as the amplitude of the echo at its center, an echo amplitude and phase dependent on the spin self-correlation function   ' -, s P r r  which denotes the conditional probability that a spin initially at r has migrated to ' r over the time  . The difference   'rr is considered as the dynamic displacement in order to clearly distinguish it from the static displacement, r . [START_REF] Callaghan | NMR Microscopy of Dynamic Displacements: K-Space and Q-Space Imaging[END_REF]] The wavevector can be presented by

  1/ 2 qg   
, q is no longer explicitly a function of time, t, but only a function of the duration of the phase-encoding period,  .

Diffusion spectrum imaging

DSI makes use of the Fourier Transform relationship between the diffusion signal and the probability distribution function (PDF). In order to visualize the diffusion function within each pixel, DSI projects the diffusion function on the sphere by integrating over the radial coordinates of the diffusion function. The result of this integration gives the orientation distribution function (ODF). [Wedeen, 2005] The model of DSI is defined by

      00 F      iqr E q E P r e dr E P r , (2.3)
where  

Eq is the diffusion signal, q is the diffusion wave vector in the direction of the corresponding acquisition,  

Pr is the probability distribution function, F denotes the Transform Fourier.

  0 0 EE 
is the signal in the absence of the applied diffusion gradient, ' r r r  is the relative spin displacement in the experimental diffusion time t . Eq and its corresponding PDF () Pr in which, bottom row has been obtained after centering the origin. We can observe that both   Eq and  

Pr are real and positive. We can also observe that, as expected,  

Eq presents less correlation since it is by essence spectral data.

The PDF is defined by

      0 0 0 , P r P x x x dx    , (2.4) 0 r x x  , (2.5)
The diffusion PDF  

Pr is related to the measured MRI diffusion signal by the Fourier relationship

    P r F E q    .
(2.6)

The orientation distribution function (ODF)   u  is defined as the projection of the diffusion probability distribution function (PDF) as follows

    0 1 u P ru dr Z     , (2.7)
where Z is a dimensionless normalization constant, ' r r r  is the relative spin displacement in the experimental diffusion time t , u is the direction of ODF.

The ODF framework is widely used in materials science to describe the orientation composition of polymers, liquid crystals, and grain composites. As it was shown in Fig. 2.3, DSI is a method for calculating ODF based on Fourier Transform and radial integration. So it requires first sampling in coordinate Cartesian and then projects the data into coordinate spherical.

Because DSI requires a lot of data, so it requires a long acquisition time and the magnetic field strength, which is not easy to achieve in practical applications. For the collection of DSI data, it must collect the data to the three-dimensional Cartesian coordinate system, which caused some difficulties in computing. Extracting the ODF from the PDF requires explicitly calculating the radial projection. The mapping between Cartesian and spherical coordinates systems may introduce Cartesian artifacts in the ODF. Cartesian coloration of the ODF may be a particular problem at the coarse Cartesian resolution typically used for DSI. Further, the radial projection is highly inefficient since the projection discards a considerable fraction of the acquired data. The efficiency of DSI is also hampered by the strong pulsed field gradients needed to satisfy the Nyquist condition for the diffusion PDF in cerebral white matter. [DS Tuch, 2004]

High angular resolution diffusion imaging (HARDI)

To deal with the sampling problem of QSI, researchers have proposed an approach based on sampling on a spherical shell (or combination of shells) in diffusion wavevector spaces. High b-values were used because at the lower b-values employed by DTI there is not sufficient contrast between the fast-diffusion component of one fiber and the slow-diffusion component of another fiber to effectively resolve the two fibers.

By high angular resolution, high b-value diffusion gradient sampling, we were capable of detecting diffusion signals with multiple, discrete maxima/minima as a function of gradient orientation, demonstrating the presence of multiple main fiber populations. [ DS Tuch, 2002] In [ DS Tuch, 2002], the Gaussian mixture formulation is proposed,

  exp( b ) D   T i k i k i k S g P g g , (2.8)
where k P is the apparent volume fraction of the voxel with diffusion tensor D k , b is diffusion sensitization factor, i g is the direction of the diffusion gradient [ DS Tuch, 2002;Kingsley, 2006].

However, this method depends on the number of fiber populations present constantly in each voxel estimated. In addition, this technique is unstable if more than two populations are present.

Q-ball Imaging

2.3.1

In order to avoid the coordinate conversion in DSI and its demand on high magnetic field strength, QBI directly performs acquisitions on a sphere, and combines such acquisitions with Funk-Radon transform, the corresponding ODF is calculated using

          , T S f x o q f x x o x q dx      . (2.9)
The QBI reconstruction is based on the FRT, also known as the spherical Radon transform or simply the Funk transform. The FRT is an extension of the planar Radon transform to the sphere. The FRT is a transform from the sphere to the sphere. Given a function on the sphere   fx, where

x is a unit direction vector, the FRT is defined as the sum over the corresponding equator, the set of points perpendicular to o . The FRT S for a direction o can be written as (2.9). It is shown in Fig. 2.4, this remarkable relationship entails that the sum of the diffusion signal over an equator approximately gives the diffusion probability in the direction normal to the plane of the equator. Consequently, to estimate the diffusion probability in a particular direction all that is needed is to sum the diffusion signal along the equator around that direction. This provides a model-free approach for estimating the diffusion probability from the spherically sampled diffusion signal.[ [START_REF] Ds Tuch | [END_REF] After collecting the diffusion signals in q-space on a sphere, several approaches can be used to reconstruct ODF. Tuch proposed a method based on Funk-Radon transform (FRT). Spherical radial basis functions (RBFs) were used as the basis in this approach. The method of calculating ODF according to FRT transform is the first linear model for calculating ODF, which establishes a direct relation between the diffusion signals acquired on the sphere in qspace and the ODF coefficient. But this model requires an important amount of computation and it is also sensitive to noise.

It samples the signals directly on a sphere. So it doesn't need conversion between the coordinate Cartesian and Spherical.

Related work about

Q-ball imaging 2.3.2
In the field of QBI, most works reported in the literature concern the research of improving the quality of orientation distribution function (ODF): some work is focused in applying new spherical basis or regularization method [ [START_REF] Ds Tuch | [END_REF]Descoteaux, 2006Descoteaux, , 2007;;[START_REF] Hess | [END_REF][START_REF] Kezele | Spherical Wavelet Transform for ODF Sharpening[END_REF]Michailovich, 2010c].

Spherical radial basis method is the original QBI method [ [START_REF] Ds Tuch | [END_REF] for the decomposition of diffusion signal, it can provide angular architecture in each voxel. Spherical wavelet based q-ball imaging is efficiency in diffusion ODF denoising and sharpening [START_REF] Kezele | Spherical Wavelet Transform for ODF Sharpening[END_REF]. The basis of spherical ridgelets has better de-noising properties with respect to HARDI signals, and as a result it is capable of better representing the directional information. [Michailovich, 2010c] Though using spherical ridgelet or spherical wavelet, the ODF quality can be improved, the spherical harmonic basis is still utilized most in q-ball imaging. Its accuracy and speed is generally as good as methods based on other basis, the more important factors are two: one is simpler to construct and compute using existing implementation; the other is Funk-Hecke theorem is used to simplify the procedure easily [Descoteaux, 2006;[START_REF] Hess | [END_REF]. Therefore, in this work, we still use spherical harmonic basis based q-ball imaging as the method for reconstructing ODFs.

Let () Eq designate the magnetic resonance (MR) diffusion signal acquired in q-space.

() Eq is in general acquired in a great number of directions (more than 60 directions).

Once the above diffusion signal defined in N directions on a sphere is obtained, we project it on a spherical harmonic basis

B of size NK          1 1 1 1 1 1 ,, ,, B      K N N K N N YY YY         , (2.10) with    1 2 / 2 K l l   
, where l represents the order of spherical harmonics. That is

[11]     1 , E    K difsig kk k q c Y  , (2.11) where   1 ,... c  T difsig K cc is a 1 K  vector. The weighting coefficients are calculated using     1 c B B L B E   difsig T T q  , (2.12)
where the KK  diagonal matrix L is a Laplace-Beltrami operator, in which each matrix element is defined as

      2 2 1 l k l k 
, and  is a regularization factor that controls the smoothing effect of the estimator.

The ODF coefficients c ODF are finally computed using Funk-Radon transform and Funk- Hecke theorem. They are given by

c Pc  ODF difsig , (2.13)
where P is a KK  Funk-Hecke diagonal matrix, whose element is defined as

() 2 (0) kk l k pP   , with () (0) lk P
designating an () lk order of Legendre function. After obtaining ODF c , we can calculate the value of the ODF in any direction using

        1 , K ODF kk k o c Y o o       , (2.14)
where o designates the direction, in which ODF is reconstructed.

There are also some work about improving q-ball imaging with new acquisition method or mathematical model, some of them suggest to exploit data from multiple q-shells [START_REF] Descoteaux | Multiple Q-Shell Diffusion Propagator Imaging[END_REF], some create ODF expression of QBI with definition constant solid angle [START_REF] Aganj | Multiple Q-Shell ODF Reconstruction in Q-Ball Imaging[END_REF]Aganj, 2010].

Spherical deconvolution 2.3.3

Spherical deconvolution method was proposed firstly for estimating directly the distribution of fiber orientations within a voxel from high angular resolution diffusion data in [Tournier, 2004].

  ( , )

i i i i i S f C R      , (2.15)
where ( , )

S  is the diffusion signal, i f is the volume fraction for the i th fiber population, and i C is the operator representing a rotation onto the direction  

, ii  . The ideal can be expressed as the convolution over the unit sphere of the response function   R  with a fiber

ODF (fODF)   , F  , as shown in (2.16),     ( , ) , ii i S F R        . (2.16)
The major limitation of spherical deconvolution is its susceptibility to noise, which often results in spurious peaks in the recovered fODF. The original method uses low-pass filtering for reducing spurious peaks, however, low-pass filtering also reduces angular resolutions.

Constrained spherical deconvolution (CSD) or super-resolved CSD [Tournier, 2008] works well even using relatively low b-values, and can neglect biased estimates of fiber orientations down to relatively small crossing angles.

Another limitation is the assumption that the diffusion characteristics of all fiber populations found in the brain are identical.

Other HARDI methods

2.3.4

Diffusion orientation transform (DOT) proposes a method for fiber orientation mapping. This method uses the Fourier transform relationship between the water displacement probabilities and diffusion attenuated MR signal expressed in spherical coordinates. The radial part of the Fourier integral is evaluated analytically under the assumption that MR signal attenuates exponentially. [Özarslan, 2006] Generalized q-sampling imaging (GQI) has been proposed in [Yeh et al., 2010] based on the Fourier transform relation between diffusion MR signals and the underlying diffusion displacement, a new relation is derived to estimate the spin distribution function (SDF) directly from diffusion MR signals.

Anisotropy measures of HARDI

The definition for the generalized anisotropy [START_REF] Ozarslan | Generalized Scalar Measures for Diffusion MRI Using Trace, Variance, and Entropy[END_REF],

(V) 1 1 1 (250 ) GA V    , (2.17)
where the exponent  

V  is defined as   1 1 1 5000 V V   
.

(2.18)

This measure was presented for scalar measures derived in this paper with the diffusivity profile. However the formulations also remain valid for other applications such as q-space imaging in which one calculates the average propagator for water molecules in an orientation distribution function is constructed as below, 2 2 2 (i)

9

(1)

N ODF i ODF c V c     , (2.19)
where ODF coefficients ODF c are computed using Funk-Radon transform and Funk-Hecke theorem [START_REF] Hess | [END_REF]Descoteaux, 2007].

GFA is a measure introduced in [DS Tuch, 2004],as an extension of the FA metric, the definition of GFA,

        2 1 2 1 () 1 ( ) n i i n i i nu std GFA rms nu            , (2.20)
where  is the mean of the ODF. Like the FA metric for DTI, the GFA metric is automatically normalized to [0, 1].

Assessment of HARDI

In [START_REF] Hess | [END_REF], the author presents that the quality of ODF can be divided into two parts: the first one is that global effects of reconstruction, such as Root Mean Square Error (RMSE), the second one is the orientation of individual peaks in ODF. It presents also in [Descoteaux, 2007] that the comparison of HARDI results can be from two parts, one is that ODF shape comparison, the other is that fiber detection and angular resolution. Also in [Seunarine, 2007], the shape of the peaks reflects the local spread of fibers and in particular, the peak anisotropy reflects local fiber bending or fanning. The shapes of the peaks of ODF by HARDI methods can provide not only information discrete fiber orientation, but also the shape of the orientation distribution function.

According to the research before, we can find that the assessment of functions of HARDI results should be from two parts, one is the information of main orientation, and the other one is the information morphological of functions, such as ODFs.

Fiber orientation assessment 2.5.1

Angular orientation error

The angular orientation error is defined as

0 180 arccos( )  ee   , (2.21)
where 0 e is the true direction and e its estimate.

Angular error  can be also defined (in degrees) based on Earth Mover's Distance (EMD) [START_REF] Rivera | Improved Diffusion Basis Functions Fitting and Metric Distance for Brain Axon Fiber Estimation[END_REF].

Most of the works used the angular orientation error as the quality criterion for comparing different reconstruction methods. Angular error is used in [START_REF] Hess | [END_REF]Descoteaux, 2007;Tournier, 2008;Michailovich, 2011]. Angular average and standard deviation errors [Descoteaux, 2007;Tournier, 2008]can be also seen as this kind of assesement.

Detection rates

The correct detection rates is defined as

  Cart( ) / Cart( ) , | (r) N(r) cc rN      , (2.22)
where r is number of the voxel,  is a set of all the fibers, c  is a set that estimated number of fiber N(r) equal to true number (r) N of fiber tracts passing through voxel r, Card() designates the cardinal number.

The percentage of success [Descoteaux, 2007;Tournier, 2008], and probability of false crossing fiber detection [Michailovich, 2011]can be also classified into detection rates.

However, all of the assessment metrics above are about discrete fiber orientation. For assessing the quality of ODF or other functions of HARDI results, the assessment methods of shape are also proposed.

Shape of HARDI results assessment 2.5.2

Qualitative assessment

Synthetic data form mixtures of Gaussians [ DS Tuch, 2002] or some other models of simulations [START_REF] Szafer | Theoretical Model for Water Diffusion in Tissues[END_REF][START_REF] Pfeuffer | Restricted Diffusion and Exchange of Intracellular Water: Theoretical Modelling and Diffusion Time Dependence of 1H NMR Measurements on Perfused Glial Cells[END_REF]] can be used as a ground truth for comparing with the quality of results obtained as fiber directions and shapes information.

In [START_REF] Descoteaux | Deterministic and Probabilistic Tractography Based on Complex Fibre Orientation Distributions[END_REF], the effect of varying the HARDI signal parameters on the ODF is qualitatively compared with synthetic simulations. However, sometimes the models of simulations simplify the biological tissues. The synthetic simulation results always compared with observations or visual results. In this paper, the effect of sharpening, the shape of peaks and the overall difference between different results are choosing as criterions qualitative.

The physical phantoms are also used for comparing different reconstruction methods or tractography methods as [START_REF] Jeurissen | Fiber Tracking on the 'Fiber Cup Phantom'using Constrained Spherical Deconvolution[END_REF][START_REF] Kezele | Spherical Wavelet Transform for ODF Sharpening[END_REF][START_REF] Fillard | Quantitative Evaluation of 10 Tractography Algorithms on a Realistic Diffusion MR Phantom[END_REF]Reisert et al., 2011]. According to the ground truth provided by physical phantoms, the methods can be clearly compared. High quality real data was also exploited as a ground truth for researching different methods. [Özarslan, 2006;Aganj, 2010;Michailovich, 2011;[START_REF] Welsh | Model-Based Reconstruction of Undersampled Diffusion Tensor K-Space Data[END_REF] For comparing different q-ball based methods, the high quality real data and synthetic simulation are used both in [ [START_REF] Ds Tuch | [END_REF][START_REF] Hess | [END_REF]Descoteaux, 2007].

Quantitative assessment

The mean square error [Descoteaux, 2007], root-mean square error (RMSE) [START_REF] Hess | [END_REF], normalized mean squared error (NMSE) [Michailovich, 2011] between the noisy ODF and the noise free ODF. These metrics are the greater, noisy ODFs is more different from noise free ODF. In [Zhan], a more global description of the ODF, called Kullback-Leibler (sKL) method coming from information theory, was proposed that uses gold standard ODFs as ground-truth to assess how accurately the diffusion profile can be reconstructed from subsampled data based on different angular sampling schemes. To assess how accurately the diffusion profiles can be reconstructed from subsampled data based on different angular sampling schemes, sKL measure the discrepancy between the reconstructed and ground truth ODFs.

              1 ( ) log / log / 2 sKL p,q p x p x q x q x q x p x dx            , (2.23)
where   px is the ODF obtained from the subsampled schemes,   qx is the noise-free ODF obtained from whole data set.

However, in real data cases, the ground-truth is difficult to obtain.

In [Aganj, 2010], the authors used an objective criterion from statistics, called Dip test [START_REF] Hartigan | The Dip Test of Unimodality[END_REF] that estimates the maximum distance between the empirical distribution function and the closest unimodal distribution function. The Dip test method does not require ground-truth ODFs. However, the Dip test method requires sampling a direction from the ODFs to be able to perform the Dip test in a one-or two-dimensional (2D) space, and such sampling should be appropriate.

The bootstrap method is proposed to assess the quality of HARDI data, as in [Cohen-Adad, 2011], the Jensen-Shannon divergence (JSD) and the angular confidence interval was derived for the first and the second ODF maxima. In [START_REF] Jiao | [END_REF], the uncertainty visualization based on ODF has been proposed. SIP functions of diffusion shapes have been researched for analyzing the underlying uncertainty. In [Seunarine, 2007], the Hessian was used for describing each peak, the shape of the peaks reflects the local spread of fibers.

Résumé en français

Les nuages tridimensionnels (3D) des points tels que la fonction d'orientation de distribution (ODF) sont largement utilisés pour résoudre les problèmes de croisement de fibres en l'imagerie de diffusion à haute résolution angulaire (high angular resolution diffusion imaging-HARDI). L'évaluation des caractéristiques ou de la qualité dans une certaine mesure des nuages de points 3D s'effectue le plus souvent en utilisant le critère visuel, bien que l'utilisation de quelques mesures objectives de la qualité soit également adressée, qui sont directement empruntées de la théorie classique de traitement du signal et de l'image, parce qu'elles sont simples et intuitives à calculer. Cependant, ces mesures ne sont pas toujours pertinentes pour caractériser les nuages de points 3D. Nous proposons un nouveau paradigme pour évaluer les caractéristiques des nuages de points 3D, basé sur la mesure des caractéristiques morphologiques de ces derniers. L'idée consiste à projeter le nuage de points 3D sur un plan angle-distance (ADM), construire une matrice angle-distance (ADMAT), et calculer des paramètres tels que le rapport de longueur, la séparabilité et l'incertitude. Trois méthodes HARDI sont comparées en utilisant les mesures proposées. Les résultats obtenus ont montré que les caractéristiques des nuages de points 3D peuvent être évaluées d'une manière relativement complète et quantitative.

Introduction

The orientation distribution function (ODF) [Wedeen, 2005] is a quantity used to describe the orientation architecture of the tissue's fibers or fiber bundles; it gives the probability of diffusion in different directions. ODF is often estimated or reconstructed from high angular resolution diffusion imaging (HARDI) such as q-ball imaging (QBI) [ [START_REF] Ds Tuch | [END_REF]] using spherical sampling. In this field, most existing works put emphasis on improving the quality of ODF using normalization and regularization [ [START_REF] Ds Tuch | [END_REF]Descoteaux, 2007], change of basis [START_REF] Hess | [END_REF][START_REF] Kezele | Spherical Wavelet Transform for ODF Sharpening[END_REF]Michailovich, 2010c], sharping deconvolution [Descoteaux, 2009], compressed sensing [Michailovich, 2011], etc. Meanwhile, other quantities were also used to describe fiber orientation or crossing, including the fiber orientation distribution (FOD) from the spherical deconvolution method [Tournier, 2004[Tournier, , 2008]], the orientation map derived from the diffusion orientation transform (DOT) based on the Fourier transform relationship between water displacement probability and diffusionattenuated magnetic resonance (MR) signal expressed in spherical coordinates (Ozarslan., 2006), and the water molecule displacement probability function [START_REF] Jian | A Novel Tensor Distribution Model for the Diffusion-Weighted MR Signal[END_REF] using the mixture of Wisharts. The above HARDI quantities share a common point, namely, they are all composed of N samplings of the sphere (N 3D points), where the direction of each sampling corresponds to one reconstruction direction of HARDI and the distance to origin of each sampling corresponds to the probability density function of water diffusion along the direction. We regard all these quantities as three-dimensional (3D) point clouds as shown in Fig. 3.1, in which each point corresponds to a vector originating from the coordination system origin. The direction and the length of the vector designate the direction and the amount of the diffusion, respectively. The difference between the 3D point cloud on the left of Fig. 3.1 and the usual ODF visualization is that the latter is a 3D surface representation of the former. In the past, a number of methods have been reported to assess the quality of 3D point clouds in HARDI. The most important metric is the minimum resolvable fiber crossing angle or angular resolution. Other angular metrics include the angular error [START_REF] Hess | [END_REF]Descoteaux, 2007;Tournier, 2008;Michailovich, 2011], percentage of success [Descoteaux, 2007;Tournier, 2008], angular average and standard deviation errors [Descoteaux, 2007;Tournier, 2008], probability of false fiber detection [START_REF] Dolui | Compressed Sensing of Diffusion MRI Data Using Spatial Regularization and Positivity Constraints[END_REF].

In addition to this basic angular information, other characteristics of 3D point clouds are also pertinent information, such as the shape and sharpness of the 3D point cloud's peaks, which reflect the local spread of fibers, and the peak anisotropy (indicating how elliptical the peak cross section) which reflects local fiber bending or fanning [Seunarine, 2007;Descoteaux, 2009]. However, the model proposed in [Seunarine, 2007] can only model the uncertainty in fiber orientation, the peak height is not used. Though multitensor reconstruction [ DS Tuch, 2002] can obtain the volume fractions of diffusion tensors, the question of how many diffusion compartments are present. Rotation invariant features describe the shape of ODFs and are applicable to any reconstruction method that represents HARDI signals in terms of an spherical harmonic basis [START_REF] Schwab | Rotation Invariant Features for HARDI[END_REF]. However, the reconstruction method that represents HARDI signals in terms of other spherical basis is also proposed in [START_REF] Kezele | Spherical Wavelet Transform for ODF Sharpening[END_REF]Michailovich, 2010cMichailovich, , 2011]].

The shapes of ODFs have been compared by some quantitative measures. The mean square error [Descoteaux, 2007], root-mean square error (RMSE) [START_REF] Hess | [END_REF], normalized mean squared error (NMSE) [Michailovich, 2011] between the noisy ODF and the noise free ODF. These metrics are the greater, noisy ODFs is more different from noise free ODF. A widely adopted assumption of these metrics is that the loss of perceptual quality is directly related to the visibility of ODF. In [Zhan] , a more complete description of the ODF, called Kullback-Leibler (sKL) method coming from information theory, was proposed that uses gold standard ODFs as ground-truth to assess how accurately the diffusion profile can be reconstructed from sub-sampled data based on different angular sampling schemes. The sKL method allows for the measurement of the discrepancy between the reconstructed and ground-truth ODFs. However, in real data cases, the ground-truth is difficult to obtain. Moreover, the abovementioned metrics for comparison of the global shape of 3D point cloud like MSE, RMSE, sKL and NMSE share a common point, namely, they are all based on computing the deviation from a ground truth 3D point cloud. They regard 3D point cloud as a set which is made up only of isolated points. Then compute the distances such as Euclidean distance or Kullback-Leibler distance between each point of ground truth and estimated 3D point cloud. In [Aganj, 2010], the authors used an objective criterion from statistics, called Dip test [Hartigan, 1985] that estimates the maximum distance between the empirical distribution function and the closest unimodal distribution function, to compare the quality of crossing fibers in various ODFs, FODs, orientation maps, and water molecule displacement probability functions. Like the sKL method, the Dip test method can quantify the quality of ODFs (or other quantities) showing the same crossing fiber. But, unlike the sKL method, the Dip test method does not require ground-truth ODFs. However, the Dip test method requires sampling a direction from the ODFs to be able to perform the Dip test in a one-or two-dimensional space, and such sampling should be appropriate.

Another type of metrics, rather different from the above ones, consists in measuring the uncertainty of ODFs or ODF parameters instead of describing the ODF itself [Cohen-Adad, 2011;[START_REF] Jiao | [END_REF]; these metrics are essentially based on the use of bootstrap analysis involving the resampling of originally acquired diffusion-weighted datasets.

The relationship between estimated ODF and actual fiber ODF is complex and unclear. [Heidi Johansen-Berg, 2009] The comparison of different ODFs should be researched for understanding this relationship more clearly. Comparison of characteristics of the ODF peaks could provide clues to distinguish the underlying fiber structure intravoxel, for example a crossing fiber arrangement can be distinguished from a kissing fiber pattern based on characteristics of the ODF peaks.[ DS Tuch et al., 2003] Therefore, we propose a novel paradigm to quantitatively and accurately assess ODF or other quantities by measuring their morphological characteristics. Our goal is then to try to gain insights into the spatial structure and shape of the 3D point in a more general and systematic manner. The motivation of our new approach is to extract reliable quantitative morphological characteristics of HARDI quantities. Our morphological characteristics is applied to find a more general and direct way to compare the shape of all HARDI quantities under noisy conditions and compare the morphological characteristics of different HARDI methods.

The idea consists in projecting a 3D point cloud onto an angle-distance map (ADM), constructing an angle-distance matrix (ADMAT), and calculating the metrics such as the length ratio, separability and uncertainty. We used our morphological metrics on simulated data for comparing with other metrics of comparing shapes of ODFs. Then we also compare three HARDI methods using the proposed metrics on both simulated and real data.

3D point cloud characteristic assessment paradigm

The flow chart for calculating the proposed morphological metrics is illustrated in Fig. 3.2. The 3D surface representation of ODFs was addressed in [ [START_REF] Ds Tuch | [END_REF]Wedeen, 2005;Descoteaux, 2007], with a variant reported in [START_REF] Aganj | Multiple Q-Shell ODF Reconstruction in Q-Ball Imaging[END_REF][Aganj, , 2010]]. It is often expressed in spherical coordinates. To better visualize and handle ODFs or more generally 3D point clouds, we express them in Cartesian coordinates by ( , , ) T q x y z qe  . Then, the 3D point cloud can be described by a set of vectors. A vector being defined by radial distance and orientation, we introduce these two parameters to characterize the 3D point clouds. If m q represents the maximal distance with m e as its direction, we will then take m q as a reference vector. From the reference vector, we construct an angle-distance map (ADM) formed of small areas . We now project the given vectors onto this ADM. To do this, we first calculate the angle between each vector q and the reference vector m q . Since the distance q of q to the origin is known (equal to its length), we then have the values of the angle and distance, which enables us to put the vectors in the corresponding area of ADM (Fig. 3.3). If the maximum distance appears in two or more directions, we can choose any of them as the reference direction. Once all the ODF points are projected on the ADM, we analyze and characterize the distribution of the projections. To this end, we first define three metrics (whose mathematical expressions will be given in Section 3.2.3): The length ratio that describes the main direction diffusivity, the separability that reflects the 3D point cloud's ability to separate main directions, and the uncertainty that indicates the width of the 3D point cloud's tine or peak. The definition of these metrics is illustrated in Fig. 3.4.(a). The values of the length ratio range from 0 to 1, those of the separability from 0 to 1 and those of the uncertainty from   0.17 10 / 180  to 1.57 90 () /180



. If there are more than 2 fibers intravoxel, we will sum the values of the length ratio, separability of every 2 fibers and sum of uncertainty of each fiber, for example, as shown in Fig. 3.8. The greater the length ratio is, the closer the lengths between the two fibers are. The greater the separability, the more the fibers can be easily separated. The smaller the uncertainty, the thinner the peaks of the 3D point cloud.

To calculate the above metrics, we will construct in the following the angle-distance matrix (ADMAT). In Fig. 3.4 (b) we illustrate an ADMAT and its relationship with the 3D point cloud. We project the 3D point cloud in Fig. 3.4 (a) onto the 2D ADMAT schematized (the details will be given in Section 3.2.2) in Fig. 3.4 (b), where the schematized ADMAT presents two rectangular bars. The distance between the two bars represent to some extend the crossing angle. Then, from the so obtained 2D map, we calculate the three metrics directly. The length ratio is defined as the ratio of shorter bar length to longer bar length using min(Length 1 ,Length 2 )/max(Length 1 ,Length 2 ), the separability as the depth from the bottom to the top of the shorter bar, and the uncertainty as the sum of the widths of the two bars divided by two.

To construct the ADMAT, we now determine how a 3D point is projected in the corresponding small area of the ADM. Since the crossing angle ranges from 0 to 90°, a 3D point is then projected in the corresponding annulus according to

    ) 1 / / / c m c d(i q i N q q i N     (3.1)
which gives the set ) d(i of the 3D points falling inside the i th annulus. At the same time, we consider the projection of the same 3D point in the fan sector using

    90( 1) / ( / ) 90 / a m m a ang j q j N arccos q q q q j N      (3.2)
which gives the set

 

ang j of the 3D points falling inside the j th fan sector. So, the set   , angd i j of 3D points falling into the small area (i, j) of the ADM is given by

      , 1 / / / ,90( 1) / ( / ) 90 / c m c a m m a angd i j q i N q q i N j N arccos q q q q j N         . (3.3)
We can count the number of the projections of 3D points on the area  

, ij using     , ( , ) N i j Card angd i j 
, where Card() designates the cardinal number. We then construct the ADMAT by calculating each of the elements in the area as

  , 1 , i i j p p M c N p j    , (3.4)
where p c denotes a coefficient which weights the influence of distance on the shape of a 3D point cloud. Because of accumulation effect in (3.4), the last line of ADMAT always has the greatest value. Equation (3.4) expresses a general idea of weighting cloud points according to their distance from the origin. Since the points near the origin do not contribute much to the shape of the 3D point cloud, we will attenuate their influence using  

( 1) , 1

,

i p ij p M a N p j     , (3.5)
where 1 a  is a constant. In the present study, we chose 10. a 

In ADMAT, the column number and the row number indicate the angular interval and the distance interval, respectively (Fig. 3.5). So, ADMAT encodes the distance and angle information of 3D point clouds. In Fig. 3.6 we give some typical examples of 3D point clouds visualized as color 3D surface, the corresponding ADMAT and the schematized ADMAT. For the anisotropic orthogonal fiber crossing cloud in Fig. 3.6 (a), ADMAT presents two regions separated by a white region, and full-zero columns or columns having very small values between two columns. For the spherical isotropic cloud shown in Fig. 3.6 (b), which represents free diffusion situation, each element of the schematized ADMAT appears as gray, there do not exist white regions, and the whole ADMAT is gray. For the anisotropic one fiber cloud shown in Fig. 3.6 (c), ADMAT presents a gray region near its first column, and the other regions are white.

It is to be underlined that ADMAT encodes richer information than 3D surface representation. Fig. 3.7 aims to illustrate such fact, in which the fiber crossing is difficult to assess in the 3D color surface representation of the point clouds, but it is clearly indicated by the schematized ADMAT due to the presence of a white region in the middle of the first lines; the schematized ADMAT indicates that the separability is very small in this 3D point cloud and the uncertainty of two fiber crossing is very large.

For 3-fiber crossing, we divide 3D point cloud into 3 groups of 3D points, each group contains one local maxima [Reisert, 2011] of the estimated ODFs and the neighbor 3D points around this local maxima. We use the methods of 2-fiber system above for analyzing each pair. As shown in Fig. 3.8.

        , i, ,1 ,2 ,3 , , i, ,1 ,2 ,3 , , 0 , , ...
, 0 , , ...

k x j k k k k l k x j k k k k g m j i k M m m m m n j i k M n n n n             
        , , 1 , ,1 ,2 ,3 , , , 1 , ,1 ,2 ,3 , 1 , , ... 1 ,
, ...

k p k x k x k k k k ls k p k x k x k k k k gs r x m m r r r r v x n n v v v v           ,1 a pN  (3.7)
where ls and gs designate the number of discontinuities in ,

    , 1 2 3 / , , ... c k p ns c k N v ns N           (3.8)
where ns is the number of lines whose , kp v is not vide.

We then define an operator () fset that takes a discontinuity as input and outputs a set of positions of zeros or non-zeros in the k th line of ADMAT

    ,1 , 1 , 1 , 1 , ,1 ,2 ,3 , , , , 1 , 2 , , , ...
) cc N lw N lw Card fset v Nv    , 1 c N   ( ( )
is the line number of ADMAT at which fiber crossing is resolved for the first time. In this case, we can further determine how well the two directions can be resolved. To do that, we now compute the above-mentioned morphological metrics: length ratio

Experiments and results

To assess the effectiveness of the proposed metrics (length ratio, separability, and uncertainty), both simulated and real 3D point clouds were used. In the case of simulated data, both noise-free and noisy (Rician noise, with different signal-to-noise ratios) data were generated. The assessment was performed in various conditions with different crossing angles (90 o , 70 o , 50 o , 30 o ) and simulation conditions (different b-values and spherical sampling densities or number of diffusion gradients ND).

Simulation results

3.3.1

The proposed morphological metrics have been tested on different 3D point clouds corresponding to different configurations of one or two fibers. For a given fiber configuration, the following three methods have been used: analytical q-ball imaging (AQBI) [Descoteaux, 2007], q-ball imaging with constant solid angle (CSA) [Aganj, 2010], and constrained spherical deconvolution (CSD) [Tournier, 2004[Tournier, , 2008]]. They generated three types of 3D point clouds we call AQBI-group (AQBI), CSA-group (CSA) and CSD-group (CSD), respectively. All the three groups were obtained using: spherical harmonic order max 8 l  , number of reconstruction points (N)=162 uniform angular distribution on a sphere, and diffusion encoding directions=81 uniform directions on a hemi-sphere. The used data concerns two fibers crossed at 90°, 70° 50°, and 30°. The simulated image is of size 10x10 pixels, and is divided into 4 regions: the upper left region corresponds to the fibers crossed at 50°, the upper right region at 70°, the lower left region at 90°, and the lower right region at 30°. Each pixel presents two main intravoxel crossed directions. The diffusion signal was simulated using the multi-tensor model  

exp T i k i k i k S g = P (-bg g )  D
, where k P is the apparent volume fraction of the voxel with diffusion tensor k D , b is the diffusion sensitization factor, and i g is the direction of the diffusion gradient [ DS Tuch, 2002]. To generate different crossing angles, we rotated the following diffusion tensor [Kingsley, 2006]. , and

( 1)

10   i p c .
In the simulation, Rician noise was used. The signal-to-noise ratio (SNR) is defined as Length ratio in Fig. 3.9 (c) 0.9 0.8 0.7 0.6 0.5 0.4

    2 1 10 2 1 20 log           K kk k K k k S SN SNR S , ( 3
In Fig. 3.9 (a), three 3D point clouds presenting a single main direction but different uncertainties are shown. In Fig. 3.9 (b) and Fig. 3.9 (c), we show different 3D point clouds which all have the same angle of maximum directions. Visually, the 3D point clouds corresponding to fiber crossing are not the same. Table 3.1 gives the quantitative measurements of this visual difference. Fig. 3.9 (b) shows the separability levels from 0.1 to 0.6 and Fig. 3.9 (c) the different length ratio levels.

In the condition of b-values=3000 s/mm 2 and spherical sampling densities (ND=81), for a 90 o crossing voxel, the analytical q-ball imaging (AQBI) was used to generate a noise free ODF and noisy ODFs at SNR=3. We compare the noise free ODF and noisy ODFs by MSE, sKL, RMSE, NMSE and our morphological metrics. NB: X denotes that no peaks are separable in the 3D point cloud.

In Fig. 3.10, the 3D point cloud of 90 o crossing is shown as ground truth for the comparison of the shape of 3D point cloud. Ten different noisy 3D point clouds are generated by adding the Rician noise to the original diffusion signal of 90 o crossing. Regarding all of the metrics, Voxel 5 is always the best. We observe that Voxel 1, 3, 4, 6, 7 and 10 shows a clearer fiber crossing than Voxel 2, 8 and 9.

However, Table 3.2 gives different quality sorting results of the 7 metrics. In terms of MSE, sKL and RMSE, Voxel 2 (3.2, 0.07 and 5.7) and Voxel 9 (3.4, 0.11 and 5.8) are regarded as high quality, because results of these metrics are smaller than Voxel 3 (4.8, 0.26 and 6.9), Voxel 4 (6.9, 0.27 and 8.3), Voxel 6 (4.1, 0.16 and 6.4) and Voxel 10 (5.1, 0.18 and 7.1) that show clear fiber crossing. In terms of NMSE, Voxel 1 (0.2) is the worst one of ten voxels, but it is clearly better than Voxel 2 (0.17). Our length ratio, separability and uncertainty metrics show Voxel 3 (0.8, 0.4 and 0.69), Voxel 4 (0.9, 0.5 and 0.69), Voxel 6 (0.8, 0.3 and 0.87) and Voxel 10 are better than Voxel 2 (0, 0 and 1.57) and Voxel 9 (0, 0, 1.57). Each morphological metric sorted the voxels in the different orders of quality in different morphological characteristics. In Fig. 3.11 are shown three noise-free 3D point cloud maps obtained using respectively three different reconstruction methods [Tournier, 2004[Tournier, , 2008;;Descoteaux, 2007;Aganj, 2010] and with b=5000s/mm 2 and ND=81. Fig. 3.12 visualizes different 3D noise-free point clouds obtained using three reconstruction methods and in different conditions (for a fixed b=3000s/mm 2 ) and each of them corresponds to a combination of spherical sampling density (ND=81 or 321) and fiber crossing angle (90 o , 70 o , 50 o and 30 o ). To get more insights into the crossing structure of the 3D point clouds corresponding to different b-values and fiber crossing angles, we give in Table 3.3 their quantitative measurements in terms of separability and uncertainty. NB: X denotes that no peaks are separable in the 3D point cloud.

In terms of separability, the separabilities of the three methods are reduced when passing from higher b-values (b=3000 s/mm 2 or 5000 s/mm 2 ) to the smaller b value (b=1000 s/mm 2 ). The smaller the b value, the smaller the separability is. Hence, the b value plays a key role in separating fiber crossings. At smaller crossing angles (50 o and 30 o ), the three methods presented smaller separabilites than at larger crossing angles (90 o and 70 o ). The smaller the crossing angle, the smaller separability is. When increasing the spherical sampling density N, AQBI-group presented more separability, and other two methods presented nearly the same separability results as before. The separability of AQBI-group is more dependent on the number of diffusion gradients. In the same experimental condition, CSD-group exhibited more separability than other two methods, except for the very small crossing angle (30 o ). CSA-group showed more separabilities at 30 o crossing.

In terms of uncertainty, increasing b value reduced uncertainty for these three methods. When the number of diffusion gradients ND increased, the uncertainty reduced. Among these three methods, CSA-group and CSD-group nearly shared the same uncertainties, and AQBIgroup presented the greatest uncertainty.

Length ratio=0.7

Length ratio=0.4 In Fig. 3.13 we show the 3D point clouds generated by combining two length ratios (0.7 and 0.4) and four crossing angles (90 o , 70 o , 50 o , and 30 o ). Visually, the 3D point clouds are different. To further quantify the difference, we give in Table III the measured length ratios of the three groups of point clouds in different crossing angle cases, with respect to the two ground truths (0.7 and 0.4). Our experiments showed that when reducing the length ratio, for smaller length ratios, fiber crossing is more difficult to resolve as can be seen in Fig. 3.13. For the same ground truth of length ratio, the smaller the crossing angle, the smaller the obtained length ratio. AQBI was the most sensitive to the reduction of length ratios. CSA presented greater length ratios than the ground truths. The length ratios of CSD were always smaller than the ground truths. NB: "X" designates denotes that no peaks are separable in the 3D point cloud.

We now consider the influence of noise on the characteristic measurement of 3D point clouds in terms of length ratio, separability, and uncertainty, as shown in Fig. 3.14 and Table 3.5. The 3D point clouds in Fig. 3.14 correspond to SNR=20.

Table 3.5 gives the results of quantitative characteristic measurement of the three groups of point clouds at different SNRs. Three metrics were computed: mean length ratio, mean separability and mean uncertainty.

In terms of mean length ratio, as SNR decreased, the length ratios of the three methods were all reduced. The length ratio of CSA-group was more sensitive to SNR in comparison with AQBI and CSD, its length ratio reducing more rapidly as SNR decreased. The decreasing of CSA-group's length ratio was also greater at smaller crossing angles than at larger angles (90 o and 70 o ), even at SNR=30; such decreasing is particularly clear compared to noise-free conditions. In contrast, CSD-group showed the most stable behavior.

In terms of separability, AQBI-group was less sensitive to noise; its separability only decreased at SNR=10. However, its peaks were not separable at small crossing angles. The separabilty of CSA-group reduced more rapidly than the others with the decrease of SNR. Note however that the separability at small angles should be interpreted with precaution since the influence of noise on small angle crossing could be too important by altering the 3D point cloud. The separability of CSD-group only decreased slightly at SNR=10.

When SNR was reduced, the uncertainties were increased. AQBI-group presented the greatest uncertainties among these three methods. The uncertainties of CSA-group were more dependent on noise level than other two methods. The uncertainties of CSD-group were stable.

We calculate the 3D point clouds by AQBI (containing 3-fiber crossed each other at 90 o in each voxel) on condition of b-value=3000s/mm 2 at noise free and SNR=10. As shown in Fig. 3.15, in noise free condition, in terms of mean length ratio, the mean length ratio, mean separability and mean uncertainty is 1.9, 1.7 and 1.04. However, in noisy condition of SNR=10, mean length ratio and mean separability are decreased to 1.52 and 0.6. The mean uncertainty is increased to 1.7. The mean separability is decreased by serious (1.1) and the uncertainty is raised by 0.48.

Results on real brain data

3.3.2

The real data comes from a fixed and excised macaque brain hemisphere. The data was acquired on a Bruker 7T scanner [Welsh, 2013]. The brain was covered with 1.5 mm slices of size 128  128. The diffusion encoding was performed in 95 directions with b-value of 5000s/mm 2 .

In Fig. 3.16 (a), we present the generalized fractional anisotropy (GFA) (Tuch, 2004) map of the brain diffusion weighted data (black: GFA=0; white: GFA=1). In Fig. 3.16 (b) are shown the 3D point clouds on a coronal slice at the level of the corpus callosum, obtained using the AQBI, CSA and CSD methods. Fig. 3.16 (c) represents the zoomed version of the fiber crossing region circled by the red box in Fig. 3.16 (b). Five 3D point clouds corresponding to five different voxels are further zoomed in Fig. 3.16 (d). For AQBI, a spherical harmonic expansion of order L=8 and Laplace-Beltrami regularization were used. For CSA and CSD, a spherical harmonic expansion of order L=4 was used. The characteristic measurements of the 3D point clouds using the proposed length ratio, separability, uncertainty are given in Table 3.6. Visually, AQBI method presents poorer separability compared to CSA and CSD methods and does not allow recovering angular information. However, the difference between 3D point clouds is not always clear visually. This is the case for example for CSA and CSD in Fig. 3.16 (d). The results in Table 3.6 give a more quantitative and complete description of the 3D point clouds in Fig. 3.16 (d). In terms of length ratio, it is CSA that presented the greatest values. In terms of separability, the AQBI 3D point clouds presented little separability with a separability value nearly null. In contrast, CSA and CSD give significantly greater (and close) separability values at the same voxel, which implies that they can allow resolving more easily fiber crossing problems. In terms of uncertainty, CSA also presented the smallest values at the five voxels. At these 5 voxels, the difference in length ratio, separability and uncertainty is not very big between CSA and CSD, this explains why it is difficult to assess their difference from their surface representation.

Let us consider a region shown in Fig. 3.16 (c) (corresponding to the boxed region in Fig. 3.16 (b)). In terms of length ratio, AQBI has a length ratio of nearly 0.9, which means that the lengths of fibers in two directions are nearly the same. CSA and CSD present almost the same length ratio (0.75 for the former and 0.74 for the latter). CSA has a little higher length ratio than CSD at 9 voxels; the sum of length ratios in the region is 38.8 for CSA and 38.7 for CSD.

In terms of separability, we found that CSA has a higher separability than CSD at 18 voxels. The sum of the separability values at all the voxels of the boxed region of Fig. 3.16 (c) was equal to 5.15 for CSA and 4.45 for CSD. More precisely, CSA has a greater separability than CSD at 35.29% (18/51) of the voxels; CSD had a greater separability than CSA at 11.76% (6/51) of the voxels; for the rest of the voxels, CSA and CSD had the same separability.

In terms of uncertainty, the sum of uncertainties in the boxed region is 144.9 for AQBI, 113.0 for CSA and 114.7 for CSD. That gives an average uncertainty at a voxel of 0.64 for AQBI, 0.50 for CSA, and 0.51 for CSD. If we only sum the voxels at which the fiber crossing is resolved, the so obtained average uncertainty is 0.6 for AQBI, 0.52 for CSA, and 0.53 for CSD. So, with uncertainty metric, we can quantitatively say that CSA and CSD resolve fiber crossing more easily than AQBI, and that between CSA and CSD, the former resolves a bit more easily than the latter.

To sum up, we have used the proposed quantitative metrics to assess the AQBI, CSA and CSD reconstruction methods. The results can be summarized as follows. Regarding length ratio, AQBI produced a length ratio close to ground truth in both noise-free and noisy cases. CSA had a length ratio greater than ground truth. CSD presented slightly smaller length ratio than ground truth. Noise can make the length ratio smaller. About separability, AQBI gave the smallest separability than the others. In simulation experiments, CSA presented a little smaller separability than CSD. However in real data experiments, CSA exhibited greater separability than the others. Noise can decrease separability, especially for CSA. Concerning uncertainty, AQBI showed the greatest uncertainty among the three methods. CSA presented a little greater uncertainty than CSD in simulation experiments, but in real data experiments, CSA presented smaller uncertainty than other two. Noise could augment uncertainty, especially for CSA.

Compared to the visual quality of 3D point clouds, the proposed morphological metrics are more complete, accurate and direct for estimating the difference between 3D point clouds.

Discussion

Characteristic assessment of 3D point clouds in HARDI is challenging to perform due to their shape and topological complexity. Although a few quantitative metrics can be found sparsely in the literature, relatively complete description of 3D point clouds remains to accomplish.

Our results of the first experiment showed that the comparison of the shapes of 3D point clouds by MSE, sKL, RMSE, and NMSE generated some unreasonable results by ignoring the morphological characteristics of 3D point clouds. Sometimes irregular 3D point clouds are decided as better quality than other 3D point clouds more cloth to ground truth in visual comparison.

Our experimental results showed that CSA can lead to length ratios greater than ground truth. This can explain why CSA has a relatively greater ability to resolve fiber crossing. CSD-group can resolve crossing problem more easily than AQBI-group and presents smaller length ratios than ground truth. That is one of the reasons why CSD cannot resolve the problem of small angle crossings.

Rician noise can also reduce the separability of 3D point clouds. When the SNR decreases, the separability of 3D point clouds is reduced. We found that this is one of the most important difficulties in separating crossing fibers. On the other hand, the separability is lower at low bvalues than at high b-values. This conforms to the initial research [Frank, 2001]. The smaller the angle between fibers, the harder it is to distinguish them [Descoteaux, 2007]. Our separability metric can furthermore quantify this. CSA and CSD always present greater separability than AQBI, which explains why they can more easily resolve fiber crossing than AQBI.

Rician noise can finally increase the uncertainty of 3D point clouds. When the SNR decreases, the uncertainty of 3D point clouds increases. Again, this induces another difficulty in separating crossing fibers. When noise increases, CSA-group can increase more uncertainties than the other two groups; CSD-group can give smaller uncertainties than AQBI-group and CSA-group. Also, in the case of resolved fiber crossing, the smaller the angle between fibers, the more uncertainty increases with noise. CSA and CSD always exhibit smaller uncertainties than AQBI. This constitutes another reason why they can more easily resolve fiber crossing.

In the proposed morphological metrics, three important parameters a N , c N and p c are to be determined. The first two parameters can be fixed, depending on the desired accuracy. In the present study, we have chosen 9  is based on the fact that the greatest power of resolving fiber crossing is greater than 20 o , as indicated in [Özarslan, 2006;[START_REF] Jian | A Mixture of Wisharts (MOW) Model for Multifiber Reconstruction[END_REF]Aganj, 2010] where the authors resolved fiber crossing with a smallest angle of 28 o . When increasing a N , the number of small areas in ADM will increase and the angle interval will be reduced. This will increase the angular precision of analysis but make the detection of maximum directions more delicate. Likewise, if a N is too great, the resulting ADM will be too coarse to perform a correct measurement of 3D point clouds. The increase of c N will increase the number of distance intervals. This can improve the distance precision of analysis, namely the tiny difference in fiber length ratio as well as in separability, but will increase computing burden. Furthermore, for most point N is too small, neighboring points will be grouped in the same set, thus reducing the precision of analysis.

Concerning reconstruction directions ( N ), its setting as 162 in simulated data or 642 in real data is based on the fact that the reconstruction directions are geometrically uniform directions on a sphere. If N is too small, the distances at some directions will not be calculated. Reducing N will increase angular error and make some rows of ADMAT be all zeros. However increasing N will increase the number of 3D points, this will make the ADMAT more accurate, but will increase computing burden. Normally, N =642 will be sufficient.

Concerning p c , its setting as ( 1) 10  i is based on the fact that the nearer the 3D point to the center of sphere, the less it contributes to the shape of 3D point clouds. It is therefore used to weight the role of points as a function of their distance to the center of sphere.

Characteristics of the ODF peaks could provide clues to distinguish the underlying fiber structure intravoxel, for example a crossing fiber arrangement can be distinguished from a kissing fiber pattern based on characteristics of the ODF peaks.[ DS Tuch, 2003] The features of HARDI can have far reaching applications in segmentation, registration, and statistical characterization of regions of interest in the brain, as in comparing features between control and diseased patients. [START_REF] Schwab | [END_REF] The shapes of the peaks of functions output by multiple-fiber reconstruction algorithms reflect the underlying distribution of fibers. This extra information can be exploited to improve the tracography of fibers. [Seunarine, 2007] The volume fractions can improve statistical analysis of derived measurements, such as fractional anisotropy (FA). [START_REF] Hao | Joint Fractional Segmentation and Multi-Tensor Estimation in Diffusion Mri[END_REF] These researches provide proof of other practical usefulness of our morphological characteristics. For example, our framework of morphological characteristics can be also used to compare features between control and diseased patients.

Characteristics of the ODF peaks could provide clues to distinguish the underlying fiber structure intravoxel, for example a crossing fiber arrangement can be distinguished from a kissing fiber pattern based on characteristics of the ODF peaks.[ DS Tuch, 2003] The features of HARDI can have far reaching applications in segmentation, registration, and statistical characterization of regions of interest in the brain, as in comparing features between control and diseased patients. [START_REF] Schwab | [END_REF] The shapes of the peaks of functions output by multiple-fiber reconstruction algorithms reflect the underlying distribution of fibers. This extra information can be exploited to improve the tracography of fibers. [Seunarine, 2007] The volume fractions can improve statistical analysis of derived measurements, such as fractional anisotropy (FA). [Hao, 2013] These researches provide proof of other practical usefulness of our morphological characteristics. For example, our framework of morphological characteristics can be also used to compare features between control and diseased patients.

The fiber ODF quantifies the fraction of fiber portions (note that fibers may vary in orientation along their length) within a voxel with each orientation. [Heidi Johansen-Berg, 2009] Though we calculate the uncertainty of 1 th  line, according to , (

)) lw Card fset v  /2 (
a N , the uncertainty could be computed at any lines of ADMAT that will provide more information in orientation along ODF's length than uncertainty estimates from methods in [Seunarine, 2007] and will provide more information to distinguish the underlying fiber structure intravoxel such as fanning or bending. Using the proposed characteristic assessment paradigm, it is possible to extract more features of 3D point clouds.

Nonetheless, several limitations of the metrics proposed here merit discussion: Firstly, projecting 3D point clouds in 2D ADMAT enables us to readily compute morphological metrics. However, in doing so, we lose rotation information in 3D space, and the projection will affect the estimation of the angle of a fiber with respect to a fixed reference axis (positive x-axis for example). Nevertheless, the proposed morphological metrics are rotationally invariant, and we are not devoting to the measurement of relative angles between two fibers. Secondly, we need to choose a reference direction (the direction of the vector having the maximal radial distance). It is possible that there are more than one maximal radial distance point in a 3D point cloud. In this case, we can select anyone of them as reference direction. This is because the proposed morphological metrics do not change with reference direction. However, the estimation of the angle of a fiber with respect to a given reference axis will change if a different direction is used as the reference direction.

Conclusion

We have proposed a novel paradigm allowing for the characteristic assessment of general 3D point clouds including the ODF in HARDI. The paradigm is based on the measurement of the morphological characteristics of 3D point clouds. In this framework, three particular quantitative metrics have been proposed. The results showed that the proposed morphological metrics are consistent with the visual quality of 3D point clouds and enable us to describe quantitatively and accurately the characteristics of the latters, which provides a new way to quantify the characteristics and potentially the quality of 3D point clouds in resolving fiber crossings. INTRODUCTION

Abstract

Q-ball imaging can resolve multiple intravoxel fiber orientations and does not require any assumptions on the diffusion process. Analytical q-ball imaging (AQBI) was proposed to model the high angular resolution diffusion imaging (HARDI) signal using spherical harmonic basis and it incorporates a regularization term which can improve the performance of original q-ball imaging. Unfortunately, AQBI still needs many acquisitions in different gradient directions to compute accurately the orientation distribution function (ODF) (souvent more than 60 directions). The acquisition of signals from such important number of gradient directions poses the problem of long acquisition time, which is known to be a major impediment to the clinical application of q-ball imaging. Therefore, reducing acquisition time is an essential concern. The success of such reduction is possible through minimization of the number of diffusion encoding gradients required for reconstructing reliable ODFs. In this perspective, we propose a technique based on Delaunay triangulation for recovering missing diffusion signals in spatial space and imposing constraints in both q-space and spatial space. The results on simulated data, physical phantom data and real brain data consistently showed that the proposed missing diffusion signal recovering approach enables us to obtain accurate ODFs with relatively fewer number of diffusion signals.

Résumé en français

L'imagerie q-ball (QBI) peut résoudre de multiples orientations de fibres dans un voxel et ne nécessite pas d'hypothèses sur le processus de diffusion. L'imagerie q-ball analytique (AQBI) a été proposée pour modéliser les signaux en imagerie de diffusion à haute résolution angulaire (HARDI) avec une base harmonique sphérique et il intègre un terme de régularisation pour améliorer la performance de QBI original. Malheureusement, AQBI nécessite toujours un nombre important d'acquisitions afin de pouvoir calculer de manière précise la fonction de distribution d'orientation (ODF) (souvent plus de 60 directions). De telles acquisitions posent le problème de long temps d'acquisition, qui est connu comme un obstacle majeur à son application en routine clinique. La préoccupation essentielle est donc de diminuer le temps d'acquisition. Le succès d'une telle réduction est possible grâce à la réduction du nombre de directions de gradients de diffusion nécessaire pour la reconstruction des ODFs fiables. Dans cette perspective, nous proposons une technique basée sur la triangulation de Delaunay pour reconstruire des signaux de diffusion manquants dans l'espace spatial et en imposant des contraintes à la fois dans l'espace q et dans l'espace spatial. Les résultats sur des données simulées, des données de fantôme physique et des données réelles du cerveau ont montré que notre approche de reconstruire des signaux de diffusion manquants permet d'obtenir des ODFs précis avec un nombre réduit de signaux de diffusion.

Introduction

Diffusion MRI [Le [START_REF] Bihan | Imagerie de Diffusion in Vivo Par Résonance Magnétique Nucléaire[END_REF] is a noninvasive tool widely used in medical imaging applications to obtain information about the neural architecture of the brain [Le [START_REF] Bihan | Diffusion Tensor Imaging: Concepts and Applications[END_REF] or myocardial fiber structure in heart [START_REF] Yang | Feature-Based Interpolation of Diffusion Tensor Fields and Application to Human Cardiac DT-MRI[END_REF]. Diffusion tensor imaging (DTI) [Le Bihan, 2001] is widely used for studying the fiber structures of the brain and the heart. However, it is well known that this technique is based on the assumption of Gaussian diffusion. High angular resolution diffusion imaging (HARDI) is proposed to make up for the insufficiency of DTI [ DS Tuch, 2002;DS Tuch, 2002;Frank, 2002]. Q-ball imaging (QBI) being one of HARDI techniques, it presents the particularity of being a linear model without any assumption on diffusion distribution. In the field of QBI, most works reported in the literature focus on improving the quality of orientation distribution function (ODF). Some researchers put the emphasis on the development of new spherical basis or regularization methods [ [START_REF] Ds Tuch | [END_REF]Descoteaux, 2006Descoteaux, , 2007;;[START_REF] Hess | [END_REF]Michailovich, 2010c]. Other researchers work more on improving q-ball imaging with new acquisition methods or mathematical models. We can cite the work that exploits data from multiple qshells [Descoteaux, 2011] and compute ODFs with constant solid angle [START_REF] Aganj | Multiple Q-Shell ODF Reconstruction in Q-Ball Imaging[END_REF][Aganj, , 2010]].

HARDI is capable of capturing complex fiber configurations such as fiber crossing by sampling the unit sphere at a much greater number of gradient directions (usually more than 60 directions) [ [START_REF] Ds Tuch | [END_REF]][Descoteaux, 2007[START_REF] Hess | [END_REF]. This also leads to the increase in acquisition time and makes it difficult to use HARDI in clinical conditions. Therefore, we give more attention to recent researches on reducing acquisitions time while resolving the problem of the reconstruction of complex fiber configurations. Meanwhile, compressed sensing (CS) [Donoho, 2006] is a recent technique of accurately reconstructing sparse signals from undersampled measurements acquired below the Shannon-Nyquist rate; it has been applied in dMRI [START_REF] Lustig | Sparse MRI: The Application of Compressed Sensing for Rapid MR Imaging[END_REF]. The CS-based HARDI reconstruction has been proposed in [START_REF] Merlet | Compressed Sensing for Accelerated EAP Recovery in Diffusion MRI[END_REF][START_REF] Merlet | Spherical Polar Fourier Eap and Odf Reconstruction via Compressed Sensing in Diffusion Mri[END_REF]Michailovich et al., 2010a;[START_REF] Cheng | Compressive Sensing Ensemble Average Propagator Estimation via l1 Spherical Polar Fourier Imaging[END_REF]. The main idea of compressed sensing used in HARDI is to exploit under-sampling signal to obtain sparsity coefficients of spherical basis to a convex optimization problem and then use sparsity coefficients to do the compression of HARDI signals. Recovering missing diffusion signal is optional solution for minimizing the number of diffusion encoding gradients required for reconstructing a reliable ODF. Instead of using compressed sensing method, another solution is also can be used to recover the missing diffusion signal by the signal acquired. Therefore, we give attention to another method in recovering new data.

In [START_REF] Lertrattanapanich | High Resolution Image Formation from Low Resolution Frames Using Delaunay Triangulation[END_REF], an algorithm based on spatial tessellation and approximation of each triangle patch in the Delaunay triangulation is advanced to construct a high resolution (HR) high quality image from a set of low resolution (LR) frames. In [Sasaki et al., 2001a[Sasaki et al., , 2001b]], the reflection signal in the inverse synthetic aperture radar is measured in the polar coordinate defined by the object rotation angle and the frequency. An interpolation method by the nearest neighbor interpolation based on Delaunay triangulation was proposed.

According to the above Delaunay triangulation based method, we propose a new idea to deal with how to obtain accurate ODFs from fewer measurements. It consists in recovering no-acquired or missing data using Delaunay triangulation interpolation. More precisely, the proposed method consists of first inverse Fourier transforming the acquired diffusion data in q-space and applying the Delaunay triangulation interpolation to recover non-acquired probability distribution function (PDF), followed by Fourier transforming the obtained PDF, then imposing in an iterative manner both data consistent principle in q-space and the realness and positivity in the spatial space, and finally computing ODFs through using spherical harmonics basis-based analysis.

The rest of the chapter is organized as follows. In Section 4.2, the proposed diffusion signal recovering method is described in detail by giving the mathematical formulation of the problem. Section 4.3 presents the experiments and results, followed by Discussion and Conclusion in Sections 4.4 and 4.5.

Methodology of recovering missing diffusion signal data

It is well known that diffusion data in q-space and probability distribution function (PDF) in spatial space constitutes a pair of Fourier transform [Wedeen, 2005;[START_REF] Poupon | Real-Time MR Diffusion Tensor and Q-Ball Imaging Using Kalman Filtering[END_REF]. Our idea is then try to recover missing diffusion data using Delaunay triangulation technique in spatial space and imposing constraints in both spaces.

In order to recover missing diffusion data, we propose the following procedure: a) Inverse Fourier transform the diffusion data acquired in saying b) Perform 3D Delaunay triangulation using the points (PDF at each direction) on the sphere. Then according to the uniform distances between the 3D points on the sphere, we regard the signal as one-dimension signal. In Fig. 4.1(a), x-axis is the gradient direction number. y-axis is the value of ODF in this direction. c) Compute nearest-neighbor interpolation using Delaunay triangulation at the uniform directions to obtain the PDF at geometrically uniform 21 directions on the hemisphere, saying additional 21 directions as shown at the middle of Fig. 4.1(b). d) Fourier transform the Delaunay triangulation-interpolated PDF data, yielding a new diffusion signal of size equal to the sum of the original directions and the additional 21 interpolated directions, in which the data samples will be called the direction sampling points or simply direction points in the following. Based on the fact that ODF should be real and positive, taking real part and absolute value of ODF data.

e) In the new diffusion signal, replace the data at the initial directions ( 29 N  in total) by the initially acquired data. We weight the interpolated signal less than original signal as shown at the right of Fig. 4.1(b).

We represent the original signal and the interpolated signal on the 3D sphere as shown in Fig. 4.1(c).

After recovering the missing diffusion data, we reconstruct the ODF using conventional methods. An example of ODFs reconstructed from the original data (29 directions in this example) and from a data set with 81 directions is shown in Fig. 4.1(d). Mathematically, given the diffusion signal () Eq ,  

,, q

x y z  , its PDF is calculated by where FD designates the constructed Delaunay triangulation.

    1 () PE   r FT q , ( 4
In Fig. 4.1(a),  

Eqis represented by the red vectors. The length of the vectors represents the strength of  

Eq and the direction of the vector represents q . The nearest neighbor interpolation based on the Delaunay triangulation is denoted as

   
,, x y z  represents, in q-space, a higher resolution surface mesh with respect to the initial surface mesh expressed by q . The directions h q form a quasi-uniform set on the sphere.

P x y z 

  ,,

E E P      h h h k q if q q q k FT real r , (4.5)
where 12  kk and we also make      hh E q E q . In this chapter 12 3 and 1  kk .

Once the above final diffusion signal defined in N directions on the sphere is obtained, we project it on a spherical harmonic basis

B of size NK          1 1 1 1 1 1 ,, ,, B      K N N K N N YY YY         , (4.6) with    1 2 / 2 K l l   
, where l represents the order of spherical harmonics. That is

[11]     1 , E    K difsig h k k k q c Y  , (4.7) where   1 ,... c  T difsig K cc is a 1 K  vector.
The weighting coefficients are calculated using

    1 c B B L B E   difsig T T h q  , (4.8)
where the KK  diagonal matrix L is a Laplace-Beltrami operator, in which each matrix element is defined as       2 2 1 l k l k  , and  is a regularization factor that controls the smoothing effect of the estimator.

The ODF coefficients c ODF are finally computed using Funk-Radon transform and Funk- Hecke theorem. They are given by c Pc  ODF difsig , (4.9)

where P is a KK  Funk-Hecke diagonal matrix, whose element is defined as

        1 , K ODF kk k o c Y o o       , (4.10)
where o designates the direction, in which ODF is reconstructed.

Experiments and Results

Simulation experiments and results

4.3.1

We used spherical harmonic order max 8 l  , reconstruction directions=642 uniform directions on a sphere, and diffusion encoding directions or sampling density (K)=81 uniform directions on a hemi-sphere. The used simulated phantom concerns two fibers crossed from 48.6 o to 90 o . The simulated image is of size 16x16 pixels, and is divided into 3 regions: the single fiber regions in the top left corner and in the lower right corner; the crossed region in the middle. The diffusion signal was simulated using the multi-tensor model

  exp D  T i k i k i k S g = P (-bg g ) ,
where k P is the apparent volume fraction of the voxel with diffusion tensor D k , b is the diffusion sensitization factor, and i g is the direction of the diffusion gradient [ DS Tuch, 2002]. The spherical harmonic based compressed sensing method [Michailovich, 2011] is used to compare our method with the original AQBI method.

Here, the angular error  is computed (in degree) based on Earth Mover's Distance (EMD) [Rivera, 2012]. The correct detection rate used in this work is defined as

  c c Cart( ) / Cart( ) , r | N(r) N(r)      , (4.11)
where r is the voxel number,  is a set of all the voxels, c  is a set that estimated number of fiber N(r) equal to the true number N(r) of the fiber passing through voxel r, and card() designates the cardinal number. Here as a performance measure, we use the average that is obtained by averaging the values of the angular error, length difference, separability and uncertainty. We observe that in the noise-free case, our method presents better results even than the complete 81 direction set case, by giving visually clearer fiber detection. With the increase of noise levels, our results are still better than original data results, but the quality of angular resolution is reduced both for original data results and for our method's results. Nevertheless, at SNR=10, our method still presents better results than original data results. The clearer the fiber crossing, the higher fiber detection capability is. In this sense, our method improves the capability of fiber detection. With the simulated phantom in Fig. 4.2(a), we compare the original AQBI, spherical harmonic basis based compressed sensing method (SH based on CS method) and our proposed method in terms of the correct detection rate, average length ratio, average separability, average uncertainty and average angular error in Fig. 4.3.

We can see from Fig. 4.3 that the best performance is reached by our proposed method, which resulted in the highest correct detection rate from Fig. 4.3(e), for all the SNRs and sampling densities. Overall, all the criterion curves demonstrate that the error values increase proportionally with the decrease of SNR, while the error values go down with the increase of the number of diffusion-encoding sampling density (K). In all cases, the curves of the correct detection rate of our method are always above other method's curves, which indicates the robustness to noise and accuracy with respect to the value of K of the proposed method.

More precisely, in terms of detection rate in Fig. 4.3(a), the original AQBI, SH based CS method and our results all improved the results by increasing sampling density. At SNR=30, the detection rate is stable for our method with sampling density=30 and the detection rate is 65%, while other two methods produce 45% detection rate at this condition. With the decrease of SNR, the detection rate of our method decreases slightly but remains higher. With the sampling density=30, this rate is reduced to 50%.

In terms of the criteria presented in Chapter 3, the results of the comparison are shown in Fig. 4.3 (b), (c) and (d). As observed, the most accurate reconstruction is obtained using our method. As expected, the length ratio decreases and the separability decreases with the decrease of SNR and increases with the increase of K. As opposed to the case of the length ration and separability criteria, the uncertainty increases as the SNR decreases and decreases when K increasing.

In terms of length ratio, our method presents larger nearly 0.2 than other two methods at any conditions. This shows that our method improves the length ratio.

In terms of separability, with sampling density=30, our method has a separability of greater than 0.04 while the other two methods have 0.02. When the SNR is 20 and 10, our method presents a separability value of 0.04 and 0.025, while the original AQBI has the separability of 0.02 and 0.01 respectively, and the SH-based CS method the 0.02 and 0.01 respectively.

In terms of uncertainty, our method presents smaller than other two methods, except at SNR=10 and sampling density=16. When decreasing the SNR, our method presents an uncertainty of about 0.66, whereas the other two methods lead to greater uncertainty, which is greater at least 0.02 than results of our method.

From the angular errors in Fig. 4.3(e), we can see that, with the increase of sampling density, the reduction of average angular errors is obvious for all the three methods. The average angular error of our method is smaller at least 2 o than other methods at SNR=30 and 20. Though SH-based CS method also presents smaller angular errors than original AQBI at SNR=10, CS method still presents more angular error than ours. When decreasing the SNR, our method always presents smaller angular errors than other methods.

Results of fiber cup phantom 4.3.2

The diffusion-weighted data of the fiber cup phantom were acquired on the 3T Tim Trio MRI system by the group [Fillard, 2011]. The b=2000s/mm 2 dataset was acquired with a field of view (FOV) =19.2 cm, a resolution of 3 mm x 3 mm x 3 mm, and a matrix size of 64 64  .

The acquisition parameters for the b=2650s/mm 2 dataset are: FOV=38.4 cm, matrix size= 64 64  , b-value=2650s/mm 2 , resolution=6 mm x 6 mm x 6 mm [Fillard, 2011]. In Fig. 4.4, we show the example of an ODF to illustrate the proposed method. For comparison, we use the ground truth taken from [Fillard, 2011], which corresponds to b=2650s /mm 2 and sampling density K=29.

On the right of To analyze in detail the quality of ODFs and the performance of the proposed method, we used the b=2000s/mm 2 phantom dataset, a set of 64 orientations as complete orientation set for comparison, and a subset of 29 orientations for reconstruction with AQBI and our method combined with AQBI.

In the fiber crossing region, we observe that analytical q-ball imaging using complete 64 directions cannot produce clearly resolved ODFs. Using 29 directions, the simple AQBI reconstruction do not allow all the fiber crossings to be brought out. Always using 29 directions, but with the proposed method, nearly all the fiver crossings come out.

Moreover, visually, the ODFs obtained with the proposed method look even better than the ODFs reconstructed from the full 64 directions. The ODFs obtained using the complete 64 directions are of medium quality. Using only 29 directions, without missing diffusion data recovering by our method, the resulting ODFs are of poor quality and cannot be accepted. For further illustration, we show in Fig. 4.6(a), the results of ODF reconstruction on another physical phantom dataset with a b-value of 2650s/mm 2 , using always the same simple AQBI method and our diffusion data recovering combined with the AQBI, and considering always the same conditions: with the complete set of 64 directions and with a subset of 29 directions.

In the crossing region, we observe nearly the same results as those in the condition of b=2000s/mm 2 , with the difference that our method exhibits even better results than in the case of b=2000s/mm 2 . Indeed, with our method, all the fiber crossings are well resolved in the ODFs for b=2650s/mm 2 .

The visual comparison leads to the results similar to the quantitative results from the simulations: simple AQBI cannot give well-defined ODFs with correct maxima when using few diffusion directions, but our method leads to better resolved ODFs.

Brain data experiments

4.3.3

We further evaluate our method on a fixed, excised brain hemisphere. The measurements were acquired on a Bruker 7T scanner [Welsh, 2013]. The brain was covered with 1.5 mm slices with 128 128  matrix size. The diffusion encoding was performed in 95 directions with a b-value of 5000s/mm 2 . From these ODF maps, we observe that the proposed method generates clearer fiber crossing ODFs in the region of fiber crossing; the ODFs have more distinguishable ODF maxima important for fiber crossing characterization.

The region inside the red box in Fig. 4.7 is likely to contain complex fiber configurations, especially fiber crossing, because it is the intersection of adjacent regions presenting different main orientations, but in each of the adjacent regions, the orientations are rather homogeneous.

The results show that the difficult crossing voxels presented by our method in Fig. 4.7 (a) ,(b) and (c). In the original ODF maps, the fibers in these adjacent regions are shown with small GFA, however, they cannot be detected with crossing directions. Comparing to pure AQBI, it is remarkable that the crossing region can be clearly shown. Especially in the marked red box, the adjacent regions are shown with more continuity, so it is more reasonable. In Fig. 4.7 (c), the more complex region was presented for comparing our method and pure qball imaging. We find that in the regions that were marked by red box, our method give a more logical results based on fiber continuity and GFA results.

Discussion

In the proposed approach, the diffusion signals were regarded as one-dimensional, according to the order of the diffusion signal directions in the whole set of geometrically uniform directions ranking from the z-axis direction to the y-axis direction.

The interpolated directions are chosen as the geometrical uniform set whose level lower than whole-set of data (For example, if the whole-set direction is 81 on hemi-sphere, the interpolated set is 21 geometrical uniform directional set.).

Regarding the weighting of recovered data and original data, the weighting coefficient is determined empirically. In the present study, we have chosen the weighting coefficient for the original diffusion data, which is three times the weighting coefficient for the recovered diffusion data. If we reduce too much the difference in weighting coefficient between the original data and recovered data, the recovered data will play a too important role in reconstruction, and the results will not be very stable. If we increase too much this difference, the recovered data will play negligible roles and will lose effect on the reconstruction of ODFs.

In some cases, our method cannot present very stable performance. This is because the reduced number of diffusion signals is not always sufficient to detect all the main directions of ODFs. Therefore, the subset of diffusion signals is an important factor. For this reason, how to choose a good subset of diffusion signals is a very interesting problem for the future. Now in this field , there have been some researches. For example, in [START_REF] Mcewen | A Novel Sampling Theorem on the Sphere[END_REF], the authors proposed a novel sampling method on the sphere.

Conclusion

We have proposed a new method of recovering missing diffusion signals for the reconstruction of high quality ODFs. The method is based on the nearest neighbor interpolation using Delaunay triangulation. Compared to the original AQBI method, the proposed method improves substantially the resolution of crossing fibers. The proposed method also provides a tool for improving ODF quality with reduced diffusion data.

Résumé en français

Le compressed sensing (CS) a été appliqué au domaine de l'imagerie de diffusion à haute résolution angulaire (HARDI) pour pallier le problème du temps d'acquisition long. Cependant, les méthodes existantes formulent le procédé actuel proposé une formulation du problème de reconstruction basé sur le CS de signaux de diffusion, dans lequel les contraintes clairsemé sont utilisées. Donc, les méthodes existantes basées sur CS représentent les signaux de diffusion par un petit nombre de coefficients, et fournissent ainsi un certain niveau de compression de données et de débruitage. Ce chapitre décrit une méthode permettant d'obtenir des fonctions de distribution d'orientation (ODF) de haute qualité en utilisant une méthode de reconstruction de signaux de diffusion manquants basée sur le CS combiné avec le ré-échantillonnage des directions de gradient. Les résultats montrent qu'avec la méthode proposée, la capacité de résolution angulaire des ODFs reconstruites à partir des données de diffusion réduites est proche de celle des ODF reconstruites à partir de l'ensemble complet de données de diffusion initiales, et que la méthode proposée permettent d'obtenir des ODFs plus précises et de résoudre le problème de croisement des fibres avec des données de diffusion réduites.

Introduction

Diffusion tensor imaging (DTI) [Le Bihan, 2001] method can only resolve a single voxel with a single fiber direction. If the diffusion of water molecules does not meet Gaussian distribution, which is the case when the environment presents complex spatial structures such as fiber crossing, DTI will not work. Diffusion spectrum imaging (DSI) [Wedeen, 2005] and high angular resolution diffusion imaging (HARDI) has been proposed to cope with the problems encountered by DTI [ DS Tuch, 2002][Wedeen, 2005] [Frank, 2002]. The DSI technique requires gradients sampling on a three-dimensional Cartesian coordinate, which is more time-intensive than HARDI [ [START_REF] Ds Tuch | [END_REF]] that uses fewer acquisitions than DSI on a sphere.

In the field of HARDI, the spherical deconvolution method proposed that the diffusionweighted signal attenuation measured over the surface of a sphere can be expressed as the convolution over the sphere of a response function (the diffusion-weighted attenuation profile for a typical fiber bundle) with the fiber orientation density function (ODF). [Tournier, 2004[Tournier, , 2008]]. It provides an estimate of the fiber orientation distribution (FOD) by assuming the diffusion weighted (DW) signal measured from any fiber bundle is adequately described by a single response function. The diffusion orientation transform (DOT) utilizes the Fourier transform relationship between water displacement probability and diffusion-attenuated signal [Ozarslan., 2006]. Orientation maps can be calculated from HARDI data by assuming that the signal attenuation along each radial line in q-space is characterized by a monoexponential decay. The mixture of Wisharts (Wisharts distribution) used the Laplace transform defined on the space of positive definite diffusion tensors by assuming each voxel is associated with an underlying probability distribution defined on the space of diffusion tensors [Jian, 2007]. Analytical q-ball imaging (AQBI) [Descoteaux, 2007] used Funk-Radon transformation and Funk-Hecke theorem to obtain the simplification of approximating the orientation distribution function (ODF). Among these HARDI methods, AQBI presents the advantages of being linear, relatively easy to program, and less affected by noise. Therefore, we use the AQBI for reconstructing ODFs in what follows.

However, the bottleneck problem encountered by the HARDI method is still how to use fewer measurements to reconstruct a relatively better quality ODF. Therefore, we give attention to compressed sensing (CS) based methods recently used in HARDI for reducing acquisitions or for resolving fiber crossing or bending problems without requiring too many acquisitions. Compressed sensing [Candes et al., 2006a[Candes et al., , 2006b;;Donoho, 2006;[START_REF] Candès | An Introduction to Compressive Sampling[END_REF] was initially proposed to reduce sampling density and shorten acquisition time and it has been applied to dMRI [Lustig, 2007]. Compressed sensing based HARDI reconstruction has been proposed in [Merlet, 2010[Merlet, , 2011;;Michailovich, 2010a;Cheng, 2011;[START_REF] Landman | Resolution of Crossing Fibers with Constrained Compressed Sensing Using Diffusion Tensor MRI[END_REF]. These methods can be classified into two groups: the first is to use compressed sensing in diffusion spectrum imaging for recovering ensemble average propagator (EAP), where CS is used in the spatial domain; the second consists of recovering HARDI signals from a much smaller number of measurements directly in the diffusion domain. In [Michailovich, 2010a], a HARDI using compressed sensing method was proposed, but the authors used spherical ridgelet basis and spherical harmonic basis of l=8 to define the compressed sensing matrix. The relationship between the quality of ODF and the spherical harmonic order is not verified.

The problem of using a much smaller number of HARDI measurements can be seen as an underdetermined problem. Overcoming such underdetermination problem will be possible by making use of the fundamental premise of the CS theory that states that an accurate estimation of spherical harmonic coefficients is attainable if the latters are sufficiently sparse and if the sampling and representation bases are sufficiently decorrelated.

In the AQBI, sampling functions on a sphere are Dirac and the most commonly used representation basis is spherical harmonic basis. The mutual coherence of the spherical harmonic basis with respect to the Dirac sampling functions is small enough.

In this chapter, we verified the method to deal with how to obtain more accurate ODFs using fewer measurements by spherical harmonic basis. The method first uses CS model to obtain the sparse coefficients of diffusion signals on spherical harmonic basis. Then based on the obtained sparse coefficients on spherical harmonic basis, a CS based recovering method [START_REF] Elad | Simultaneous Cartoon and Texture Image Inpainting Using Morphological Component Analysis (MCA)[END_REF] is employed to obtain missing or non-acquired diffusion signals in the diffusion domain. The recovered diffusion signals are finally used as inputs of AQBI to improve fiber detection.

The rest of the chapter is organized as follows. In Section 5.2, the subsets of diffusion signals and the proposed CS-based method are presented. Section 5.3 presents the experiments and results, followed by Discussion and Conclusion in Sections 5.4 and 5.5.

Methodology of recovering missing diffusion signals

As shown in Fig. 5.1, to recover the missing diffusion signals, in the experiments, the independent (here, Dirac sampling functions and the spherical harmonic representation basis) and uniform sampling method is needed on the sphere. We first use a method based on the electrostatic sampling and tessellation icosahedron (geometrically uniform) sampling to obtain randomly generated quasi-uniform subsets of the directions of diffusion signals. Then, we use the 1 l norm constraint to obtain the sparsity expression (coefficients c) of the subset of diffusion signals represented by spherical harmonic basis, which is used to recover missing or non-acquired diffusion signals. Finally, the AQBI method is applied to the recovered diffusion signals to reconstruct ODFs and we compare the results with the reconstruction results from the original complete set of diffusion signals.

Tessellation icosahedron (geometrically uniform sampling)

The icosahedron (20 faces, 12 vertices) is a polyhedral with the highest possible number of identical faces that inscribe a sphere. Therefore, it is used as the basis of subdivision for obtaining uniformly distributed spherical points. [Teanby, 2006] The subdivision of the icosahedron proceeds as follows:  Define each vertex such that the center of the icosahedron is at the origin.  Project each vertex onto the unit sphere by normalizing the length of each vertex vector to unity. This projection ensures that the angular separations obtained from successive subdivisions are equal.  Define an array of triangles, where each triangle is defined by three adjacent vertices.  For each triangle, bisect each side to give a new vertex point, by summing the position vectors of the two end vertices. Project this vector onto the unit sphere by normalizing to unit length.  There are now four triangles for each original triangle.  Reject repeated vertices. The vertices of the resulting subdivided spherical icosahedron provide an evenly distributed grid of M points.

The number of times the subdivision is performed refers to as the "level" of the grid, where level 0 corresponds to the initial icosahedron having a number of vertices M=12, level 1 to the M=42 grid, and level 2 to the M=162 grid.

Thomson problem (electrostatic uniform sampling)

Coulomb's law states that the magnitude of the electrostatic force of interaction between two point charges is directly proportional to the scalar multiplication of the magnitude of charges and inversely proportional to the square of the distance between them. For distinct particles in the point set

  1 NS i i Ep e  
, the energy is defined as

2 11 ()       NN ij i j i
Energy e e Ep .

(5.1)

A configuration of spherical points is considered optimal if it minimizes the Coulomb potential. If two points are close, then the reciprocal of their distance will be large, so minimizing the energy is a way of uniformly spreading out N points on a sphere [START_REF] Erber | Equilibrium Configurations of N Equal Charges on a Sphere[END_REF].

Quasi-uniform subsets from the geometrically uniform complete set

To simulate the reduction of diffusion directions, we denote the gradient direction of a diffusion signal by q . We can then represent the complete set of diffusion directions and its

subsets by   1 M m m q  and     1 1 NS M pm m p qq   
respectively, here, M and NS is the number of directions in a complete set and a chosen subset respectively. Two techniques of uniformly sampling a sphere are used: the electrostatic sampling and the tessellation icosahedron (geometrically uniform) sampling. It is well known that diffusion signals are symmetric, namely

  () E q E q  (  
Eq is the diffusion signal in the diffusion direction q ). Therefore, we use the signals only on a hemisphere.

Then we exploit the electrostatic technique to obtain a subset of NS quasi-uniform directions on the hemisphere, here

NS M  . The quasi-uniform subset 1 {} NS pp q  of   1 M m m
q  can be obtained as follows.

Let

  1

NS i i

Ep e   denote a quasi-uniform direction set randomly generated by electrostatic sampling technique using (5.1). To generate subsets of diffusion directions, we defined the formulation below,

  2 | min( ), [1 ], [1, ] m m i q q e m ,M i NS    . (5.2)
If meeting the condition in (5.2), we will have pm qq  . To obtain a quasi-uniform subset

  1 NS p p q  from a quasi-uniform complete set   1 M m m q  , for each i e
, when the distance between m q from the complete set   1 M m m q  and i e from the randomly generated quasi-uniform set

  1 NS i i e  is minimum. Note that p q is necessarily a direction of the set   1 M m m q  .
As shown in Fig. 5.

2 (a),   1 M m m
q  is generated by tessellation icosahedron method, M=81, the points are quasi-uniform on a hemisphere. In Fig. 5.2 (b), the NS=21, 29 and 39 (from left to right) subsets from the complete set in (a) are represented by red points. ) of NS directions on a hemisphere is obtained, we project it on a spherical harmonic basis sub B of size

 NS K         1 1 1 1 1 1 ,, ,, B      K sub NS NS K NS NS YY YY         , (5.3)
with

   1 2 / 2    K l l
, where l represents the order of spherical harmonics, here we use

6  l . Spherical harmonic decomposition of diffusion signal is,     1 , K kk k q c Y     E , (5.4)
where

  1 ,... c  T sub K cc is a 1  K
vector. [Descoteaux, 2007] We now try to combine the compressed sensing technique with missing signal recovering technique to obtain more accurate ODFs.

Compressed sensing was used to reconstruct sparse vector of diffusion signals in the directions of

1 {}  NS pp q
. The basic model is:

1 arg min cc  sub c subject to   1 Bc  p Eq  , (5.5)
where B is the compressed sensing matrix defined by spherical harmonics basis of the order l (in what follows, we will call it spherical harmonic sensing matrix), c is defined as a vector of spherical harmonic coefficients of the subset diffusion signals with its sparsity vector expressed as c sub , and   E p q is the subset signals in q-space.

If assuming

 

Bc e  p Eq

, where e is an error vector 2  e  , then the theory of CS are useful to find the sparsity coefficients c sub . (5.5) is a convex optimization problem, we use 1 lmagic software to resolve (5.5). [START_REF] Cand | l1-Magic : Recovery of Sparse Signals via Convex Programming[END_REF] Then we use compressed sensing based missing data recovering model,

        1 1 1 1 1 1 ,, ,, B      K RS RS RS K RS RS YY YY         , (5.6)   RS E q B c  sub RS , (5.7)
where the sparse coefficient vector c sub is obtained by means of (5.5). The re-sampling directions are the uniform directions on the hemisphere. Here, we choose geometrically uniform set at level 1 where the  KK diagonal matrix L is a Laplace-Beltrami operator, in which each matrix element is defined as

      2 2 1  l k l k
, and  is a regularization factor that controls the smoothing effect of the estimator.. where P is a  KK Funk-Hecke diagonal matrix, whose element is defined as

() 2 (0)  kk l k pP 
, with () (0) lk P designating an () lk order of Legendre function. After obtaining c ODF , we can calculate the value of the ODF in any direction using where o designates the directions in the set 1 {} 

RS nn q

, in which ODF is reconstructed.

        1 , K ODF kk k o c Y o o       , (5.10)

Experiments

Two fiber crossed at 70 o was simulated for comparing the proposed method at different subsets of directions and different spherical harmonic sensing matrix with different spherical harmonic order l. b-value 3000 s/mm 2 is used in the simulation. The diffusion signal was simulated using the multi-tensor model  

exp D  T i k i k i k S g = P (-bg g ) ,
where k P is the apparent volume fraction of the voxel with diffusion tensor D k , b is the diffusion sensitization factor, and i g is the direction of the diffusion gradient [ DS Tuch, 2002]. For generating the crossing fibers in this simulation, we use the diffusion tensors SN is the signal with rician noise. [Pierpaoli et al., 1996][Kingsley, 2006] The used data in Fig. 5.5 concerns two fibers crossed at 90°, 70° 50°, and 30°. The simulated image is of size 10x10 pixels, and is divided into 4 regions: the upper left region corresponds to the fibers crossed at 50°, the upper right region at 70°, the lower left region at 90°, and the lower right region at 30°. Each pixel presents two main intravoxel crossed directions.

In simulated phantom experiments in Fig. 5.7, we used spherical harmonic order max 8 l  , reconstruction directions=162 uniform directions on a sphere, and diffusion directions=81 uniform directions on a hemisphere. The used data concerns two fibers crossed from 48.6 o to 90 o . The simulated image is of size 16x16 pixels, and is divided into 3 regions: the single fiber region in the top left corner and in the lower right corner; the crossed region in the middle. To generate different crossing angles, we have rotated the diffusion tensor [Kingsley, 2006]. We do 10 times for the same experiments using this simulated phtantom.

In real data experiments, to test the performance of the proposed method, fully sampled spin echo data (96 encoding directions, b=5000 s/mm 2 , 128 128  matrix size, 4.5 cm FOV, 1.5 mm slice thickness, TE=30 ms, TR=1000ms) was acquired on a fixed, excised macaque brain hemisphere using a Bruker 7T scanner. Images reconstructed by the subsets of encoded directions of 19, 25, 31 quasi-unifrom directions.

Results on simulated data of different subset gradient directions 5.3.1 The upper subplot of Fig. 5.5 shows the subset diffusion signals that we randomly chose for reconstructing ODF from different whole-sets and the ODFs reconstructed by the subset signals. At the top of Fig. 5.5(a), we chose from a geometrically uniform whole-set of 81 directions by the rule of sampling rate of 5 (81/5). By the subset of diffusion signals, analytical q-ball imaging (AQBI) is unable to directly provide an estimation with fiber directions. The ODF obtained using our proposed method represents clearly the 90 o crossing at all of the orders of spherical harmonic basis used above. The ODF obtained using our proposed method represents clearly the 70 o crossing at only spherical harmonic order l =6 for all of these three subsets.

The middle subset of 18 quasi-uniform directions of diffusion signals from geometric uniform 81 directions shows the results after processing by our proposed method. The results of middle subset shows a better result closer to ground truth in Fig. 5.4 than the results of other two subsets. The right subset of 18 quasi-uniform directions from electrostatic uniform 81 directions shows nearly the same result of the middle one, but the result by middle subset is better for representing the crossing structure at spherical harmonic order l=4. For all of these three subsets, at spherical harmonic order l=6, the results present clearer fiber crossing structures. Therefore, in experiments below, we only used our proposed method based on spherical harmonic sensing matrix order l=6. The upper subplot of Fig. 5.6 shows the subset diffusion signals that we randomly chose for reconstructing ODF from different whole-sets and the ODFs reconstructed by the subset signals. At the top of Fig. 5.6 (a), we chose from a geometrically uniform whole-set of 81 directions by the rule of sampling rate of 5 (81/5). By the subset of diffusion signals, analytical q-ball imaging (AQBI) is unable to directly provide an estimation with fiber directions. The ODF obtained using our proposed method represents clearly the 70 o crossing. The middle subset of 18 quasi-uniform directions of diffusion signals from geometric uniform 81 directions shows the results of directly AQBI (left) and after processing by our proposed method (right). Our method also shows a better result closer to ground truth than the results of AQBI directly. The right subset of 18 quasi-uniform directions from electrostatic uniform 81 directions shows nearly the same result of the middle one by our method, but the result by AQBI is always not clear for representing the crossing structure.

The results of different subsets also show that our method always present the better results than directly using AQBI method. This shows that after processing by our recovering missing data method, the less diffusion signals can be used for reconstructing ODFs more accurate.

Results on simulated data of different sampling density and noise levels 5.3. Here angular error  is also computed (in degrees) based on Earth Mover's Distance (EMD) [Rivera, 2012]. The correct detection rates used in this work is defined as

  ( ) / ( ) , | (r) N(r) cc Cart Cart r N      , (5.11)
where r is number of the voxel,  is a set of all the fibers, c  is a set that estimated number of fiber N(r) equal to true number (r) N of fiber passing through voxel r . Here as a performance measure, we use an average for all of the measures, which is obtained by averaging the values of angular error, length difference, separability and uncertainty.

With the simulation phantom in Fig. 5.7 (a), we contrast the performance of directly using AQBI, spherical harmonic basis based compressed sensing method (SH based CS method without processing by recovering missing data) and our proposed method in terms of the correct detection rate, average length ratio, average separability, average uncertainty and average angular error. We can see that the best performance here is attained by our proposed method, which results in the best correct detection rate, all the tested values of SNR and sampling density. Overall, all the criterions curves demonstrate that the error values are increasing proportionally with the decrease of SNR, while going down with an increase in the number of diffusion-encoding sampling density (K). However, as comparing with other two methods, the curves of the correct detection rate of our method are always above other methods, which indicate the more accurate results at noisy condition and the different values of K.

In terms of detection rate, with reduction of SNR, our method presents a detection rate more than 50%. This is higher than the results of others that are less than 40%.

Our approach has been also compared by morphological metrics. The results of this comparison are summarized in Fig. 5.7 (c), (d) and (e), which again indicates that the most accurate ODF is obtained by using our method. As expected, the length ratio and the separability decreases as SNR decreases and grows with K increasing. As opposed to the case of these two criterions, the uncertainty increases as SNR decreases and grows with K increasing.

In terms of average length ratio, with reduction of SNR, the length ratio of three methods is reduced, but our method present higher nearly 0.2 than other methods.

In terms of separability, the result of our method is always higher than other methods. With the reduction of SNR, our method present average separability higher than other two methods nearly 0.04 at each sampling density.

In terms of uncertainty, with the reduction of SNR, our method always present lower uncertainty, comparing with other two methods, our result is at least 0.04 lower than results of other methods for sampling density more than 20.

The angular error in Fig. 5.7(f) is also used for comparing different methods. We can see that, in any cases here, our method could yield a smallest average angular error. As shown in Fig. 5.8 (a), the b0 image of a coronal cross-section of the brain and its reference field of ODFs corresponding to the region indicated by the red rectangular in (b). The whole-set of sampling density is equal to 96 diffusion gradients. We reconstructed ODF map using this whole-set directly by AQBI, as shown in Fig. 5.8 (b). The ODFs obtained in (b) are used as a reference against which different subset reconstruction results were compared.

Three different subsets of 19, 23, and 29 gradient directions are chosen from the original set of 96 diffusion gradients. In each of these subsets, their corresponding directions were chosen as a quasi-uniform coverage of the hemisphere.

According to the visual evaluation with assumption of fiber continuity, it is very possible that the regions in red box are the regions crossing.

As expected, we can see that our reconstruction results are more accurate than AQBI, according to the clearer crossing directions detected. We can see from the ODFs indicating with red box in our chosen region, the ODFs by our method is clearer for showing fiber crossing than original ODFs. The most accurate results in the regions in red box is even better than original reference results.

Discussion

There is an infinite number of solutions which would fit the constraint 1 arg min cc  sub c in (5.5), then a useful sparsity estimate c sub can be obtained as a solution to the convex optimization problem. However, the constraint here could also be changed to l 0 norm [START_REF] Daducci | Sparse Regularization for Fiber ODF Reconstruction: From the Suboptimality of ℓ2 and ℓ1 Priors to ℓ0[END_REF] or TV norm [Michailovich, 2011]. The representation basis here is spherical harmonic, it can be also replaced by other basis, such as the spherical ridgelet basis and spherical wavelet basis in [START_REF] Kezele | Spherical Wavelet Transform for ODF Sharpening[END_REF]Michailovich, 2010c;[START_REF] Leistedt | Exact Wavelets on the Ball[END_REF] In [Michailovich, 2010a], a HARDI using compressed sensing method was proposed, but the authors used spherical ridgelet basis and compared with spherical harmonic basis of l=8.

Here we used compressed sensing method in HARDI of spherical harmonic basis of l=6. We show that it is more stable than l=8 in [Michailovich, 2010a]. We also show ODFs reconstructed by the different subsets of diffusion directions from different complete sets. The morphological metrics are used to compare different ODFs reconstructed from different subsets in different noise level cases.

Comparing with the method that we introduced in Chapter 4, this method is more stable for the changing of the subset of diffusion signals. However, the ODF obtained by this method is less angular resolution than the method presented in Chapter 4.

The best subsets for qualified ODFs are an interesting research problem. There have been some researches, for example [McEwen, 2011], they proposed a novel sampling method on the sphere.

Conclusion

We have proposed a method for recovering missing or non-acquired diffusion data. The method presents the particularity of combining the geometrically uniform sampling and random electrostatic sampling to obtain any quasi-uniform sampling, using the CS technique to sparsify the spherical harmonic representation of diffusion signals, and increasing resampling density to recover missing diffusion signals in new diffusion gradient directions. We show that spherical harmonic order l of compressed sensing matrix and subsets of diffusion directions can affect the quality of ODF using CS-based method. Spherical harmonic order l=6 used in compressed sensing matrix and quasi uniform subset of diffusion directions can improve the quality of ODF maps.

The results on simulated phantom and real brain data consistently showed that the proposed method allows main crossing directions in ODFs to be more clearly brought out, which is interesting for reducing acquisition time and for clinical application of q-ball imaging.

Abstract

Diffusion tensor imaging (DTI) and high angular resolution diffusion imaging (HARDI) are often used to analyze the fiber complexity of tissues. In these imaging techniques, the most commonly calculated metric is anisotropy, such as fractional anisotropy (FA), generalized anisotropy (GA), and generalized fractional anisotropy (GFA). The basic idea underlying these metrics is to compute the deviation from free or spherical diffusion. However, in many cases, the question is not really to know whether it concerns spherical diffusion. Instead, the main concern is fiber complexity such as fiber crossing in a voxel. In this context, it would be more direct and effective to compute the deviation from a single fiber bundle instead of a sphere. We propose a new metric, called PEAM (PEAnut Metric), which is based on computing the deviation of orientation diffusion functions (ODFs) from a single fiber bundle ODF represented by a peanut. As an example, the proposed PEAM metric is used to classify intravoxel fiber configurations. The results on simulated data, physical phantom data and real brain data consistently showed that the proposed PEAM provides greater precision than FA, GA and GFA and enables parallel and complex fibers to be better distinguished.

Résumé en français

L'imagerie du tenseur de diffusion (DTI) et l'imagerie de diffusion à haute résolution angulaire (HARDI) sont souvent utilisées pour analyser la complexité de fibres dans un voxel. Dans ces techniques d'imagerie, la mesure la plus utilisée est l'anisotropie, comme l'anisotropie fractionnelle (FA), l'anisotropie généralisée (GA), et l'anisotropie fractionnelle généralisée (GFA). L'idée de base qui sous-tend ces mesures est de calculer l'écart par rapport à la diffusion libre ou sphérique. Cependant, dans de nombreux cas, la question n'est pas vraiment de savoir s'il s'agit d'une diffusion sphérique. Au contraire, la principale préoccupation est la complexité des fibres, comme le croisement de fibres dans un voxel. Dans ce contexte, il serait plus direct et plus efficace de calculer la déviation par rapport à un faisceau de fibres au lieu d'une sphère. Nous proposons une nouvelle mesure, appelé PEAM (PEAnut Metric), qui est basée sur le calcul de l'écart des fonctions d'orientation distribution (ODF) par rapport à l'ODF d'un faisceau de fibres, qui est caractérisée par une forme cacahuète (peanut). Comme un exemple, la mesure PEAM proposée est utilisée pour classifier des configurations de fibres dans un voxel. Les résultats sur des données simulées, des données de fantôme physique et des données réelles du cerveau ont montré que la PEAM proposée offre une plus grande précision que FA, GA et GFA et permet de mieux distinguer les fibres parallèles et des fibres complexes.

Introduction

Diffusion magnetic resonance imaging (dMRI) [Le Bihan, 1985] including diffusion tensor imaging (DTI) [Le Bihan, 2001] and high angular resolution diffusion imaging (HARDI) [DS Tuch, 2002, 2004;Tournier, 2004] is a promising technique allowing us to get insights into the fiber structures of white matter in the brain [Basser, 1996;DS Tuch, 2002] and myocardium in the heart [START_REF] Geerts | Characterization of the Normal Cardiac Myofiber Field in Goat Measured with MR-Diffusion Tensor Imaging[END_REF]Yang, 2012;[START_REF] Wei | Assessment of Cardiac Motion Effects on the Fiber Architecture of the Human Heart in Vivo[END_REF]. In this technique, diffusion anisotropy metrics are often computed to quantify intravoxel fiber complexity [Basser, 1996;[START_REF] Ds Tuch | [END_REF]Descoteaux, 2006;Yeh et al., 2013].

In DTI, diffusion direction information is encoded in diffusion tensors and diffusion anisotropy is typically quantified through fractional anisotropy (FA) metric [Basser, 1995]. However, the main limitation of DTI is that this technique can only resolve a single fiber or fiber bundle direction within each voxel [ DS Tuch, 2002]. To overcome this limitation, HARDI was proposed in which the generalized anisotropy (GA) metric [Ozarslan, 2005] based on high-order diffusion tensor was defined. In [ [START_REF] Ds Tuch | [END_REF], based on q-ball imaging, scalar metrics were derived from orientation distribution function (ODF) to define tissue contrast, especially generalized fractional anisotropy (GFA) was defined as an extension of the FA metric.

The above-mentioned metrics like FA, GA and GFA share a common point, namely, they are all based on computing the deviation from free or spherical diffusion. However, in many cases, the question is not really to know whether it concerns spherical diffusion. Instead, the main concern is fiber complexity such as fiber crossing in a voxel. For example, recently the prevalence of complex fiber configurations in white matter tissue with HARDI is investigated [START_REF] Haroon | Using the Model-Based Residual Bootstrap to Quantify Uncertainty in Fiber Orientations From-Ball Analysis[END_REF][START_REF] Jeurissen | Investigating the Prevalence of Complex Fiber Configurations in White Matter Tissue with Diffusion Magnetic Resonance Imaging[END_REF]. Though complexity of fiber configuration may be reflected by number of maxima, sometimes the number of maxima detected is not clear (typically at small angle differences). Therefore, the new metrics are needed to interpret HARDI results. [START_REF] Jeurissen | Investigating the Prevalence of Complex Fiber Configurations in White Matter Tissue with Diffusion Magnetic Resonance Imaging[END_REF] More recently a measure called multi-directional anisotropy (MDA) metric was proposed in [START_REF] Tan | Multi-Directional Anisotropy from Diffusion Orientation Distribution Functions[END_REF]. In multi-fiber ODFs, MDA may alternatively be defined for one of the many possible ODF maxima, however MDA requires the prior knowledge of the number of maxima detected. Also in [START_REF] Schwab | [END_REF], the invariant scalar features can be extracted from any spherical function written in an SH basis. However, the reconstruction method that represents HARDI signals in terms of other spherical basis is also proposed in [START_REF] Kezele | Spherical Wavelet Transform for ODF Sharpening[END_REF]][Michailovich et al., 2010b]. These measures will bring additional information on the angular complexity of diffusion functions, complementary to the indices already in use, based on tensor or on diffusion anisotropy.

In this context, it would be more direct and effective to compute the deviation from a single fiber bundle. The motivation of this paper is to propose a new metric describe fiber complexity, which can be used to provide more specific and more readily results for analyzing fiber complexity intravoxel. With a fast visual inspection, our method allows the quantification of fiber complexity in cerebral regions from real data. The prevalence of complex fiber configurations in white matter tissue with HARDI will be better investigated by our new metric. This can potentially be of great help to clinicians and neurosurgeons [Descoteaux, 2006].

The idea is to compute the deviation from a peanut representing a single fiber or fiber bundle in HARDI by comparing the ODFs under investigation with their ideal single fiber (for simplicity, we will call a fiber bundle a fiber) ODFs, thus leading to a so-called PEAM (PEAnut Metric). More precisely, the proposed method consists of first projecting ODF on a two-dimensional (2D) angle-distance plane, and then computing the degree of fitting of the angle-distance point map to a curve representing a single fiber.

The rest of the paper is organized as follows. After a brief recall of existing anisotropy quantities, we formulate in detail the proposed PEAM metric in section 6.2. Section 0 firstly a brief introduction of the methods of fiber configuration classification except anisotropy measures, then presents the application of complexity metrics in classifying fiber configurations and results, followed by Discussion and Conclusion in Sections 6.4 and 0, respectively.

Methodology

To derive the PEAM metric, we first project ODF on a 2D angle-distance plane, then compute the deviation of such projection from a curve representing a single fiber by calculating the sum of squared errors, and finally take this deviation reflecting the degree of fitting as PEAM metric.

Brief Recall of Existing Anisotropy Metrics 6.2.1 FA is defined by [Basser, 1996]

      2 2 2 1 2 3 2 2 2 1 2 3 ˆˆ3 2 FA                 , (6.1)
where  designates tensor eigenvalue and  the mean diffusivity of the diffusion tensor with  

1 2 3 ˆ/3       
. FA ranges from 0 to 1. When FA=0, the diffusion tensor takes the shape of a sphere with 1 2 3

    and it concerns an isotropic diffusion. When FA>0, the diffusion tensor is an ellipsoid and it concerns an anisotropic diffusion. So, FA is a metric that quantifies the deviation from spherical diffusion. The definition of GA in HARDI is given by [Ozarslan, 2005] (

V) 1 1 1 (250 ) GA V    , (6.2)
where the exponent  

V  is defined as   1 1 1 5000 V V    (6.3) with 2 2 2 (i) 9 (1) N ODF i ODF c V c     , (6.4)
where N denotes the number of ODF reconstruction directions, ODF coefficients ODF c are computed using Funk-Radon transform and Funk-Hecke theorem [START_REF] Hess | [END_REF]Descoteaux, 2007]. The coefficient (1) ODF c in the denominator of V is that of the spherical harmonic function of degree 0 l  . When GA=0, the ODF takes the shape of a sphere with 2 2 (i) 0

N ODF i c   
and it concerns an isotropic diffusion. When GA>0, it concerns an anisotropic diffusion. So, GA is a metric that quantifies the deviation from spherical diffusion.

GFA is expressed as [START_REF] Ds Tuch | [END_REF]]

        2 1 2 1 () 1 ( ) n i i n i i nu std GFA rms nu            , (6.5)
where  is the mean of the ODF. Like FA in DTI, GFA is normalized to [0, 1]. When GFA=0, the ODF takes the shape of a sphere with () i u   , indicating an isotropic diffusion. When GFA>0, it concerns an anisotropic diffusion. So, GFA is also a metric that quantifies the deviation from spherical diffusion.

Projection of ODFs on the Angle-Distance plane 6.2.2

The 3D surface representation of ODFs was addressed in [ [START_REF] Ds Tuch | [END_REF]Wedeen, 2005;Descoteaux, 2007], with a variant reported in [Aganj, 2010]. It is often expressed in spherical coordinates. To better visualize and handle the ODFs, we regard them as point clouds, and each ODF is composed of a cloud of points. We then express each point in Cartesian coordinates as a vector ( , , ) T d x y z du  . Thus, the 3D point cloud can be described by a set of vectors. A vector being defined by radial distance and orientation, we use them to characterize 3D point clouds. If m d represents the max distance with m u its direction, we will then take m d as a reference vector. Using the max distance ( m d ) for normalizing d , we project the 3D point cloud onto an angle-distance plane (Fig. 6.1). The cloud points can then be determined by their angle with the reference vector and distance to the origin. To achieve such projection, we first calculate the angle between each vector d and the reference vector In the following subsection, we will derive the analytical expression for the projections of a single fiber ODF in the angle-distance plane.

Analytical expression of single fiber ODF in the Angle-Distance plane 6.2.3 An ODF for Gaussian diffusion is initially defined by [ [START_REF] Ds Tuch | [END_REF]

1 1 ( ) / T u u u Z     D , (6.6)
where u is the unit vector denoting the direction of ODF,  the diffusion time, Z the normalization constant, and D the apparent diffusion tensor. () u  denotes the value of ODF in the direction u , which also represents the length q of the vector q . It can be geometrically represented by a peanut (Fig. 6  designate the three eigenvalues of the diffusion tensor D,  is the angle with the reference vector m q (without loss of generality, we take the direction m q as the direction of z-axis), according to the crossing angle varying from 0 o to 90 o , we normalized  to a range from 0 to pi/2 (1.57),  is the angle between the positive x-axis and the projection in the x-y plane of the vector q .

If we use the uniaxial model [START_REF] Wiegell | Fiber Crossing in Human Brain Depicted with Diffusion Tensor MR Imaging[END_REF] to represent unidirectional or single fibers, (6.7) reduces to

zz xx yy     [

   

). In (6.9), ()

f  represents the radial distance to the origin. So, (6.9) depicts a curve in the angle- distance plane.

Intuitively, if the projections of an ODF in the angle-distance plane coincide with the curve () f  , it will concern a single fiber. The more the point cloud projections deviate from the single fiber curve () f  , the more it would concern complex fiber configurations. We can then calculate a PEAM metric that reflects such derivation by fitting the projections to the single fiver curve via a non-linear least square curve fitting, where

T is normalized to [0,1], N denotes the number of ODF reconstruction directions, i T the distance to the origin of ith 3D point, and i  the angle with the reference vector of ith 3D point. The PEAM value ranges from 0 to N , if i T can perfectly fit to () i f  , PEAM will be reduced to 0. For a fixed number of ODF reconstruction directions (N), the maximum of PEAM will be smaller than N at the condition of no i T can be fitted to () i f  and all of i T is equal to maximum of T . The value of PEAM is therefore the squared error of curve fitting, and reflects the degree of fitting. In other words, PEAM can be seen as the difference between an ODF and its ideal single fiber model.

Application of PEAM to the classification of intravoxel fiber configurations

In what follows, we apply the above-described PEAM to classify intravoxel fiber orientation configurations into two classes: complex configuration (multiple orientations such as fiber crossing, which are characterized by non-Gaussian diffusion) and parallel configuration (single orientation characterized by Gaussian diffusion). To do that, the ODFs are first reconstructed from HARDI acquisitions. Then the metrics such as PEAM are calculated from the ODFs, thus generating parameter maps whose value at a given voxel is equal to the value of the given metric at that voxel. Finally, applying the automatic Otsu thresholding method [Otsu, 1979] to the parameter map leads to classifying intravoxel fiber configurations into the two classes.

The investigation of intravoxel fiber configurations is an important issue in the study of microstructures of tissues in dMRI. A number of methods were reported in the literature on this matter. The authors of [Westin et al., 1997] used a linear diffusion metric defined by where the three eigenvalues are assumed to satisfy 1 2 3    .

l C describes how close the tensor is to the generic cases of line, and hence can be used for classification of the diffusion tensor according to its geometry [Westin et al., 2002]. In [START_REF] Hosey | Inference of Multiple Fiber Orientations in High Angular Resolution Diffusion Imaging[END_REF], a Markov chain method was introduced to determine the number of fibers within a voxel. The authors of [START_REF] Behrens | Probabilistic Diffusion Tractography with Multiple Fibre Orientations: What Can We Gain?[END_REF] used Bayesian modeling to estimate complex architecture that is used to perform subsequent tractography. The above-mentioned methods are by essence only the methods for fiber classification. In contrast, the proposed PEAM can be used both for classifying diffusion complexity and for distinguishing complex configuration and parallel configuration.

To evaluate the proposed PEAM metric, simulated data, physical phantom data and real brain data were used. In the experiments on simulated data, the PEAM was compared with FA, C l , GA and GFA while in the experiments on real phantom and brain data, the PEAM was compared with GA and GFA.

The simulated diffusion signal was obtained using the multi-tensor model,

  exp( bg ) T i k i k i k S g P g   D , (6.12)
where k P is the apparent volume fraction of the voxel with diffusion tensor k D , b the diffusion sensitization factor, i g the direction of the diffusion gradient [ DS Tuch, 2002;Kingsley, 2006]. To generate different crossing angles, we rotated the following diffusion tensor [Kingsley, 2006] The simulated ODFs are obtained from the diffusion signals simulated by (6.12). Different representative fiber configurations were simulated: parallel fibers and crossed fibers at the angles of 90°, 70° 50° and 30°. Each pixel presents two main intravoxel crossing directions. 3 fibers are crossing at 90 o . 2 fibers are crossing at 90 o , but the ratio of volume fractions is 1:2. Different types of diffusion (spherical diffusion, oblate diffusion, single direction diffusion, etc.) signals with different b-values were also simulated and simulations with different tensors (diag([1.7,1,1]), diag([1.7,1,0.9]), diag ([1.7,1,0.7]), and diag ([1.7,1,0.5] with "diag" indicating a diagonal diffusion tensor) are also performed.

The proposed metric PEAM was also evaluated on simulation data at different noise levels. The signal-to-noise ratio (SNR) in the simulations was defined as the ratio of non diffusion weighted signal intensity to noise standard deviation.

To assess true positive rate (TPR) and false positive rate (FPR) of the classification, we simulated 500 diffusion signals of parallel configurations and other 500 diffusion signals of complex fiber configurations that are 2-fiber crossing from 20 o to 90 o at different noise levels. The receiver operating characteristic (ROC) is then calculated, and the performance of the metrics is measured by the area under the curve (AUC) of ROC. The greater AUC is, the better performance is [Fawcett, 2006].

The classification accuracy is measured using the following formula where c N represents the number of correctly classified voxels and t N the total number of voxels. The greater the value of Accuracy, the better the classification.

The physical phantom data concerns the diffusion-weighted images of the physical "fiber cup" phantom from [Fillard, 2011]. The data has a spatial resolution of 3 mm x 3 mm x 3 mm, a b-value of 2000s/mm 2 and a matrix size of 64 64

 . The real brain data concerns the brain hemisphere of a macaque. The data was acquired on a Bruker 7T scanner [Welsh, 2013]. The brain was covered with 1.5 mm slices of size 128 128 

. The diffusion encoding was performed in 95 directions with b-value of 5000s/mm 2 . Since we do not know ground truth in real brain data cases, to still be able to get further insights into intravoxel fiber configurations, we used two ODF reconstruction methods, analytical q-ball imaging (AQBI) [Descoteaux, 2007] and q-ball imaging within constant solid angle (CSA) [Aganj, 2010]. With AQBI, we used: spherical harmonic order max 8 l  , diffusion directions=81 uniformly on a hemi-sphere, number of reconstruction points n=642 uniformly on a sphere. With CSA, we used the spherical harmonic order max 4 l  and the same number of reconstruction points as in the case of CSA.

Results on simulated data 6.3.1 Fig. 6.1 shows the ODFs of a typical one single fiber (1-fiber) and two fibers (2-fiber) crossed at 90°, their 3D point clouds, their projections (in blue) in the angle-distance plane, and the corresponding fitting curves (in red). We observe that the projection points of the 1fiber ODF (i.e. a peanut shaped ODF) coincide almost perfectly (fitting error=0.02, so approximately zero) with the fitting curve while the projection points of the 2-fiber ODF are more scattered and deviate clearly (fitting error=16.57) from the fitting curve in red (Fig. 6.1, bottom row). By calculating the deviation of a given ODF's projections (in the angle-distance plane) from that of a 1-fiber ODF's projections, we can then obtain the PEAM value as expressed by (6.10). Fig. 6.2 shows the ODFs of different 2-fiber systems (crossed respectively at 30°, 50° and 70°) as well as their PEAM values reflecting the fitting error between peanut ODF (1-fiber) and 2-fiber ODF. For the 30 o crossing ODF, the PEAM is 0.37; it is 18 times greater than the PEAM value of 1-fiber ODF. For 50 o crossing, the PEAM increases to 1.82. For 70 o and 90 o crossings, the PEAM further increases to 14.95 and 16.57 respectively. Therefore, the greater crossing angle is, the greater the PEAM is. In Fig. 6.3, a 3-fiber ODF crossing at 90 o is shown. The PEAM becomes greater than 2fiber ODF crossing at 90 o . The result shows that because there are more points that cannot fit to single fiber curve, the 3-fiber ODF crossing at 90 o is more complex than 2-fiber ODF crossing at 90 o . In Fig. 6.4, an ODF crossing at 90 o with the ratio of volume fractions=1:2 is shown. The PEAM becomes smaller than normal 2-fiber ODF crossing at 90 o with equal volume fractions. The result show that because the points can be fitted to single fiber curve better, the 2-fiber ODF crossing at 90 o is less complex than 2-fiber ODF crossing at 90 o . It should be noted that the spherical ODF is a special case with its fitting curve that is a line parallel to the horizontal angle axis (Fig. 6.5). Since the fitting coefficients b and c in (6.9) are equal, the fitting error is 0. But this does not mean that it concerns a single fiber voxel. Four 1-fiber ODFs simulated with different b-values are given in Fig. 6.7(a). The results show that the angle-distance points (Fig. 6.7(b)) can be always fitted to the fitting curve (Fig. 6.7(c)) with nearly the same graph, but the distance value tends to decrease when increasing b-values. A greater b-value will make PEAM greater, especially when b-value is greater than 1000s/mm 2 . However, the PEAMs of these four b-values are all small, which implies that PEAM should be small when ground truth is 1-fiber ODF and will not change much when varying b-values. Fig. 6.9 shows the ROC curves obtained when using the five metrics FA, C l , GA, GFA and PEAM to classify whether an ODF has complex fiber configuration or not, and at different SNR levels. At SNR=16, the GA and the proposed PEAM classifiers have nearly the same greatest area under the curve (i.e. AUC=0.98) and therefore the best average performance. The FA and C l classifiers get the smallest AUC (0.91 and 0.92) and therefore the worst average performance. The GFA has an AUC of 0.95 and therefore a medium performance. However, when the SNR is reduced to 10, the proposed PEAM got the highest AUC (0.95) and therefore the best performance. The AUCs for GA and GFA at this noise level are respectively 0.93 and 0.92. FA and C l exhibit the worst performance (0.85 and 0.83, respectively). As SNR continues to be reduced to 5, the proposed PEAM still leads to the highest AUC (0.72) and therefore always the best performance, in comparison with the GA (0.68), GFA (0.64), FA (0.65) and C l (0.65). In terms of the Accuracy criterion, Table 6.1 gives different classification accuracy values of the 5 metrics in different noise level cases. Classification accuracy decreases with the increase of noise for all the 5 metrics. The GA and GFA always lead to greater accuracy than FA and Cl. But, GA is more sensitive to noise; its accuracy decreases by 0.07 with a decrease in SNR from 16 to 10, more than GFA (0.04), FA (0.05), C l (0.04) and PEAM (0.04). As SNR decreases to 5, all the metrics have an accuracy less than 0.7. However, our PEAM always gives the greatest accuracy value whatever the noise level. This means that our metric is better than other metrics not only in accuracy but also in robustness to noise. Results on physical phantom 6.3.2 Since the simulations in the preceding subsection have qualitatively and quantitatively demonstrated that the proposed PEAM metric produces better results than FA, C l , GA and GFA metrics, for the physical phantom data in the present subsection and the real brain data in the next subsection, we will compare the proposed PEAM only with the GA and GFA metrics, which are more HARDI-based metrics.

Fig. 6.10(a) shows the ground truth fiber configuration of the physical phantom (second slice), with the numbers indicating different fiber configuration regions. In Fig. 6.10(b) are shown the results of quantification of intravoxel complexities using GA (left), GFA (middle) and the proposed (right) metrics. For the proposed PEAM, the lighter the voxel is, the less complex the intravoxel fiber configuration. For the other two metrics, darker voxels indicate greater intravoxel complexity. To have a simpler view of the obtained intravoxel complexity maps, we thresholded them to classify the voxels into homogenerous orientation (e.g. 1-fiber bundle) class and heterogeneous orientation (e.g. crossing fibers) class using the Otsu thresholding method [19], which calculates the optimum threshold separating the two classes so that their intra-variance is minimal. The thresholds given by Otsu's method are respectively 0.15, 0.51 and 6.27 for GA, GFA and PEAM. The so classified fiber configurations are given in Fig. 6.10(c).

We observe that in region 1, the classification is the best for the proposed PEAM, because the results are closer to the ground truth (in view of the larger number of complex voxels and their positions). GFA is better than GA in region 1. In regions 4 and 5 of Fig. 6.10(c), neither of GA, GFA and PEAM gave very good results, but GFA and PEAM are more reasonable than GA.

For 1-fiber configurations such as regions 2 and 3 of Fig. 6.10(c), the proposed metric generated the best results by classifying the fewest voxels as complex ones (i.e. non-Gaussian ones). GFA made more misclassifications than the proposed PEAM. GA tended to consider 1fiber configurations as complex ones, which is the case for regions 2 and 3. Results on real brain data 6.3.3 Fig. 6.11 shows two typical ODFs of real brain data, their angle-distance point maps and fitting curves. Visually, the left simple ODF has a typical peanut shape, thus giving strong indication on the presence of parallel fibers in the voxel. Its smaller PEAM value (0.56) tends to support such inference. In contrast, the right ODF has visually much more complex shape. This is confirmed by its much more widely scattered angle-distance points and much greater PEAM value (17.53), which suggests the presence of complex fiber configurations. In Fig. 6.12(a), we show the b0 image of a coronal slice. The ODF maps obtained using AQBI and CSA are respectively given in Fig. 6.12(b) and (c), whose three columns correspond to the three metrics, GA, GFA, and PEAM. With the AQBI method (Fig. 6.12(b)), visually, GA leads to slightly greater parallel fiber region surrounded by complex fiber regions and misclassifies a complex fiber region into the parallel one (top left corner). GFA underestimates the central parallel fiber region, but overestimates the parallel fiber region in the bottom right corner. PEAM gives the closest results to visual inspection in all the regions.

With the CSA method (Fig. 6.12(c)), GA classifies almost all the voxels as the parallel fiber class, which is clearly unreasonable. GFA gives rise to scattered parallel fiber regions; many simple voxels are taken as complex ones and vice versa. PEAM produces almost the same result as in Fig. 6.12(b).

So, the proposed PEAM gives the best classification results in the sense that the results conform more to visual results and generates nearly the same results as those obtained with AQBI and CSA, and this despite the fact that AQBI is less sensitive to fiber crossing. This implies that the use of PEAM can give more accurate classification of fiber configurations.

To better apprehend the classification results given by the three metrics, we come back on ODFs. We zoom the 6 ODFs in the three rectangles in Fig. 6.12 (b) and (c) (each red rectangle contains 2 voxels), and show them in Fig. 6.13. Visually, AQBI shows the following fiber configurations (from left to right): Rectangle 1-complex, parallel; Rectangle 2-parallel, complex; Rectangle 3-complex, complex. CSA leads to the same results as those with AQBI, but with sharper ODFs that facilitate visual inspection particularly at fiber crossing voxels. The GA, GFA and PEAM values of the ODFs at these 6 voxels, calculated using AQBI and CSA, are given in Table 6.2. The corresponding results of voxel classification using the GA, GFA and PEAM metrics as classifiers are provided in Table 6.3. NB: "P" designates the parallel diffusion, "C" the complex diffusion, and "Visual observation" the results visually observed as in Fig. 6.13. The red letters indicate that the metric-based classifications are not consistent with visual observations.

As shown in Table 6.3, the classification results based on PEAM are always consistent to the visual inspection results while GA and GFA either misclassified the parallel fiber voxels into complex fiber voxels, or inversely, and this irrespective of the used AQBI or CSA.

As observed, using the GA as classifier may produce misclassifications (parallel fiber voxels classified into complex fiber voxels, or inversely), but this seems independent of the used ODF reconstruction methods (AQBI or CSA). In contrast, the GFA produces better intravoxel fiber classification for AQBI ODFs than for CSA ODFs. In all cases, the proposed PEAM did not lead to any misclassifications regardless of the ODF reconstruction methods used. Since the ODF is sharper in CSA than in AQBI, this implies that the proposed PEAM can still generate correct classifications even when the maximal values of the ODF are not easy to be distinguished. In other words, the proposed PEAM can more accurately account for change in intravoxel fiber configurations of voxels.

Discussion

The number of ODF reconstruction directions was 642 in the present study. If using less reconstruction directions, the number of projected points in the angle-distance plane will be reduced, and the fitting accuracy may be decreased. If using more reconstruction directions, the number of projected angle-distance points will be increased, and the fitting could be more accurate but will be more time-consuming. Fortunately, in practice, one always use 642 or more reconstruction directions to obtain ODFs [D Tuch, 2004][Descoteaux, 2007]. Therefore, the proposed PEAM metric will remain accurate.

The unidirectional or single fiber model involved in the calculation of PEAM does not depend on acquisition parameters; the fitting (and consequently the PEAM value) changes only when the ODF under investigation changes. In other words, for a given ODF field, we are not fixing a unique single fiber model for all the ODFs, each ODF (at a voxel) is fitted to its own ideal single fiber model defined by the parameters and bc in (6.9).

Although focused on parallel and complex fiber configurations, we can also recognize free diffusion by checking whether the fitting error is 0 and the fitting coefficients are equal.

Any HARDI results other than those from q-ball imaging can also be used as inputs of the proposed metric-based intravoxel classification.

As an illustration, we have used the Otsu thresholding method for automatically classifying intravoxel fiber configurations. Other more sophisticated classification methods can be applied to the metric maps. Furthermore, the use of PEAM for other applications such as segmentation, registration, statistical characterization of regions of interest in the tissue, and comparison between control and diseased patients can also be envisaged.

Conclusion

We have proposed a new metric to analyze and quantify the complexity of ODFs. The proposed PEAM metric is based on computing the deviation of ODFs from a peanut-shaped ODF, thus allowing us more accurately account for changes in shape of ODFs. The results on simulation, physical phantom and real brain data consistently showed that the PEAM is more effective than existing metrics such as FA, C l , GA and GFA for analyzing intravoxel fiber complexity. When applied to the classification of intravoxel fibers, the PEAM provides better differentiation between parallel and complex intravoxel fiber configurations, and the resulting classification is more consistent with visual inspections in comparison with existing metrics.

Part III General

Conclusions

The fiber architecture of tissues is important for understanding the neural function of the brain or the myocardial function of the heart. High angular resolution diffusion imaging (HARDI) signals are widely used to reconstruct intravoxel fiber crossing. In this thesis, orientation distribution functions (ODFs) in HARDI are regarded as 3D point clouds and a novel paradigm is then proposed to assess their characteristics. The proposed paradigm of assessing the characteristics accounts for the morphological characteristics of ODFs. More precisely, three morphological metrics (length ratio, separability and uncertainty) were proposed. The results show that the metrics are consistent with the visual quality of 3D point clouds and enable us to describe quantitatively and accurately the characteristics of the latters. The paradigm provides a new way to quantify the characteristics of the 3D point clouds in resolving fiber crossing, which can be potentially used to quantify the quality of the 3D point clouds. Several HARDI methods have been assessed and their characteristics have been quantified on both simulated and real brain data.

Long acquisition time is the main limitation of HARDI methods. We demonstrated that recovering missing or non-acquired data is possible through interpolation based on Delaney triangulation. When compared with the original AQBI or CS-based method, our approach improves the detections of fiber crossing using a few number of acquisitions. We have also proposed another method based on the CS. The idea of first compressing diffusion signals through sparsifying the coefficients of spherical harmonic decomposition and then resampling gradient directions is useful for recovering diffusion signals in non-acquired directions. The results on both simulated and real signals show that the capability of fiber detection is improved in the situation of fewer diffusion data.

The proposed PEAM metric based on computing the deviation of ODFs from a peanutshaped ODF can be used to analyze and quantify the complexity of ODFs. It accounts for more accurately changes in shape of ODFs. The results show that PEAM is more effective than other existing metrics for the analysis of intravoxel fiber complexity. The classification between parallel and complex intravoxel fiber configurations is better achieved using the proposed PEAM, and the results obtained with our PEAM are more consistent with visual inspections.

Conclusions en français:

L'architecture des fibres des tissus est importante pour la compréhension de la fonction neuronale du cerveau ou de la fonction myocardique du coeur. Les signaux d'imagerie de diffusion à haute résolution angulaire (HARDI) sont largement utilisés pour reconstruire le croisement des fibres à l'intérieur d'un voxel (intravoxel). Dans cette thèse, les fonctions de distribution d'orientation (ODFs) en HARDI sont considérées comme des nuages des points 3D. Une méthodologie est ensuite proposée pour évaluer leurs caractéristiques. La nouvelle méthode prend en compte les caractéristiques morphologiques de ces nuages de points 3D. Plus précisément, trois mesures morphologiques (ratio de longueurs, séparabilité et incertitude) sont proposées. Ces mesures sont compatibles avec l'appréciation visuelle des nuages de points 3D et nous permettent de décrire quantitativement et précisément les caractéristiques de ces derniers. La méthode offre une nouvelle façon de quantifier les caractéristiques des nuages de points 3D pour la résolution du croisement de fibres, qui peut être potentiellement utilisée pour quantifier la qualité des nuages de points 3D. Plusieurs méthodes HARDI ont été évaluées et leurs caractéristiques ont été quantifiées à l'aide des données aussi bien simulées que réelles.

Le temps d'acquisition long est une limitation principale des méthodes HARDI. Nous démontrons que la reconstruction des données manquantes ou non-acquises est possible par interpolation basée sur la triangulation de Delaunay. En comparaison avec la méthode initiale AQBI ou AQBI basée sur le CS, notre méthode améliore la détection des croisements de fibres en utilisant un nombre réduit d'acquisitions. Nous avons également développé une autre méthode basée sur le CS. L'idée de comprimer d'abord les signaux de diffusion en rendant les coefficients de décomposition harmonique sphérique parcimonieux et de rééchantillonner ensuite les directions de gradients s'avère efficace pour reconstruire des signaux de diffusion dans des directions non-acquises. Les résultats sur des données simulées et réelles montrent que la capacité de détection de fibres est améliorée dans le cas des données de diffusion réduites.

La mesure PEAM proposée est basée sur le calcul de l'écart par rapport à des ODFs en forme de cacahuète et peut être utilisée pour l'analyse et la quantification de la complexité des ODFs. Elle prend en compte de manière précise les changements en forme des ODFs. Les résultats montrent que PEAM est plus efficace que d'autres mesures existantes pour l'analyse de la complexité des fibres introvoxel. La classification entre la configuration de fibres intravoxel parallèles et la configuration de fibres intravoxel complexes est meilleure avec notre PEAM, et les résultats obtenus avec notre PEAM sont plus conformes à des inspections visuelles.

Chapter 7 Conclusions and perspectives Changyu SUN 117

Perspectives

In this thesis, we have addressed intravoxel fiber crossing by assessing the characteristics of ODFs, reducing gradients directions, and describing intravoxel fiber complexity. In the future, our goal would apply the proposed methods to the myocardium. However, a main question remains about whether complex fiber configurations (such as fiber crossing) exist in the myocardium. To answer this question, we should make experiments to answer for example the following questions in the context of the heart organ:

-What does the fiber crossing mean? -What is the difference between fiber crossing in the myocardium and that in the brain? -What is the influence of the observation scale on perceptible fiber crossing? -What are appropriates b-values for better reconstruction of fiber configurations in the heart? -How many directions of diffusion gradients are necessary for reconstructing sufficiently high quality ODFs (to image eventual fiber crossing)? -What are the best (if exist) HARDI methods for fiber reconstruction in the heart? -What is the relationship between crossing voxels and neighboring voxels in the directions of fibers? Ce travail de thèse porte sur la reconstruction et la description des fonctions de distribution d'orientation (ODF) en imagerie de diffusion à haute résolution angulaire (HARDI) telle que l'imagerie par q-ball (QBI). Dans ce domaine, la fonction de distribution d'orientation (ODF) en QBI est largement utilisée pour étudier le problème de configuration complexe des fibres. Toutefois, jusqu'à présent, l'évaluation des caractéristiques ou de la qualité des ODFs reste essentiellement visuelle et qualitative, bien que l'utilisation de quelques mesures objectives de qualité ait également été reportée dans la littérature, qui sont directement empruntées de la théorie classique de traitement du signal et de l'image. En même temps, l'utilisation appropriée de ces mesures pour la classification des configurations des fibres reste toujours un problème. D'autre part, le QBI a souvent besoin d'un nombre important d'acquisitions pour calculer avec précision les ODFs. Ainsi, la réduction du temps d'acquisition des données QBI est un véritable défi.

Dans ce contexte, nous avons abordé les problèmes de comment reconstruire des ODFs de haute qualité et évaluer leurs caractéristiques. Nous avons proposé un nouveau paradigme permettant de décrire les caractéristiques des ODFs de manière plus quantitative. Il consiste à regarder un ODF comme un nuage général de points tridimensionnels (3D), projeter ce nuage de points 3D sur un plan angle-distance (ADM), construire une matrice angle-distance (ADMAT), et calculer des caractéristiques morphologiques de l'ODF telles que le rapport de longueurs, la séparabilité et l'incertitude. En particulier, une nouvelle métrique, appelé PEAM (PEAnut Metric) et qui est basée sur le calcul de l'écart des ODFs par rapport à l'ODF (représenté par une forme arachide) d'une seule fibre, a été proposée et utilisée pour classifier des configurations intravoxel des fibres. Plusieurs méthodes de reconstruction des ODFs ont également été comparées en utilisant les paramètres proposés. Les résultats ont montré que les caractéristiques du nuage de points 3D peuvent être évaluées d'une manière relativement complète et quantitative. En ce qui concerne la reconstruction de l'ODF de haute qualité avec des données réduites, nous avons proposé deux méthodes. La première est basée sur une interpolation par triangulation de Delaunay et sur des contraintes imposées à la fois dans l'espace-q et dans l'espace spatial. La deuxième méthode combine l'échantillonnage aléatoire des directions de gradient de diffusion, le compressed sensing, l'augmentation de la densité de ré-échantillonnage, et la reconstruction des signaux de diffusion manquants. Les résultats ont montré que les approches de reconstruction des signaux de diffusion manquants proposées nous permettent d'obtenir des ODFs précis à partir d'un nombre relativement faible de signaux de diffusion.

MOTS-CLES : HARDI, q-ball imaging, ODF, nuage de points, PEAM, reconstruction.

  ................................................................................................................... 6.5 CONCLUSION ................................................................................................................... ............................................................................................................... 7.2 PERSPECTIVES ................................................................................................................. AUTHOR'S PUBLICATIONS ............................................................................................................. ..........................................................................................................

Fig. 1

 1 Fig. 1.1 A cut-away view of a scanner with a subject inside the scanner bore, illustrating the major features of a medical grade, closed-MRI scanner. ......................................................................... 12 Fig. 1.2 No external magnetic field. [Edelman et al., 1993] ................................................................... 12 Fig. 1.3 With external magnetic field. [Edelman, 1993] ........................................................................ 13 Fig. 1.4 A 90 o pulse moves the spin into the transverse plane with the flip angle which M has rotated from the z axis. .............................................................................................................................. 13 Fig. 1.5 Left: After the pulse knocks M into the transverse plane, it rotates about the z axis. Right: the receiver coil detects the signal, which decreases over time. .......................................................... 14 Fig. 1.6 To excite a slice of spins in the xy plane, a gradient G z in the z direction is used.[Directors et al., 2006] ........................................................................................................................................ 14 Fig. 1.7 The experiment for frequency and phase encoding can be diagramed as above.[Directors, 2006] ....................................................................................................................................................... 15 Fig. 1.8 Illustration of Spin-echo pulse sequence [Jung et al., 2013] ..................................................... 16 Fig. 1.9 Diffusion. Over time, molecules within gases or liquids will move freely though the medium.This motion is known as diffusion. Shown here are sample random paths that can be taken by molecules within a medium that allows isotropic (the same is in every direction) diffusion. As time passes, the net distance traveled by a molecule increases.[START_REF] Huettel | Functional Magnetic Resonance Imaging[END_REF] ................... 16 Fig. 1.10 Diagram shows the diffusion-driven random trajectory (red line) of a single water molecule during diffusion. The dotted white line (vector r ) represents the molecular displacement during the diffusion time interval, between 1 t0  and 2 t .[Hagmann et al., 2006] ............................ 17

Fig. 1 .

 1 Fig. 1.11 Histogram shows a typical displacement distribution due to diffusion in a one-dimensional model. For each displacement distance r, there is a corresponding probability n/N, which is the proportion of molecules within a voxel that distance within a time interval  . ............................ 17 Fig. 1.12 Left: diffusion in isotropic sample. Right: diffusion in anisotropic sample. In the isotropic case the diffusion is similar in all directions. In the anisotropic case is the diffusion larger in one direction than the others. ................................................................................................................ 18 Fig. 1.13 The Stejskal-Tanner PFG sequence.[Stejskal, 1965[Stejskal, , 1965] ] Diffusion-weighting gradients G are inserted either side of the 180 o pulse of the spin echo sequence to sensitize the protocol to diffusion processes. ........................................................................................................................ 18 Fig. 1.14 Principle of diffusion-weighting preparation: static spins (thick line) are dephased by the first diffusion gradient according to their respective position. 180 o pulse and a second diffusion gradient rephrase those spins. Yet moving spins do not have the same position when the second diffusion gradient is applied, there at echo time all spins are not perfectly in phase, resulting in an echo of smaller magnitude. [Chabert et al., 2007] ......................................................................... 18 Fig. 1.15 Schematic representation of diffusion displacement distribution for diffusion tensor. Ellipsoid are used to represent diffusion displacements. The diffusion is highly anisotropic in fibrous tissues such as white matter and the direction of greatest diffusivity is generally assumed to be parallel to the local direction of white matter. [Jing, 2011] ........................................................... 21 Fig. 2.1 Narrow-pulse approximation PGSE. [Callaghan, 1993] ........................................................... 25 Fig. 2.2 Left: ()

  Fig. 1.11 Histogram shows a typical displacement distribution due to diffusion in a one-dimensional model. For each displacement distance r, there is a corresponding probability n/N, which is the proportion of molecules within a voxel that distance within a time interval  . ............................ 17 Fig. 1.12 Left: diffusion in isotropic sample. Right: diffusion in anisotropic sample. In the isotropic case the diffusion is similar in all directions. In the anisotropic case is the diffusion larger in one direction than the others. ................................................................................................................ 18 Fig. 1.13 The Stejskal-Tanner PFG sequence.[Stejskal, 1965[Stejskal, , 1965] ] Diffusion-weighting gradients G are inserted either side of the 180 o pulse of the spin echo sequence to sensitize the protocol to diffusion processes. ........................................................................................................................ 18 Fig. 1.14 Principle of diffusion-weighting preparation: static spins (thick line) are dephased by the first diffusion gradient according to their respective position. 180 o pulse and a second diffusion gradient rephrase those spins. Yet moving spins do not have the same position when the second diffusion gradient is applied, there at echo time all spins are not perfectly in phase, resulting in an echo of smaller magnitude. [Chabert et al., 2007] ......................................................................... 18 Fig. 1.15 Schematic representation of diffusion displacement distribution for diffusion tensor. Ellipsoid are used to represent diffusion displacements. The diffusion is highly anisotropic in fibrous tissues such as white matter and the direction of greatest diffusivity is generally assumed to be parallel to the local direction of white matter. [Jing, 2011] ........................................................... 21 Fig. 2.1 Narrow-pulse approximation PGSE. [Callaghan, 1993] ........................................................... 25 Fig. 2.2 Left: ()

  Fig. 1.11 Histogram shows a typical displacement distribution due to diffusion in a one-dimensional model. For each displacement distance r, there is a corresponding probability n/N, which is the proportion of molecules within a voxel that distance within a time interval  . ............................ 17 Fig. 1.12 Left: diffusion in isotropic sample. Right: diffusion in anisotropic sample. In the isotropic case the diffusion is similar in all directions. In the anisotropic case is the diffusion larger in one direction than the others. ................................................................................................................ 18 Fig. 1.13 The Stejskal-Tanner PFG sequence.[Stejskal, 1965[Stejskal, , 1965] ] Diffusion-weighting gradients G are inserted either side of the 180 o pulse of the spin echo sequence to sensitize the protocol to diffusion processes. ........................................................................................................................ 18 Fig. 1.14 Principle of diffusion-weighting preparation: static spins (thick line) are dephased by the first diffusion gradient according to their respective position. 180 o pulse and a second diffusion gradient rephrase those spins. Yet moving spins do not have the same position when the second diffusion gradient is applied, there at echo time all spins are not perfectly in phase, resulting in an echo of smaller magnitude. [Chabert et al., 2007] ......................................................................... 18 Fig. 1.15 Schematic representation of diffusion displacement distribution for diffusion tensor. Ellipsoid are used to represent diffusion displacements. The diffusion is highly anisotropic in fibrous tissues such as white matter and the direction of greatest diffusivity is generally assumed to be parallel to the local direction of white matter. [Jing, 2011] ........................................................... 21 Fig. 2.1 Narrow-pulse approximation PGSE. [Callaghan, 1993] ........................................................... 25 Fig. 2.2 Left: ()

  Fig. 1.11 Histogram shows a typical displacement distribution due to diffusion in a one-dimensional model. For each displacement distance r, there is a corresponding probability n/N, which is the proportion of molecules within a voxel that distance within a time interval  . ............................ 17 Fig. 1.12 Left: diffusion in isotropic sample. Right: diffusion in anisotropic sample. In the isotropic case the diffusion is similar in all directions. In the anisotropic case is the diffusion larger in one direction than the others. ................................................................................................................ 18 Fig. 1.13 The Stejskal-Tanner PFG sequence.[Stejskal, 1965[Stejskal, , 1965] ] Diffusion-weighting gradients G are inserted either side of the 180 o pulse of the spin echo sequence to sensitize the protocol to diffusion processes. ........................................................................................................................ 18 Fig. 1.14 Principle of diffusion-weighting preparation: static spins (thick line) are dephased by the first diffusion gradient according to their respective position. 180 o pulse and a second diffusion gradient rephrase those spins. Yet moving spins do not have the same position when the second diffusion gradient is applied, there at echo time all spins are not perfectly in phase, resulting in an echo of smaller magnitude. [Chabert et al., 2007] ......................................................................... 18 Fig. 1.15 Schematic representation of diffusion displacement distribution for diffusion tensor. Ellipsoidare used to represent diffusion displacements. The diffusion is highly anisotropic in fibrous tissues such as white matter and the direction of greatest diffusivity is generally assumed to be parallel to the local direction of white matter.[START_REF] Jing | [END_REF] ........................................................... 21 Fig. 2.1 Narrow-pulse approximation PGSE. [Callaghan, 1993] ........................................................... 25 Fig. 2.2 Left: () Pr , Right:   Eq ........................................................................................................... 26 Fig. 2.3 The procedure of DSI. [Heidi Johansen-Berg, 2009] ................................................................ 27 Fig. 2.4 The procesure of Tuch's QBI .................................................................................................... 28

  Fig.3.1 Examples of 3D point cloud representation. Left: ODF from q-ball imaging[Descoteaux, 2007] Middle: ODF from q-ball imaging within constant solid angle[Aganj, 2010] Right: FOD from Constrained spherical deconvolution[Tournier, 2007]. ................................................................ 38 Fig. 3.2 Flow chart for calculating morphological measures. ................................................................. 40 Fig. 3.3 Projection of 3D points on the angle-distance map. From left to right: ODF, 3D point cloud with the red point indicating the reference vector, illustraton of 3D points in two annuluses, and projection of the 3D points on the ADM. ...................................................................................... 41 Fig. 3.4 Illustration of the definition of the length ratio, separability, and uncertainty metrics in the 3D point cloud (a) and in the schematized ADMAT (b). .................................................................... 42 Fig. 3.5 Definition of the ADMAT. ........................................................................................................ 43 Fig. 3.6 3D point cloud and schematized ADMAT of (a) 90 o crossing fiber, (b) free diffusion, and (c) single fiber. .................................................................................................................................... 44 Fig. 3.7 Example where fiber crossing is difficult to assess with 3D color surface representation (left), but much easier to see with the schematized ADMAT (right). ...................................................... 44 Fig. 3.8 Example where 3-fiber crossing is paired off, and the method of 2-fiber crossing can be used on each pair for calculating length ratio, separability and uncertainty. ......................................... 44 Fig. 3.9(a) Three 3D point clouds presenting different uncertainties. (b) Six 3D point clouds presenting different separabilities. (c) Six 3D point clouds presenting different length ratios. ....................... 48 Fig. 3.10 Examples of one 3D point cloud at noise free (Ground truth) and ten different 3D point clouds under noisy condition at SNR=3. ................................................................................................... 48 Fig. 3.11 Three noise-free 3D point cloud maps obtained using respectively three different reconstruction methods and with b=5000s/mm2 and N=81. (a) AQBI-group 3D point clouds. (b) CSA-group 3D point clouds. (c) CSD-group 3D point clouds. ..................................................... 49 Fig. 3.12 Visualization of different 3D noise-free point clouds obtained using three reconstruction methods and in different conditions. (a) AQBI-group 3D point clouds. (b) CSA-group 3D point clouds. (c) CSD-group 3D point clouds. ........................................................................................ 50 Fig. 3.13 3D point clouds having different length ratios and crossing angles. Top line: AQBI-group 3D point clouds. Middle line: CSA-group 3D point clouds. Bottom line: CSD-group 3D point clouds. b-value=5000s/mm2. ..................................................................................................................... 51 Fig. 3.14 3D point clouds computed from data having SNR=20 (b=5000s/mm2, N=81) using the three reconstruction methods. (a) AQBI-group. (b) CSA-group. (c) CSD-group. .................................. 52 Fig. 3.15 3D point clouds (containing 3 fibers) computed from data noise free (left) and having SNR=10 (right) (b=3000s/mm2, ND=81) using AQBI. ................................................................ 53 Fig. 3.16 3D point clouds of real brain data. Left: AQBI 3D point cloud map. Middle: CSA 3D point cloud maps. Right: CSD 3D point cloud map. The 3D point clouds at each voxel is superimposed on a grayscale background modulated by the GFA at that voxel (black: GFA=0; white: GFA=1). ....................................................................................................................................................... 54 Fig. 4.1 Procedure of recovering missing diffusion signal using the proposed Delaunay triangulation method ........................................................................................................................................... 64 Fig. 4.2 (a) ODF map reconstructed using the complete set of original noise-free diffusion data (K=81). (b) K=30, the ODF map reconstructed in the noise-free condition. Left: the map reconstructed using the original data (K=30). Right: the map reconstructed from the data recovered using our method. (c), (d) and (e) are the ODFs maps reconstructed from the original diffusion signal (left) and after processing by our proposed method (right), at SNR=30, 20 and 10 respectively. .......... 67 Fig. 4.3 Comparison of the original AQBI, SH-based CS and the proposed methods on simulated one single-and two-fiber diffusion data (b=3000s/mm 2 ) in terms of (a) the correct detection rate, (b)

  Fig.3.1 Examples of 3D point cloud representation. Left: ODF from q-ball imaging[Descoteaux, 2007] Middle: ODF from q-ball imaging within constant solid angle[Aganj, 2010] Right: FOD from Constrained spherical deconvolution[Tournier, 2007]. ................................................................ 38 Fig. 3.2 Flow chart for calculating morphological measures. ................................................................. 40 Fig. 3.3 Projection of 3D points on the angle-distance map. From left to right: ODF, 3D point cloud with the red point indicating the reference vector, illustraton of 3D points in two annuluses, and projection of the 3D points on the ADM. ...................................................................................... 41 Fig. 3.4 Illustration of the definition of the length ratio, separability, and uncertainty metrics in the 3D point cloud (a) and in the schematized ADMAT (b). .................................................................... 42 Fig. 3.5 Definition of the ADMAT. ........................................................................................................ 43 Fig. 3.6 3D point cloud and schematized ADMAT of (a) 90 o crossing fiber, (b) free diffusion, and (c) single fiber. .................................................................................................................................... 44 Fig. 3.7 Example where fiber crossing is difficult to assess with 3D color surface representation (left), but much easier to see with the schematized ADMAT (right). ...................................................... 44 Fig. 3.8 Example where 3-fiber crossing is paired off, and the method of 2-fiber crossing can be used on each pair for calculating length ratio, separability and uncertainty. ......................................... 44 Fig. 3.9(a) Three 3D point clouds presenting different uncertainties. (b) Six 3D point clouds presenting different separabilities. (c) Six 3D point clouds presenting different length ratios. ....................... 48 Fig. 3.10 Examples of one 3D point cloud at noise free (Ground truth) and ten different 3D point clouds under noisy condition at SNR=3. ................................................................................................... 48 Fig. 3.11 Three noise-free 3D point cloud maps obtained using respectively three different reconstruction methods and with b=5000s/mm2 and N=81. (a) AQBI-group 3D point clouds. (b) CSA-group 3D point clouds. (c) CSD-group 3D point clouds. ..................................................... 49 Fig. 3.12 Visualization of different 3D noise-free point clouds obtained using three reconstruction methods and in different conditions. (a) AQBI-group 3D point clouds. (b) CSA-group 3D point clouds. (c) CSD-group 3D point clouds. ........................................................................................ 50 Fig. 3.13 3D point clouds having different length ratios and crossing angles. Top line: AQBI-group 3D point clouds. Middle line: CSA-group 3D point clouds. Bottom line: CSD-group 3D point clouds. b-value=5000s/mm2. ..................................................................................................................... 51 Fig. 3.14 3D point clouds computed from data having SNR=20 (b=5000s/mm2, N=81) using the three reconstruction methods. (a) AQBI-group. (b) CSA-group. (c) CSD-group. .................................. 52 Fig. 3.15 3D point clouds (containing 3 fibers) computed from data noise free (left) and having SNR=10 (right) (b=3000s/mm2, ND=81) using AQBI. ................................................................ 53 Fig. 3.16 3D point clouds of real brain data. Left: AQBI 3D point cloud map. Middle: CSA 3D point cloud maps. Right: CSD 3D point cloud map. The 3D point clouds at each voxel is superimposed on a grayscale background modulated by the GFA at that voxel (black: GFA=0; white: GFA=1). ....................................................................................................................................................... 54 Fig. 4.1 Procedure of recovering missing diffusion signal using the proposed Delaunay triangulation method ........................................................................................................................................... 64 Fig. 4.2 (a) ODF map reconstructed using the complete set of original noise-free diffusion data (K=81). (b) K=30, the ODF map reconstructed in the noise-free condition. Left: the map reconstructed using the original data (K=30). Right: the map reconstructed from the data recovered using our method. (c), (d) and (e) are the ODFs maps reconstructed from the original diffusion signal (left) and after processing by our proposed method (right), at SNR=30, 20 and 10 respectively. .......... 67 Fig. 4.3 Comparison of the original AQBI, SH-based CS and the proposed methods on simulated one single-and two-fiber diffusion data (b=3000s/mm 2 ) in terms of (a) the correct detection rate, (b)

  ......................................................................................................... 70 Fig. 4.5 ODF maps of the circled region of the b=2000s/mm 2 phantom dataset, obtained using our method and analytical q-ball imaging. (a) Non-diffusion weighted (b0) image of the physical phantom from the b=2000s/mm 2 dataset. The rectangle indicates the fiber crossing region. (b) ODF maps of the circled region obtained using analytical q-ball imaging with complete 64 directions (left), analytical q-ball imaging with 29 directions (middle), and our method with 29 directions (right). ........................................................................................................................... 71 Fig.

  ........................................................................................................................................ 73 Fig. 5.1 Flow chart for recovering diffusion signal by compressed sensing. .......................................... 80 Fig. 5.2(a) The complete set of 81 geometrically uniform directions on a hemisphere. (b) Left: 21 quasi-uniform directions from 81 geometrically uniform directions. Middle: 29 quasi-uniform directions. Right: 39 quasi-uniform directions. ............................................................................. 82 Fig. 5.3 From left to right: the directions of subset diffusion signals original in red and the directions of recovering signals in black. ........................................................................................................... 83 Fig. 5.4 The top: The 81 directions uniform complete set. The bottom: ODFs obtained by the diffusion signals in the directions of its top subset by AQBI. ....................................................................... 85 Fig. 5.5 The top: The different subsets of directions. The bottom: ODFs obtained by the diffusion signals in the directions of its top subset by our method (from top to below: spherical harmonic order used in compressed sensing matrix l=4,6,8,12 respectively)(a) no-uniform 16 directions come from 81 geometrically uniform directions, (b) quasi-uniform 18 directions come from 81 geometrically uniform directions, (c) simulation diffusion image corresponding to quasi-uniform 18 directions come from electrostatic quasi-uniform 81 directions. .............................................. 86 Fig. 5.6 The top: The different subsets of directions. The bottom: ODFs obtained by the diffusion signals in the directions of its top subset by AQBI directly (Left) and our method with spherical harmonic sensing matrix order l=6 (Right). (a) no-uniform 16 directions come from 81 geometrically uniform directions, (b) quasi-uniform 18 directions come from 81 geometrically uniform directions, (c) simulation diffusion image corresponding to quasi-uniform 18 directions come from electrostatic quasi-uniform 81 directions. .................................................................... 87 Fig. 5.7(a) ODFs map of simulation phantom by whole 81 gradient directions reconstructed directly by

  Fig. 6.9 ROC curves of five metrics-based classifications. From left to right: SNR=16, 10 and 5. ..... 105 Fig. 6.10 Quantification and classification of the intravoxel complexities of a physical Fiber Cup phantom image. Each circled region designates a constant complexity region where each pixel has close complexity. (a) Original phantom image with fiber ground-truth indicating various crossing, splitting and kissing fiber configurations that are numbered. For example, the number 1 represents a crossing region, and the numbers 2 and 3 the parallel fibers regions. (b) Results of intravoxel complexity quantification of the image in (a) using GA (left), GFA (middle) and the proposed (right) metrics. (c) Results of classifying the voxels into two classes: homogenerous orientations (e.g. parallel fibers) and heterogeneous orientations (e.g. crossed fibers). In (b) for the GA and GFA images, darker pixels indicate greater complexity and for PEAM image , lighter pixels indicate greater complexity .In (c), for the GA , GFA and PEAM images, darker pixels indicate greater complexity. In all images, the black regions represent background. ............................... 107 Fig. 6.11 Typical ODFs of real brain data, their angle-distance point maps and fitting curves. Columns 1 and 2 show the feature of parallel fiber configuration, and Columns 3 and 4 the feature of complex fiber configuration. The PEAM values are 0.56 and 17.53 respectively. ...................... 108

  la méthode de reconstruction des signaux de diffusion basée sur le compressed sensing (CS). Sur la gauche : Un nombre réduit des données HARDI acquises (points rouges sur la sphère), leur projection sur une basis harmonique sphérique et la reconstruction de l'ODF. Afin de ré-échantillonner les signaux de diffusion, nous utilisons d'abord le modèle CS et on obtient les coefficients * PD c en base harmonique sphérique. Puis (sur la droite), en multipliant ces coefficients avec une nouvelle matrice de base harmonique sphérique à haute résolution angulaire (RS>P), on obtient le rééchantillonnage des signaux de diffusion (points noirs sur la sphère) et par voie de conséquence une ODF de meilleure qualité (en haut à droite). Chapitre 6 : Représentation et description quantitatives de la complexité des fibres à l'intérieur d'un voxel Dans ce chapitre, deux configurations spatiales des fibres sont considérées : configuration simple avec des fibres en parallèle (à direction unique, décrite par diffusion gaussienne) et configuration complexe des fibres (décrite par des diffusions non-gaussiennes). Nous proposons un indice sur la base de l'ODF, mais pas le relais sur le nombre de maxima d'ODF. En utilisant les informations de l'ODF, les classifications simples ou complexes fibre peuvent être déterminées. La nouvelle idée est que l'on compare l'ODF de fibre unique avec tous ODFs par notre méthode. La première partie : Définition du point carte angle-distance. De gauche à droite : ODF, nuage de points 3D avec le point rouge indiquant le vecteur de référence, angle-distance carte de points en bleu et son ajustement de la courbe en rouge, l'axe désigne l'angle entre la vecteur cible et de référence et un axis dénote la distance entre le point et point de l'origine. L'erreur d'ODF croisement en 90 o est 16.57.

Fig. 1 . 2

 12 Fig. 1.2 No external magnetic field. [Edelmanet al., 1993]

  /www.cis.rit.edu/htbooks/mri/inside.htm While the spins are aligned with the B 0 field, no signal can be detected. A radio frequency (RF) pulse is used to the spins. The net magnetization M can be changed by application, via a coil, of a radiofrequency (RF) pulse applied at the Larmour frequency 0  . The degree of rotation of the axis M is called the flip angle. This angle is determined by the strength and duration of the RF pulse. A 90 o pulse moves the direction of M from the z plane into the transverse (xy) plane as shown in Fig. 1.4.

Fig. 1

 1 Fig. 1.4 A 90 o pulse moves the spin into the transverse plane with the flip angle which M has rotated from the z axis.

Fig. 1 . 5

 15 Fig. 1.5 Left: After the pulse knocks M into the transverse plane, it rotates about the z axis. Right: the receiver coil detects the signal, which decreases over time. http://wikidoc.org/index.php/Basic_MRI_Physics K-space 1.1.2

Fig. 1 . 6

 16 Fig. 1.6 To excite a slice of spins in the xy plane, a gradient G z in the z direction is used.[Directors et al., 2006]

Fig. 1 . 7

 17 Fig. 1.7 The experiment for frequency and phase encoding can be diagramed as above.[Directors, 2006] After slice-selective, as shown in Fig. 1.7, the identified signal   St that multiplied by

Fig. 1 . 8

 18 Fig. 1.8 Illustration of Spin-echo pulse sequence [Jung et al., 2013]

Fig. 1 .

 1 Fig. 1.9 Diffusion. Over time, molecules within gases or liquids will move freely though the medium. This motion is known as diffusion. Shown here are sample random paths that can be taken by molecules within a medium that allows isotropic (the same is in every direction) diffusion. As time passes, the net distance traveled by a molecule increases. [Huettel et al., 2004]

Fig. 1 .

 1 Fig. 1.10 Diagram shows the diffusion-driven random trajectory (red line) of a single water molecule during diffusion. The dotted white line (vector r ) represents the molecular displacement during the diffusion time interval, between 1 t0  and 2 t .[Hagmann et al., 2006]

Fig. 1 .

 1 Fig. 1.11 Histogram shows a typical displacement distribution due to diffusion in a one-dimensional model.For each displacement distance r, there is a corresponding probability n/N, which is the proportion of molecules within a voxel that distance within a time interval  .

Fig. 1 .

 1 Fig. 1.12 Left: diffusion in isotropic sample. Right: diffusion in anisotropic sample. In the isotropic case the diffusion is similar in all directions. In the anisotropic case is the diffusion larger in one direction than the others. http://bmia.bmt.tue.nl/Research/MVIAV/DTI/20020501-DTI/BER03_DTI.htm

Fig. 1 .

 1 Fig. 1.13 The Stejskal-Tanner PFG sequence. [Stejskal, 1965, 1965] Diffusion-weighting gradients G are inserted either side of the 180 o pulse of the spin echo sequence to sensitize the protocol to diffusion processes.

Fig. 1 .

 1 Fig. 1.14 Principle of diffusion-weighting preparation: static spins (thick line) are dephased by the first diffusion gradient according to their respective position. 180 o pulse and a second diffusion gradient rephrase those spins. Yet moving spins do not have the same position when the second diffusion gradient is applied, there at echo time all spins are not perfectly in phase, resulting in an echo of smaller magnitude. [Chabert et al., 2007]

S

  is the diffusion weighted signal observed at one of the diffusion direction, tensor is a symmetric 33  matrix with six unique elements,

Fig. 1 .

 1 Fig. 1.15 Schematic representation of diffusion displacement distribution for diffusion tensor. Ellipsoid are used to represent diffusion displacements. The diffusion is highly anisotropic in fibrous tissues such as whitematter and the direction of greatest diffusivity is generally assumed to be parallel to the local direction of white matter.[START_REF] Jing | [END_REF]

  the gyromagnetic ratio and G the imaging field gradient. In (2.1),   Sk is the time domain signal and   r is the spin density.

Fig. 2 .

 2 Fig.2.1 Narrow-pulse approximation PGSE.[Callaghan, 1993]

Fig. 2 . 2

 22 Fig. 2.2 Left: () Pr , Right:   Eq Fig. 2.2 show an example of diffusion signal  

Fig. 2 .

 2 Fig. 2.3 The procedure of DSI. [Heidi Johansen-Berg, 2009]

Fig. 2 .

 2 Fig. 2.4 The procesure of Tuch's QBI

Fig. 3 . 1

 31 Fig. 3.1 Examples of 3D point cloud representation. Left: ODF from q-ball imaging [Descoteaux, 2007] Middle: ODF from q-ball imaging within constant solid angle [Aganj, 2010] Right: FOD from Constrained spherical deconvolution [Tournier, 2007].

Fig. 3 . 2

 32 Fig. 3.2 Flow chart for calculating morphological measures. Construction of the Angle-Distance Map (ADM) 3.2.1

N

  designates the number of partitions in the radial direction and a N the number of partitions of the angle range. A small area is then the intersection of an annulus and a fan sector. The annulus is determined by  

Fig. 3 . 3

 33 Fig. 3.3 Projection of 3D points on the angle-distance map. From left to right: ODF, 3D point cloud with the red point indicating the reference vector, illustraton of 3D points in two annuluses, and projection of the 3D points on the ADM.

  of the definition of the length ratio, separability, and uncertainty metrics in the 3D point cloud (a) and in the schematized ADMAT (b). Construction of the Angle-Distance Matrix (ADMAT) 3.2.2

Fig. 3 . 5

 35 Fig. 3.5 Definition of the ADMAT.

  Fig. 3.6 3D point cloud and schematized ADMAT of (a) 90 o crossing fiber, (b) free diffusion, and (c) single fiber.

Fig. 3 . 8

 38 Fig. 3.8 Example where 3-fiber crossing is paired off, and the method of 2-fiber crossing can be used on each pair for calculating length ratio, separability and uncertainty.

  at k th line, , kp v  , this means that either the 3D point cloud contains only one main direction or it cannot distinguish different directions. Otherwise, at k th line, the 3D point cloud contains   k Card v main directions. For example, in case there are two directions, we can determine the number c N   of lines at which two directions can be resolved

S

  is the pure signal without noise and kSN is the signal with rician noise. To analyze 3D point clouds, we used the length ratio, separability and uncertainty. The experiments were performed in various conditions: with and without noise, 4 different crossing angles, 3 b-values, 2 spherical sampling densities (ND), and 3 SNRs. In the noisy cases, we calculated the mean value of the metrics (length ratio, separability, uncertainty), which was obtained by averaging 25 measures over a neighborhood of size 5x5, in which the 3D point clouds have the same crossing angle.

  Fig. 3.10 Examples of one 3D point cloud at noise free (Ground truth) and ten different 3D point clouds under noisy condition at SNR=3. Table 3.2 Measurement of different shape comparison metrics of the 3D point clouds in Fig. 3.10 and the quality of 3D point cloud sorted by metrics.

Fig. 3 .

 3 Fig. 3.11 Three noise-free 3D point cloud maps obtained using respectively three different reconstruction methods and with b=5000s/mm2 and N=81. (a) AQBI-group 3D point clouds. (b) CSA-group 3D point clouds. (c) CSD-group 3D point clouds.

Fig. 3 .

 3 Fig. 3.13 3D point clouds having different length ratios and crossing angles. Top line: AQBI-group 3D point clouds. Middle line: CSA-group 3D point clouds. Bottom line: CSD-group 3D point clouds. b-value=5000s/mm2.

  point clouds computed from data having SNR=20 (b=5000s/mm2, N=81) using the three reconstruction methods. (a) AQBI-group. (b) CSA-group. (c) CSD-group.

Fig. 3 .

 3 Fig. 3.15 3D point clouds (containing 3 fibers) computed from data noise free (left) and having SNR=10 (right) (b=3000s/mm2, ND=81) using AQBI.

  Fig. 3.16 3D point clouds of real brain data. Left: AQBI 3D point cloud map. Middle: CSA 3D point cloud maps. Right: CSD 3D point cloud map. The 3D point clouds at each voxel is superimposed on a grayscale background modulated by the GFA at that voxel (black: GFA=0; white: GFA=1).

  an angular resolution of 10 o

  the left ofFig. 4.1(b). Based on the fact that PDF should be real and positive, taking real part and absolute value of PDF data.

  : x-axis is the direction number of original diffusion signal, y-axis is the ODF value at each direction Middle: inverse Fourier transform of a diffusion signal in 1-D. Right: the Delaunay triangulation constructed by the direction of original diffusion signal and PDF(the length of the red line denotes the strength of the original signal). (b) Left: the red points represent the original PDF data (obtained after inverse Fourier transform) and the blue points the interpolated PDF data points using our method. Right column: Fourier transform of the entire signal composed of original and interpolated PDF data points). (c) Left: the red points denote the original diffusion data and blue ones are the recovered signal (weighting coefficient of the recovered signal is 1/3). Right: representation of the original and recovered diffusion data on a sphere. (d) Left: ODF obtained by analytical q-ball imaging. (Original 29 directions) Middle: ODF obtained after using our proposed method. (Interpolating to 50 directions) Right: ODF obtained by analytical q-ball imaging with 81 orientations as ground truth. Fig. 4.1 Procedure of recovering missing diffusion signal using the proposed Delaunay triangulation method

  triangulation. Here, we have used the nearest-neighbor interpolation based on Delaunay triangulation.Fourier transforming the absolute value of the real part of the PDF()

  Fig. 4.2 (a) ODF map reconstructed using the complete set of original noise-free diffusion data (K=81). (b) K=30, the ODF map reconstructed in the noise-free condition. Left: the map reconstructed using the original data (K=30). Right: the map reconstructed from the data recovered using our method. (c), (d) and (e) are the ODFs maps reconstructed from the original diffusion signal (left) and after processing by our proposed method (right), at SNR=30, 20 and 10 respectively. In Fig. 4.2, the ODFs reconstructed respectively directly by AQBI and after processing by our method are shown. The ODFs with 81 orientations of AQBI are used as the ground truth (Fig. 4.2(a)) for the comparison with the results in Fig. 4.2(b) to 2(e) where a subset of 30 directions was used and different noise levels were considered.We observe that in the noise-free case, our method presents better results even than the complete 81 direction set case, by giving visually clearer fiber detection. With the increase of noise levels, our results are still better than original data results, but the quality of angular resolution is reduced both for original data results and for our method's results. Nevertheless, at SNR=10, our method still presents better results than original data results. The clearer the fiber crossing, the higher fiber detection capability is. In this sense, our method improves the capability of fiber detection.

Fig. 4 . 3

 43 Fig. 4.3 Comparison of the original AQBI, SH-based CS and the proposed methods on simulated one singleand two-fiber diffusion data (b=3000s/mm 2 ) in terms of (a) the correct detection rate, (b) the average length ratio, (c) the average separability, (d) the average uncertainty, and (e) the average angular error. From left to right columns: SNR=30, 20 and 10 dB.

Fig. 4 . 4

 44 Fig. 4.4 Overview of the diffusion signal recovering method based on Delaunay triangulation. A reduced HARDI dataset is used to reconstruct the phantom ("fiber cup") data. ODF on the left: constructed using the original q-ball imaging. ODF on the right: constructed after processing using the proposed approach.

  Fig. 4.4 are shown the recovered diffusion-weighted images. Using the original signals and the recovered signals together, the reconstructed ODF (bottom right of Fig. 4.4) is sharper than the original ODF obtained using directly AQBI. This implies that the ODF reconstructed after recovering missing diffusion data can allow us to better resolve fiber crossing.

  Fig. 4.5(a), we show a b0 image with a circled region containing fiber crossings. To analyze in detail the quality of ODFs and the performance of the proposed method, we used the b=2000s/mm 2 phantom dataset, a set of 64 orientations as complete orientation set for comparison, and a subset of 29 orientations for reconstruction with AQBI and our method combined with AQBI.In the fiber crossing region, we observe that analytical q-ball imaging using complete 64 directions cannot produce clearly resolved ODFs. Using 29 directions, the simple AQBI reconstruction do not allow all the fiber crossings to be brought out. Always using 29 directions, but with the proposed method, nearly all the fiver crossings come out.Moreover, visually, the ODFs obtained with the proposed method look even better than the ODFs reconstructed from the full 64 directions. The ODFs obtained using the complete 64 directions are of medium quality. Using only 29 directions, without missing diffusion data recovering by our method, the resulting ODFs are of poor quality and cannot be accepted.

  maps of the circled fiber crossing region of the b=2650s/mm 2 phantom dataset, obtained using our method and analytical q-ball imaging. (a) b0 image of the physical phantom from the b=2650s/mm 2 dataset. The circled region is a fiber crossing region. (b) ODF maps of the circled region obtained using analytical q-ball imaging with complete 64 directions (left), analytical q-ball imaging with 29 directions (middle), and our method with 29 directions (right).

Fig. 4

 4 Fig. 4.7 ODF maps of real brain data, obtained using analytical q-ball imaging and our diffusion data recovering combined with analytical q-ball imaging. Left: the analytical q-ball imaging results. Middle: b0 image of the brain with the red box indicating the region where ODFs are reconstructed. Right: our method's results. The ODFs at each voxel is superimposed on a grayscale background modulated by the GFA [DSTuch, 2004] at that voxel (black: GFA=0; white: GFA=1).

Fig. 4 .

 4 Fig. 4.7 shows the ODF maps of the real brain data, obtained using the AQBI (left column) and our diffusion data recovering combined with AQBI (right column).From these ODF maps, we observe that the proposed method generates clearer fiber crossing ODFs in the region of fiber crossing; the ODFs have more distinguishable ODF maxima important for fiber crossing characterization.

Fig. 5 .

 5 Fig. 5.1 Flow chart for recovering diffusion signal by compressed sensing. Subsets of diffusion signals 5.2.1

  complete set of 81 geometrically uniform directions on a hemisphere. (b) Left: 21 quasiuniform directions from 81 geometrically uniform directions. Middle: 29 quasi-uniform directions. Right: 39 quasi-uniform directions. Recovering missing diffusion signals 5.2.2Once the diffusion signal defined in the above randomly generated subset (

Fig. 5 . 3

 53 Fig. 5.3 From left to right: the directions of subset diffusion signals original in red and the directions of recovering signals in black.The new coefficient representation of diffusion signals obtained by our method are calculated using  

  The ODF coefficients ODF c are finally computed using Funk-Radon transform and Funk-Hecke theorem. They are calculated by

Fig. 5 .Fig. 5 . 5

 555 Fig. 5.4 The top: The 81 directions uniform complete set. The bottom: ODFs obtained by the diffusion signals in the directions of its top subset by AQBI.

Fig. 5 . 6

 56 Fig. 5.6 The top: The different subsets of directions. The bottom: ODFs obtained by the diffusion signals in the directions of its top subset by AQBI directly (Left) and our method with spherical harmonic sensing matrix order l=6 (Right). (a) no-uniform 16 directions come from 81 geometrically uniform directions, (b) quasi-uniform 18 directions come from 81 geometrically uniform directions, (c) simulation diffusion image corresponding to quasi-uniform 18 directions come from electrostatic quasi-uniform 81 directions.

2

 2

  Fig. 5.7(a) ODFs map of simulation phantom by whole 81 gradient directions reconstructed directly by AQBI, (b) the correct detection rates of a single and two-fiber diffusion pattern by our proposed method in red, the directly using q-ball imaging method in blue and the spherical harmonic based CS method in green

Fig. 5 .

 5 Fig. 5.8(a) Left: b0 image of a coronal view of the macaque brain, (b) the ODF field of the indicated region recovered by analytical q-ball imaging with sampling density 96, (c) (Upper row of subplots) the ODF maps of the indicated region recovered by analytical q-ball imaging reconstructed by a subset of 19, 23 and 29 encoding directions. (Lower row of subplots) estimated ODF maps obtained using our proposed method on the same conditions. Examples of recovered ODFs marked by red box.

md.

  Since the distance d to the origin is known (i.e. the length of the vector d ), we then know the values of the angle and distance, which enables us to put the cloud points onto the corresponding angle-distance plane. If the maximum distance appears in two or more directions, we can choose any of them as the reference direction.

Fig. 6 . 1

 61 Fig. 6.1 Typical examples of one single fiber ODF and orthogonal crossing fiber ODF. From left to right: ODF, 3D point cloud with the red point indicating the reference vector, angle-distance point map in blue and its fitting curve in red. The axis x denotes the angle between the target vector and reference one and the aixs y the distance between the target point and origin. The fitting errors (i.e. PEAM values) for the single fiber ODF and 90 o crossing fiber ODF are respectively 0.02 and 16.57.

Fig. 6 .

 6 Fig. 6.2 2-fiber ODFs with 30 o , 50 o and 70 o crossing angles, their angle-distance point maps in blue, and the corresponding fitting curves in red. The PEAM values of the three ODFs are respectively 0.37, 1.82 and 14.95.

Fig. 6

 6 Fig. 6.3. 3-fiber 90 o crossing ODF (left), its angle-distance point map in blue and fitting curve in red. The fitting error of the 3-fiber ODF is 21.31.

Fig. 6

 6 Fig. 6.4. 2-fiber 90 o crossing ODF with the ratio of volume fractions=1:2(left), its angle-distance point map in blue and fitting curve in red. The fitting error of this 2-fiber ODF is 2.92.

Fig. 6 . 5

 65 Fig. 6.5 Spherical ODF (left), its angle-distance point map in blue and fitting curve in red. The fitting error of the spherical ODF is 0 because the two coefficients of fitting are equal. In Fig. 6.6, four oblate ODFs corresponding to four b-values are shown. When increasing b-values, the fitting curve approaches more the vertical distance axis, and the PEAM becomes greater.

Fig. 6 .

 6 Fig. 6.6 Oblate ODFs with different b-values (top row and from left to right: b=500, 1000, 3000, 5000s/mm 2 ), and their angle-distance point maps in blue and fitting curves in red (bottom). Their PEAM values are respectively 22.89, 24.48, 30.49 and 32.35.

  Fig. 6.7 (a) 1-fiber ODFs with (from left to right) b=500, 1000, 3000, 5000s/mm2. (b) The corresponding angle-distance point maps and (c) fitting curves. Their PEAM values are respectively 0.10, 0.007, 0.02 and 0.05.The ODFs reconstructed from diffusion signals simulated by different tensors are shown in Fig.6.8(a), together with their angle-distance point maps and fitting curves in Fig.6.8(b) and Fig.6.8(c). When the ODFs become fatter, the angle-distance points become more dispersed. As a result, the fitting curve exhibits greater distance values, and the PEAM value becomes greater. This implies that the proposed PEAM is sensitive to the change in shape (thin or fat) of ODFs.

  Fig. 6.8(a) ODFs simulated by four different tensors. (b) The angle-distance point map and (c) the fitting curves of 1-fiber ODF simulated by four different tensors, from left to right, the tensor is 1 2 3    . The fitting errors of the four ODFs are respectively 0.10, 0.17, 2.75 and 5.33.

Fig. 6

 6 Fig. 6.9 ROC curves of five metrics-based classifications. From left to right: SNR=16, 10 and 5.

  and classification of the intravoxel complexities of a physical Fiber Cup phantom image. Each circled region designates a constant complexity region where each pixel has close complexity. (a) Original phantom image with fiber ground-truth indicating various crossing, splitting and kissing fiber configurations that are numbered. For example, the number 1 represents a crossing region, and the numbers 2 and 3 the parallel fibers regions. (b) Results of intravoxel complexity quantification of the image in (a) using GA (left), GFA (middle) and the proposed (right) metrics. (c) Results of classifying the voxels into two classes: homogenerous orientations (e.g. parallel fibers) and heterogeneous orientations (e.g. crossed fibers). In (b) for the GA and GFA images, darker pixels indicate greater complexity and for PEAM image , lighter pixels indicate greater complexity .In (c), for the GA , GFA and PEAM images, darker pixels indicate greater complexity. In all images, the black regions represent background.

Fig. 6 .

 6 Fig. 6.11 Typical ODFs of real brain data, their angle-distance point maps and fitting curves. Columns 1 and 2 show the feature of parallel fiber configuration, and Columns 3 and 4 the feature of complex fiber configuration. The PEAM values are 0.56 and 17.53 respectively.

  Fig. 6.13 Enlarged view of the six ODFs at the voxels indicated by the red box in Fig. 6.12 (b) and (c). The three numbered rows (1, 2, and 3) correspond respectively to Rectangle 1, Rectangle 2 and Rectangle 3.

  accessible à l'adresse : http://theses.insa-lyon.fr/publication/2014ISAL0119/these.pdf © [C. Sun], [2014], INSA de Lyon, tous droits réservés ............................................................................................................... 7.2 PERSPECTIVES ................................................................................................................. AUTHOR'S PUBLICATIONS .............................................................................................................

  Reconstruction et description des fonctions de distribution d'orientation en imagerie de diffusion à haute résolution angulaire

  Jean-Philippe THIRAN, François ROUSSEAU, Rachid DERICHE , Su RUAN, Nicole VINCENT, Patrick CLARYSSE, Yue-Min ZHU

  INTRODUCTION .................................................................................................................. 6.2 METHODOLOGY ................................................................................................................. Brief Recall of Existing Anisotropy Metrics ................................................................. 6.2.1 Projection of ODFs on the Angle-Distance plane ........................................................ 6.2.2 Analytical expression of single fiber ODF in the Angle-Distance plane ......................

	6.2.3
	6.3	APPLICATION

OF PEAM TO THE CLASSIFICATION OF INTRAVOXEL FIBER CONFIGURATIONS 100 Results on simulated data ........................................................................................... 6.3.1 Results on physical phantom ...................................................................................... 6.3.2 Results on real brain data .......................................................................................... 6.3.

Chapter 1 1 Magnetic Resonance Imaging (MRI) and Diffusion Weighted Magnetic Resonance Imaging (DW-MRI) Contents

  MRI) .......................................................................... 12 Basis of MRI ................................................................................................................. 12 1.1.1 K-space ........................................................................................................................Slice selective technique.................................................................................................... 14 1.1.2.

	1.1 1.1.2	MAGNETIC RESONANCE IMAGING (. 14
	1.1.2.1

2 Frequency and phase encoding technique ......................................................................... 15

Spin

echo ......................................................................................................................

15 1.1.3 1.2 DIFFUSION WEIGHTED IMAGING (DWI) ............................................................................ 16 Diffusion ....................................................................................................................... 16 1.2.1 Principles for DW-MRI ................................................................................................ 17 1.2.2 Apparent diffusion coefficient (ADC) ........................................................................... 19 1.2.3 Diffusion tensor imaging (DTI) .................................................................................... 19 1.2.4 1.2.4.1 DTI model ......................................................................................................................... 20 1.2.4.2 Extraction of diffusion characteristic................................................................................. 20

Chapter 2 2 Fiber Crossing Reconstruction Techniques Contents RÉSUMÉ

  EN FRANÇAIS .................................................................................................................... SPACE ............................................................................................................................ 2.2 DIFFUSION SPECTRUM IMAGING ........................................................................................ ANISOTROPY MEASURES OF HARDI ................................................................................. 2.5 ASSESSMENT OF HARDI ................................................................................................... Angular orientation error ................................................................................................... 2.5.1.2 Detection rates ...............................................

	2.1 Q-2.3 HIGH ANGULAR RESOLUTION DIFFUSION IMAGING (HARDI) ............................................
	2.3.1	Q-ball Imaging .............................................................................................................
	2.3.2	Related work about Q-ball imaging .............................................................................
	2.3.3	Spherical deconvolution ...............................................................................................
	2.3.4	Other HARDI methods .................................................................................................
	2.4
	2.5.1	Fiber orientation assessment ........................................................................................
	2.5.1.1

.................................................................... Shape of HARDI results assessment ............................................................................. 2.5.2 2.5.2.1 Qualitative assessment ...................................................................................................... 2.5.2.2 Quantitative assessment ....................................................................................................

  INTRODUCTION .................................................................................................................. 38

	3.2 3.2.3	3D 45
	3.3	EXPERIMENTS

POINT CLOUD CHARACTERISTIC ASSESSMENT PARADIGM............................................. 40 Construction of the Angle-Distance Map (ADM) .

........................................................ 40 3.2

.1 Construction of the Angle-Distance Matrix (ADMAT) ................................................. 42 3.2.2 Calculation of morphological criteria from ADMAT ................................................... AND RESULTS ............................................................................................. 46 Simulation results ......................................................................................................... 46 3.3.1 Results on real brain data ............................................................................................ 53 3.3.2 3.4 DISCUSSION ....................................................................................................................... 55 3.5 CONCLUSION ..................................................................................................................... 58

Table 3 .3 Quantitative measurements of 3D point clouds using separability and uncertanty in various conditions (length ratio=1).

 3

Table 3 .4 Length ratios of the three groups of 3D point clouds with respect to the ground-truth.

 3

	Crossing angle	90 o		70 o		50 o		30 o
	Ground-truths	0.7	0.4	0.7	0.4	0.7	0.4	0.7	0.4
	AQBI-group	0.6	0.4	0.6	X	X	X	X	X
	CSA-group	0.8	0.7	0.8	0.6	0.8	0.6	0.7	0.6
	CSD-group	0.6	0.4	0.6	0.3	0.5	0.3	X	X

NB: "X" designates denotes that no peaks are separable in the 3D point cloud.

Table 3 .5 Characteristic assessment of the three groups (A,B and C) of point clouds at different SNRs, in terms of length ratio, separability and uncertainty.

 3

	Metric	Angle estimation	AQBI	SNR=30 CSA	CSD	AQBI	SNR=20 CSA	CSD	AQBI	SNR=10 CSA	CSD
		90 o	0.90	0.90	0.90	0.90	0.85	0.90	0.87	0.80	0.88
	Length ratio	70 o 50 o	0.90 X	0.90 0.86	0.90 0.90	0.90 X	0.83 0.83	0.90 0.90	0.84 X	0.80 0.80	0.85 0.83
		30 o	X	0.76	X	X	0.61	X	X	0.60	X
		90 o	0.5	0.7	0.9	0.5	0.6	0.9	0.4	0.4	0.8
	Separability	70 o 50 o	0.2 X	0.6 0.5	0.8 0.6	0.2 X	0.5 0.5	0.8 0.6	0.2 X	0.4 0.3	0.7 0.6
		30 o	X	0.3	X	X	0.3	X	X	0.2	X
		90 o	0.26	0.26	0.26	0.26	0.28	0.26	0.28	0.29	0.26
	Uncertainty	70 o 50 o	0.30 0.70	0.26 0.28	0.26 0.26	0.30 0.70	0.28 0.28	0.26 0.26	0.37 0.70	0.34 0.36	0.26 0.28
		30 o	0.52	0.35	0.52	0.52	0.37	0.48	0.52	0.37	0.44

Table 3 .6 Measurements of length ratio, separability and uncertainty of the point clouds obtained with the three reconstruction methods.

 3

	Methods	Metric	Voxel 1	Voxel 2	Voxel 3	Voxel 4	Voxel 5
		Length ratio	0.90	X	0.90	X	X
	AQBI	Separability Uncertainty	0.15 0.53	X 0.70	0.10 0.52	X 1.57	X 1.22
		Length ratio	0.90	0.75	0.95	0.9	0.95
	CSA	Separability Uncertainty	0.45 0.35	0.30 0.44	0.45 0.26	0.25 0.35	0.20 0.26
		Length ratio	0.85	0.75	0.95	0.85	0.9
	CSD	Separability Uncertainty	0.40 0.44	0.25 0.44	0.40 0.35	0.20 0.44	0.10 0.52
		NB: X denotes that no peaks are separable in the 3D point cloud.

  .................................................................................................................. 62 4.2 METHODOLOGY OF RECOVERING MISSING DIFFUSION SIGNAL DATA ................................. 63

	4.3 4.3.3	EXPERIMENTS . 72
	4.4 DISCUSSION ..Changyu SUN
		60

AND RESULTS ............................................................................................. 66 Simulation experiments and results .

............................................................................. 66 4.3

.1 Results of fiber cup phantom ........................................................................................ 70 4.3.2 Brain data experiments .................................................................................................................................................................................................................... 74 4.5 CONCLUSION ..................................................................................................................... 74 Chapter 4 High Quality ODF Reconstructed by Interpolating Probability Distribution Function

  If the image is contaminated by noise, the type of noise is always Rician noise. The definition of the signal to noise ratio is as

	1 D	3    3 3 1.7 10 0 0 0 0.3 10 0 0 0 0.3 10       	2 / mm s .
	with FA=0.81, MD=	3 0.76 10  	s/mm 2 and 12 0.5 PP 	.
	For reconstruction of ODF, we used spherical harmonic order max 8  l	, reconstruction
	directions=162 uniform directions on a sphere, and diffusion directions=81 uniform directions
	on a hemi-sphere.
	10 SNR 10 log  	 11    22     10 22 11 20 log KK k k k k kk KK kk kk S SN S SN SS               	,
	where

k S is the pure signal without noise and k

  INTRODUCTION .................................................................................................................. 6.2 METHODOLOGY ................................................................................................................. Brief Recall of Existing Anisotropy Metrics ................................................................. 6.2.1 Projection of ODFs on the Angle-Distance plane ........................................................ 6.2.2 Analytical expression of single fiber ODF in the Angle-Distance plane ...................... Results on simulated data ........................................................................................... 6.3.1 Results on physical phantom ...................................................................................... 6.3.2 Results on real brain data .......................................................................................... 6.3.3 6.4 DISCUSSION ..................................................................................................................... 6.5 CONCLUSION ................................................................................................................... Chapter 6 Quantitative Representation and Description of Intravoxel Fiber Complexity in HARDI

	6.2.3
	6.3	APPLICATION OF PEAM TO THE CLASSIFICATION OF INTRAVOXEL FIBER CONFIGURATIONS
		100
	Changyu SUN
		94

   as a fitting function, and put it under the form Chapter 6 Quantitative Representation and Description of Intravoxel Fiber Complexity in HARDI

		f	( ) 		1 m qZ		/ ( sin b 	22 cos ) c    	,	(6.9)
	where	b and c are two constants with c smaller than b (because
										zz	xx	yy
		( )  	1  2 1 sin zz    1 / ( xx  	2 cos ) 	,	(6.8)
					qZ
					m
	We now take ()
	Changyu SUN
										100

Table 6 .1 Classification accuracies of five metrics for the simulation data at SNR=16, 10 and 5.

 6

	Metric	GA	GFA	FA	Cl	PEAM
	SNR
	16	0.93 0.89 0.85 0.85	0.95
	10	0.86 0.85 0.80 0.81	0.91
	5	0.65 0.61 0.63 0.62	0.67

Table 6 .2 The GA, GFA and PEAM values of the ODFs at different voxels and calculated using AQBI and CSA.

 6

	AQBI

Table 6 .

 6 3 gives the fiber classifications of the same voxels as those in Table6.2, obtained using the Otsu thresholding method.

Table 6 .3 Classification results of the same voxels as in Table 6.2 both visually and using the GA, GFA and PEAM metrics

 6

			AQBI
	metrics	Rectangle 1 Rectangle 2 Rectangle 3
	GA	C	C	C	C	P	C
	GFA	C	P	P	C	P	C
	PEAM	C	P	P	C	C	C
	Visual observation	C	P	P	C	C	C
			CSA
	Metrics	Rectangle 1 Rectangle 2 Rectangle 3
	GA	C	P	P	P	P	P
	GFA	C	C	P	P	C	P
	PEAM	C	P	P	C	C	C
	Visual observation	C	P	P	C	C	C

Cette thèse est accessible à l'adresse : http://theses.insa-lyon.fr/publication/2014ISAL0119/these.pdf © [C. Sun], [2014], INSA de Lyon, tous droits réservés

Changyu SUN

(a) (b) (c) (d) (e) Cette thèse est accessible à l'adresse : http://theses.insa-lyon.fr/publication/2014ISAL0119/these.pdf © [C. Sun], [2014], INSA de Lyon, tous droits réservés

(a) (b) (c)