PROCEDE DE TRAITEMENT D'IMAGE BASE SUR LA TECHNIQUE DES ELEMENTS FINIS POUR LA RESOLUTION DIRECTE DE PROBLEMES INVERSES EN MECANIQUE DES STRUCTURES DOMAINE TECHNIQUE La présente invention concerne le domaine technique général du traitement d'image en mécanique des structures, et notamment du traitement d'image pour la résolution des problèmes inverses en élasticité linéaire. Un exemple non limitatif d'application de la présente invention concerne le domaine de la reconstruction de l'élasticité de lésions athérosclérotiques pour l'estimation d'un risque de rupture d'une plaque d'athérome. Toutefois, il est bien évident pour l'homme du métier que la présente invention ne se limite pas à cette application particulière. PRESENTATION GENERALE DE L'ART ANTERIEUR 1. Conséquence d'une rupture d'une plaque d'athérome 20 L'athérome ou athérosclérose correspond à un remaniement de l'intima des artères de gros et moyen calibre (aorte, artères coronaires, carotides, artères cérébrales, artères des membres inférieurs, etc.) par accumulation segmentaire de lipides, glucides complexes, sang et produits sanguins, tissus adipeux, dépôts calcaires et autres minéraux au sein de l'intima (qui est la couche interne du vaisseau). 25 Cette pathologie vasculaire est généralement à progression lente (sur des décennies). Elle peut se stabiliser et ne pas présenter de danger notable pour le patient. Mais elle peut aussi dégénérer en une forme instable menant à une rupture de la plaque et, en quelques jours, provoquer des accidents cardiovasculaires ou cérébraux (ACV), mortels 30 ou morbides. En effet, la rupture d'une plaque met son contenu en contact avec le sang véhiculé par l'artère, ce qui peut entraîner la formation d'un thrombus. Celui-ci perturbe la circulation sanguine dans l'artère atteinte. Il peut aussi se détacher et être transporté par le flux 35 sanguin, et, dans les cas les plus sévères, bloquer complètement la lumière de l'artère, stopper l'irrigation de la région post lésion et provoquer son ischémie. 15 La caractérisation des propriétés mécaniques des tissus (c'est à dire de leur élasticité) présente un intérêt fondamental en diagnostic médical, notamment pour l'estimation d'un risque de rupture d'une plaque d'athérome. 2. Mesure de l'épaisseur de la chape fibreuse pour diagnostiquer les plaques d'athérosclérose à risque Pour permettre l'estimation d'un risque de rupture d'une plaque d'athérome, plusieurs techniques d'imagerie intravasculaire ont été développées. Il s'agit notamment de : - l'échographie intravasculaire (IVUS), de - la tomographie par cohérence optique (OCT), et de - l'imagerie par résonance magnétique (IV-I RM). Ces différentes techniques d'imagerie permettent d'obtenir une image de la plaque d'athérome, et de mesurer l'épaisseur de la chape fibreuse de celle-ci.
Toutefois, la seule connaissance de l'épaisseur de la chape fibreuse n'est pas un indicateur suffisant de stabilité de la plaque pour estimer un risque de rupture d'une plaque d'athérome. 3. Estimation de l'amplitude de la contrainte de la chape fibreuse pour détecter un risque de rupture Des travaux antérieurs ont permis d'établir que l'amplitude de la contrainte maximale au niveau de la chape fibreuse (ou « PCS », sigle de l'expression anglo-saxonne « Peak Cap Stress ») est un bon indicateur biomécanique de la vulnérabilité d'une plaque d'athérome. Toutefois, quantifier l'amplitude de la contrainte maximale au niveau de la chape fibreuse (PCS) in vivo reste un défi puisqu'une telle contrainte mécanique dans la chape fibreuse dépend non seulement de la morphologie de la plaque d'athérome, mais aussi des propriétés mécaniques des composants de la plaque, et notamment des inclusions contenues dans celle-ci. Bien que plusieurs méthodes aient été développées pour extraire les distributions des déformations dans une plaque d'athérome, la complexité des géométries des plaques d'athérosclérose ne permet pas d'en déduire directement les propriétés mécaniques de la plaque à partir de la connaissance des déformations.
Depuis les vingt dernières années, une nouvelle méthode d'imagerie médicale a été développée : il s'agit de l'élastographie ultrasonore. Basée sur les mêmes principes que la palpation, l'élastographie étudie localement le comportement élastique d'un milieu sous l'action d'une contrainte. Cette étude repose sur l'analyse de signaux ultrasonores radiofréquences acquis avant et après application d'une contrainte, ou acquis pour différents niveaux de contrainte. 4. Art antérieur connu dans le domaine de l'élastographie 4.1 Elastographie L'élastographie ultrasonore consiste à estimer l'élasticité d'un tissu mou à partir de techniques de traitements de séquences d'images échographiques rétrodiffusées par celui-ci alors qu'il se déforme. Le traitement de ces images permet de déterminer la répartition de la rigidité/élasticité de ce tissu. 4.2 Problème direct et problème inverse en élastographie En calcul de structure classique, le problème direct consiste à prédire les contraintes au sein d'un tissu, connaissant la distribution de son élasticité. Plus précisément, 20 connaissant : - la structure du tissu (i.e. sa géométrie en tout point), - la rigidité/l'élasticité des différents matériaux constituant le tissu, et - le chargement du tissu (i.e. le champ des contraintes externes exercé sur le tissu), il est possible de déterminer un champ de déplacement en tout point du tissu et ainsi de 25 remonter à la distribution spatiale des contraintes connaissant les lois de comportement élastique des milieux. À l'opposé, le problème inverse consiste à déterminer l'élasticité d'un tissu, étant donné les déplacements issus des signaux ultrasonores qu'il émet à différents états de 30 compression. Plus précisément, connaissant : - la structure du tissu (obtenue en utilisant des méthodes de traitement d'images échographiques), - le chargement évolutif du tissu (obtenu par mesure des contraintes externes exercées sur le tissu), et 35 - le champ de déplacement ou des déformations en tout point du tissu (obtenu par comparaison en tout point d'images ultrasonore représentatives des signaux ultrasonores acquis pour différents niveaux de chargement), 15 il est possible de déterminer la rigidité/l'élasticité des différents matériaux constituant le tissu. 4.3 Méthode connues pour déterminer la rigidité/l'élasticité d'un tissu Plusieurs procédés basés sur la technique des éléments finis ont été développées pour estimer une carte d'élasticité vasculaire à partir de l'estimation du champ de déformation à l'intérieur d'une lésion athéroscléreuse obtenue par imagerie intravasculaire.
La technique des éléments finis en mécanique consiste à réaliser un maillage d'une image d'un corps à étudier. Le maillage permet une discrétisation spatiale du corps. Ce maillage est composé de cellules élémentaires appelées aussi "mailles" (triangulaires, polygonales, carrées, etc.) adjacentes constituées chacune de noeuds et de bords délimitant une portion du corps.
Zhu et col. Zhu et col. (2003) ont développé un procédé direct de reconstruction basé sur la technique des éléments finis.
Ce procédé utilise des hypothèses sur les propriétés mécaniques de sorte à obtenir une relation directe entre le champ de déformation et le module de Young. Notamment, il est supposé dans ce procédé que les propriétés mécaniques de la plaque sont constantes pour chaque cellule élémentaire.
Malgré sa robustesse, ce procédé présente un temps de calcul très long du fait du nombre important de cellules élémentaires lors de l'examen d'une plaque athéroscléreuses très hétérogène. Oberai et col. (2003) Pour surmonter cette limitation, Oberai et col. (2003) ont proposé un procédé de reconstruction d'une carte d'élasticité (pour tous les cas de problèmes inverses en mécanique des structures) dans lequel non seulement la force et les déplacements sont considérés comme des variables nodales des cellules, mais également le module d'Young et le coefficient de Poisson. Ladite approche a pour nom "Propriétés Mécaniques Nodales" et pour acronyme « PMN ». Ainsi, pour une cellule élémentaire donnée, Oberai et col. proposent de considérer les variables d'élasticité aux noeuds de la cellule, et de déduire les élasticités à l'intérieur de la cellule élémentaire par interpolation en utilisant une fonction d'interpolation polynomiale. Un inconvénient de ce procédé est qu'il ne permet pas de prendre en compte les discontinuités de matériau dans la plaque d'athérome. En effet, les fonctions d'interpolation polynomiales nécessitent une continuité de la variable interpolée. Ainsi, les fonctions de forme spatiale utilisées ne permettent pas de caractériser une plaque d'athérome hétérogène présentant des discontinuités élevées de matériaux entre les différentes inclusions de la plaque d'athérome.
LeFloc'h et col. (2009) LeFloc'h et col. (2009) ont développé un procédé itératif de reconstruction connu sous le nom de « IMOD » (acronyme de l'expression anglo-saxonne « Imaging MODulography »).
Les procédés itératifs de reconstruction consistent généralement à : - assigner des valeurs d'élasticité aux différents composants/inclusions d'une plaque d'athérome, à - estimer un champ de déformation à partir des valeurs d'élasticité assignées, à - comparer le champ de déformation calculé au champ de déformation mesuré par imagerie intravasculaire, et à - réitérer les étapes précédentes jusqu'à ce que la différence entre le champ de déformation calculé et le champ de déformation mesuré soit inférieure à un seuil. La méthode IMOD (proposée par LeFloc'h et col.) comprend une étape de pré- conditionnement utilisant des critères de détection d'hétérogénéités de sorte à identifier automatiquement les contours de tous les composants de la plaque d'athérome. L'utilisation de tels critères est effectuée en employant un modèle paramétrique piloté par une segmentation par ligne de partage des eaux.
Malgré l'efficacité et la robustesse de l'approche IMOD, cette méthode ne permet pas de véritable reconstruction de l'élasticité en temps réel. En effet, cette méthode - comme toutes les méthodes itératives de reconstruction nécessite un temps de calcul très important (de l'ordre de plusieurs minutes) du fait du nombre important de réitération des étapes de calcul (assignation, estimation, comparaison) à mettre en oeuvre pour obtenir une carte d'élasticité vasculaire représentant la répartition locale du module de Young de la plaques d'athérosclérose. 5. But de l'invention Un but de la présente invention est de proposer un procédé de reconstruction permettant de pallier au moins l'un des inconvénients précités. Notamment, un but de la présente invention est de proposer un procédé de reconstruction en temps réel d'une image d'élasticité d'un corps - tel qu'une plaque d'athérome - à partir de signaux ultrasonores radiofréquences acquis avant et après application d'un chargement, ou acquis pour différents niveaux de chargement. PRESENTATION DE L'INVENTION A cet effet, l'invention propose un procédé de traitement d'image permettant la formation 15 d'une image d'élasticité d'un corps en fonction du ou des matériaux constituant ledit corps, caractérisé en ce que le procédé comprend les étapes consistant à : - Recevoir une image de déformation illustrant un champ de déplacement ou de déformation des points du corps en fonction d'une différence de pression dans le corps, 20 - Définir un maillage de l'image de déformation, ladite étape de maillage mettant en oeuvre une méthode des éléments finis, le maillage étant composé d'une pluralité de cellules élémentaires (ou mailles) délimitant chacune un même matériau et comportant chacune au moins trois noeuds, chaque noeud appartenant à une ou plusieurs cellules adjacentes du maillage, 25 - Assigner à chaque noeud i du maillage un couple de variables nodales inconnues (correspondant aux coefficients de Lamé (X' 1.0 ou au module d'Young et au coefficient de Poisson (E' v,)), chaque couple étant représentatif de propriétés élastiques à déterminer, - Détecter au moins un noeud commun à au moins deux cellules élémentaires 30 adjacentes délimitant des matériaux différents, - Enrichir ledit et au moins un noeud commun détecté, ladite étape d'enrichissement consistant à remplacer le couple de variable inconnues assigné audit noeud commun détecté par au moins deux couples de variables inconnues d'enrichissement, chaque couple de variables inconnues d'enrichissement du 35 noeud commun détecté étant attribué à une cellule élémentaire respective parmi lesdites au moins deux cellules élémentaires adjacentes délimitant des matériaux différents, 10 - Calculer une matrice élémentaire pour chaque cellule élémentaire du maillage de sorte à obtenir une pluralité de matrices élémentaires, la matrice élémentaire d'une cellule élémentaire étant calculée en tenant compte des couples de variables inconnues assignés à ladite cellule élémentaire, et des couples de variables inconnues d'enrichissement attribués à ladite cellule élémentaire, - Assembler les matrices élémentaires de la pluralité de matrices élémentaires pour obtenir une matrice de structure, et - Calculer une image d'élasticité du corps à partir de l'image de déformation et de la matrice de structure.
Ainsi, la présente invention propose une nouvelle formulation directe par éléments finis pour la résolution du problème inverse en élasticité linéaire. Cette formulation s'appuie sur l'introduction de nouvelles fonctions d'interpolation qui permet de tenir compte des discontinuités de matière à l'intérieur du corps étudié. Ce nouveau procédé peut avantageusement être mis en oeuvre dans n'importe quel logiciel de calcul de structure par éléments finis en élasticité linéaire. Des aspects préférés mais non limitatifs du procédé selon l'invention sont les suivants : - Le procédé comprend en outre, préalablement à l'étape consistant à définir un maillage, une étape consistant à segmenter l'image de déformation pour obtenir une image de déformation segmentée composée d'une pluralité de régions représentatives de zones du corps supposées être dans des matériaux différents ; - l'étape consistant à définir un maillage est mise en oeuvre sur l'image de déformation segmentée, les dimensions et positions des cellules du maillage étant déterminées en fonction des positions et dimensions des régions de l'image de déformation segmentée ; - l'étape consistant à calculer une image d'élasticité comprend la résolution de l'équation matricielle suivante : {R}= ([Q']T [Q1)-1 [Q']T {F'} où : o [Q'] représente la matrice de structure réduite, et [Qïr sa transposée, o {R} représente la matrice du champ des élasticités du corps aux noeuds, o {F'} représente la matrice réduite du champ de forces appliquées aux noeuds ; - l'étape consistant à enrichir un noeud commun détecté comprend le remplacement du couple de variables inconnues assigné audit noeud commun détecté par n couples de variables inconnues d'enrichissement, n correspondant au nombre de matériaux différents à la frontière desquels est situé ledit noeud commun. L'invention concerne également un système de traitement d'image permettant la formation d'une image d'élasticité d'un corps en fonction du ou des matériaux constituant ledit corps, caractérisé en ce que le système comprend : - Un récepteur pour recevoir une image de déformation illustrant un champ de déplacement des points du corps en fonction d'une différence de pression dans le corps, - Un processeur configuré pour : o Définir un maillage de l'image de déformation en mettant en oeuvre une méthode des éléments finis, le maillage étant composé d'une pluralité de cellules élémentaires délimitant chacune un même matériau et comportant chacune au moins trois noeuds, chaque noeud appartenant à une ou plusieurs cellules adjacentes du maillage, o Assigner à chaque noeud i du maillage un couple de variables nodales inconnues (X' III chaque couple étant représentatif de propriétés élastiques à déterminer, o Détecter au moins un noeud commun à au moins deux cellules élémentaires adjacentes délimitant des matériaux différents, o Enrichir ledit et au moins un noeud commun détecté, l'enrichissement consistant à remplacer le couple de variable inconnues assignées audit noeud commun détecté par au moins deux couples de variables inconnues d'enrichissement, chaque couple de variables inconnues d'enrichissement du noeud commun détecté étant attribué à une cellule élémentaire respective parmi lesdites au moins deux cellules élémentaires adjacentes délimitant des matériaux différents, o Calculer une matrice élémentaire pour chaque cellule élémentaire du maillage de sorte à obtenir une pluralité de matrices élémentaires, la matrice élémentaire d'une cellule élémentaire étant calculée en tenant compte des couples de variables inconnues assignés à ladite cellule élémentaire, et des couples de variables inconnues attribués à ladite cellule élémentaire, o Assembler les matrices élémentaires de la pluralité de matrices élémentaires pour obtenir une matrice de structure, et o Calculer une image d'élasticité du corps à partir de l'image de déformation et de la matrice de structure.
Des aspects préférés mais non limitatifs du système selon l'invention sont les suivants : - le système comprend en outre un afficheur pour afficher l'image d'élasticité calculée ; - le processeur est configuré pour mettre en oeuvre les étapes du procédé décrit ci- dessus. L'invention concerne également un produit programme d'ordinateur comportant un code programme enregistré sur un support de données lisible par un ordinateur pour exécuter le procédé décrit ci-dessus lorsque le programme d'ordinateur est appliqué à un ordinateur pour être exécuté. PRESENTATION DES FIGURES D'autres caractéristiques et avantages de l'invention ressortiront encore de la description qui suit, laquelle est purement illustrative et non limitative et doit être lue en regard des dessins annexés, sur lesquels : - la figure 1 illustre une plaque d'athérome, - la figure 2 illustre les étapes du procédé selon l'invention, - les figures 3A et 3B illustrent des cellules élémentaires enrichies (figure 3B) et non enrichies (figure 3A). DESCRIPTION DE MODES DE REALISATION DE L'INVENTION 1. Introduction Comme indiqué dans la partie présentation de l'art antérieur, Oberai et col. ont proposé un procédé de reconstruction par éléments finis pour résoudre un problème inverse en mécanique des structures.
Ce procédé - qui permet de reconstruire une image représentant les différentes élasticités d'un corps - est original en ce non seulement la force et les déplacements sont considérées comme des variables nodales, mais aussi le module d'Young et le coefficient de Poisson (où ce qui revient au même les deux coefficients de Lamé).
Toutefois, un inconvénient de ce procédé concerne l'utilisation de fonctions de forme spatiale continues. Ceci ne permet pas de caractériser les corps présentant des discontinuités élevés des matières.
Par ailleurs, ce procédé de reconstruction étant itératif, il ne permet pas de reconstruire une image représentant les différentes élasticités d'un corps en temps réel. La présente invention a donc été développée pour permettre une reconstruction en temps réel d'une image représentant les différentes élasticités d'un corps composé d'au moins deux zones présentant une discontinuité de matière. 2. Domaine d'application du procédé selon l'invention On va maintenant décrire plus en détail un exemple de procédé selon l'invention. Ce procédé est notamment applicable à la détermination d'une carte d'élasticité d'un corps incluant au moins une inclusion et au moins une cavité, le ou les différents matériaux constituant le corps étant sensiblement incompressibles.
Dans la suite, la présente invention sera décrite en référence à la fourniture - en temps réel - d'une information relative à l'élasticité des différentes inclusions contenues dans un corps. Plus précisément, la présente invention sera décrite en référence à un procédé de génération d'une image d'élasticité de module d'Young d'une structure deformable hétérogène entourant une cavité. Notamment, le procédé sera décrit en référence à l'étude d'un vaisseau avec plaque d'athérome. Toutefois, il est bien évident pour l'homme du métier que le procédé selon l'invention peut être appliqué : - à d'autres types de corps incluant ou non une cavité, tels que le coeur, le pancréas, le foie, la prostate, le sein, ... ou - a tout type de structure industrielle contenant ou non une cavité et ou la déformation doit être quantifiée. afin de reconstruire simultanément les cartographies du module d'Young E et du coefficient de Poisson y (ou ce qui revient au même, aux cartographies des deux coefficients de Lamé X et p) 3. Principe du procédé selon l'invention Le nouveau procédé décrit ci-après est basé sur une méthode d'analyse aux éléments finis (MEF) permettant d'obtenir un maillage du corps étudié. Le maillage est composé d'une pluralité de cellules, les contours de chaque cellule délimitant un matériau avec ses propres propriétés mécaniques.
Ce procédé a été développé en utilisant l'approche des propriétés mécaniques nodales pour le matériau (PMN). Pour modéliser les discontinuités entre les cellules du maillage, la méthode aux éléments finis étendue (ou « XFEM » pour « Extended Finite Element Method ») a été adaptée à l'approche PMN (PMN-XFEM). La méthode des XFEM est aujourd'hui la plus fiable pour traiter la fissuration en éléments finis et gérer les discontinuités de déplacements. La présence de singularités (fissures, perforations, etc.) dégrade fortement la convergence de la méthode aux éléments finis (MEF), et donc il ne suffit pas de raffiner fortement le maillage à proximité des singularités pour obtenir une bonne solution. Différentes approches ont été proposées pour pallier ce problème, la plupart reposant sur l'introduction au niveau des déplacements de fonctions de forme capables de représenter le comportement d'un corps fissuré au niveau de la fissure. La méthode des XFEM quant à elle propose un enrichissement des noeuds (i.e. dédoublement des noeuds) situés sur la fissure. L'adaptation de cette méthode à l'approche PMN est nouvelle et permet de gérer les discontinuités des matériaux au sein du corps étudié.
Le procédé de reconstruction directe de l'élasticité d'une plaque d'athérome selon l'invention a été appliqué à six lésions coronariennes de patients imagées in vivo par échographie intravasculaire (IVUS). 4. Généralités relatives au procédé selon l'invention Un exemple de plaque d'athérome est illustré à la figure 1. Cette plaque est une portion d'artère en coupe. Elle comprend une paroi intérieure 10 définissant et délimitant une cavité à l'intérieure de laquelle circule le sang, et une paroi extérieure 20 définissant la surface externe de l'artère. La plaque comprend également un espace 30 remplit d'une accumulation segmentaire de lipides (appelé aussi core nécrotique). La plaque d'athérome peut également comprendre d'autres inclusions telles qu'une région calcifiée 40. Ainsi, une plaque d'athérome comprend généralement une ou plusieurs inclusions (lipidique, calcique, etc.) dans des matériaux différents.
En référence à la figure 2, le procédé selon l'invention comprend les étapes suivantes : - Réception (étape 100) d'une image de déformation illustrant un champ de déplacement ou de déformation des points du corps en fonction d'une différence de pression dans le corps, - Optionnellement, segmentation (étape 200) de l'image de déformation, - Définition d'un maillage (étape 300) de l'image de déformation, - Attribution à chaque noeud i du maillage d'un couple de variables inconnues (X' 1.0 (où (E' v,)) représentatif de propriétés élastiques du corps, - Détection d'un (ou plusieurs) noeud(s) localisé(s) à la frontière entre deux (ou plusieurs) matériaux différents du corps, - Enrichissement (étape 400) du (ou des) noeud(s) détecté(s), - Calcul d'une matrice élémentaire (étape 500) pour chaque cellule élémentaire du maillage, - Assemblage (étape 600) des matrices élémentaires pour obtenir une matrice de structure globale, - Calcul (étape 700) d'une image d'élasticité du corps à partir de la matrice de structure globale. 4.1. Réception d'une image de cartographie de déformation/déplacement des points constituant le corps à étudier Comme indiqué ci-dessus, le procédé comprend une étape de réception 100 d'une image de déformation - dite « élastogramme de déformation» - illustrant un champ de déformation des points du corps en fonction d'une différence de pression dans le corps. L'élastogramme de déformation représente les déformations internes qui résultent de la compression du corps analysé - ici un tissu vasculaire sanguin - en fonction de la pression sanguine. L'élastogramme de déformation peut être obtenu par toute méthode connue de l'homme du métier. Par exemple, l'élastogramme de déformation peut être obtenu de la façon suivante. Une sonde ultrasonore est introduite dans une artère d'un patient au niveau d'une portion d'artère que l'on souhaite analyser. Une séquence d'images ultrasonore de la portion de l'artère est acquise pendant que le tissu artériel se comprime/dilate sous l'effet de la pulsation cardiaque. L'élastogramme de déformation est obtenu en étudiant la cinétique de la séquence temporelle d'images ultrasonore. Notamment, l'étude cinétique peut comprendre les étapes consistant à : - comparer deux images différentes de la portion d'artère acquises à deux instants distincts d'un cycle cardiaque, lesdites images différentes étant acquises pour deux chargements différents exercés sur la portion d'artère, - déterminer une carte des déplacements en comparant la position de chacun des points ou pixels entre les deux images différentes, - déterminer une différence de pression artérielle entre les deux images différentes en soustrayant les pressions de chargement artériel mesurées lors de l'acquisition desdites images différentes. On obtient ainsi une image de déformations - ou élastogramme de déformation- de la portion d'artère d'une part, et la différence de pression ayant permis d'obtenir cette image de déformation d'autre part.
Ces deux informations (i.e. cartographie des déplacements/déformations du corps analysé, et différence de pression associées aux cartographies des déplacements/déformations) sont utiles à la mise en oeuvre d'une étape ultérieure de détermination de l'élasticité de la paroi interne de la cavité du corps. 4.2. Segmentation L'étape de segmentation de l'image de déformation permet l'obtention d'une image de déformation segmentée comprenant une pluralité de régions délimitées par leur contour. 20 Chaque région de la pluralité de région est représentative d'une zone du corps à étudier, et plus précisément d'une zone du corps supposée être dans un matériau différent de la (ou des) zone(s) qui l'entoure(nt). Une première région délimite par exemple une inclusion lipidique de la plaque d'athérome, une deuxième région délimite une inclusion calcique, 25 etc. Plus précisément, l'étape de segmentation permet de détecter les contours des différentes inclusions de la plaque d'athérome, et notamment : - le contour intérieur de la portion d'artère, et représentatif de la paroi intérieure 10 définissant une cavité à l'intérieure de laquelle circule le sang, - le contour extérieur de la portion d'artère, et représentatif de la paroi extérieure 20 30 définissant la surface externe de l'artère, - les contours extérieurs des différentes inclusions 30, 40 de la plaque d'athérome. L'étape de segmentation permet de réduire de manière significative la quantité de données contenue dans l'image de déformation, et élimine les informations non 35 pertinentes pour la mise en oeuvre des étapes suivantes du procédé, tout en préservant les propriétés structurelles importantes de l'image de déformation.
L'étape de segmentation peut utiliser une technique de segmentation par ligne de partage des eaux, ou toute autre technique de segmentation connue de l'homme du métier. 4.3. Maillage de l'image de déformation L'étape de maillage met en oeuvre une méthode des éléments finis. Le maillage permet une discrétisation spatiale du corps. Le maillage consiste à quadriller l'image de déformation (éventuellement segmentée) par l'intermédiaire d'une pluralité de cellules élémentaires. Chaque cellule élémentaire délimite un milieu avec ses propres propriétés mécaniques. Chaque cellule élémentaire comporte des noeuds et des bords. Un noeud peut appartenir à plusieurs cellules élémentaires adjacentes. De même, un bord ne peut appartenir qu'à deux cellules élémentaires adjacentes. 4.4. Assignation de couples de variables inconnues aux noeuds des cellules Pour mémoire, le procédé permet de déterminer une image d'élasticité d'un corps en fonction du ou des matériaux constituant ledit corps. Les élasticités des différents matériaux (inclusions, etc.) constituant le corps sont donc inconnues en tout point de l'image de déformation.
Plutôt que de rechercher l'élasticité en toutes cellules élémentaires du corps, un couple de variables inconnues (connues sous le nom de coefficients de Lamé) (X' 1.0 représentatif des propriétés élastiques à déterminer est assigné à chaque noeud i des cellules élémentaires. Notons que l'on peut assigner à chaque noeud i soient les coefficients de Lamé (X' 1.0 soient le module d'Young et le coefficient de Poisson (E' v,).
Après avoir déterminé les valeurs associées à ces couples de variables inconnues au niveau des noeuds des différentes cellules élémentaires, il est possible de calculer les propriétés élastiques à l'intérieur des cellules élémentaires à l'aide d'une fonction d'interpolation polynomiale.
Plus précisément pour une cellule élémentaire donnée, lorsque les valeurs des couples de variables inconnues ont été calculées aux noeuds de la cellule, il est possible de calculer des couples de variables inconnues en tout point à l'intérieur de la cellule par interpolation à partir des couples de variables inconnues aux noeuds de la cellule. Le fait d'assigner les couples de variables inconnues aux noeuds de chaque cellule élémentaire autorise la variation des propriétés mécaniques au sein de la cellule et va permettre de reconstruire plus finement les cartographies d'élasticité. 4.5. Détection et enrichissement Une fois les couples de variables inconnues assignés aux noeuds des cellules élémentaires, le procédé comprend une étape de détection des noeuds situés à la frontière entre deux (ou plusieurs) régions représentatives de zones en matériaux différents du corps.
Un enrichissement des noeuds détectés est ensuite effectué. Pour chaque noeud détecté, l'enrichissement consiste à remplacer le couple de variables nodales inconnues attribué audit noeud commun détecté par n couples d'enrichissement, « n » correspondant au nombre de régions représentatives de zones de matériaux différents.
Par exemple, si un noeud détecté est situé à la frontière entre de deux régions, alors les deux variables inconnues assignées au noeud (un couple) sont remplacées par quatre (2 x 2) variables inconnues (deux couples). Si le noeud détecté est situé à la frontière de trois régions, alors les deux variables inconnues assignées au noeud (un couple) sont remplacées par six (2 x 3) variables inconnues (trois couples). Si le noeud détecté est situé à la frontière de quatre régions, alors les deux variables inconnues assignées au noeud (un couple) sont remplacées par huit (2 x 4) variables inconnues (quatre couples), etc. Les cas d'enrichissement lorsque le noeud est situé à la frontière de deux régions et trois régions sont illustrés sur les Figures 3A et 3B.
Chaque couple de variables inconnues d'enrichissement du noeud commun détecté est associée à une cellule élémentaire respective parmi les cellules élémentaires adjacentes délimitant des matériaux différents. Par exemple, dans le cas d'un noeud détecté N2 situé à la frontière entre deux régions représentatives de zones de matériaux différents M1, M2, alors le premier couple de variables d'enrichissement est attribué au noeud N10 de la cellule élémentaire délimitant une portion de la première région M1, et le deuxième couple de variables d'enrichissement est attribué au noeud N11 de la cellule élémentaire délimitant une portion de la deuxième région M2. Dans le cas d'un noeud détecté N5 situé à la frontière entre trois régions représentatives de zones de matériaux différents M1, M2, M3, alors un premier couple de variables d'enrichissement est attribué au noeud N14 de la cellule élémentaire délimitant une portion de la première région M1, un deuxième couple de variables d'enrichissement est attribué au noeud N15 de la cellule élémentaire délimitant une portion de la deuxième région M2, et un troisième couple de variables d'enrichissement est attribué au noeud commun N16 des cellules élémentaires délimitant une portion de la troisième région M3. L'étape d'enrichissement (basée sur un principe similaire à la méthode XMEF mais ramenée aux propriétés mécaniques nodales des cellules élémentaires) permet de gérer les discontinuités entre les matériaux constituant le corps étudié (dans le cas d'espèce, les inclusions de la plaque d'athérome). 4.6. Calcul des matrices élémentaires des cellules et assemblage Une fois les noeuds communs à plusieurs régions enrichis, le procédé comprend une étape de calcul d'une matrice élémentaire pour chaque cellule élémentaire. Cette matrice élémentaire est calculée en tenant compte des coordonnées des noeuds, du champ de déplacements (ou des déformations) aux noeuds, et des fonctions de forme utilisées pour les constantes nodales de matériaux et de déplacements.
Deux exemples de calculs de matrices élémentaires sont donnés en Annexe pour une cellule élémentaire triangulaire à 3 noeuds et une cellule élémentaire rectangulaire à 4 noeuds.
Lorsque les matrices élémentaires ont été calculées pour toutes les cellules élémentaires, celles-ci sont assemblées selon une méthode classique d'assemblage utilisée dans la méthode aux éléments finis connue de l'homme du métier. On obtient ainsi une matrice de structure globale.
Cette matrice globale et le champ de forces appliquées aux noeuds sont utilisés pour calculer les élasticités aux noeuds des différentes cellules du maillage.
On obtient ainsi une carte des élasticités des différents matériaux constituant la plaque d'athérome. Cette carte d'élasticité est ensuite affichée (étape 800) sur des moyens d'affichage tel qu'un écran d'ordinateur pour permettre à un utilisateur d'établir un diagnostic, et plus précisément d'identifier les plaques d'athérome susceptibles de se rompre. Avantageusement, le procédé décrit ci-dessus peut être implémenté sous forme de code programme pour permettre sa mise en oeuvre sur un ordinateur comprenant un processeur (configuré pour mettre en oeuvre les étapes du procédé décrit ci-dessus). Ce procédé permet à partir de signaux ultrasonore acquis par l'utilisateur par exemple par IVUS, d'afficher une carte d'élasticité d'une plaque d'athérome en temps réel sur un afficheur de l'ordinateur. On va maintenant décrire plus en détail des aspects théoriques associés au procédé selon l'invention. 5. Théorie relative à l'invention Nouvelle formulation de calcul des structures en éléments finis -prédictif : Afin d'introduire la notion d'enrichissement de matériau discontinu, on considère dans la suite le cas simple d'une plaque composée de trois matériaux distincts et maillée en 25 utilisant un maillage composé de quatre cellules élémentaires comportant chacune quatre noeuds (l'illustration de ce cas est représentée sur les figures 3A et 3B). La représentation par éléments finis de Galerkin des équations de l'élasticité appliquées à cette structure est donnée sous forme matricielle par : 30 [K(X, p.)] {U} = {F} (1) où : - {F} représente les forces nodales, 35 - {U} représente les vecteurs de déplacement, - X et p sont les coefficients de Lamé associés aux propriétés des matériaux des cellules, et - [K] est la matrice de rigidité symétrique globale. 20 Les approximations par éléments finis des déplacements associées à un maillage régulier tel qu'illustré sur les figures 3A et 3B sont données par: N9 N9 "()-C>)= 49.()-c>) et v()-e) = 1v, (p, (x) (2a, b) i=N1 i=N1 où : - u et y sont les deux composantes du vecteur de déplacement, - u et 111 sont les déplacements nodaux au noeud i, - (pi est la fonction de forme pour les déplacements associée au noeud i et - z est le vecteur de position.
Les principaux enjeux de cette technique sont la sélection des noeuds appropriés à enrichir, et la forme des fonctions de forme d'enrichissement associés. La méthode d'enrichissement nodal permet de satisfaire la contrainte de discontinuité des matériaux constituant la plaque d'athérome en introduisant des noeuds supplémentaires qui augmentent les degrés de liberté des noeuds situés à la frontière entre deux (ou plusieurs) matériaux différents de la plaque d'athérome. Il a été considéré dans le procédé selon l'invention qu'un noeud doit être enrichi s'il est situé à la frontière de discontinuités matérielles.
Les figures 3A et 3B illustrent l'application de cette règle sur l'exemple simple de la plaque. Les noeuds N2, N4, N5 et N6 ont été enrichis. En conséquence, l'approximation par éléments finis des deux constantes matérielles (c'est à dire les deux coefficients de Lamé X et p des cellules) associées au maillage dans les figures 3A et 3B est : 4)-e>) =1A, Oi()-c>)+IA, 0,()c) (3a) lei jeJ et p()-c>) = 1mi 01()-c>) /I, ()-c>) (3b) lei jeJ avec 0(x) = 02 (>) isi )- 4k ) si )-é G el (3c-e) é *i( = E el 10.2()) )cEe2 0: °4( (3f-h) )) 2())=10 sinon 01*> 1°5 ( )-é 0 sinon 0 sinon ) si e2 04 ( S )é e3 0 sinon 44,[*3()) = 01*.4())=°e si )-e> E el 0 sinon 0 sinon 4,6 ()_c) 105 ()-é) si )-é G e3, e4 OA si )-CEe2 le) si YeGe4 (3i-k) 0 sinon 07.7()) = 0 sinon 411*8())= 0 sinon OÙ : - 0, est la fonction de forme pour les propriétés mécaniques associée au noeud i et - e, (i = 1 à 4) est le nombre cellules, - I est l'ensemble de tous les noeuds non enrichis (i.e. I = {noeuds N1, N3, N7, N8 et N9}), - J est l'ensemble de tous les noeuds considérés à l'emplacement des noeuds enrichis (i.e. J = {noeuds N10 à N18}) et - (k E I,J) sont les composants des vecteurs matières nodaux {X} et {p} respectivement. Des fonctions de forme distinctes 'A et (pi ont été utilisées pour les propriétés des matériaux et les champs de déplacement. Par exemple une même géométrie de cellule élémentaire peut être considérée à 3 noeuds pour les propriétés mécaniques et à 6 noeuds pour les déplacements. Toutefois, il est bien évident pour l'homme du métier que les fonctions de forme peuvent-être prises identiques pour les propriétés des matériaux et les champs de déplacement.
Nouvelle formulation de calcul des structures en éléments finis inverse : Les équations de la méthode par éléments finis prédictive PMN-XFEM ont été réécrites en tenant compte des fonctions d'interpolation liées aux déplacements nodaux (i.e. déplacement au noeud) et des fonctions d'interpolation liées aux propriétés mécaniques nodales des matériaux. Puis ces mêmes équations de la méthode par éléments finis prédictive PMN-XFEM ont été reformulées de façon à isoler les inconnues du problème d'élastographie inverse qui sont les propriétés mécaniques nodales {R}. Ainsi, une matrice de déplacement [Q] a été extraite et utilisée pour déterminer les propriétés mécaniques aux noeuds {R}. On obtient ainsi le jeu suivant d'équations d'équilibre linéaires pour la nouvelle méthodologie de calcul des structures en éléments finis inverse : [Q(uk , vk)] {R} = {F} avec {R} = (4) L'expression de la matrice de déplacement élémentaire [Qe] pour une maille triangulaire à 3 noeuds et pour une maille rectangulaire à 4 noeuds sont données lorsque l'on fait l'hypothèse des déformations planes (c. à d. lorsque les déplacements de tous les points de la plaque reste sur le plan de coupe) (cf Annexe ci-dessous). La matrice de déplacement élémentaire [Qe] est carrée mais pas symétrique. Pour obtenir la matrice de déplacement globale [Q], les matrices de déplacement des différentes cellules élémentaires [Qe] sont assemblées de manière connue de l'homme du métier pour en déduire la matrice globale de rigidité [K].
Résolution numérique: La matrice globale de déplacement [Q] est rectangulaire tandis que la matrice de rigidité [K] est carrée.
Plus précisément, les deux matrices [Q] et [K] ont un même nombre de lignes, mais leur nombre de colonnes est différent. Le nombre de colonnes de la matrice globale de déplacement [Q] est égal à la somme de son nombre de lignes et de deux fois le nombre de noeuds supplémentaires au niveau des noeuds enrichis. Un système similaire peut être obtenu en utilisant des déformations nodales au lieu de variables de déplacement nodales.
Considérant le problème inverse, la matrice [Q] est connue puisque tous les déplacements nodaux uk et vk sont obtenus à partir de l'image de déformation illustrant le champ de déplacement. Le problème inverse est mal conditionné parce que les forces nodales ne sont pas connues sur les noeuds où les déplacements ont été imposés. Le vecteur de force global {F} a par conséquent, certaines composantes inconnues. En supprimant les équations impliquant ces composantes de force inconnues du système par éléments finis global, on obtient le système réduit d'équations linéaires suivant : [Q'] {R} = {F'} (5) où : - {F'} représente le vecteur force nodal réduit, et - [Q'] représente la matrice de déplacement réduite. Pour résoudre ce système linéaire rectangulaire, on utilise des algorithmes mathématiques connus de l'homme du métier. Par exemple, on effectue un produit matriciel des termes situés de part et d'autre de l'égalité par la transposée de la matrice de déplacement réduite ; on obtient ainsi un système linéaire carré : [Q']T {R} = [Q']T {F'} (6) On résout ensuite ce système pour déterminer les propriétés mécaniques aux noeuds {R} : {R}= ([Q']T [a])-1 [Q']T {F'} (7) Le procédé décrit ci-dessus permet donc de fournir une information temps réel exploitable à l'utilisateur lui permettant de prédire les risques de rupture d'une plaque d'athérome. 6. Annexe 6.1. Expression de la matrice élémentaire de déplacement [Qe] lorsque l'on considère une même cellule élémentaire triangulaire à 3 noeuds pour les propriétés mécaniques et les déplacements On fournit ici l'expression de la matrice de déplacement élémentaire [Qe] (qui est une matrice carrée et non symétrique) pour une cellule élémentaire triangulaire à 3 noeuds (1, 2 et 3) lorsque l'on fait l'hypothèse des déformations planes (c. à d. lorsque les déplacements de tous les points de la plaque restent sur le plan de la structure). Soient respectivement (xi, yi) et (ui, vi) les coordonnées et les déplacements du noeud i (pour i= 1 à 3). Les 36 composantes Que (i: indice de la ligne; j: indice de la colonne) de la matrice carrée (6 lignes x 6 colonnes) de déplacement non symétrique [Qe] sont données ci-dessous en fonction de xij, yij (où xij = xj-xi , yij = yj-yi avec i, j= 1 à 3) et où A est la surface d'une cellule élémentaire triangulaire (A= (x21y31- y21x31)/2): Ligne 1 Qu Q12 22 1)23 - 0,2 .1/31 -1123 + - 3- 1l12-1)23 r1 32 lu) - 3 Y1-2.123 x2 i) c:3-T:32 Y23 - U- 2 Y31 Y2:3 - 03- 1121122 + U. r1 -'221/23 ' - - Y'S2 .Y1') Y23 U.2.!)31 Y23 2.1112/)2:3 ± -1'32 4131 (2Y12 .Y23 + .r2i) + r3-1".32.Y12 Q13 Q14 Q15 Q16 )1 ..Y23 21 - 1:12:3 Ligne 2 Q21 1 T22,2:3 + ' `>0 ') - 3') .Y31 + = U -i. ..,--/- =,.-2 ..T.!/2'3 23 - '2 - + . 32 .r1:?.. 1)23 .1.1?..1) - 0:3 -121 .Y-').?, - ..r.?..2.7-21 + 912 Y23 ) ti1 T22.!)23 - - - .02- :3') .Y31 1 + 03.-912.. + c:3-1"21 . = tt1.....2.!./23 7r -2 ..Y 2 - Y 3 + ' 3-2 r 1 -:. + 923 .!». 0 :3 .-1"21 .Y23 - til. ..T2 ..Y2:3 - 11.32 ., --T3..) . Y:31 + r2» ui...1.2 . Y23 - '') + V2-1, 3-2 ..-r1... -I- Y23 .1:hl 3-21922 - 32-1'21 + 1112.922) 11.c.732.b(12 + Y12. Y23 32 ..Y12 2 3 Ligne 3 V3 .17.) + î Y31 + U2.91 + V2 113 .Y31 - 3. Y12 Y31 Q2201 .1,423 .Y31 3 32 Y13 ) "13 + u2-( 2 + 2 ..T13 931 ( 1112.931 .x21) + = rl - x3') Y31 + U2.1J1 + - x13 .Y31 - 3 -Y12 g31 - )1 = - xi-3) 9231/31 2"13 ) -.1"13-Y2 + + r2 93. ( Y12 Y:31 ..T21) + r:3 Y12 . 2 ll 1 ...Y:31 4-U2 Y31 a 13 ..Y31 - :3-912.931 (2923 .Y31 - + - 1"13-1/722 + 02 1 + Ét3.(2/112. + + C:3 -T13 Yr? Vy: 21 Ligne 4 Y23 Q42 i-.TU. Y:31 X . = U T13.j23 - 4 U1 :-/-2 .Y31 3' V. -T1:3 Y23 - Ét - + V2. - té.? Y:31 t12 Y12 + ) .) - thl ..T.)1.!/31 - X3.) - té.) ...T13 .Y31 - )2-/11 1 Y31 912 .Y31 ) -X3 té.? .X13 . Y31 - - Y12 4- + L1'2. /)j v3.(2.x13/21 ± 912- 2/131) Y12 4- Ligne 5 'té 1.1/12 - Y2:3 + Y12 4- él 2-/Jif112 12..Ti./Ji2 - ) - Y12 U.1(2 t123 Y1"' -T 32.a'21 ) :1 - X.,1 - u2.':Y.:31.y1 2 4- .T 1.3 1) + V2 .X21.93/ + 11-3 yi2 ± .121 UI Y12 -V23 ') Y12 ti2 el.1112 y12 u { 292:3 Yl") -.) X21 ) u1.21.2;=, ee2 - 431.1/12 1.2.3'21.im + 37,:221) - 111 ./112.5J2i + .1.112 + .Y12 - + 1..9r2 + Y21 ../J12 ui-{e'"Lj - 132 "T 21 ) °X') -Y-2;=., 4- U2. L Y:31 - Y12 4- v2-x-21.9;n. + U2 ,L +x 912 Ligne 6 QI = 21 2 = -T21 )') ,T 21 -//31 ri:3 -21 ±Le:?..-X-)1 -Y1.2 U1 - U23.Y12) - .Y12 112 . - 931. Y12 ) Y12 - U1 .X21 .Y2:3 - 1 -X21 -X3') - 1.12 Y21.131 - M. 21 ek X') - /112 y - T.«.
21 Y )3. Y12 ) V2 - - ..931 .Y12 U3 .X21 .Y1', et1 .X-21 -Y23 - 1 - I''1 -X32 - 112 X21 .Y31 -- U1 ..TS2 - .X21 + Y)3-912) - -) - Y12 . -X21 931- Y12 U:3-X21 .Y12 + + .Y12 + 6.2. Expression de la matrice élémentaire de déplacement [Qe] lorsque l'on considère une même cellule élémentaire rectangulaire à 4 noeuds pour les propriétés mécaniques et les déplacements On fournit ici l'expression de la matrice de déplacement élémentaire [Cr] (qui est une matrice carrée et non symétrique) pour une cellule élémentaire rectangulaire à 4 noeuds (1, 2, 3 et 4 de coordonnées respectives (-a,-b), (a,-b), (a,b) et (-a,b)) lorsque l'on fait l'hypothèse des déformations planes (c. à d. lorsque les déplacements de tous les points de la plaque restent sur le plan de la structure). Soit (ui, vi) les déplacements du noeud i (pour i= 1 à 4). Les 64 composantes Que (i: indice de la ligne; j: indice de la colonne) de la matrice carrée (8 lignes x 8 colonnes) de déplacement non symétrique [Cr] sont données ci-dessous. Ligne 1 20 Ligne 2 Ligne 3 Ligne 4 Ligne 5 Ligne 6 Ligne 7 Ligne 8 REFERENCES BIBLIOGRAPHIQUES Doyley MM. Model-based elastography: a survey of approaches to the inverse elasticity problem. Phys Med Biol 2012;57:R35-73. Le Floc'h S, Ohayon J, Tracqui P, Finet G, Gharib AM, Maurice RL, Cloutier G, Pettigrew RI. Vulnerable atherosclerotic plaque elasticity reconstruction based on a segmentation-driven optimization procedure using strain measurements: theoretical framework. IEEE Trans Med Imaging 2009;28:1126-37. Moés N, Dolbow J, Belytschko T. A finite element methods for crack growth without remeshing. Int J Numer Meth Engng 1999;46:131-50 Oberai AA, Gokhale NH, Feijoo GR. Solution of inverse problems in elasticity imaging using the adjoint method. Inverse Problems 2003;19:297-313. Zhu Y, Hall TJ, Jiang J. A finite-element approach for Young's modulus reconstruction. IEEE Trans Med Imaging 2003;22:890-901.