FR 3036189 A1 20161118 FR 1554358 A 20150513 1 « Procédé et dispositif d'imagerie par résonance magnétique » Domaine technique La présente invention concerne un procédé d'imagerie par résonance magnétique (IRM). Elle concerne aussi un dispositif mettant en oeuvre ce procédé. Un tel procédé ou dispositif peut permettre par exemple de préférence à un utilisateur la génération et l'exploitation d'images co-localisées, de même résolution spatiale, de même temps d'écho (TE) et de répétition (TR) pour l'aimantation nucléaire, les temps de relaxation longitudinal (Ti) et transverse (T2), le coefficient de diffusion (D), de l'amplitude du champ radiofréquence (B1) et des inhomogénéités du champ magnétique principal BO. Etat de la technique antérieure L'IRM est une modalité d'imagerie clinique et préclinique utilisée en routine dans le diagnostic médical radiologique ou pour l'évaluation de modèles animaux. Son aspect non invasif et les multiples contrastes disponibles la rendent indispensable dans l'évaluation de nombreuses maladies. Un examen typique d'IRM contient un protocole d'acquisition de plusieurs séquences d'imagerie ayant chacune un objectif différent. Les séquences dites 'en densité de proton' vont renseigner sur la quantité d'eau libre dans le tissu. Les séquences dites 'pondérées-T 1' servent en général à l'évaluation anatomique car elles génèrent un bon contraste entre les tissus biologiques. Les séquences 'pondérées-T1' servent également à visualiser des agents de contraste qui peuvent être injectés en angiographie ou pour la caractérisation tumorale. Les séquences dites `pondérées-T2' ont globalement la propriété de permettre de bien distinguer les contenus liquidiens et renseignent par exemple sur la présence d'oedèmes. Là encore, des agents de contraste dit T2 peuvent être utilisés pour moduler la pondération T2 pour certaines applications. La diffusion, qui influence également le contraste des images générées avec des séquences particulières dite 'pondérées en diffusion' sont également utilisées pour pousser la caractérisation tissulaire, par exemple pour distinguer les zones viables des zones nécrotiques dans les tumeurs, une distinction pas toujours possible avec les autres types de pondération, ou pour faire de la tractographie de fibres cérébrales basée sur la direction privilégiée de la diffusion de l'eau. En fonction de la maladie étudiée et de la zone anatomique, plusieurs séquences sont appliquées lors d'un examen. Celles-ci sont appliquées de manière séquentielle, ce qui nécessité déjà un temps de paramétrage des manipulateurs en électroradiologie pour adapter chaque séquence à la zone étudiée. Pour des contraintes de temps d'acquisition, de type de pondération et localisation, les différentes informations ne sont en général pas acquises avec les mêmes résolutions spatiales, ne couvrent pas les mêmes champs de vue et peuvent même souffrir de déformations spatiales différentes en fonction notamment des temps d'écho et des 3036189 2 bandes passantes d'acquisition différents. On citera par exemple la nécessité dans certains cas d'avoir recourt à une correction des images a posteriori ou à des techniques de recalage d'image pour compenser partiellement ces effets. Ces aspects concourent à rendre l'IRM soit limitée à une interprétation exclusivement basée sur des comparaisons visuelles, soit à introduire des biais ou des imprécisions 5 résultant des algorithmes de correction. Pour améliorer le diagnostic, Il y a donc un besoin de pouvoir co- localiser exactement à l'échelle des mêmes voxels d'imagerie les différentes informations. Par ailleurs, le choix des séquences appliquées et le paramétrage des séquences incombent à chaque centre de radiologie, et les implémentations des séquences peuvent varier entre les constructeurs. Il n'existe pas de standard de séquence et de paramètres associés applicables dans tous les cas. Ces éléments rendent difficiles la 10 comparaison multi-centriques des résultats d'imagerie et ralentissent, voire empêchent, l'établissement de consensus sur les séquences et sur leur paramétrage. Pour fiabiliser le diagnostic, la quantification directe des paramètres physiques tels que l'aimantation nucléaire, les Ti et T2, et le coefficient de diffusion est une possibilité permettant de standardiser les résultats d'imagerie. Là encore, pour quantifier les différents paramètres Ti, T2 et le coefficient de diffusion, de 15 multiples séquences existent, chacune adaptée au paramètre à quantifier. On citera par exemple la séquence d'inversion-récupération avec plusieurs temps d'inversion pour quantifier T 1, la séquence d'écho de spin avec plusieurs temps d'écho de spin pour quantifier T2, la séquence d'écho de spin pondéré en diffusion pour quantifier le coefficient de diffusion. Ces différents motifs de sensibilisation peuvent ensuite être combinés avec différentes séquences permettant la localisation avec tous les 20 désavantages précités de déformation spatiale entre les acquisitions, de champ de vue ou zone couvertes différents, réduisant la possibilité de coupler voxel à voxel les paramètres extraits. Une séquence standard rapide d'acquisition est basée sur l'écho de gradient rapide. Cette séquence d'imagerie consiste à répéter tous les TR un motif incluant une excitation d'amplitude fixe, éventuellement combinée à un gradient dit 'de sélection' pour réaliser une sélection spatiale, suivie de 25 gradients d'imagerie (de lecture et de codages de phase). De multiples variantes de cette séquence ont été proposées. Pour synthétiser les différences entre les variantes, on notera les éléments suivants : la compensation totale des aires de gradient entre les répétitions (séquence dite balancée), ou au contraire la non-compensation de cette aire de gradient (séquence dite non-balancée). Cette dernière variante non-balancée créée des ordres de cohérence différents qui peuvent 30 éventuellement être exploités. Un ordre de cohérence est noté k. C'est un nombre entier qui correspond à l'ordre de la série de Fourier discrète qui peut être utilisée pour décrire les états de configuration de l'aimantation, comme proposé et utilisé dans l'art. 3036189 3 Le cyclage particulier de la phase 0, de l'impulsion radiofréquence à la répétition n sans changer son amplitude avec 0n+1 = On + (n + 1) 4 . Les variantes travaillent soit à phase fixe en fonction des répétitions (4 = 0 ), soit avec un cyclage quadratique de la phase ( 0 ), équivalent à un incrément A constant et non nul d'une répétition à l'autre. A est bien entendu défini modulo 27r 5 La mesure entre les impulsions de 1 ou plusieurs échos, notamment correspondant aux ordres de cohérence de phase (k=0, en général le seul acquis, et/ou k=-1 et/ou k=1 qui peuvent l'être également) L'ajout ou non de gradients supplémentaires entre les impulsions pour sensibiliser de manière différente à la diffusion 10 L'acquisition du signal pendant la phase transitoire ou pendant la phase stationnaire Pour la phase transitoire, l'aimantation de départ peut être l'état d'équilibre thermique ou un état préparé à l'aide de motifs de préparation et l'information utilisée pour l'analyse est la variation du signal en fonction des répétitions. L'acquisition stationnaire est, au contraire, basée sur l'acquisition d'un équilibre dynamique atteint à la suite de plusieurs répétitions. 15 Rappelons ici que l'aimantation est décomposable en une composante alignée sur le champ magnétique principal BO, appelé composante longitudinale, et d'une composante transverse, cette dernière pouvant être détectée par induction dans un détecteur radiofréquence et pouvant être modulée spatialement avec les gradients. Rappelons que l'effet d'une impulsion radiofréquence appliquée à résonance et d'amplitude Bl(t), de phase phi (par rapport à une référence interne des synthétiseurs RF) et 20 appliquée pendant une durée T est de basculer l'aimantation d'un angle a avec BO, l'angle étant proportionnel à l'intégrale du champ RF pendant la durée T, l'axe de rotation est lui déphasé de 0 = phi + zr / 2 radians par rapport au référentiel tournant à la résonance. Rappelons que les composantes longitudinale et transverse relaxent avec un temps caractéristique T1 et T2, respectivement, temps dépendant de l'environnement local et éventuellement influencés par le transfert d'aimantation. 25 Rappelons enfin que les atomes ou molécules porteurs de l'aimantation nucléaire diffusent au cours du temps, ce qui est macroscopiquement décrit par les équations de Bloch-Torrey. Une complexité importante dans les séquences rapides (avec en général TR k+1, Figure 2c), le complexe conjugué est une symétrie centrale k-*-k (Figure 2d, les coefficients sont réels). La convolution avec S, 'penche' les coefficients tel que j-j-k.v (Figure 2e). 5 Pour calculer les signaux transitoires et stationnaires pour Ny x Nz états de configuration et des séries d'amplitude et de phase arbitraires, un calcul par récurrence peut être réalisé à partir de l'équation Eq.9 sur des grilles 2D pour Pns et Qns avec j e[-Ny12+1, Ny/21 et k e[-I\Ç /2 +1, N /2] initialisés à l'équilibre thermique (Pos=1, Q0s=0). Pour obtenir l'ordre k de l'aimantation effectivement mesuré, une transformée de Fourier de Qns peut être réalisée selon la direction Y-1 ce qui fournit directement les 10 valeurs pour tous les incréments A. Notons que la discrétisation selon Fl (le fait de ne prendre que Ny pas) peut être traitée en implémentant le repliement dans l'opérateur de décalage S'. En revanche, la limitation du nombre de pas dans la direction Zi nécessite d'estimer quand les coefficients des différents polynômes deviennent négligeables. Pour estimer cette valeur Nz de pas nécessaire, on peut considérer la relaxation transversale et la diffusion car leur action tend à atténuer la composante transversale. Lors du transfert 15 d'état avec le terme Z-1(W2* as) , en considérant k>0 avec un transfert vers l'état k+1, l'atténuation des coefficients est de l'ordre de E2E/ , si bien que l'atténuation cumulée depuis l'ordre zéro jusqu'à k est de l'ordre de E2kEdk3/3. On peut alors définir une limite sur l'atténuation cumulée E,;.,2 des coefficients, ce qui permet d'estimer l'ordre k qui produirait cette atténuation via la relaxation : k2 = f (E 2) , oum via la diffusion : k- .131n(E,',,, yin (Ed ) (arrondis à l'entier le plus proche). Le minimum 20 N = 2x min (k2,kdiff) fournit alors une estimation du nombre d'ordre à considérer dans le calcul. Si ka'ff > k2 , cela implique également un effet limité de la diffusion par rapport à la relaxation transversale. Bien qu'applicable à des séries d'impulsions arbitraires, le cas de séries périodiques permet d'atteindre un état stationnaire Le calcul de cet état nécessite un nombre d'itérations suffisantes neq. Si l'on considère Ti>T2, la relaxation longitudinale apparaît comme le facteur limitant pour atteindre l'équilibre. De 25 manière similaire, définir une limite pour l'atténuation cumulée à Eh no fournit un ordre de grandeur du nombre de répétitions nécessaires pour atteindre l'état stationnaire ri n n ( - eq hm1 (E1) Numériquement, un critère tel que la variation relative peut être utilisé pour arrêter le calcul itératif. On notera que le cas particulier du signal SSFP atteint avec cyclage quadratique de la phase et angle constant (an=a, u, 1) peut être calculé plus efficacement. Cet aspect a été proposé par ailleurs avec de la relaxation 30 (33). Pour prendre en compte des polynômes de pondération W1 et W2 quelconque, notamment incluant la 3036189 29 diffusion, une dérivation similaire basée sur la formulation d'un système linéaire tridiagonal est donnée dans l'annexe 3 pour calculer efficacement l'état stationnaire avec cyclage quadratique de phase et angle constant. 5 Effets hors-résonance et de flux Le procédé de calcul proposé est de préférence mis en oeuvre avec une fréquence de résonance qui est exactement celle de l'impulsion RF d'excitation. En pratique, le champ statique Bo n'est pas parfaitement homogène ce qui produit des effets sur le signal mesuré. Les effets hors-résonance tout comme les déphasages induits par le mouvement produisent une modulation de phase supplémentaire 10 q = yAB,TR +7r leen, , où AB0 représente les inhomogénéités de champ ou éventuellement le décalage chimique, et ou Veo, représente le codage de vitesse lié au moment d'ordre 1 de la forme du gradient. En introduisant l'exponentielle complexe z, = exp (-400), l'action de ce déphasage supplémentaire selon Z-1 dans la représentation 2D en état de configuration (9) est exprimée de manière similaire en considérant un polynôme déphasé en Y et (ZoZ)-1 pour lequel tous les coefficients restent réels comme dans le cas 15 précédemment considéré. Entre la fin de l'impulsion RF et le temps d'écho, le signal est atténué par les effets de la relaxation transversale apparente E; = exp (-TE/T2) et déphasé d'un facteur ZTE = exp ((pTE) avec (pTE = (pi + yABoTE + (pv qui inclus un éventuel déphasage (pi entre l'émission et la réception, les effets de l'inhomogénéité du champ Bo et de codage de vitesse éventuel (p v . 20 Variante du premier mode de réalisation d'un procédé et d'un dispositif selon l'invention : imagerie multiparamétrique ; Implémentation d'un problème inverse basé sur l'acquisition de l'état k=0 Dans une variante (dite « d'imagerie multiparamétrique ») du premier mode de réalisation de procédé selon l'invention mis en oeuvre par le premier mode de réalisation de dispositif selon l'invention (et décrite uniquement par rapport à ses différences ou particularités par rapport à ce premier mode de 25 réalisation), le procédé selon l'invention comprend : - une acquisition d'un signal (de préférence pour l'ordre de cohérence égal à zéro) pour des états d'équilibres correspondants à des amplitudes différentes de l'impulsion radiofréquence et/ou des A différents et/ou différentes valeurs de la suite un_,1 = --tin périodique constante ou non, les moyens de commande 5 étant agencés pour ou programmés pour commander les moyens d'acquisition 6 pour 30 effectuer une telle acquisition, et - pour au moins un point dans l'échantillon, une détermination (par les moyens de calcul 7, qui sont agencés pour ou programmés pour cela) de plusieurs données (de préférence au moins 3, de préférence 3036189 30 toutes) parmi une aimantation nucléaire, un angle de bascule de l'aimantation, un coefficient de diffusion, un taux ou temps de relaxation longitudinal R1 ou Tl, et un taux ou temps de relaxation transversal R2 ou T2, à partir du signal acquis pour ces états d'équilibres. Comme nous le verrons plus en détails par la suite, la détermination est de préférence effectuée, 5 (par les moyens de calcul 7, qui sont agencés pour ou programmés pour cela) soit par comparaison avec un dictionnaire pré-calculé soit par estimation itérative, selon une minimisation d'une norme de la différence entre le signal acquis exprimé sous forme complexe avec une partie réelle et une partie imaginaire et un modèle du signal exprimé sous forme complexe avec une partie réelle et une partie imaginaire. Dans l'exemple, la minimisation comprend un algorithme de minimisation au sens des 10 moindres carrés, de préférence en utilisant l'algorithme de Gauss-Newton appliqué aux problèmes non- linéaires. On peut se placer, pour cette variante d'imagerie multiparamétrique, dans le cas dans lequel on exploite des états d'équilibre séquentiels tels que précédemment décrit avec plusieurs ensembles d'applications de séquence. Dans ce cas, cette variante d'imagerie multiparamétrique comprend une 15 acquisition d'un signal (de préférence pour l'ordre de cohérence égal à zéro) pour des états d'équilibres (de préférence avec une suite un+1 = vn+1 -vn constante pour chaque ensemble) correspondants à des différents (les moyens de commande 5 étant agencés pour ou programmés pour commander les moyens d'acquisition 6 pour effectuer une telle acquisition), dont : - plusieurs points pour une valeur absolue de 4 proches de 0°, entre 4=0° inclus et 4=32° 20 inclus modulo 360°, et/ou - plusieurs points pour une valeur absolue de 4 proches de 180°, entre 4=164° inclus et 4= 196° inclus modulo 360°, et/ou - plusieurs points pour une valeur absolue de 4 proches de 120°, entre 4=112° inclus et 4=128° inclus modulo 360°, et/ou 25 - plusieurs points, pour une valeur absolue de 4 parmi 45°, 60°, 72°, 90°, et 144° modulo 360° pour un total de préférence d'au moins cinq points. On voit sur les figures 6 et 7 l'intérêt de ces points, en particulier proches de 0°, 120° et 180°. Dans cette variante d'imagerie multiparamétrique, l'amplitude de l'impulsion radiofréquence 30 correspond toujours à un angle de bascule de l'aimantation : - supérieur à l'angle de Ernst pour un temps de relaxation longitudinal T1 égal à 2000 millisecondes (ce qui laisse la liberté sur le TR et sur le type de tissu ciblé), et - inférieur ou égal à 90° les moyens de commande 5 étant agencés pour ou programmés pour cela. 3036189 31 Dans cette variante d'imagerie multiparamétrique, la détermination peut comprendre une condition de continuité spatiale de chacune des données déterminées parmi l'aimantation nucléaire, de l'angle de bascule de l'aimantation, du coefficient de diffusion, du taux ou temps de relaxation longitudinal R1 ou T 1, et du taux ou temps de relaxation transversal R2 ou T2 entre différents points de 5 l'échantillon (les moyens de commande 5 étant agencés pour ou programmés pour cela), de préférence l'angle de bascule. Nous allons à présent expliquer plus en détail cette variante d'imagerie multiparamétrique par un premier exemple détaillé. De multiples procédés d'imagerie paramétrique ont été proposées basées sur l'acquisition de plusieurs états k (28-31), mais beaucoup moins basées sur différents cyclage de phase 10 (33-35). Pour illustrer le calcul direct selon la variante « d'imagerie multiparamétrique » selon l'invention et le valider en présence de relaxation et de diffusion, un problème inverse est proposé basé sur l'acquisition de l'état k=0. La description faite ici considère le cas du signal à l'état stationnaire pour un angle constant en amplitude a et avec des incréments quadratiques de phases différents, mais l'expression est générale et peut être étendue aux autres cas d'équilibres stationnaires, ou même au cas de 15 phases transitoires (hors équilibre dynamique). Pour chaque voxel de l'espace image, les signaux mesurés sont placés dans un vecteur noté Q0,','. Le problème inverse est exprimé comme une minimisation au sens des moindres carrés : H2 min 11Q0,' -M x Q0 (El, E2, Ed, a .a Ma = Dans cette expression, un certain nombre de paramètres sont inconnus : M qui représente 20 l'amplitude du signal au temps d'écho, les atténuations par relaxation et par diffusion, ainsi que l'amplitude de l'angle effectivement réalisé. De multiples procédés peuvent être utilisés pour résoudre ce problème des moindres carrés non-linéaires. Une implémentation basée sur la combinaison d'une initialisation par comparaison à un dictionnaire, suivie par l'application de l'algorithme de Gauss-Newton est proposée ici. Tout d'abord, on peut noter que le problème est linéaire en M, de sorte que si le vecteur 25 modèle Q0 (E1, E2, E d , a) est connu, alors la meilleure solution pour M au sens des moindres carrés est : M = Q01 I QI° (16) La notation en exposant H représente la transposée Hermitienne du vecteur Q0. Pour initialiser les paramètres à (Mrru, Eure, E2,ioit, amit), on peut calculer Ed,init en considérant le coefficient de diffusion libre de l'eau et choisir am, comme étant l'angle prescrit. Cela réduit le calcul de 30 Qo à deux dimensions (selon E1 et E2) et différentes valeurs pour E1 et E2 peuvent être choisies dans la gamme de 0 à 1. L'estimation linéaire de M pour chaque couple de valeurs (El, E2) en utilisant l'équation Eq.16 permet de choisir des valeurs (M;rut, Eure, E2,'111) qui minimise la norme résiduelle donnée par (15) 3036189 32 l'équation Eq.15. Par la suite, dans le calcul non-linéaire itératif, le signal est normalisé par l'estimation de M et, pour s'assurer que des valeurs physiques d'atténuation sont obtenues, on effectue un changement de variables avec 4= exp(-11x12), E2/4 = exp(-11x2n2), Ed=exp(-11xd2). La matrice Jacobienne est alors donnée par : 5 J = real (a ) real (Q ° real imag initial ( real ( real fflo (17) et peut être calculée jour est ensuite obtenu imag (Q, ) imag (ffl ° imag Q'`) imag de mise-à- numériquement au par : aie a \,a2/1 '/ \, d '/ fflo' \ 1 ) ( \ 1 \ \ a Q'" \ 1 ) \,a2/1 ) \ d ) . Un vecteur vecteur [1, xunit, [dr,a5c1, cbc211,c1xcl,daf =(-11141 JH real(Q0) - real(Q,m' /M) imag(Q0) imag (Q0,m' /M) (18) Les fonctions real et imag représentent respectivement les parties réelles et imaginaires. 10 L'addition de ce vecteur mise-à-jour à l'estimation précédente affine la solution. M peut ensuite être ré- estimé comme M x (1+ dr) et utilisé pour normaliser le signal mesuré pour la prochaine itération. La procédure est répétée jusqu'à ce que la variation relative de la norme résiduelle donnée par l'équation Eq.15 soit suffisamment petite. La matrice JHJ est supposée inversible. Si l'on considère la nature Gaussienne du bruit en IRM, caractérisé par sa variance 62 identique sur les parties réelles et 15 imaginaires, la matrice d'information de Fisher (36) est égale à F =( JHJ. Si elle est inversible, son inverse F-1 = (01 M1)2 -1 représente la matrice de covariance du bruit à partir de laquelle les bornes inférieures de Cramer-Rao peuvent être extraites pour estimer la précision sur les paramètres ajustés si un estimateur non biaisé est obtenu à l'aide de l'algorithme proposé. Si la matrice de Fisher n'est pas inversible ou est mal conditionnée, la pseudo-inverse est calculée pour assurer la stabilité 20 numérique. De plus, des contraintes peuvent être ajoutées sur les paramètres pour éviter d'obtenir des valeurs aberrantes et réduire les instabilités numériques : M est positif, l'angle de bascule reste dans une gamme déterminée (de 2° à 90° par exemple), les termes d'atténuation sont bornés entre exp (-1) et exp (-101 . Cette procédure d'ajustement générique peut-être adaptée aussi bien à l'ajustement du signal transitoire, qu'inclure des acquisitions obtenues de manière entrelacée ou séquentielle en faisant varier les 25 séries d'impulsions RF. 3036189 33 Si l'exemple de cette variante d'imagerie multiparamétrique venant d'être détaillé se place dans le contexte d'équilibres séquentiels avec plusieurs ensembles n'ayant chacun qu'un seul équilibre dynamique, cette variante d'imagerie multiparamétrique peut aussi être mis en oeuvre selon l'invention avec plusieurs ensembles en exploitant des états d'équilibres entrelacés (de préférence au moins 2 états 5 d'équilibre entrelacés), ainsi : dans un second exemple, on utilisera plutôt une suite un+1 = vn+1 - vn constante, 4 bien évidemment constant pour chaque ensemble, et l'amplitude de l'impulsion radiofréquence qui suit une suite périodique non constante (avec donc une période au moins égale à deux répétitions avec deux amplitudes différentes suffisamment différentes entre 0° et 90°, de 10 manière à obtenir les au moins deux états d'équilibres entrelacés différents), comme montré sur les figures 6c et 6d. Pour chaque point d'intérêt de l'échantillon, les données (aimantation nucléaire, angle de bascule de l'aimantation, coefficient de diffusion, taux ou temps de relaxation longitudinal R1 ou Tl, et taux ou temps de relaxation transversal R2 ou T2) sont alors déterminées avec l'ajustement proposé, comme montré sur les figures 7c et 7d. On 15 notera que l'utilisation d'amplitudes différentes n'ajoute pas d'inconnue au problème inverse car, s'il peut y avoir une incertitude sur l'angle de basculement réalisé avec une amplitude donnée, les angles sont proportionnels à la série d'amplitude. On adaptera donc par exemple l'ajustement en considérant l'angle le plus grand de la série d'amplitude comme une inconnue à déterminer, les autres angles étant directement proportionnels à cet angle d'un 20 facteur connu déduit de la suite d'amplitude imposée. Dans un troisième exemple, on utilisera plutôt une amplitude de l'angle constante et une série un+1 = v.+1 périodique non constante avec 4 constant pour chaque ensemble, mais prioritairement proche de 0° et/ou proche de 180°, comme montré sur la figure 6b. Pour chaque point d'intérêt de l'échantillon, les données (aimantation nucléaire, angle de bascule 25 de l'aimantation, coefficient de diffusion, taux ou temps de relaxation longitudinal R1 ou Ti, et taux ou temps de relaxation transversal R2 ou T2) sont alors déterminées avec l'ajustement proposé. Evidemment, on peut imaginer cette variante « d'imagerie multiparamétrique » en mélangeant l'acquisition séquentielle d'états d'équilibres entrelacés générés par la périodicité de l'amplitude et de la 30 suite un+1 = vn+1 -vn , par exemple en mélangeant le second exemple et le troisième exemple. Le cas extrême de l'acquisition d'un seul ensemble d'équilibres entrelacés avec des suites périodiques non constante d'amplitude et/ou un_,1 = vn+1 --tin non constante est également possible. Pour cela, on choisira de préférence : 3036189 34 - Dans le cas où l'amplitude est constante, un+1 Vn+1 lin non constante aura une périodicité de préférence d'au moins cinq points et A sera constant d'une valeur de préférence de 180°. La suite u'_,1= v ' pourra être choisie arbitrairement. - Dans le cas où la suite d'amplitude sera périodique non constante de préférence d'une 5 période d'au moins cinq répétitions, la suite u'_,1= v ',1-v ' pourra être choisi soit constante soit non-constante avec la même période. On choisira A, qui est alors unique et constant pour ce seul ensemble, de préférence d'une valeur de 0° ou de 180° ou proche de ces valeurs. La suite d'amplitude pourra être choisie arbitrairement entre -90° et 90° avec des valeurs suffisamment différentes entre elles, par exemple réparties uniformément entre 0° et 90°. La 10 suite un+1=v,,,1-- T, pourra être choisie arbitrairement, de préférence avec une moyenne nulle sur la période. Pour chaque point d'intérêt de l'échantillon, les données (aimantation nucléaire, angle de bascule de l'aimantation, coefficient de diffusion, taux ou temps de relaxation longitudinal R1 ou T 1, et taux ou temps de relaxation transversal R2 ou T2) sont alors déterminées avec l'ajustement proposé. De manière 15 similaire, on adaptera l'ajustement dans le cas d'amplitude périodique non-constante par exemple en considérant l'angle le plus grand de la série d'amplitude comme une inconnue à déterminer, les autres angles étant directement proportionnels à cet angle d'un facteur connu déduit de la suite d'amplitude imposée. 20 Expériences sur fantôme Pour illustrer l'extension en 2D du formalisme des états de configuration, des expériences ont été réalisées avec différentes solutions d'agents de contraste dilués dans l'eau. L'imagerie a été faite à 1.5 T (Achieva; Philips Healthcare, Best, The Netherlands) en utilisant une antenne tête en quadrature pour la réception. Des chélates de Gadolinium (10, 5, 2.5, 1.25, 0 mM, Dotarem; Guerbet, Villepinte, France) et 25 des nanoparticules d'oxydes de fer (0.32, 0.16, 0.08, 0.04, 0 mM, CL-30Q02-2, Molday ION, Biopal, Worcester, MA, USA) ont été préparés et placés dans des contenants cylindriques de 15 mL alignés au champ magnétique Bo. Une séquence en écho de gradient 3D non-sélective spatialement a été modifiée pour permettre de multiples valeurs de Fi séquentiellement en commençant par ri = 1 (pas de brouillage RF, 4=0°) jusqu'à couvrir tout le cercle unité. Des impulsions RF rectangulaires ont été utilisées (durée 30 de 150 jts pour 45° d'angle prescrit). Un paramètre a également été ajouté pour permettre l'entrelacement de l'amplitude des angles tel que cc2 = a et 6c2.+1 = a/2 (en divisant par deux en temps réel l'amplitude du champ B1 d'émission). Le gradient de brouillage a été placé sur la direction de la lecture tel que son 3036189 aire corresponde à a = 1/Akz égal à la taille du pixel d'imagerie. Les paramètres d'acquisition étaient soit TRITE=9.214 ms, matrice d'acquisition 168x84x9 et taille du voxel 0.5x1x8 mm, bande passante par pixel 217 Hz/pix, temps d'acquisition Tacq = 7 s par volume (que l'on notera par la suite 'séquence l', avec une sensibilité à la diffusion de b=1.45 s mm-2) ou TR/TE=4.8/2.4 ms, matrice d'acquisition 84x84x9 et taille 5 de voxel 1x1x8 mm, bande passante par pixel 434 Hz/pix, temps d'acquisition Tacq = 3.5 s par volume (que l'on notera par la suite 'séquence 2', avec une sensibilité à la diffusion de b=0.19 s mm 2) La taille des voxels de reconstruction était de 0.47x0.47x8 mm Plusieurs angles prescrits ont été testés séquentiellement (entre 7.5° to 75°), ainsi que de manière entrelacée. De multiples pas Ny ont été effectués (par pas de 2° pour A). Etant donné ces paramètres d'acquisition, pour chaque pas A, le centre de l'espace 10 de Fourier a été échantillonné après quelques secondes. Il n'y avait pas de délai entre les acquisitions des différents pas A, et 2 à 4 volumes de préparation ont été réalisés avant de stocker le premier volume de données si bien qu'on peut considérer que l'état stationnaire était atteint. La reconstruction de chaque volume a été effectuée par le logiciel constructeur en retirant les corrections de phase présentes par défaut et en utilisant une échelle de stockage absolue du signal RMN dans les fichiers DICOM. Les images 15 d'amplitude et de phase ont été sauvegardées. Analyse des données Les images DICOM ont été analysées hors ligne en utilisant Matlab (2011b version; The Mathworks, Natick, MA) sur un ordinateur portable récent (Intel core i7, RAM 16 Go, 2.7 GHz) exploité sous 20 Microsoft Windows 8. Seule la coupe centrale du volume 3D a été analysée. Une région d'intérêt a été dessinée sur chaque tube pour réaliser l'analyse du signal moyen par tube. Une reconstruction pixel à pixel a également été réalisée. Les images ont été normalisées par l'écart type du bruit a qui a été estimé sur une région sans signal. Pour garantir que le signal ait une symétrie Hermitienne par rapport à Y-1 ( la phase du premier pas A=0° a été utilisée pour estimer Z1F qui fut ensuite Qo,mes (y -1) = Q0,mes (Y ).' j 25 retiré des autres pas. Un léger décalage temporel a été remarqué lors de l'acquisition des différents pas, et corrigé. Ce décalage a été attribué à une dérive thermique des gradients et il a été considéré que cette dérive variait linéairement dans le temps. Pour estimer cette dérive, l'acquisition correspondant à A=0° a été répétée au début et à la fin de tous les autres pas, ou alors l'acquisition obtenue pour A=2° (au début) et pour A=358° (à la fin) ont été comparées. Ces différents pas sont supposés être complexes conjuguées 30 si bien que leurs produits doivent être réels. La phase de ce produit est donc caractéristique à de la dérive de phase temporelle. Une fois la dérive caractérisée, le signal mesuré pour chaque pas a été rephasé pour fournir le signal Q0. sur lequel le problème inverse proposé a été appliqué. Lors de l'estimation itérative des paramètres dans le problème inverse, l'algorithme de calcul de l'état stationnaire proposée 3036189 36 dans l'annexe 3 a été utilisé pour les cas où un angle constant en amplitude et en décalage apparent de fréquence a été appliqué. Le nombre d'ordre Nz a été choisi en fixant Ehm,2=0.01. Dans tous les autres cas, le calcul itératif a été effectué en fixant le nombre maximum de répétitions à max(neq, no) où no correspondait au nombre de répétitions appliquées pour atteindre le centre du plan de Fourier en imagerie, 5 et nec, obtenu pour Eh',,/=0.01. L'initialisation a été effectuée sur 32 valeurs pour E1 et E2/Ei. Les itérations ont été arrêtées quand la norme du résidu normalisée par la norme du signal variait de moins de 10-6. Résultats Simulations 10 Le calcul de l'aimantation en utilisant l'extension en 2D du formalisme des états de configuration est illustré sur la Figure 3. Les coefficients de l'aimantation transversale q,,k,n, qui ont été initialisés à l'équilibre thermique, sont affichés pour n=1 (avec seulement un coefficient central à j=k=0, Figure 3-a). Progressivement, des coefficients non-nuls apparaissent, jusqu'à atteindre un état stationnaire après un large nombre d'impulsions. Les coefficients de l'état stationnaire sont très fortement influencés par la 15 diffusion (Figure 3-e et f). Chaque itération a duré -1 ms avec Ny= 360, et le nombre d'itérations nécessaires pour atteindre l'état stationnaire était de l'ordre de 5 Ti/TR, quelle que soit la série d'angle appliquée. En ce qui concerne le calcul de l'état stationnaire à la suite du cyclage quadratique de phase, les deux approches (itérative basée sur l'application successive d'impulsion et directe dérivée de l'algorithme tridiagonal) donnaient les mêmes états stationnaires L'algorithme rapide prenait -10 ms 20 pour fournir l'état stationnaire sur le même nombre de points Ny = 360. L'algorithme proposé fournissait les mêmes résultats que publiés précédemment lorsque les effets de la diffusion étaient négligés (Figure 4-a). L'état stationnaire est très fortement dépendant de A avec la création de multiples échos stimulés. Des pics spécifiques sont visibles lorsque A est un nombre rationnel de 360° (0°, 180°, 120°, 90°, ...). La relaxation transverse et la diffusion réduisent toutes deux la formation d'échos stimulés ce qui rapproche 25 le signal d'un équilibre de Ernst pur (Figure 4-b-d). L'aimantation transversale varie entre des valeurs plus grandes et plus petites, et n'atteint jamais de l'équilibre de Ernst qui est purement réel, même s'il existe des valeurs de A pour lesquelles le module de l'aimantation lui est égal. La dépendance de l'état stationnaire obtenu à la suite d'une série constante en amplitude et en décalage apparent de fréquence est complexe et des contrastes mélangés en T1, T2 et diffusion sont obtenus. Néanmoins, des caractéristiques 30 générales peuvent être extraites. En effet, les états stationnaires sont contenus dans un plan qui contient lui-même l'équilibre de Ernst et qui est basculé d'un angle effectif oce = arctan ((E1 -c)/s)(Figure 5) par rapport à l'axe du champ principal, et ce quel que soit T2 ou le coefficient de diffusion. Ceci implique une dépendance particulière par rapport à l'angle de Ernst. Pour des angles d'amplitude plus petites que 3036189 37 l'angle de Ernst (Figure 5-a-b), l'angle effectif est négatif et l'aimantation varie entre des 'pics' qui sont plus petits que l'équilibre de Ernst et des valeurs 'brouillées' qui sont plus grandes. Pour des angles plus grands que l'angle de Ernst (Figure 5-e-f), l'inverse est observé. Entre les deux, si l'angle est exactement égal à l'angle de Ernst (Figure 5-c-d), l'angle effectif est nul : seule la partie imaginaire de la composante 5 transversale de l'aimantation varie avec A et la partie réelle est toujours égale à l'équilibre de Ernst. L'algorithme proposé est non seulement capable de simuler le cas du cyclage quadratique de phase standard (amplitude et décalage apparent de fréquence constants), mais aussi de calculer les états stationnaires entrelacés qui peuvent être obtenus avec des séries périodiques d'amplitude et de décalage apparent de fréquence. Ceci est démontré sur la Figure 6 sur laquelle un angle constant de 75° avec 10 décalage apparent de fréquences constant est comparé à une série ayant un angle constant de 75° avec une oscillation du décalage apparent de fréquence (+A, -A), ainsi qu'une série qui oscille en amplitude (75°/37.5°) et est constante en décalage apparent de fréquence. Les états stationnaires obtenus sont très différents. Le cas de l'amplitude alternée donne naissance à une dépendance très complexe avec A. Le cas d'un angle d'amplitude constante mais avec une alternance du décalage apparent de fréquence (u2 1, 15 /12.+1--1) donne naissance à deux états stationnaires complexes conjugués. Pour ce dernier cas, une transition lente est obtenue entre A=0 et A=n (ces deux états sont respectivement les mêmes que ceux obtenus pour un décalage apparent de fréquence constant avec A=0 et A=n). Cette alternance de décalage apparent de fréquence est obtenue par exemple avec une phase effective des impulsions RF O qui alterne entre +A/2 et -412. 20 Expériences sur fantôme Les signaux moyens sur les tubes ont pu être modélisés efficacement en utilisant les approches proposées sur la gamme des angles de basculement testés, aussi bien en utilisant une série constante d'amplitude et de décalage apparent de fréquence qu'avec des séries périodiques plus complexes. Après l'étape 25 d'initialisation, le nombre moyen d'itérations de l'algorithme non linéaire de minimisation au sens des moindres carrés était de -35. Pour toutes les séries d'angle testées, les valeurs de relaxivité R1 et R2 extraites de l'analyse sur le signal moyen des tubes étaient linéaires avec la concentration pour les solutions de Gd et d'ION, les angles ajustés ainsi que la valeur de l'atténuation Ed étaient consistants avec respectivement l'angle prescrit et le coefficient de diffusion libre de l'eau. Pour le cas d'amplitude 30 constant à 45° obtenu avec la séquence 1, les relaxivités molaires étaient ri/r2 = 3.9/4.6 mM-is-1 pour l'agent à base de Gd et 86.4/12.7 mM-1 s-1 pour l'agent à base d'ION, l'angle ajusté était de 44.3±0.3°, l'aimantation nucléaire Mo était 471±54 (moyenne ± écart-type sur les tubes) et le coefficient d'atténuation par diffusion était Ed = 0.9967±0.0003, correspondant à D = 2.29±0.17 10' mes', ce qui est en accord avec la théorie. Des signaux mesurés représentatifs ainsi que les paramètres ajustés 3036189 38 correspondants sont présentés sur la Figure 7 pour les tubes de 1.25 mM Gd et de 160 jtM ION dans les cas 45° et 45/22.5° entrelacés. Comme attendu, le cas de l'amplitude de l'angle entrelacé produit un équilibre complexe. Pour démontrer la possibilité de réaliser de l'imagerie multiparamétrique à partir de l'approche proposée, les signaux obtenus avec la séquence 1 avec 45° ont été ajustés pixel à pixel. Les 5 cartes paramétriques d'aimantation, d'angle de basculement, de R1, R2 et de coefficient de diffusion sont présentées sur la Figure 8. Les valeurs moyennes ainsi que les écart-types calculés sur ces images reconstruites sont donnés dans la Table 1 ci-dessous, et associés à la valeur théorique calculée à partir de l'élément diagonal de la matrice de covariance du bruit. Tous les paramètres à ajuster sont en très bon accord avec les valeurs attendues. R1 et R2 sont linéaires avec la concentration (Figure 9) et permettent 10 d'estimer des relaxivités molaires en accord avec celles mesurées sur les signaux moyens des tubes. L'amplitude de l'angle de bascule est très proche de la valeur prescrite, et le coefficient de diffusion indique une diffusion libre dans tous les tubes. Comme on peut le voir, l'estimation des paramètres est influencée par la propagation du bruit, et la précision dépend d'une manière complexe des paramètres. Les éléments diagonaux de la matrice de covariance du bruit, alors qu'ils n'incluent qu'une partie de la 15 propagation du bruit, fournissent déjà une bonne estimation de la précision sur R1, R2 et D, mais sont moins efficaces pour prédire la propagation du bruit sur M et sur l'amplitude de l'angle, notamment parce que les estimations de ces deux paramètres sont plus fortement corrélées. 20 25 30 3036189 39 Eau Gd Gd Gd Gd ION ION ION ION 1.25 mM 2.5 mM 5 mM 10 mM 40 jtM 80 jtM 160 jtM 320 jtM M/6 797 588 570 548 511 620 623 635 604 ±217 ±114 ±113 ±112 ±117 ±107 ±89 ±89 ±84 (±115) (±6) (±4) (±3) (±4) (±36) (±33) (±29) (±22) a (°) 45.53 44.06 44.01 44.13 44.17 43.90 44.28 44.21 44.23 ±3.54 ±0.26 ±0.21 ±0.23 ±0.51 ±1.87 ±0.85 ±0.52 ±0.41 (±0.01) (±0.05) (±0.10) (±0.13) (±0.29) (±0.11) (±0.05) (±0.11) (0.24) Ri (s-1) 0.28 5.54 10.05 20.14 39.57 0.78 1.26 2.19 4.24 ±0.07 ±0.13 ±0.20 ±0.42 ±1.98 ±0.12 ±0.15 ±0.18 ±0.29 (±0.04) (±0.06) (±0.10) (±0.19) (±0.57) (±0.05) (±0.07) (±0.11) (±0.18) R2 (s-1) 0.73 6.47 11.69 23.82 46.78 3.10 6.43 13.47 27.90 ±0.27 ±0.20 ±0.28 ±0.55 ±1.79 ±0.39 ±0.49 ±0.56 ±0.88 (±0.17) (±0.10) (±0.14) (±0.25) (±0.67) (±0.16) (±0.22) (±0.33) (±0.55) D 1.96 2.23 2.27 2.36 2.72 2.27 2.21 2.19 2.15 10- ±0.28 ±0.14 ±0.18 ±0.52 ±2.53 ±0.27 ±0.27 ±0.40 ±0.78 91112S-1 (±0.14) (±0.08) (±0.12) (±0.32) (±1.66) (±0.17) (±0.20) (±0.29) (±0.63) La Table 1 ci-dessus présente les résultats obtenus après analyse des images multiparamétriques reconstruites à partir des données obtenues avec la séquence 1 avec un angle prescrit de 45° et Ny=91 (A entre 0° to 180° par pas de 2°). Valeur moyenne ± écart-type mesurés sur les tubes pour les différents paramètres ajustés. Le nombre entre parenthèse correspond à la borne de Cramer-Rao issue des éléments 5 diagonaux de l'inverse de la matrice de Fisher estimée aux valeurs moyennes ajustées. Bien que toutes les sources de bruit ne soient pas présentes dans les éléments diagonaux, ce nombre théorique fournit déjà une bonne estimation de l'écart-type, avec des valeurs proches de l'écart-type mesuré sur les tubes, pour Ri, R2 et D. 10 Discussion L'application répétée d'impulsions RF entrelacés d'une aire constante de gradient a été analysée à l'aide du formalisme des états de configuration étendu pour prendre en compte la modulation de la phase des 3036189 impulsions RF. Des algorithmes efficaces de calcul de l'aimantation incluant la relaxation et la diffusion ont été déduits. L'imagerie multiparamétrique basée sur l'ajustement de signaux acquit à l'ordre Ic=0 avec plusieurs décalages de fréquence apparente ainsi que des angles alternés en amplitude a été démontré pour une large gamme de R1, R2 et rapport R2/Ri. L'aimantation nucléaire, l'angle de basculement, le 5 coefficient de diffusion libre et les paramètres de relaxation ont pu être mesurés, validant les approches directe et inverse. Le formalisme des états de configuration est équivalent à l'approche du graphe de phase étendue (6-9), dans laquelle l'aimantation est décomposée en coefficients de Fourier évoluant dans le temps, qui peut 10 donc être étendu avec une représentation graphique similaire (8) de l'action d'un angle de basculement d'amplitude et de phase données, et d'atténuation par relaxation et diffusion. L'approche est également similaire à la décomposition en polynôme effectuée dans le cadre de la conception d'impulsions RF (10,11) étendue en 2D, suggérant la possibilité d'utiliser des techniques de filtrage numérique pour la conception d'impulsions RF potentiellement sensibles à la relaxation et à la diffusion. On notera que la 15 dérivation proposée peut être adaptée à des TR et des aires de gradient variables (9) en adaptant les polynômes de pondération et les décalages selon Fi et Zi pour calculer efficacement la réponse à de telles séquences. Dans le cas d'un TR et d'une aire de gradient constants comme traité ici, le cyclage quadratique de phase mène bien à l'état stationnaire très étudié (3-5). L'ordre k= 0 à l'état stationnaire dans ce cas est contenus dans un plan qui passe par l'équilibre de Ernst et qui est penché d'un angle 20 effectif déterminé menant à une dépendance particulière qui pourrait être utilisée pour rapidement estimer si l'angle prescrit est plus grand ou plus petit que l'angle de Ernst. Ceci pourrait être fait en comparant les signaux de quelques valeurs de A (comme par exemple 0° et un autre incrément de phase). Une décroissance de la partie réelle indique un angle plus grand que l'angle de Ernst, et réciproquement. La dérivation d'un algorithme rapide basée sur l'expression d'un système linéaire tri-diagonal pour calculer 25 l'état stationnaire dans ce cas, proposée précédemment avec une dérivation différente (33), inclut également ici les effets de la diffusion et peut être généralisée à des polynômes de pondérations plus complexes. De plus, la description selon les états de configuration permet de démontrer que des séries périodiques, modulées en amplitude et en phase, produisent des états stationnaires entrelacés, étendant donc les possibilités de manipulation du contraste dans ce type de séquences. Comme les paramètres 30 physiques (taux de relaxation, coefficient de diffusion) et de séquences (TR, amplitude de l'angle) sont pris en compte dans le calcul du problème direct, le problème inverse est basé sur l'acquisition de différentes séries d'angles. Ce type de problème inverse a déjà été proposé en utilisant différents états k (28-31) ou en utilisant du brouillage partiel autour de 0° avec cyclage quadratique de la phase pour quantifier T2 (33-35). Il a été choisi ici de n'utiliser que l'état k= 0 et de multiples décalages apparent de 3036189 41 fréquence. Cela a permis de valider expérimentalement, sur un système clinique standard, le problème direct, ainsi que de démontrer la possibilité de réaliser de l'imagerie multiparamétrique. Le problème inverse peut être étendu pour inclure différents états k et de multiples séries d'angle. La détermination des états k et des séries optimaux à mesurer est un problème complexe qui nécessiterait l'étude approfondie 5 de la matrice de covariance du bruit. Pour la question de la sélection d'un nombre minimum de mesures pour estimer de manière optimale un paramètre d'intérêt (relaxation, diffusion, ...) tout en maintenant une durée totale d'acquisition réalisable en pratique, des éléments peuvent être précisés. Dans un premier temps, l'aimantation à une symétrie Hermitienne de sorte qu'il n'y a pas besoin d'acquérir des valeurs d'incrément de -A (modulo 360°) si la valeur A est acquise. D'autre part, si on se restreint à l'acquisition 10 de l'état k = 0 avec un décalage apparent de fréquence, comme proposé avec du brouillage partielle autour de 0° (33-35), les régions autour de fraction de 360° semblent être les plus informatives (0°, 180°, 120°, ...) de sorte que le brouillage partiel, définit comme un échantillonnage autour de ces valeurs, pourrait être suffisant. De plus, comme on peut voir sur la Fig. 7a-b, les deux premières ellipses correspondent au brouillage partiel autour de 0° et de 180°, et la vitesse de variation augmente autour de 180° par rapport à 15 0° (visible en comparant la position des signaux obtenus pour 0-2-4° et 180-182-184°, par exemple), ce qui suggère que 1) une évolution géométrique du type 0-2-4-8-16° autour de 0° par exemple, et 2) une réduction des pas telles que 180-181-182-184-188° autour de 180° peut être suffisant pour moduler l'état stationnaire suffisamment pour estimer les paramètres. 20 Le signal mesuré a été pré-traité pour retirer la phase au temps d'écho, ainsi qu'une dérive dans le temps. Ces paramètres ont été estimés et corrigés avant l'application du problème inverse mais auraient pu également être inclus comme paramètres à ajuster. La phase du signal, toujours disponible dans une séquence IRM, reste sensible aux inhomogeneités du champ Bo, permettant de les cartographier, ou de quantifier les déphasages dus au flux. Pour réduire les effets de phases dus au flux (9) qui ne restent pas 25 constants pendant l'acquisition et donc n'atteignent pas forcément un état stationnaire, des formes de gradients compensées en flux peuvent être utilisées. En ce qui concerne la dérive attribuée à des effets thermiques des gradients, la réduction du nombre de pas échantillonnés limiterait ces effets. Une calibration est toujours possible en acquérant plusieurs fois des incréments de 0°, ou en comparant des images obtenues avec des incréments A et -A acquises à différents moments. 30 L'implémentation proposée pour résoudre le problème inverse possède la limite du temps de calcul important lorsqu'il est appliqué sur la base d'un ajustement pixel à pixel. Néanmoins, cet algorithme est fortement parallélélisable, et des processeurs de carte graphique peuvent être utilisés pour accélérer la reconstruction. L'angle de basculement a été considéré comme un paramètre ajustable fournissant une 3036189 42 estimation moins précise des autres paramètres par propagation du bruit. Pour améliorer la précision sur les autres paramètres, on peut fixer l'angle de basculement à la valeur de l'angle prescrit ou ajouter de la régularité spatiale sur celui-ci. Les valeurs d'angle utilisées ici (entre 7.5 et 75°) semblent fournir une variation suffisante du signal mesuré pour ajuster les différents paramètres. Pour des petits angles (Fig. 5 5a-b), l'aimantation varie moins que pour des angles plus grand (Fig.5c-f). Pour des angles encore plus grand (non présenté ici), l'aimantation souffre de phénomènes de saturation ce qui indique qu'un angle maximisant les variations du signal en fonction de A se situe dans la gamme d'amplitude de 0° à 90°. Cet angle optimal dépend de la relaxation et de la diffusion. La possibilité de cartographier l'angle de basculement, la relaxation et la diffusion dans la gamme des petits angles REVENDICATIONS1. Procédé d'imagerie par résonance magnétique, comprenant : - une application continue d'un champ magnétique principal BO selon un axe Z sur un échantillon, et - au moins un ensemble d'applications répétées, sur l'échantillon et selon une période TR, d'une séquence, ladite séquence comprenant : o une impulsion radiofréquence éventuellement d'amplitude et/ou de phase variable à chaque répétition, et o après l'impulsion radiofréquence de la séquence, un gradient spatial de la composante selon l'axe Z du champ magnétique, caractérisé en ce que, au cours des applications répétées de l'impulsion radiofréquence et du gradient spatial de champ magnétique de la séquence d'un même ensemble : - l'impulsion radiofréquence suit, entre ses différentes applications répétées, une suite périodique : o pour son amplitude et o pour un+1 = y',1 -va , où n est un entier supérieur ou égal à 1 représentant le numéro de la répétition de la séquence pour cet ensemble, et où un FI est une suite de nombres entiers relatifs qui permet de définir la suite y' à une constante arbitraire près y1 tel que yin = v' x A et pn=9,-9, i avec A qui est un nombre réel constant pour toutes les applications répétées de la séquence de cet ensemble et cp, qui est l'incrément entre la phase 0, de l'impulsion radiofréquence à sa n'ème répétition dans cet ensemble et la phase 0, 1 de l'impulsion radiofréquence à sa (n-1)1èmC répétition dans cet ensemble, avec 00 une valeur arbitraire, et - chaque application répétée du gradient spatial de champ magnétique de la séquence a, selon une direction spatiale de gradient dite direction de codage identique pour chaque application de gradient de cet ensemble, une intégrale temporelle égale à A non nulle et identique pour chaque application de gradient de cet ensemble, le procédé comprenant en outre une acquisition, pendant au moins une des répétitions de la séquence, d'au moins un signal de résonance magnétique nucléaire. 2. Procédé selon la revendication 1, caractérisé en ce que, au cours des applications répétées de l'impulsion radiofréquence et du gradient spatial de champ magnétique de la séquence d'un même ensemble, la suite un+1 = y',1 -v' est une suite périodique non constante. 3036189 57 3. Procédé selon la revendication 1 ou 2, caractérisé en ce que, au cours des applications répétées de l'impulsion radiofréquence et du gradient spatial de champ magnétique de la séquence d'un même ensemble, l'amplitude de l'impulsion radiofréquence suit une suite périodique non 5 constante. 4. Procédé selon l'une quelconque des revendications précédentes, caractérisé en ce que l'on acquiert le signal de résonance magnétique nucléaire à au moins un état d'équilibre dynamique de l'aimantation de l'échantillon. 10 5. Procédé selon la revendication 4 considérée comme dépendante de la revendication 2 ou 3, caractérisé en ce qu'au cours des applications répétées de l'impulsion radiofréquence et du gradient spatial de champ magnétique de la séquence d'un même ensemble : - la suite un+1 = vn+1 -vn est une suite périodique non constante, et/ou 15 - l'amplitude de l'impulsion radiofréquence suit une suite périodique non constante, de manière à obtenir différents états d'équilibre de l'aimantation de l'échantillon entrelacés au sein de cet ensemble respectivement pour différentes valeurs périodiques de la suite 1'1)1+1 vn+1 lin et/ou de l'amplitude de l'impulsion radiofréquence, et en ce que l'on acquiert le signal de résonance magnétique nucléaire à ces différents états d'équilibre entrelacés. 20 6. Procédé selon la revendication 4 ou 5, caractérisé en ce que l'on émet séquentiellement plusieurs ensembles, avec : - A dont la valeur diffère entre les différents ensembles, et/ou - l'amplitude de l'impulsion radiofréquence qui suit une suite périodique qui diffère entre les 25 différents ensembles, de manière à obtenir différents états d'équilibres de l'aimantation de l'échantillon séquentiels pour différentes valeurs de A et/ou de l'amplitude dont au moins un état d'équilibre par ensemble, et en ce que l'on acquiert le signal de résonance magnétique nucléaire à ces différents états d'équilibre séquentiels. 30 7. Procédé selon la revendication 6, caractérisé en ce qu'il comprend : - une acquisition d'un signal pour deux états d'équilibres correspondants à deux amplitudes de l'impulsion radiofréquence identiques mais deux A différents, 3036189 58 - une comparaison du signal pour ces deux états d'équilibres, et - une déduction, à partir de cette comparaison, de si l'amplitude de l'impulsion radiofréquence produit un angle de bascule effectif plus grand que ou plus petit que ou égal à l'angle de Ernst pour l'échantillon. 5 8. Procédé selon la revendication 7, caractérisé en ce qu'un des deux 4 différents est égal à zéro. 9. Procédé selon l'une quelconque des revendications 4 à 6, caractérisé en ce qu'il comprend : - une acquisition d'un signal pour des états d'équilibres correspondants à des amplitudes différentes de 10 l'impulsion radiofréquence et/ou des 4 différents et/ou différentes valeurs de la suite un_,1 = vn+1 - périodique, et - pour au moins un point dans l'échantillon, une détermination d'une aimantation nucléaire, d'un angle de bascule de l'aimantation, d'un coefficient de diffusion, d'un taux ou temps de relaxation longitudinal R1 ou T 1, et d'un taux ou temps de relaxation transversal R2 ou T2, à partir du signal acquis pour ces états 15 d'équilibres. 10. Procédé selon la revendication 9, caractérisé en ce que la détermination est effectuée, soit par comparaison avec un dictionnaire pré-calculé soit par estimation itérative, selon une minimisation d'une norme de la différence entre le signal acquis exprimé sous forme complexe 20 avec une partie réelle et une partie imaginaire et un modèle du signal exprimé sous forme complexe avec une partie réelle et une partie imaginaire. 11. Procédé selon la revendication 10, caractérisé en ce que la minimisation comprend un algorithme de minimisation au sens des moindres carrés, de préférence en utilisant l'algorithme de Gauss- 25 Newton appliqué aux problèmes non-linéaires. 12. Procédé selon l'une quelconque des revendications 9 à 11, caractérisé en ce qu'il comprend une acquisition d'un signal pour des états d'équilibres correspondants à des 4 différents, dont : - plusieurs points pour une valeur absolue de 4 entre 4=0° inclus et 4=32° inclus modulo 30 360°, et/ou - plusieurs points pour une valeur absolue de 4 entre 4=164° inclus et 4= 196° inclus modulo 360°, et/ou - plusieurs points pour une valeur absolue de 4 entre 4=112° inclus et 4=128° inclus modulo 360°, et/ou 3036189 59 - plusieurs points, pour une valeur absolue de A parmi 45°, 60°, 72°, 90°, et 144° modulo 360°. 13. Procédé selon l'une quelconque des revendications 9 à 12, caractérisé en ce que l'amplitude de 5 l'impulsion radiofréquence correspond toujours à un angle de bascule de l'aimantation : - supérieur à l'angle de Ernst pour un temps de relaxation longitudinal T1 égal à 2000 millisecondes, et - inférieur ou égal à 90°. 10 14. Procédé selon l'une quelconque des revendications 9 à 13, caractérisé en ce que la détermination comprend une condition de continuité spatiale de l'aimantation nucléaire, de l'angle de bascule de l'aimantation, du coefficient de diffusion, du taux ou temps de relaxation longitudinal R1 ou Ti, et du taux ou temps de relaxation transversal R2 ou T2 entre différents points de 1' échantillon. 15 15. Procédé selon l'une quelconque des revendications 1 à 14, caractérisé en ce que l'on émet séquentiellement plusieurs ensembles, et en ce que : - la direction de codage diffère entre différents ensembles, et/ou - la valeur de A diffère entre différents ensembles, et/ou 20 - la forme du gradient de champ magnétique selon la direction de codage diffère entre différents ensembles. 16. Procédé selon la revendication 15, caractérisé en ce qu'il comprend en outre une quantification d'un coefficient de diffusion dans l'échantillon ou une détermination d'une direction privilégiée 25 de diffusion dans l'échantillon ou une pondération de la diffusion dans l'échantillon par exploitation de la différence de A et/ou direction de codage et/ou de forme du gradient entre les différents ensembles. 17. Dispositif (1) d'imagerie par résonance magnétique, comprenant - 30 - des moyens (2) pour appliquer, de manière continue, un champ magnétique principal BO selon un axe Z sur une zone d'échantillon (8), et - des moyens (3) pour émettre un gradient de champ magnétique et des moyens (4) pour émettre une impulsion radiofréquence, et des moyens de commande (5) agencés pour commander les moyens pour émettre un gradient de champ magnétique et les moyens pour émettre une impulsion 3036189 60 radiofréquence, les moyens de commande étant agencés pour ou programmés pour réaliser au moins un ensemble d'applications répétées, sur la zone d'échantillon et selon une période TR, d'une séquence, ladite séquence comprenant : o une impulsion radiofréquence éventuellement d'amplitude et/ou de phase variable à chaque 5 répétition, et o après l'impulsion radiofréquence de la séquence, un gradient spatial de la composante selon l'axe Z du champ magnétique, caractérisé en ce que les moyens de commande sont agencés pour ou programmés pour que, au cours des applications répétées de l'impulsion radiofréquence et du gradient spatial de champ magnétique de la 10 séquence d'un même ensemble : - l'impulsion radiofréquence suive, entre ses différentes applications répétées, une suite périodique : o pour son amplitude et o pour un+1 = y',1 -y' , où n est un entier supérieur ou égal à 1 représentant le numéro de 15 la répétition de la séquence pour cet ensemble, et où un ,1 est une suite de nombres entiers relatifs qui permet de définir la suite y' à une constante arbitraire près y1 tel que yin = v' x A et pn=9,-9, i avec A qui est un nombre réel constant pour toutes les applications répétées de la séquence de cet ensemble et cp, qui est l'incrément entre la phase 0, de l'impulsion radiofréquence à sa n'ème répétition dans cet ensemble et la phase 20 0'_1 de l'impulsion radiofréquence à sa (n-1)1èmC répétition dans cet ensemble, avec 00 une valeur arbitraire, et - chaque application répétée du gradient spatial de champ magnétique de la séquence a, selon une direction spatiale de gradient dite direction de codage identique pour chaque application de gradient de cet ensemble, une intégrale temporelle égale à A non nulle et identique pour chaque 25 application de gradient de cet ensemble, le dispositif comprenant en outre des moyens (6) pour acquérir, pendant au moins une des répétitions de la séquence, au moins un signal de résonance magnétique nucléaire. 30