International audienceThis article contains two contributions. First, we describe the state-of-the-art theories in compressed sensing for Magnetic ResonanceImaging (MRI). This allows us to bring out important principles that should guide the generation of sampling patterns. Second, wedescribe an original methodology to design efficient sampling schemes. It consists of projecting a sampling density on the space of feasiblemeasures for MRI. We end up by proposing comparisons to current sampling strategies on simulated data. This illustrates the well-foundednessof our approach.Cet article a deux finalités. Premièrement, nous proposons un état de l'art des théories d'échantillonnage compressé pour l'Imagerie par résonance magnétique (IRM). Ceci nous permet de dégager quelques grands principes à suivre pour générer des schémas performants en termes de temps d'acquisition et de qualité de reconstruction. Dans une deuxième partie, nous proposons une méthodologie originale de conception de schémas qui repose sur des algorithmes de projection de densités sur des espaces de mesures. Nous proposons finalement des comparaisons avec des stratégies actuelles d'échantillonnage sur des simulations et montrons ainsi le bien-fondé de notre approche. Abstract – This article contains two contributions. First, we describe the state-of-the-art theories in compressed sensing for Magnetic Resonance Imaging (MRI). This allows us to bring out important principles that should guide the generation of sampling patterns. Second, we describe an original methodology to design efficient sampling schemes. It consists of projecting a sampling density on the space of feasible measures for MRI. We end up by proposing comparisons to current sampling strategies on simulated data. This illustrates the well-foundedness of our approach

Boyer, Claire

Chauffert, Nicolas

Ciuciu, Philippe

Kahn, Jonas

Weiss, Pierre

HAL-INSA Toulouse

Sur la génération de schémas d'échantillonnage compressé en IRM

INRIA a CCSD electronic archive server

HAL Id: hal-01255325https://hal.inria.fr/hal-01255325Submitted on 13 Jan 2016HAL is a multi-disciplinary open accessarchive for the deposit and dissemination of sci-entific research documents, whether they are pub-lished or not. The documents may come fromteaching and research institutions in France orabroad, or from public or private research centers.L’archive ouverte pluridisciplinaire HAL, estdestinée au dépôt et à la diffusion de documentsscientifiques de niveau recherche, publiés ou non,émanant des établissements d’enseignement et derecherche français ou étrangers, des laboratoirespublics ou privés.Sur la génération de schémas d’échantillonnagecompressé en IRMPierre Weiss, Nicolas Chauffert, Claire Boyer, Jonas Kahn, Philippe CiuciuTo cite this version:Pierre Weiss, Nicolas Chauffert, Claire Boyer, Jonas Kahn, Philippe Ciuciu. Sur la génération deschémas d’échantillonnage compressé en IRM. GRETSI, Patrice Abry; Paulo Gonçalves, Sep 2015,Lyon, France. pp.4. ￿hal-01255325￿Sur la génération de schémas d’échantillonnage compressé en IRMPierre WEISS1, Nicolas CHAUFFERT2, Claire BOYER3, Jonas KAHN4, Philippe CIUCIU21ITAV (USR3505) et IMT (UMR5219), Université de Toulouse, CNRS.2CEA Neurospin & INRIA Saclay, Parietal team.3IMT (UMR5219), Université de Toulouse, CNRS.4LPP (UMR8524), Université de Lille, CNRS.pierre.weiss@itav.fr, nicolas.chauffert@inria.fr, Claire.Boyer@math.univ-toulouse.fr,jonas.kahn@math.univ-lille1.fr, philippe.ciuciu@cea.frRésumé – Cet article a deux finalités. Premièrement, nous proposons un état de l’art des théories d’échantillonnage compressé pour l’Imageriepar résonance magnétique (IRM). Ceci nous permet de dégager quelques grands principes à suivre pour générer des schémas performantsen termes de temps d’acquisition et de qualité de reconstruction. Dans une deuxième partie, nous proposons une méthodologie originale deconception de schémas qui repose sur des algorithmes de projection de densités sur des espaces de mesures. Nous proposons finalement descomparaisons avec des stratégies actuelles d’échantillonnage sur des simulations et montrons ainsi le bien-fondé de notre approche.Abstract – This article contains two contributions. First, we describe the state-of-the-art theories in compressed sensing for Magnetic Res-onance Imaging (MRI). This allows us to bring out important principles that should guide the generation of sampling patterns. Second, wedescribe an original methodology to design efficient sampling schemes. It consists of projecting a sampling density on the space of feasiblemeasures for MRI. We end up by proposing comparisons to current sampling strategies on simulated data. This illustrates the well-foundednessof our approach.1 IntroductionLa génération de schémas d’échantillonnage compressé effi-caces en Imagerie par Résonance Magnétique (IRM) est unequestion qui reste aujourd’hui largement ouverte d’un pointde vue théorique et applicatif. A notre connaissance, la majo-rité des méthodes implémentées actuellement sur des imageursrepose sur le bon sens, mais n’a pas de fondements mathé-matiques solides. Ceci s’explique par la grande diversité desmodalités d’IRM (anatomique, fonctionnelle, parallèle,...), lacomplexité des contraintes d’acquisition et des images, ainsique les nombreuses difficultés mathématiques pour produireune théorie unifiée capable de prendre en compte ces différentsaspects. Il nous semble qu’il est encore possible d’obtenir desgains significatifs en terme de temps d’échantillonnage ou derésolution spatio-temporelle. Ceci aurait des conséquences trèsimportantes pour la médecine et la biologie.Dans une série de travaux récents [2, 4, 8, 10, 3], nous avonsproposé de nouvelles idées pour combler une partie de ces la-cunes. Bien que beaucoup de questions restent en suspens, cestravaux nous ont permis d’acquérir une bonne intuition deslignes à suivre pour générer des schémas performants. L’ob-jectif de cet article est i) de résumer les conclusions principalesde ces travaux et ii) de proposer une technique de générationde schémas d’échantillonnage originale suivant ces principes.2 NotationsDans ce travail, nous considérons une image u ∈ L2(Ω) àvaleurs dans C. Nous supposons pour simplifier la discussionque Ω = [0, 1]d avec d = 2 ou d = 3. La transformée de Fou-rier d’une image u est notée û. Les images discrétisées sontnotées en gras u ∈ Cn, où n représente le nombre de pixelsou de voxels. La matrice de la transformée de Fourier discrèteest notée F ∈ Cn×n. On suppose de plus que les images sontcompressibles avec une certaine transformée Ψ ∈ Cn×n. Dansles expériences numériques, la transformée Ψ est une transfor-mée en ondelettes orthogonale. La compressibilité signifie quel’essentiel de l’énergie de u est concentré dans un petit nombrede coefficients de la transformée Ψ∗u. L’espace M désignel’ensemble des mesures de probabilité sur Ω. L’espace ∆n estle simplexe de Rn.3 Etat de l’art3.1 Echantillonnage en IRML’IRM permet de mesurer un ensemble de valeurs de la trans-formée de Fourier d’une image u. 1 Les valeurs de la transfor-mée de Fourier sont acquises le long d’une courbe paramétrée1. Ceci est vrai de façon conventionnelle, cependant, on peut modifier lasignification des mesures en utilisant des pulses radio-fréquence particulierss : [0, T ] → Rd qui doit appartenir à un certain ensemble decontraintes de faisabilité ST [10]. Par exemple, dans le cas oùl’acquisition est obtenue à partir d’une seule excitation radio-fréquence, la courbe s doit appartenir à l’ensemble suivant :ST = {s : [0, T ]→ Rd, ‖ṡ‖∞ ≤ α, ‖s̈‖∞ ≤ β}. (1)Les paramètres α et β dépendent de l’imageur. Si l’acquisitionse fait à partir de plusieurs excitations radio-fréquences, s ap-partient à des produits cartésiens de l’ensemble ST .Après une acquisition en IRM, on a donc accès à l’ensembleE de mesures :E ={û(s(i∆t)), 1 ≤ i ≤⌊T∆t⌋}(2)où ∆t > 0 décrit la période d’échantillonnage de la courbe s.Dans ce papier, nous négligeons les dégradations de mesuresliées au bruit ou aux inhomogénéité des champs magnétiques.Maintenant que les possibilités d’une IRM sont posées, nouspouvons présenter les questions typiques qui se posent en échan-tillonnage pour l’IRM. A T fixé comment trouver une courbes ∈ ST optimale, au sens où elle permet d’obtenir la recons-truction la plus fidèle ? Quel est le temps d’acquisition T mi-nimal qui permette d’obtenir une précision de reconstructiondonnée ? Comment reconstruire les images connaissant E ?Ces questions semblent extrêmement complexes mathémati-quement, ce qui explique la prédominance d’heuristiques [13]dans la littérature. Dans la suite, nous présentons quelques ré-sultats théoriques qui permettent tout de même de guider laconception d’une trajectoire.3.2 Echantillonnage compressé usuelDans cette partie, nous décrivons la théorie de l’échantillon-nage compressé telle qu’elle a été posée dans les papiers fon-dateurs [7] et plus récemment [6]. Les auteurs considèrent unematrice orthogonale A0 = [a∗1; . . . ; a∗n] où les vecteurs a∗i sontles lignes du dictionnaire (le séparateur ; indique une conca-ténation verticale). Ils proposent de construire une matrice demesure aléatoire A = [a∗J1 ; . . . ; a∗Jm] où les quantités Jk ∈{1, . . . , n} sont des variables aléatoires i.i.d. uniformes. Connais-sant y = Ax, les auteurs proposent de reconstruire une estiméex̄ de x en résolvant le problème `1 suivant :x̄ ∈ ArgminAx=y‖x‖1Le type de théorème qu’ils obtiennent dans [6] est de laforme suivante 2.Théorème 1. Supposons que x est s-parcimonieux, i.e. qu’il necontient que s composantes non nulles parmi n. Si le nombre[11], l’IRM parallèle [12] ou encore en modifiant l’utilisation des correcteursde champs magnétiques [14].2. Notez que nous ne parlons à aucun moment de Restricted Isometry Pro-perty (RIP). Ce type de notion est très répandu car il permet la reconstructionde tous les signaux s-parcimonieux. Cette propriété - bien que très plaisante -est extrêmement exigeante. Il nous semble clair qu’elle doit être abandonnée àl’avenir pour obtenir des résultats théoriques convaincants en IRM.de mesures m satisfait :m ≥ Cs(n max1≤k≤n‖ak‖2∞)log(nε),où C > 0 est une constante universelle, alors x̄ = x avecprobabilité 1− ε.Ce théorème et tous ceux annoncés dans le papier peuventde plus être étendus au cas de données bruitées en relaxant lacontrainte Ax = y. La cohérence [6], κ(A0) = n max1≤k≤n‖ak‖2∞est comprise entre 1 et n. En particulier, κ(F) = 1 et κ(Id) =n. Dans le cas favorable de la transformée de Fourier, ce théo-rème indique donc qu’il suffit d’environ s log(nε)mesures pourobtenir une reconstruction exacte d’un signal s-parcimonieux.Bien que ce type de théorème ait eu un impact énorme dansla littérature, il n’est pas applicable en IRM. En effet, la trans-formée naturelle en IRM est A0 = F∗Ψ, i.e. le produit d’unetransformée de Fourier et d’une transformée en ondelettes. Onpeut montrer que κ(F∗Ψ) = O(n).3.3 Les avancées récentesDans cette partie, nous indiquons les avancées récentes enéchantillonnage compressé qui nous semblent les plus impor-tantes. Des omissions sont possibles car une telle entreprisecontient nécessairement une bonne dose de subjectivité.3.3.1 Echantillonnage à densité variableDans la majorité des applications réelles, les transforméesà utiliser sont cohérentes. C’est notamment le cas de l’IRM.Une technique possible pour briser la cohérence consiste à tirerles échantillons cohérents plus souvent. Précisons cette idée.Soit π ∈ ∆n la distribution des variables aléatoires Jk. Lethéorème suivant [8] justifie l’utilisation d’échantillonnage àdensité variable. Notons que cette technique est la base (nonjustifiée) de l’article fondateur [13].Théorème 2. Supposons que x est s-parcimonieux, i.e. qu’il necontient que s composantes non nulles parmi n. Posons πk =‖ak‖2∞∑ni=1 ‖ai‖2∞. Si le nombre de mesures m satisfait :m ≥ Cs(n∑i=1‖ai‖2∞)log(nε),où C > 0 est une constante universelle, alors x̄ = x avecprobabilité 1− ε.On peut montrer que pour le cas de l’IRM (i.e. A0 = F∗Ψ),∑ni=1 ‖ai‖2∞ ∝ log(n). Ainsi, il suffit de O(s log(n)2) coef-ficients de Fourier pour pouvoir reconstruire exactement uneimage s-parcimonieuse.3.3.2 Parcimonie structuréeBien que le théorème (2) soit plutôt attrayant théoriquement,des expériences numériques [1] montrent qu’il est encore in-suffisant pour une application en IRM : les constantes apparais-sant dans lesO sont en effet grandes, et même le terme log(n)2n’est pas négligeable en pratique. Une avancée récente qui noussemble importante est le papier [1]. Les auteurs montrent qu’ilest possible d’exploiter une parcimonie structurée (e.g. plus decoefficients d’ondelettes non nuls dans les sous-bandes corres-pondant aux basses fréquences) pour obtenir de meilleures ga-ranties théoriques. Nous ne décrivons pas les théorèmes résul-tants par manque de place. Ils permettent cependant d’obte-nir deux conclusions importantes. Premièrement, il faut utili-ser des schémas d’échantillonnage à densité variable commedans le paragraphe précédent. Les densités doivent dépendre àla fois de la cohérence locale ‖ak‖∞ et de l’a priori qu’on sedonne sur le support de x. Deuxièmement, comme la parcimo-nie est asymptotique (elle n’apparait que dans les sous-bandeshautes-fréquences), l’acquisition comprimée en IRM peut dif-ficilement être utilisée comme une technique pour gagner dutemps. Elle permet par contre - à temps d’acquisition constant- d’obtenir de meilleures résolutions spatiales.3.3.3 Acquisition structuréeTous les résultats donnés jusqu’à présent ne contiennent au-cune notion de structure dans l’acquisition. Ils ne sont donc pasapplicables en pratique. A notre connaissance, les seuls travauxcommençant à proposer des solutions d’échantillonnage plau-sibles et justifiées théoriquement sont [8, 3].Dans [8] nous avons d’abord proposé une définition géné-rique précise de ce qu’est un échantillonneur à densité variable :c’est un processus dont la mesure empirique convergence versune densité cible quand le temps tend vers l’infini. Nous avonsanalysé l’utilisation de tels échantillonneurs théoriquement eten avons déduit deux propriétés essentielles pour qu’ils per-mettent des reconstruction précises :1. La mesure empirique du processus doit correspondre àune mesure cible, i.e. il faut effectuer un échantillonnageà densité variable.2. Le processus ne sera efficace que si son temps de mé-lange (i.e. sa rapidité à converger vers la densité cible)est faible.Ces résultats indiquent qu’un bon échantillonneur doit i) cou-vrir l’espace rapidement (temps de mélange) ii) avoir une den-sité cible fixée et iii) appartenir à un ensemble admissible pourl’IRM. Ces conclusions sont très proches des travaux plus heu-ristiques menés par les leaders mondiaux en échantillonnagecompressé pour l’IRM [15]. Dans la suite de ce travail, nousproposons un algorithme original permettant de générer de telsschémas automatiquement.4 Génération automatique de schémasd’échantillonnageNous supposons connue une densité cible π ∈ M. Celle-cipeut être obtenue par des considérations heuristiques [13] ou enutilisant des théories de l’échantillonnage compressé [1, 8, 3].Nous fixons aussi le temps total d’acquisition T et cherchons àtrouver une courbe d’échantillonnage s ∈ ST efficace.À une courbe s ∈ ST , nous pouvons associer une mesureimage µs ∈M(ST ) définie pour tout Borélien B ⊆ Ω parµs(B) = γ(s−1(B))où γ est la mesure de Lebesgue normalisée par T . Ainsi µs(B)mesure le temps relatif que s passe dans l’ensembleB. On sou-haite que µs ' π et il faut donc définir une notion de distanceentre mesures. Dans un travail récent [9], nous avons proposéde définir une distance entre mesures comme suit :dist(µs, π) = ‖h ? (µs − π)‖2,où h est un noyau de convolution continu sur Ω. Si la transfor-mée de Fourier de h est réelle et non nulle partout, on a montréque dist est effectivement une distance et qu’elle métrise laconvergence faible entre mesures. Elle est de plus majorée parla distance de Wasserstein si h est Lipschitz continue.L’idée principale de ce papier est de définir s comme lacourbe paramétrisant la solution du problème variationnel sui-vant :minµs∈M(ST )12‖h ? (µs − π)‖22. (3)Ce problème de projection consiste à chercher la courbe s ad-missible dont la mesure approche au mieux la densité cible π.Nous avons montré dans [9] que des algorithmes de pro-grammation non-linéaire permettent de trouver des point cri-tiques de la fonctionnelle (3). De plus, ces algorithmes donnentdes résultats de projection d’une bonne qualité quelque soitl’initialisation. Nous renvoyons le lecteur intéressé à ce papierpour plus de détails.5 RésultatsPour évaluer la pertinence, de cette approche, nous travaillonssur un fantôme couramment utilisé en IRM de taille 10242, voirFig. 1 (a). Nous projetons une densité cible avec décroissanceradiale sur : des sommes de N Diracs, Fig. 2 (b), les mesuresempiriques de courbes dans l’ensemble ST défini en (1), Fig.2 (c), la mesure empirique d’un ensemble de 200 segments delongueurs arbitraires, Fig. 2 (d). L’exemple de la Fig. 2 (a) cor-respond à des tirages i.i.d. Il sert uniquement de référence car iln’est pas implémentable, mais correspond aux approches stan-dards utilisées en simulation. Nous indiquons uniquement lesSNR des images reconstruites par manque de place. On peutdonner plusieurs conclusions. i) Utiliser une répulsion entrepoints (telle que le Poisson disk sampling [5]) améliore trèslégèrement les résultats par rapport à des tirages i.i.d. ii) Leschéma de projection permet de retrouver des schémas corres-pondant à l’état de l’art (e.g. spirales bruitées) automatique-ment, Fig. 2 (c). iii) Le schéma de projection permet aussi detrouver de nouveaux schémas inédits (e.g. lignes de longueurvariables, Fig. 2 (d)). iv) Ajouter des contraintes cinématiquessur les trajectoires semble réduire la qualité de reconstruction,mais la baisse est raisonnable (~1dB).Les schémas proposés sont en cours d’implémentation surles scanners pré-cliniques Bruker 7T et 17T du projet Neuros-pin.(a) (b)FIGURE 1 – (a) image à reconstruire de taille 1024x1024. (b)Image reconstruite avec le schéma de Fig. 2 (a), SNR=19.6dB.Références[1] B. Adcock, A. C. Hansen, C. Poon, and B. Roman. Brea-king the coherence barrier : asymptotic incoherence andasymptotic sparsity in compressed sensing. Preprint,2013.[2] J. Bigot, C. Boyer, and P. Weiss. An analysis of blocksampling strategies in compressed sensing. arXiv preprintarXiv :1305.4446, 2013.[3] C. Boyer, J. Bigot, and P. Weiss. Compressed sensingwith structured sparsity and block structured acquisition.arXiv preprint, 2015.[4] C. Boyer, P. Weiss, and J. Bigot. An algorithm for va-riable density sampling with block-constrained acquisi-tion. SIAM Journal on Imaging Sciences, 7(2) :1080–1107, 2014.[5] R. Bridson. Fast poisson disk sampling in arbitrary di-mensions. In ACM SIGGRAPH, volume 2007, page 5,2007.[6] E. Candès and Y. Plan. A probabilistic and ripless theoryof compressed sensing. Information Theory, IEEE Tran-sactions on, 57(11) :7235–7254, 2011.[7] E. Candès, J. Romberg, and T. Tao. Stable signal recoveryfrom incomplete and inaccurate measurements. Commu-nications on pure and applied mathematics, 59(8) :1207–1223, 2006.[8] N. Chauffert, P. Ciuciu, J. Kahn, and P. Weiss. Variabledensity sampling with continuous sampling trajectories.SIAM Journal on Imaging Sciences, 2014.[9] N. Chauffert, P. Ciuciu, J. Kahn, and P. Weiss. A projec-tion method on measures sets. submitted, 2015.(a) (b)(c) (d)FIGURE 2 – Exemples de schémas d’échantillonnage gé-nérés avec 50000 échantillons, soit 4.7% des mesures. (a)Tirages i.i.d. (infaisable - SNR=19.6dB) (b) projection surles mesures discrètes (infaisable - SNR=19.8dB) (c) projec-tion sur les mesures admissibles pour une trajectoire d’IRM(SNR=18.6dB). (d) projection sur 200 segments de longueurvariable (SNR=18.9dB).[10] N. Chauffert, P. Weiss, J. Kahn, and P. Ciuciu. Gradientwaveform design for variable density sampling in magne-tic resonance imaging. arXiv preprint arXiv :1412.4621,2014.[11] J. P. Haldar, D. Hernando, and Z.-P. Liang. Compressed-sensing MRI with random encoding. Medical Imaging,IEEE Transactions on, 30(4) :893–903, 2011.[12] D. J. Larkman and R. G. Nunes. Parallel magneticresonance imaging. Physics in medicine and biology,52(7) :R15, 2007.[13] M. Lustig, D. Donoho, and J. M. Pauly. Sparse MRI :The application of compressed sensing for rapid mr ima-ging. Magnetic resonance in medicine, 58(6) :1182–1195, 2007.[14] G. Puy, J. P. Marques, R. Gruetter, J. Thiran, D. VanDe Ville, P. Vandergheynst, and Y. Wiaux. Spread spec-trum magnetic resonance imaging. Medical Imaging,IEEE Transactions on, 31(3) :586–598, 2012.[15] S. Vasanawala, M. Murphy, M. T. Alley, P. Lai, K. Keut-zer, J. M. Pauly, and M. Lustig. Practical parallel imagingcompressed sensing MRI : Summary of two years of ex-perience in accelerating body MRI of pediatric patients.In Proc. ISBI, pages 1039–1043, 2011.

HAL Descartes

Scientific Publications of the University of Toulouse II Le Mirail

Sur la ge´ne´ration de sche´mas d’e´chantillonnagecompresse´ en IRMPierre Weiss, Nicolas Chauffert, Claire Boyer, Jonas Kahn, Philippe CiuciuTo cite this version:Pierre Weiss, Nicolas Chauffert, Claire Boyer, Jonas Kahn, Philippe Ciuciu. Sur la ge´ne´rationde sche´mas d’e´chantillonnage compresse´ en IRM. GRETSI, Sep 2015, Lyon, France. pp.4, 2015,25 e`me actes du colloque GRETSI. <hal-01255325>HAL Id: hal-01255325https://hal.inria.fr/hal-01255325Submitted on 13 Jan 2016HAL is a multi-disciplinary open accessarchive for the deposit and dissemination of sci-entific research documents, whether they are pub-lished or not. The documents may come fromteaching and research institutions in France orabroad, or from public or private research centers.L’archive ouverte pluridisciplinaire HAL, estdestine´e au de´poˆt et a` la diffusion de documentsscientifiques de niveau recherche, publie´s ou non,e´manant des e´tablissements d’enseignement et derecherche franc¸ais ou e´trangers, des laboratoirespublics ou prive´s.Sur la génération de schémas d’échantillonnage compressé en IRMPierre WEISS1, Nicolas CHAUFFERT2, Claire BOYER3, Jonas KAHN4, Philippe CIUCIU21ITAV (USR3505) et IMT (UMR5219), Université de Toulouse, CNRS.2CEA Neurospin & INRIA Saclay, Parietal team.3IMT (UMR5219), Université de Toulouse, CNRS.4LPP (UMR8524), Université de Lille, CNRS.pierre.weiss@itav.fr, nicolas.chauffert@inria.fr, Claire.Boyer@math.univ-toulouse.fr,jonas.kahn@math.univ-lille1.fr, philippe.ciuciu@cea.frRésumé – Cet article a deux finalités. Premièrement, nous proposons un état de l’art des théories d’échantillonnage compressé pour l’Imageriepar résonance magnétique (IRM). Ceci nous permet de dégager quelques grands principes à suivre pour générer des schémas performantsen termes de temps d’acquisition et de qualité de reconstruction. Dans une deuxième partie, nous proposons une méthodologie originale deconception de schémas qui repose sur des algorithmes de projection de densités sur des espaces de mesures. Nous proposons finalement descomparaisons avec des stratégies actuelles d’échantillonnage sur des simulations et montrons ainsi le bien-fondé de notre approche.Abstract – This article contains two contributions. First, we describe the state-of-the-art theories in compressed sensing for Magnetic Res-onance Imaging (MRI). This allows us to bring out important principles that should guide the generation of sampling patterns. Second, wedescribe an original methodology to design efficient sampling schemes. It consists of projecting a sampling density on the space of feasiblemeasures for MRI. We end up by proposing comparisons to current sampling strategies on simulated data. This illustrates the well-foundednessof our approach.1 IntroductionLa génération de schémas d’échantillonnage compressé effi-caces en Imagerie par Résonance Magnétique (IRM) est unequestion qui reste aujourd’hui largement ouverte d’un pointde vue théorique et applicatif. A notre connaissance, la majo-rité des méthodes implémentées actuellement sur des imageursrepose sur le bon sens, mais n’a pas de fondements mathé-matiques solides. Ceci s’explique par la grande diversité desmodalités d’IRM (anatomique, fonctionnelle, parallèle,...), lacomplexité des contraintes d’acquisition et des images, ainsique les nombreuses difficultés mathématiques pour produireune théorie unifiée capable de prendre en compte ces différentsaspects. Il nous semble qu’il est encore possible d’obtenir desgains significatifs en terme de temps d’échantillonnage ou derésolution spatio-temporelle. Ceci aurait des conséquences trèsimportantes pour la médecine et la biologie.Dans une série de travaux récents [2, 4, 8, 10, 3], nous avonsproposé de nouvelles idées pour combler une partie de ces la-cunes. Bien que beaucoup de questions restent en suspens, cestravaux nous ont permis d’acquérir une bonne intuition deslignes à suivre pour générer des schémas performants. L’ob-jectif de cet article est i) de résumer les conclusions principalesde ces travaux et ii) de proposer une technique de générationde schémas d’échantillonnage originale suivant ces principes.2 NotationsDans ce travail, nous considérons une image u ∈ L2(Ω) àvaleurs dans C. Nous supposons pour simplifier la discussionque Ω = [0, 1]d avec d = 2 ou d = 3. La transformée de Fou-rier d’une image u est notée uˆ. Les images discrétisées sontnotées en gras u ∈ Cn, où n représente le nombre de pixelsou de voxels. La matrice de la transformée de Fourier discrèteest notée F ∈ Cn×n. On suppose de plus que les images sontcompressibles avec une certaine transformée Ψ ∈ Cn×n. Dansles expériences numériques, la transformée Ψ est une transfor-mée en ondelettes orthogonale. La compressibilité signifie quel’essentiel de l’énergie de u est concentré dans un petit nombrede coefficients de la transformée Ψ∗u. L’espace M désignel’ensemble des mesures de probabilité sur Ω. L’espace ∆n estle simplexe de Rn.3 Etat de l’art3.1 Echantillonnage en IRML’IRM permet de mesurer un ensemble de valeurs de la trans-formée de Fourier d’une image u. 1 Les valeurs de la transfor-mée de Fourier sont acquises le long d’une courbe paramétrée1. Ceci est vrai de façon conventionnelle, cependant, on peut modifier lasignification des mesures en utilisant des pulses radio-fréquence particulierss : [0, T ] → Rd qui doit appartenir à un certain ensemble decontraintes de faisabilité ST [10]. Par exemple, dans le cas oùl’acquisition est obtenue à partir d’une seule excitation radio-fréquence, la courbe s doit appartenir à l’ensemble suivant :ST = {s : [0, T ]→ Rd, ‖s˙‖∞ ≤ α, ‖s¨‖∞ ≤ β}. (1)Les paramètres α et β dépendent de l’imageur. Si l’acquisitionse fait à partir de plusieurs excitations radio-fréquences, s ap-partient à des produits cartésiens de l’ensemble ST .Après une acquisition en IRM, on a donc accès à l’ensembleE de mesures :E ={uˆ(s(i∆t)), 1 ≤ i ≤⌊T∆t⌋}(2)où ∆t > 0 décrit la période d’échantillonnage de la courbe s.Dans ce papier, nous négligeons les dégradations de mesuresliées au bruit ou aux inhomogénéité des champs magnétiques.Maintenant que les possibilités d’une IRM sont posées, nouspouvons présenter les questions typiques qui se posent en échan-tillonnage pour l’IRM. A T fixé comment trouver une courbes ∈ ST optimale, au sens où elle permet d’obtenir la recons-truction la plus fidèle ? Quel est le temps d’acquisition T mi-nimal qui permette d’obtenir une précision de reconstructiondonnée ? Comment reconstruire les images connaissant E ?Ces questions semblent extrêmement complexes mathémati-quement, ce qui explique la prédominance d’heuristiques [13]dans la littérature. Dans la suite, nous présentons quelques ré-sultats théoriques qui permettent tout de même de guider laconception d’une trajectoire.3.2 Echantillonnage compressé usuelDans cette partie, nous décrivons la théorie de l’échantillon-nage compressé telle qu’elle a été posée dans les papiers fon-dateurs [7] et plus récemment [6]. Les auteurs considèrent unematrice orthogonale A0 = [a∗1; . . . ; a∗n] où les vecteurs a∗i sontles lignes du dictionnaire (le séparateur ; indique une conca-ténation verticale). Ils proposent de construire une matrice demesure aléatoire A = [a∗J1 ; . . . ; a∗Jm] où les quantités Jk ∈{1, . . . , n} sont des variables aléatoires i.i.d. uniformes. Connais-sant y = Ax, les auteurs proposent de reconstruire une estiméex¯ de x en résolvant le problème `1 suivant :x¯ ∈ ArgminAx=y‖x‖1Le type de théorème qu’ils obtiennent dans [6] est de laforme suivante 2.Théorème 1. Supposons que x est s-parcimonieux, i.e. qu’il necontient que s composantes non nulles parmi n. Si le nombre[11], l’IRM parallèle [12] ou encore en modifiant l’utilisation des correcteursde champs magnétiques [14].2. Notez que nous ne parlons à aucun moment de Restricted Isometry Pro-perty (RIP). Ce type de notion est très répandu car il permet la reconstructionde tous les signaux s-parcimonieux. Cette propriété - bien que très plaisante -est extrêmement exigeante. Il nous semble clair qu’elle doit être abandonnée àl’avenir pour obtenir des résultats théoriques convaincants en IRM.de mesuresm satisfait :m ≥ Cs(n max1≤k≤n‖ak‖2∞)log(n),où C > 0 est une constante universelle, alors x¯ = x avecprobabilité 1− .Ce théorème et tous ceux annoncés dans le papier peuventde plus être étendus au cas de données bruitées en relaxant lacontrainte Ax = y. La cohérence [6], κ(A0) = n max1≤k≤n‖ak‖2∞est comprise entre 1 et n. En particulier, κ(F) = 1 et κ(Id) =n. Dans le cas favorable de la transformée de Fourier, ce théo-rème indique donc qu’il suffit d’environ s log(n)mesures pourobtenir une reconstruction exacte d’un signal s-parcimonieux.Bien que ce type de théorème ait eu un impact énorme dansla littérature, il n’est pas applicable en IRM. En effet, la trans-formée naturelle en IRM est A0 = F∗Ψ, i.e. le produit d’unetransformée de Fourier et d’une transformée en ondelettes. Onpeut montrer que κ(F∗Ψ) = O(n).3.3 Les avancées récentesDans cette partie, nous indiquons les avancées récentes enéchantillonnage compressé qui nous semblent les plus impor-tantes. Des omissions sont possibles car une telle entreprisecontient nécessairement une bonne dose de subjectivité.3.3.1 Echantillonnage à densité variableDans la majorité des applications réelles, les transforméesà utiliser sont cohérentes. C’est notamment le cas de l’IRM.Une technique possible pour briser la cohérence consiste à tirerles échantillons cohérents plus souvent. Précisons cette idée.Soit pi ∈ ∆n la distribution des variables aléatoires Jk. Lethéorème suivant [8] justifie l’utilisation d’échantillonnage àdensité variable. Notons que cette technique est la base (nonjustifiée) de l’article fondateur [13].Théorème 2. Supposons que x est s-parcimonieux, i.e. qu’il necontient que s composantes non nulles parmi n. Posons pik =‖ak‖2∞∑ni=1 ‖ai‖2∞ . Si le nombre de mesuresm satisfait :m ≥ Cs(n∑i=1‖ai‖2∞)log(n),où C > 0 est une constante universelle, alors x¯ = x avecprobabilité 1− .On peut montrer que pour le cas de l’IRM (i.e. A0 = F∗Ψ),∑ni=1 ‖ai‖2∞ ∝ log(n). Ainsi, il suffit de O(s log(n)2) coef-ficients de Fourier pour pouvoir reconstruire exactement uneimage s-parcimonieuse.3.3.2 Parcimonie structuréeBien que le théorème (2) soit plutôt attrayant théoriquement,des expériences numériques [1] montrent qu’il est encore in-suffisant pour une application en IRM : les constantes apparais-sant dans lesO sont en effet grandes, et même le terme log(n)2n’est pas négligeable en pratique. Une avancée récente qui noussemble importante est le papier [1]. Les auteurs montrent qu’ilest possible d’exploiter une parcimonie structurée (e.g. plus decoefficients d’ondelettes non nuls dans les sous-bandes corres-pondant aux basses fréquences) pour obtenir de meilleures ga-ranties théoriques. Nous ne décrivons pas les théorèmes résul-tants par manque de place. Ils permettent cependant d’obte-nir deux conclusions importantes. Premièrement, il faut utili-ser des schémas d’échantillonnage à densité variable commedans le paragraphe précédent. Les densités doivent dépendre àla fois de la cohérence locale ‖ak‖∞ et de l’a priori qu’on sedonne sur le support de x. Deuxièmement, comme la parcimo-nie est asymptotique (elle n’apparait que dans les sous-bandeshautes-fréquences), l’acquisition comprimée en IRM peut dif-ficilement être utilisée comme une technique pour gagner dutemps. Elle permet par contre - à temps d’acquisition constant- d’obtenir de meilleures résolutions spatiales.3.3.3 Acquisition structuréeTous les résultats donnés jusqu’à présent ne contiennent au-cune notion de structure dans l’acquisition. Ils ne sont donc pasapplicables en pratique. A notre connaissance, les seuls travauxcommençant à proposer des solutions d’échantillonnage plau-sibles et justifiées théoriquement sont [8, 3].Dans [8] nous avons d’abord proposé une définition géné-rique précise de ce qu’est un échantillonneur à densité variable :c’est un processus dont la mesure empirique convergence versune densité cible quand le temps tend vers l’infini. Nous avonsanalysé l’utilisation de tels échantillonneurs théoriquement eten avons déduit deux propriétés essentielles pour qu’ils per-mettent des reconstruction précises :1. La mesure empirique du processus doit correspondre àune mesure cible, i.e. il faut effectuer un échantillonnageà densité variable.2. Le processus ne sera efficace que si son temps de mé-lange (i.e. sa rapidité à converger vers la densité cible)est faible.Ces résultats indiquent qu’un bon échantillonneur doit i) cou-vrir l’espace rapidement (temps de mélange) ii) avoir une den-sité cible fixée et iii) appartenir à un ensemble admissible pourl’IRM. Ces conclusions sont très proches des travaux plus heu-ristiques menés par les leaders mondiaux en échantillonnagecompressé pour l’IRM [15]. Dans la suite de ce travail, nousproposons un algorithme original permettant de générer de telsschémas automatiquement.4 Génération automatique de schémasd’échantillonnageNous supposons connue une densité cible pi ∈ M. Celle-cipeut être obtenue par des considérations heuristiques [13] ou enutilisant des théories de l’échantillonnage compressé [1, 8, 3].Nous fixons aussi le temps total d’acquisition T et cherchons àtrouver une courbe d’échantillonnage s ∈ ST efficace.À une courbe s ∈ ST , nous pouvons associer une mesureimage µs ∈M(ST ) définie pour tout Borélien B ⊆ Ω parµs(B) = γ(s−1(B))où γ est la mesure de Lebesgue normalisée par T . Ainsi µs(B)mesure le temps relatif que s passe dans l’ensembleB. On sou-haite que µs ' pi et il faut donc définir une notion de distanceentre mesures. Dans un travail récent [9], nous avons proposéde définir une distance entre mesures comme suit :dist(µs, pi) = ‖h ? (µs − pi)‖2,où h est un noyau de convolution continu sur Ω. Si la transfor-mée de Fourier de h est réelle et non nulle partout, on a montréque dist est effectivement une distance et qu’elle métrise laconvergence faible entre mesures. Elle est de plus majorée parla distance de Wasserstein si h est Lipschitz continue.L’idée principale de ce papier est de définir s comme lacourbe paramétrisant la solution du problème variationnel sui-vant :minµs∈M(ST )12‖h ? (µs − pi)‖22. (3)Ce problème de projection consiste à chercher la courbe s ad-missible dont la mesure approche au mieux la densité cible pi.Nous avons montré dans [9] que des algorithmes de pro-grammation non-linéaire permettent de trouver des point cri-tiques de la fonctionnelle (3). De plus, ces algorithmes donnentdes résultats de projection d’une bonne qualité quelque soitl’initialisation. Nous renvoyons le lecteur intéressé à ce papierpour plus de détails.5 RésultatsPour évaluer la pertinence, de cette approche, nous travaillonssur un fantôme couramment utilisé en IRM de taille 10242, voirFig. 1 (a). Nous projetons une densité cible avec décroissanceradiale sur : des sommes de N Diracs, Fig. 2 (b), les mesuresempiriques de courbes dans l’ensemble ST défini en (1), Fig.2 (c), la mesure empirique d’un ensemble de 200 segments delongueurs arbitraires, Fig. 2 (d). L’exemple de la Fig. 2 (a) cor-respond à des tirages i.i.d. Il sert uniquement de référence car iln’est pas implémentable, mais correspond aux approches stan-dards utilisées en simulation. Nous indiquons uniquement lesSNR des images reconstruites par manque de place. On peutdonner plusieurs conclusions. i) Utiliser une répulsion entrepoints (telle que le Poisson disk sampling [5]) améliore trèslégèrement les résultats par rapport à des tirages i.i.d. ii) Leschéma de projection permet de retrouver des schémas corres-pondant à l’état de l’art (e.g. spirales bruitées) automatique-ment, Fig. 2 (c). iii) Le schéma de projection permet aussi detrouver de nouveaux schémas inédits (e.g. lignes de longueurvariables, Fig. 2 (d)). iv) Ajouter des contraintes cinématiquessur les trajectoires semble réduire la qualité de reconstruction,mais la baisse est raisonnable (~1dB).Les schémas proposés sont en cours d’implémentation surles scanners pré-cliniques Bruker 7T et 17T du projet Neuros-pin.(a) (b)FIGURE 1 – (a) image à reconstruire de taille 1024x1024. (b)Image reconstruite avec le schéma de Fig. 2 (a), SNR=19.6dB.Références[1] B. Adcock, A. C. Hansen, C. Poon, and B. Roman. Brea-king the coherence barrier : asymptotic incoherence andasymptotic sparsity in compressed sensing. Preprint,2013.[2] J. Bigot, C. Boyer, and P. Weiss. An analysis of blocksampling strategies in compressed sensing. arXiv preprintarXiv :1305.4446, 2013.[3] C. Boyer, J. Bigot, and P. Weiss. Compressed sensingwith structured sparsity and block structured acquisition.arXiv preprint, 2015.[4] C. Boyer, P. Weiss, and J. Bigot. An algorithm for va-riable density sampling with block-constrained acquisi-tion. SIAM Journal on Imaging Sciences, 7(2) :1080–1107, 2014.[5] R. Bridson. Fast poisson disk sampling in arbitrary di-mensions. In ACM SIGGRAPH, volume 2007, page 5,2007.[6] E. Candès and Y. Plan. A probabilistic and ripless theoryof compressed sensing. Information Theory, IEEE Tran-sactions on, 57(11) :7235–7254, 2011.[7] E. Candès, J. Romberg, and T. Tao. Stable signal recoveryfrom incomplete and inaccurate measurements. Commu-nications on pure and applied mathematics, 59(8) :1207–1223, 2006.[8] N. Chauffert, P. Ciuciu, J. Kahn, and P. Weiss. Variabledensity sampling with continuous sampling trajectories.SIAM Journal on Imaging Sciences, 2014.[9] N. Chauffert, P. Ciuciu, J. Kahn, and P. Weiss. A projec-tion method on measures sets. submitted, 2015.(a) (b)(c) (d)FIGURE 2 – Exemples de schémas d’échantillonnage gé-nérés avec 50000 échantillons, soit 4.7% des mesures. (a)Tirages i.i.d. (infaisable - SNR=19.6dB) (b) projection surles mesures discrètes (infaisable - SNR=19.8dB) (c) projec-tion sur les mesures admissibles pour une trajectoire d’IRM(SNR=18.6dB). (d) projection sur 200 segments de longueurvariable (SNR=18.9dB).[10] N. Chauffert, P. Weiss, J. Kahn, and P. Ciuciu. Gradientwaveform design for variable density sampling in magne-tic resonance imaging. arXiv preprint arXiv :1412.4621,2014.[11] J. P. Haldar, D. Hernando, and Z.-P. Liang. Compressed-sensing MRI with random encoding. Medical Imaging,IEEE Transactions on, 30(4) :893–903, 2011.[12] D. J. Larkman and R. G. Nunes. Parallel magneticresonance imaging. Physics in medicine and biology,52(7) :R15, 2007.[13] M. Lustig, D. Donoho, and J. M. Pauly. Sparse MRI :The application of compressed sensing for rapid mr ima-ging. Magnetic resonance in medicine, 58(6) :1182–1195, 2007.[14] G. Puy, J. P. Marques, R. Gruetter, J. Thiran, D. VanDe Ville, P. Vandergheynst, and Y. Wiaux. Spread spec-trum magnetic resonance imaging. Medical Imaging,IEEE Transactions on, 31(3) :586–598, 2012.[15] S. Vasanawala, M. Murphy, M. T. Alley, P. Lai, K. Keut-zer, J. M. Pauly, and M. Lustig. Practical parallel imagingcompressed sensing MRI : Summary of two years of ex-perience in accelerating body MRI of pediatric patients.In Proc. ISBI, pages 1039–1043, 2011.

Sur la génération de schémas d'échantillonnage compressé en IRM

Abstract

Similar works

Full text

Available Versions

HAL-INSA Toulouse

INRIA a CCSD electronic archive server

HAL Descartes

Scientific Publications of the University of Toulouse II Le Mirail

HAL-CEA