Un Programme ANR 2008-2011: Optimised Computational Functional Imaging for Arteries

OCFIA (Optimised Computational Functional Imaging for Arteries) est une méthode de modélisation qui part de l’imagerie médicale et fournit des informations sur la biomécanique de la zone explorée.

http://www.ocfia.org/

une chaîne de traitement des images médicales a été conçue et optimisée afin de répondre à un cahier de charges incluant les éléments suivants :

•     Compatible avec l’imagerie cardiovasculaire IRM et Scanographie
•     Mouvements des structures au plus près de la réalité observée dans les images, calcul explicite afin d’éviter des erreurs sur la rhéologie des parois.
•     Temps de calcul < 2H par dossier patient
•     Possibilité de faire des études cliniques multicentriques afin d’évaluer de nouveaux facteurs prédictifs pour les pathologies cardiovasculaires (exemple : Anévrisme de l’aorte abdominale, Dissection aortique, Cardiopathies, ..)

 

Post-processing

Contents

1 Présentations des Logicels
 1.1 MRIcro et le format ANALYZE
 1.2 MATLAB 7.5
 1.3 AMIRA 4.1 et son interface 3D
2 Filtrage des images
 2.1 Filtre Gaussien
 2.2 Filtre de Flou Selectif
  2.2.1 Modèle théorique
  2.2.2 Algorithme
  2.2.3 Résultats
 2.3 Filtre par diffusion anisotrope
  2.3.1 Modèle théorique
  2.3.2 Algorithme
  2.3.3 Résultats
3 Extraction de la géométrie "native"
 3.1 Méthode de segmentation par Level Set
  3.1.1 Méthode Level Set et Algorithme
  3.1.2 Résultats 3D
 3.2 Méthode de segmentation par seuillage et proche voisin
  3.2.1 Exemple d’une dissection de l’aorte
4 Déplacement de la géométrie
 4.1 Principe de déformation
 4.2 Méthode de déformations non-linéaires
  4.2.1 Principe d’optimisation
  4.2.2 Algorithme
  4.2.3 Probabilité de trouver la "bonne" déformation
  4.2.4 Inverser un champ de déformation
  4.2.5 Exemple 3D avec un Cas test : Modèle "dalle Carrée"/ "dalle Ronde"
  4.2.6 Exemple 3D avec un modèle de fantôme d’anévrisme
5 Flux artériels, ajustement et correction
Conclusion

1 Présentations des Logicels

Les programmes utilisés pour développer le concept de chaîne décrit plus haut ont la caractéristique d’avoir des fonctionnalités bien diverses, mais tournent autour d’un ou plusieurs formats compatibles. Ils ont été choisis parmi ceux qui étaient à notre disposition et parmi ceux qui se présentaient avec un meilleur rapport entre convivialité des interfaces utilisateur et fonctionnalités. Il est bien accepté qu’un nouveau code puisse prendre la place d’un de ceux décrits ci-dessous, car l’approche modulaire de notre chaîne de traitement le permet.

1.1 MRIcro et le format ANALYZE

MRIcro (www.mricro.com) est un logiciel libre qui permet de visualiser facilement les images DICOM provenantes de l’imageur, au format ANALYZE-7 (Mayo Clinic, Rochester, USA., www.mayo.edu/bir/). Il est basé sur l’organisation des données dans deux fichiers, au lieu d’en avoir un par image. Ceci est pratique quand le volume de données peut atteindre 1000 images. Alors il est possible de transformer ce volume de données en deux fichiers de types différents : un fichier d’entête (*.hdr) et un fichier de données brutes (*.img) qui portent le même nom.

• Le fichier image *.img : Le fichier image est une suite ininterrompue de valeurs correspondant aux données des voxels dans un paquet d’images DICOM (selon les cas, entiers non signés, entiers signés, virgule flottante ou double précision). Chaque fichier *.img est généralement associé à un fichier d’entête contenant des informations sur l’organisation de ces données image.
• Le fichier d’entête *.hdr : ANALYZE adopte un fichier de 348 octets qui contient les informations nécessaires au traitement géométrique. Il est divisé en différents champs :
  • taille de l’image en voxels (x,y et z)
  • taille des voxels en mm (x,y et z)
  • codage des valeurs des voxels
  • coefficient d’échelle.
  • offset des valeurs des voxels dans *.img (en octets).
  • coordonnée du point d’origine, le centrage de la séquence par exemple.
  • description (un string).

Comme j’ai choisi de travailler de façon modulaire, j’utilise le format analyze pour passer d’une tâche de post-traitement à une autre sans avoir des soucis de format. L’autre avantage est qu’il me permet d’accéder à tout moment aux résultats intermédiaires. Un module slice Viewer permet de parcourir un ou plusieurs volumes par l’intermédiaire d’une glissière. Le volume peut être observé en différents points de vue (coronal, sagittal, axial, ou les trois).

1.2 MATLAB 7.5

MATLAB (www.mathworks.fr) est un langage de calcul scientifique de haut niveau. C’est un environnement interactif pour le développement d’algorithmes, la visualisation et l’analyse de données, ou encore pour le calcul numérique. En utilisant MATLAB, nous pouvons résoudre des problèmes de calcul scientifique et réaliser des prototypes logiciels presque aussi performants qu’avec les langages de programmation traditionnels, tels que C, C++ et Fortran.

MATLAB offre un certain nombre de fonctionnalités pour la documentation et le partage du travail. On peut aussi intégrer le code MATLAB avec d’autres langages et applications, et distribuer les algorithmes et applications MATLAB.

Ses principales fonctionnalités sont :

• Langage de haut niveau pour le calcul scientifique.
• Environnement de développement pour la gestion du code, des fichiers et des données.
• Outils interactifs pour l’exploration itérative, la conception et la résolution de problèmes.
• Fonctions mathématiques pour l’algèbre linéaire, les statistiques, l’analyse de Fourier, le filtrage, l’optimisation et l’intégration numérique.
• Fonctions graphiques 2-D et 3-D pour la visualisation des données.
• Outils pour la construction d’interfaces graphiques personnalisées.
• Fonctions pour l’intégration d’algorithmes développés en langage MATLAB, dans des applications et langages externes, tels que C/C++, Fortran, Java, COM et Microsoft Excel.

Avec nos algorithmes développés in house, nous avons intégré la boîte à outils SPM qui est distribuée gratuitement dans la communauté scientifique avec une licence libre GNU (General Public Licence http://www.fil.ion.ucl.ac.uk/spm). Cette intégration nous a servi de point de départ pour le calcul des déformations de la paroi vasculaire, sur la base des volumes dynamiques obtenus par IRM.

Développé initialement en 1994 par John Ashburner et Karl Friston, membres et collaborateurs du Welcome Department of Cognitive Neurology (Institute of neurology, University College of London), la boîte à outils SPM représente l’implémentation des concepts du traitement statistique de séries d’images réunissant toutes les fonctions nécessaires au traitement des différentes modalités d’imagerie d’activation cérébrale (imagerie par résonance magnétique fonctionnelle et tomographie par émission de positrons). Il propose toutes les fonctions nécessaires aux prétraitements, telles que la normalisation spatiale, le recalage et l’échantillonnage entre les images de modalités différentes.

Principales tâches développées sous l’environnement MATLAB :

• Filtrage des images.
• Extraction de la géométrie vasculaire par la méthode Level Set.
• Calcul des déformations de la paroi basé sur les images d’IRM dynamiques.
• Écriture des maillages mobiles et quelques opérations sur le maillage natif.
• Préparation et correction des CL pour les débits d’entrée / sortie.

Nous avons utilisé des fichiers au format ANALYZE pour les manipulations géométriques, les maillages ont été écrits au format *.cas (FLUENT/UNS) et les CL au format *.txt (Texte).

1.3 AMIRA 4.1 et son interface 3D

AMIRA est un système de modélisation et visualisation 3D, développé conjointement par Konrad-Suze-Zentrum et Visual Concepts GmbH. Il vient du Open Inventor Toolkit proposé par TGS Template Graphics Software, Inc. Le logiciel Amira 4 est un produit de Mercury Computer System mis à jour sur sa version 4 en 2006. Il s’agit d’un système multifonctionnel de visualisation et modélisation de données. Il permet de manipuler virtuellement tout type de données scientifiques à partir de domaines d’applications aussi variés que la médecine, la chimie, la physique ou l’ingénierie. Son principe général de fonctionnement repose sur la représentation de toute base de données dans un espace de gestion appelé "object pool", où l’on peut associer différentes opérations. Ce logiciel autorise une représentation numérique fidèle d’une image y compris en 3 dimensions, avec un outil 3D de navigation et de visualisation de cette image. Ce logiciel fonctionne sur la base de modules et de données-objets. Les modules permettent de visualiser les données-objets ou de réaliser des calculs sur ces données-objets. Une fenêtre représente en tridimensionnel l’objet étudié sur demande, ce qui permet de vérifier en permanence la correspondance des modifications apportées à la réalité anatomique 3-D. La correspondance entre la géométrie discrétisée et son volume d’imagerie associé peut être validée visuellement par l’affichage des deux modules dans le même espace virtuel. Un exemple de résultat de segmentation d’un AAA par la méthode level set en 3D est donné à la figure.1 où une coupe morphologique nous sert de référence pour valider l’extraction de la forme anantomique.

---

Figure 1: Corrélation entre la représentation surfacique d’un anévrisme de l’aorte abdominale (AAA) avec son volume d’imagerie morphologique. Sert à contrôler le processus d’extraction de la géométrie vasculaire native

---

Les structures en 3D peuvent être représentées par des maillages convenables aux simulations numériques (CFD), à savor les surfaces triangulaires (éléments finis) et les volumes tétraédriques (volumes finis). Amira propose des méthodes pour générer de tels maillages à partir des données de segmentation manuelles ou automatiques.

2 Filtrage des images

Les images obtenues sont en niveaux de gris et possèdent un rapport signal/bruit plus ou moins important suivant le type d’acquisition. En effet, des coupes fines d’IRM (<  3mm) ont un rapport signal/bruit plus faible que les coupes plus épaisses (> 5mm). En réalité c’est le volume du voxel d’acquisition qui détermine l’intensité de signal utile sur le bruit de fond. Il est important d’ôter ce bruit afin de faciliter l’opération de segmentation qui sera développée plus loin dans ce manuscrit. L’objectif des filtres est d’homogénéiser les niveaux de gris des voxels qui correspondent à la lumière vasculaire et d’accentuer ceux qui marquent la limite avec la paroi artérielle. Nous avons à disposition deux filtres qui ont pour objectif l’amélioration des contours, un filtre de diffusion anisotropique basé sur les équations de diffusion, et un filtre de flou sélectif basé sur le rapport signal/bruit. Mais nous avons aussi un filtre qui a un autre objectif : le filtre gaussien part d’un niveau de filtrage très élevé pour ajuster la déformation de la paroi par un processus de filtrage dégressif.

Par la suite, nous allons décrire ces filtres qui sont utilisés dans la chaîne de traitement.

2.1 Filtre Gaussien

Un filtrage gaussien est appliqué aux volumes d’images avant l’analyse statistique. L’algorithme qui estime les déformations de la paroi par rapport au volume natif a en effet besoin de ce pré-traitement pour deux raisons :

• Les données d’imagerie présentent des auto-corrélations spatiales. En effet, le signal acquis dans un voxel n’est pas indépendant de celui des voxels voisins. Rappelons que l’image réelle acquise est assimilable à la convolution d’une image idéale par une fonction de dispersion ponctuelle liée à la résolution spatiale intrinsèque de la méthode d’acquisition. La connaissance de ces auto-corrélations est nécessaire par exemple pour faire des corrections dans les images. Cependant, leurs caractéristiques sont difficiles à estimer et elles ne sont de toute façon jamais parfaitement connues par les techniques de restauration d’image (déconvolution). Ici, un filtrage passe-bas (gaussien) est imposé dont les caractéristiques de corrélation spatiale sont, elles, parfaitement connues. Il est utile dans la recherche d’une déformation optimale par filtrage dégressif des volumes natif et cible.
• La déformation spatiale n’est jamais parfaite. L’application d’un filtrage gaussien a aussi l’intérêt de réduire l’effet des imperfections locales de déformation spatiale.

---

(a) Gaussienne

(b) LMH = 20mm

(c) LMH = 10mm

(d) LMH = 5mm

(e) Image Brute

Figure 2: Exemples du filtrage gaussien sur une coupe provenante de la segmentation d’un volume d’un fantôme d’anévrisme. Plusieurs niveaux de filtrage sont appliqués pour arriver à optimiser la transformation entre deux volumes.

---

(1)

(2)

Le filtrage gaussien est réalisé par convolution discrète du volume natif et du volume cible par un noyau gaussien Gj dont l’amplitude à j unités du centre est définie par la relation suivante. Où LMH est la largeur en mm à mi-hauteur de la Gaussienne. Un exemple de filtrage dégressif est exposé à la figure 2(a).

2.2 Filtre de Flou Selectif

L’utilisation de ce filtre est inspirée des travaux de Soltanian-Zadech et al [1, 2]. Ce filtrage utilise les informations contenues dans les pixels voisins du point considéré. Il a été développé en 2D et le masque (ensemble des pixels pris en compte) est de taille M × M et est appelé G.

2.2.1 Modèle théorique

On obtient l’image filtrée après convolution de l’image originale par le masque G :

(3)

(4)

Avec i,j les éléments du masque de convolution (lignes, colonnes). Les différents coefficients de G sont calculés avec normalisation par l’équation :

(5)

Où g_kl est le poids attribué aux différents pixels voisins et avec 0 < g_ij < 1, 1 ≤ i ≤ M et 1 ≤ j ≤ M. G est donc calculé suivant l’environnement du point considéré comme dans le filtre de diffusion anisotropique décrit précédemment. Ce filtre tient compte du SD (bruit de l’appareil d’imagerie), la déviation standard ou l’écart type du signal. En posant S(X,Y ) l’intensité du pixel considéré et S(x,y) l’intensité des pixels voisins, les coefficients de pondération du filtre proviennent des équations suivantes :

(6)

avec,

(7)

On utilise un rapport signal/bruit connu de l’appareil d’imagerie (IRM Philips Intera 1.5T) qui est entre dix et vingt, suivant la résolution spatiale, le nombre d’acquisitions et autres paramètres de la séquence. Ce filtre permet de fixer un seuil statistique de telle sorte que seuls les pixels qui atteignent ce seuil seront filtrés. Quand la variation de l’intensité des pixels contenus dans le filtre est proche des valeurs de fluctuation du bruit, le poids attribué est maximal et égal à 1 mais si ces variations ne sont pas du niveau des fluctuations du bruit, et donc proviennent d’un détail important de l’image, le poids est limité pour conserver cette information locale.

2.2.2 Algorithme

On calcule la nouvelle valeur du point central en chargeant les valeurs des pixels voisins, leurs différents poids (calculés en fonction de l’écart type donné et des différences entre le point central et ces voisins) chacun dans une vecteur 1D.

% On remplit la matrice de gradient  
mat_gradient = abs(matrice - mat_pt_central);  
 
for num_coeff = 1 : nb_lignes * nb_colonnes  
    if mat_gradient(num_coeff) < 2 * SD  
        mat_poids(num_coeff) = 1;  
    else  
        mat_poids(num_coeff) = 0; % 1 / mat_gradient(num_coeff);  
    end  
end  
 
% On calcule la matrice resultats  
for num_coeff = 1 : nb_lignes * nb_colonnes  
    mat_resultats(num_coeff) = matrice(num_coeff) * mat_poids(num_coeff);  
end  
 
% Calcul de la somme des elements de la matrice contenant les resultats  
somme_mat_resultats = sum(mat_resultats);  
 
% Calcul de la somme des elements de la matrice contenant les poids  
somme_mat_poids = sum(mat_poids);  
 
% Détermination de la valeur du point central  
point_central = somme_mat_resultats / somme_mat_poids;  
point_central = uint16(point_central);

2.2.3 Résultats

Concernant la taille du masque, nous nous sommes basés sur les travaux de D. Gensanne [3] sur les tissus adipeux et nous avons implémenté un filtre 9x9. L’exemple présenté ici illustre l’influence du filtre sur une image de l’aorte thoracique traité par Stent-Graft.

---

(a) Image brute

(b) SD=14

(c) SD=20

(d) SD=50

Figure 3: Coupe Sagittale Oblique d’une séquence injectée. Opérations de filtrage avec plusieurs niveaux de SD pour le filtre de flou sélectif

---

Les paramètres de séquence de l’image étudiée sont assez bons car l’injection de produit de contraste lui donne un excellent contraste (Fig.3(a)), mais il faut savoir aussi que la finesse de ses voxels (0.88 × 0.88 × 1.6mm³) donne un rapport signal sur bruit qui est estimé à 14 [3] pour l’appareil d’imagerie en question. Nous pouvons observer les résultats du filtrage pour SD = 14 (Fig.3(b)), ce qui équivaut à une image théoriquement "non bruitée". Nous sommes allés plus loin dans le filtrage en forçant sur le paramètre SD, pour SD = 20 (Fig.3(c)) et même SD = 50 (Fig.3(d)). Les résultats sont visiblement intéressants pour le traitement numérique des images. Le seul inconvénient est que ce filtre est actuellement fait pour travailler avec des images 2D et un développement vers l’analyse volumique est assez difficile à implémenter. Nous l’utilisons alors en mode multicoupes 2D.

Nous allons voir par la suite un filtre qui utilise l’information des trois directions de l’espace pour filtrer un volume d’imagerie. La possibilité de combiner deux filtres nous a paru intéressante. La segmentation des images est très sensible aux problèmes d’irrégularités de signal entre coupes voisines, surtout quand elles sont épaisses (< 4mm) par rapport à la direction de la coupe (~ 1mm). Utiliser les deux filtres nous donne la possibilité de coupler les fonctionnalités des deux méthodes.

2.3 Filtre par diffusion anisotrope

C’est un filtre que nous avons trouvé approprié pour réduire le bruit et homogénéiser les zones de niveaux de gris voisins. Il est décrit par P. Perona et J. Malik pour le filtrage des images en 2D [4]. Ce filtrage s’inscrit dans la lignée des filtres EPSF (Edge Preserving Smoothing Filtering), c’est-à-dire que théoriquement il préserve les contours des objets tout en filtrant les zones homogènes.

2.3.1 Modèle théorique

L’équation de diffusion de la chaleur est isotrope (8a) mais ce filtre modifie cette condition par l’utilisation de l’équation de diffusion anisotrope (8b) pour que les contours ne soient pas traités de la même manière :

(8a)

(8b)

Ici h(f,t) est une fonction généralement du module de la dérivée de f (signal), maximale pour une dérivée nulle et décroissante lorsque le gradient croît ; on choisit généralement deux options :

• Option 1 : favorise les forts contrastes sur les faibles contrastes (9a)
• Option 2 : favorise les régions étendues sur les petites régions (9b)

(9a)

(9b)

Si κ est petit, les petits gradients arrivent à bloquer la conduction. Une grande valeur réduit l’influence des gradients d’intensité sur la conduction.

Ce filtre, initialement utilisé en 2D, a été développé au CHU de Rangueil pour traiter des volumes d’images représentant des morphologies en 3D. L’algorithme codé en langage MATLAB tient compte des 26 voisins du voxel traité. La longueur de son code source, est compensée par son efficacité et ses possiblités de parallélisation. En effet, l’algorithme itératif utilise des opérations matricielles simples et indépendantes pour calculer massivement le coefficient de diffusion (eq.9a) ou (eq.9b), pour ensuite l’appliquer massivement aux voxels filtrés suivant la formule (eq.8b).

2.3.2 Algorithme

---

Figure 4: Calcul du gradient dans la direction y-1

---

Pour l’implémentation discrète du filtre, le gradient de l’image est calculé pour chaque direction en regardant le pixel voisin, ici l’exemple pour les colonnes :

(10)

Cette écriture est simplifiée par la suite :

(11)

En posant,

(12)

Nous pouvons définir la contribution du voxel voisin à droite :

(13)

où c est le coefficient de diffusion et Δx = 1 voxel.

Ensuite, par analogie, nous pouvons faire la même chose pour les contributions des voxels à gauche, devant, derrière, en haut, en bas et ainsi tous les 26 voisins. Nous pouvons alors compléter l’équation avec toutes les contributions :

(14)

Le codage se fait par une méthode matricielle pour laquelle MATLAB répond de façon très efficace. Pour ce faire, au volume initial on ajoute un voxel dans toutes les directions de l’espace. Si par exemple nous avons un volume à filtrer contenant 512 × 512 × 50 voxels, après modification nous aurons un volume de 514 × 514 × 52 voxels. Un exemple en 2D est donné à la Fig.4. Le calcul du gradient I(x,y - Δy,z) - I(x,y,z) revient à soustraire la matrice initiale avec son homologue déplacée de dy. Le code pour le calcul des gradients se présente finalement de la façon suivante :

Le coefficient de diffusion étant calculé à partir des gradients précédents le codage est simple et rapide. Un bon filtrage nécessite plusieurs itérations. A chaque itération l’image sera modifiée, un coefficient d’atténuation λ doit être ajouté pour ne pas perdre le contrôle de l’algorithme. Dans la relation 15 les itérations sont identifiées par le facteur temps :

(15)

avec I₁ le volume image de départ qui est non filtré.

2.3.3 Résultats

Suivant l’intensité du coefficient κ et le nombre d’itérations, les résultats peuvent être très différents. Ce qu’il faut garder en tête c’est qu’il ne faut pas détruire l’information utile de l’image initiale,tout en essayant de lisser le plus possible les régions qui sont morphologiquement homogènes, facilitant par la suite le calcul de la segmentation. Ci-dessous, nous présentons les résultats obtenus avec différents paramètres de filtre sur une image du volume traité.

La figure 5(a) représente une image du volume original, avant le traitement par filtrage. Un zoom sur la région des vaisseaux du cou 5(b) permet d’observer l’effet "salt & pepper", ce qui confirme l’existence d’un bruit de fond dans l’image. La figure 5(c) illustre le résultat du filtrage par diffusion anisotrope utilisant l’option 1 (eq.9a) et les paramètres suivants : κ = 40, nombre d’itérations=3, λ = 0.06. Le bruit est fortement diminué et l’image résultat conserve bien les informations "vraies". La figure 5(d) illustre le résultat du filtrage par diffusion anisotrope utilisant l’option 2 (eq.9b) et les paramètres suivants : κ = 30, nombre d’itérations=3, λ = 0.06. Le bruit est enlevé mais les structures géométriques sont mieux préservées, la précision de la segmentation en sera améliorée. C’est l’option qui sera choisie pour la segmentation par la méthode Level Set.

---

(a) Image brute

(b) Zoom

(c) filtrage option 1

(d) filtrage option 2

Figure 5: Coupe Sagittale Oblique d’une séquence injectée. Opérations de filtrage avec les deux options du filtre de diffusion anisotrope

---

3 Extraction de la géométrie "native"

De nos jours, la segmentation des structures anatomiques est un sujet très important en imagerie médicale. En effet, c’est un problème difficile à résoudre dans bien des situations puisque les images sont généralement très bruitées et leur niveau moyen de signal n’est pas constant, en raison d’une hétérogénéité du champ magnétique et d’une collecte non linéaire du signal au niveau des antennes (section ??). La majorité des méthodes est basée sur deux grands principes : une séparation par seuils d’intensité observés sur l’histogramme du volume (segmentation par seuillage) ou bien avec des fonctions "supérieures" qui utilisent les gradients de l’image pour déplacer des hyper-surfaces fermées. Mais jusqu’à présent, aucune technique véritablement efficace et robuste n’a vu le jour, aucune qui puisse être complètement automatique aux yeux du clinicien.

Dans cette section nous allons exposer deux techniques d’extraction mises à preuve au Laboratoire sur des volumes d’IRM à effet angiographique (séquence injectée). La méthode Level Set, implémentée sous un algorithme MATLAB 7 et la méthode par seuils d’intensité, pratiquée dans l’environnement 3D de AMIRA 4.

3.1 Méthode de segmentation par Level Set

Les méthodes Level Set, développées par Osher et Sethian en 1988 [5] sont aujourd’hui utilisées avec succès dans de nombreux domaines. Quelques applications sont la reconnaissance des formes en traitement d’images, la robotique, l’optimisation des trajectoires, la réalisation d’images de synthèse et la mécanique des fluides numérique. Le principal avantage de l’approche Level Set est de pouvoir représenter implicitement une interface statique sur un maillage eulérien conventionnel. La plupart des applications utilisant cette méthode sont décrites dans un livre édité en 1999 par Cambridge University Press à Berkeley [6]. En imagerie médicale, la méthode Level Set est très répandue et l’on retrouve plusieurs applications combinant l’IRM, la mécanique des fluides numérique et la méthode Level Set. Mais elles sont encore difficilement applicables en clinique [7, 8, 9, 10].

Dans ce qui va suivre, je vais exposer le principe de la méthode Level Set et de son implémentation dans un algorithme MATLAB codé au CHU de Rangueil pendant mon travail de thèse. Dans une deuxième partie j’exposerai les résultats obtenus sur quelques exemples concrets d’extraction de géométries vasculaires.

3.1.1 Méthode Level Set et Algorithme

La méthode standard de segmentation bouge des marqueurs sur une courbe, tandis que la méthode Level Set crée un lien entre l’évolution de la courbe et les lois de conservation de l’énergie en physique. Nous allons nous intéresser au principe et aux équations de la méthode level set rapide [11].

---

Figure 6: Fonction signée de distance ϕ définie pour un rectangle.

---

L’équation pour la propagation d’un front avec une vitesse F, fonction de la courbure C, est mise en relation avec une loi hyperbolique de conservation pour la propagation de ce front. L’idée centrale est de suivre l’évolution d’une fonction ϕ de distance pour laquelle ses valeurs zero correspondent toujours à la position du front propagé (zero-distance)(Fig.6). Le déplacement de cette fonction d’évolution est déterminé par une équation différentielle partielle d’un ordre supérieur qui autorise des bords tranchants, des coins aigus et des changements de topologie dans l’ensemble des zéros qui décrivent l’interface.

Cette approche générale attire un grand nombre d’applications à trois dimensions, comme c’est le cas des problèmes qui sont sensibles à la courbure (tension de surface), ou bien ceux qui subissent des changements complexes de topologie. Pour une interface à une dimension se déplaçant dans un espace à deux dimensions (image 2D), l’algorithme Level Set est une méthode qui consomme 0(n²) par pas de temps, où n est le nombre de points dans l’espace 2D (nombre de pixels de l’image). C’est le point coûteux de la méthode. Et si nous cherchons à résoudre un problème à trois dimensions, alors ça devient très coûteux pour le calculateur avec 0(n³) opérations par pas de temps. Pour l’exemple, un volume d’image peut faire 512 × 512 × 50 = 13107200 calculs par pas de temps !

Mais une solution pour accélérer le calcul existe déjà et consiste à utiliser uniquement les points qui sont près de l’interface. Cette technique appelée Narrow Band (NB) permet de faire face au grand nombre de données présentes dans un volume 3D, tout en gardant les fonctionnalités de la méthode initiale. L’algorithme itératif fait alors les étapes suivantes :

Avec une courbe initiale,

• (1) l’algorithme calcule les distances eulériennes signées sur le "tube", de part et d’autre de la courbe d’interface. Intérieur négatif, extérieur positif.
• (2) l’algorithme met à jour la matrice du domaine et fait évoluer ϕ dans le "tube".
• (3) si la courbe est près du bord du domaine, l’algorithme ré-initialise la surface en passant par (1) et refait un nouveau maillage rectangulaire pour continuer les calculs (Fig.7).
• (4) sinon, il va sur (2).

---

Figure 7: Fonction signée de distance ϕ définie pour une zone (NB) proche à l’interface.

---

Si on définit γ₀ comme le front initial fermée et ϕ(x,t),xϵR² une fonction scalaire telle que au temps t, l’ensemble de niveau-zéro de ϕ(x,t) est la courbe γ_t. Ensuite les distances eulériennes entre la courbe γ₀ et tout point de ϕ sont définies par ϕ(x,0)  =   ± d(x). Maintenant si on considère le déplacement vers un niveau défini par ϕ(x(t),t) = C, avec x(t) la trajectoire de la particule, nous pouvons écrire que sa vitesse instantanée dans la direction normale au front de propagation est donnée par la fonction de vitesse F :

(16)

(17)

Naturellement, si l’algorithme se base uniquement sur la courbure du front pour définir la vitesse de ce dernier, notre interface ne va pas s’arrêter au niveau de la paroi vasculaire. Si nous voulons trouver la position de la paroi, nous devons alors intégrer le gradient du volume d’imagerie (gradients de niveaux de gris en 3D) dans F. La vitesse de propagation sera donc fonction de l’inverse du gradient, ainsi, à fort gradient, la vitesse sera nulle et le front s’arrêtera. Pour des problèmes d’instabilité, un gradient d’ordre 2 est aussi intégré dans l’expression de la vitesse. Il sert à piéger le front au niveau de la paroi et contrôle les débordements.

La forme générale de l’équation d’ensemble de niveaux de Sethian correspond à la relation (18). Si nous faisons intervenir les gradients de l’image et le terme attractif (gradient d’ordre 2) nous avons la relation (20). Le terme de courbure est représenté par la relation (19).

(18)

(19)

Finalement, ϕ est solution de l’équation différentielle partielle (20) avec la condition ϕ(x,t = 0) = ±d(x) sur les distances.

(20)

Avec k_I(1 - εk) la vitesse de propagation du front, k_I fonction des gradients du volume d’imagerie et k la courbure moyenne de poids ε. Le terme β∇P correspond au terme attractif où P est fonction du gradient d’ordre deux de poids β.

3.1.2 Résultats 3D

En pratique, la mise en place de l’algorithme de segmentation demande beaucoup de ressources mémoire. La taille des multiples volumes image est importante : matrice 512 × 512 × 50 voxels; avec un nombre d’images fonction de l’épaisseur de coupe, la mémoire vive d’un ordinateur de bureau est vite saturée. Plusieurs solutions sont donc envisagées pour remédier à ce problème :

• Définition d’une bande de calcul autour du front (Fig.7).
• Diminution de la taille de la matrice (Fig.8).

Dans une image acquise, beaucoup d’informations sont inutiles autour du volume à explorer. Il faut donc découper le volume à travers la sélection d’une région d’intérêt. Par exemple n’en garder que l’aorte abdominale si l’objectif est l’étude du comportement du flux sanguin et de sa paroi au niveau d’un AAA. Cette simplification peut être faite sous MATLAB (Fig.8(a)) avant l’opération de filtrage (gain de temps) ou bien sous AMIRA (Fig.8(b)). Le volume obtenu sera soumis au traitement d’extraction de la géométrie native et représente une base morphologique, pour le premier maillage statique et ses homologues dynamiques. Alors pour être cohérent en termes de centrage et taille des images, nous avons que cette opération de simplification sera identique et automatique pour les volumes dynamiques.

---

(a) Matlab

(b) Amira

Figure 8: Coupe sagittale sur un AAA, simplification du volume avec les interfaces Matlab et Amira. L’opération de filtrage se fait automatiquement après la simplification avec MATLAB

---

Pour diminuer la taille de la matrice, nous pouvons encore partager le volume image en plusieurs sous-volumes. La segmentation peut-être lancée en parallèle (à partir de plusieurs points de départ) sur chacun d’entre eux. Les fronts se propagent et se retrouvent aux frontières pour faire un seul volume.

Les premiers tests ont été élaborés sur des formes simples acquises dans des conditions "idéales", telle une bouteille remplie d’un mélange eau et liquide de contraste (iode). Comme il y a un grand nombre d’images (près de 350 coupes au scanner), il est nécessaire de partager le volume en plusieurs petits volumes et de les segmenter bout à bout. L’algorithme effectue sa tâche. Le front se propage sans problème dans les coupes fines (environ 1 mm d’épaisseur) et le contraste est excellent. Il n’y a pas de problèmes d’arrêt du front. La forme de la bouteille est retrouvée aisément (Fig.9). Le premier résultat est mis en forme sous Matlab grâce à la fonction "isosurface", on ne voit que le haut de la bouteille et on peut remarquer la présence de bulles d’air sur le côté (la bouteille était couchée et remplie d’eau lors de l’acquisition des images).

---

Figure 9: Premier résultat par la méthode Level Set sur une bouteille couchée remplie d’eau, affichage du volume avec l’interface de MATLAB

---

Un deuxième résultat est obtenu sur un fantôme d’anévrisme de l’aorte thoracique (Fig.10(a)). Pour cette expérimentation nous avons utilisé un montage à écoulement pulsé grâce à un ventricule de prothèse de Thoratec. Nous avons directement obtenu les images dynamiques (20 phases) avec une séquence IRM du type balanced-SSFP (Fig.10(b)). Étant donné que les images étaient de bonne qualité, nous avons fait plusieurs tests à différents moments du cycle cardiaque, ce qui nous a permis de reconstruire la géométrie dynamique pour ce fantôme (Fig.10(c)).

---

(a) fantôme d’AAT

(b) Images b-SSFP

(c) Extraction Level Set multi phases

Figure 10: Résultats d’extraction par la méthode Level Set à partir d’une séquence dynamique b-SSFP sur un fantôme d’AAT. Affichage avec l’interface AMIRA 4

---

3.2 Méthode de segmentation par seuillage et proche voisin

Les travaux réalisés par deux étudiants en médecine (un en 4^ième année et l’autre, interne en chirurgie vasculaire) dans le cadre d’unités d’enseignement, ont été l’occasion d’évaluer les techniques d’extraction basées sur deux méthodes. La méthode Level Set (automatique mais lente) que nous venons de voir et la méthode par seuillage et proche voisin (plus manuelle mais rapide). L’application de la méthode de filtrage et d’extraction des géométries vasculaires a été bien intégrée par les deux étudiants. Mais ils se sont familiarisés beaucoup plus vite avec AMIRA qui a une interface intuitive en 3D. Les images provenant de séquences à effet angiographique (Section ??) ont été importées au format ANALYZE avec MRIcro (Sec.1.1) pour être ensuite filtrées avec MATLAB et traitées dans l’interface AMIRA 4. Ce travail a eu lieu au Laboratoire d’Imagerie Interventionnelle Vasculaire du Service du Prof. H. Rousseau, au CHU de Rangueil.

Aucune difficulté pour la prise en main du logiciel. Pour ce travail nous avons utilisé une méthode de segmentation de l’histogramme du volume, combinée à une contrainte de connectivité entre les voxels voisins. A travers une interface 3D, l’utilisateur a pu sélectionner un seuil sur l’histogramme associé à l’étude. En général, pour ce type de séquences, l’histogramme est divisé en deux pics très nets, alors le placement du seuil est assez facile. Ensuite, tout le travail consiste à vérifier et corriger les voxels qui ont été pris comme étant à l’intérieur de la lumière artérielle. L’interface AMIRA permet aussi de bien définir les limites de la zone à explorer. Une fois définis, le programme AMIRA réalisera une extraction de la surface par dessus ces voxels, et ceci par une méthode de triangulation sur la surface. Le rôle du stagiaire était d’arriver à cette reconstruction de surface en accord avec les images IRM et tout autre matériel à disposition. En cas de doute, certains patients avaient d’autres types d’examens, tels qu’un scanner ou une angiographie.

3.2.1 Exemple d’une dissection de l’aorte

Cette pathologie a été analysée grâce à l’outil informatique. L’aorte était celle d’un patient souffrant d’une pathologie vasculaire complexe : il avait développé une dissection aortique dont le traitement avait été médical 2 années plus tôt. Aucun geste thérapeutique n’ayant été réalisé, l’aorte comportait au lieu d’une seule lumière, un vrai et un faux chenal. Le contrôle évolutif avait permis d’objectiver la survenue d’un anévrisme sacciforme de l’aorte au niveau de la crosse. L’attitude médicale a consisté en la pose d’un Stent-Graft localisée au site de l’anévrisme.

Pour l’extraction de la géométrie avant la mise en place du Stent-Graft, une séquence IRM à effet angiographique a permis d’identifier le vrai et le faux chenal de la dissection, ainsi que l’anévrisme sacciforme développé au niveau de l’arche aortique. Nous avons décidé de garder l’aorte dans toute sa longueur pour observer d’un côté les portes d’entrée et de réentrée de la dissection chronique, et de l’autre, l’anévrisme au niveau de l’arche, raison de cet examen diagnostique avec produit de contraste.

Pour ce dossier, il a été demandé de préparer plusieurs reconstructions.

• une première avec la morphologie réaliste,
• la deuxième avec quelques simplifications (collatérales),
• la troisième avec l’exclusion virtuelle du sac anévrismal
• et la quatrième avec la fermeture virtuelle d’une porte d’entrée au niveau de la dissection.

---

(a) Gestions des objets

(b) Editeur de segmentation

Figure 11: Gestion des objets dans le "pool" et Interface de l’éditeur de segmentation dans AMIRA 4. Segmentation par seuillage de l’histogramme et contrainte de voxels voisins

---

A travers les fenêtres de l’éditeur de segmentation (Fig.11(b)), les images sont analysées et les structures anatomiques sont identifiées afin de décider de la démarche à suivre. L’histogramme donne un point de départ car il choisit des voxels qui sont au dessus du seuil. La contrainte des voxels voisins impose que tous les voxels ayant une intensité supérieure au seuil et étant connectés au premier voxel sélectionné, sont dans la lumière artérielle. Ceci permet d’exclure certaines zones comme le cœur, et quelques collatérales, en effaçant leur connexion. La figure 13(a) et 13(b) illustrent l’aorte avec ses principales collatérales, tandis que la figure 13(c) et 13(d) montre déjà une simplification. Ces manipulations sont nécessaires pour la préparation du maillage destiné aux calculs CFD. L’entrée et les sorties seront par la suite des "patch" bien définis avec des CL hémodynamiques. Ils doivent correspondre à la position des coupes de vélocimétrie par contraste de phase effectuées par IRM. Actuellement il est très difficile de faire des mesures de flux assez précises au niveau des coronaires, du tronc cœliaque et des artères intercostales. Alors elles sont naturellement exclues virtuellement, sachant que la perte de flux par ces collatérales n’est pas grande face au volume d’ejection systolique. Toutes les décisions sont basées sur d’autres examens s’il y en a, tel que l’angiographie ou le scanner, la figure 3.2.1, montre des images de scanner centrées sur une porte d’entrée d’une dissection aortique.

---

Figure 12: Images scanner suivant différentes orientation sur une porte d’entrée d’une dissection aortique, reconstruction MRP

---

D’autres exclusions virtuelles ont été possibles, la figure 14(a) illustre une sélection de voxels à effacer lors de l’exclusion de l’AAT (Fig.14(b). La figure 14(d) expose la possibilité de fermer virtuellement une porte d’entrée dans une dissection de l’aorte. A partir de ces reconstructions surfaciques, des maillages statiques décrivant chaque cas peuvent être créés dans la même interface.

---

(a) Initiale devant

(b) Initiale derrière

(c) Simplification

(d) Derrière

Figure 13: Première reconstruction et simplification par l’exclusion des principales collatérales

---

---

(a) Identification virtuelle

(b) Exclusion de l’AAT

(c) Diss 2 portes

(d) Diss 1 porte

Figure 14: Manipulations virtuelles, exclusion numérique de l’anévrisme et fermeture numérique d’une porte d’entrée pour la dissection aortique

---

4 Déplacement de la géométrie

Dans le chapitre précédent nous avons présenté les difficultés à obtenir des acquisitions dynamiques de bonne qualité au niveau de l’aorte thoracique. Cette difficulté ne se présente pas au niveau de la bifurcation carotidienne, car l’accès est plus facile, le patient n’a pas besoin de garder des apnées successives et parce que l’antenne Synergy Head-Neck a un excellent rapport signal sur bruit. Les antennes utilisées, avec les paramètres de séquence décrits à la section ?? pour les explorations cardiaques, ne permettent pas actuellement d’obtenir des images morphologiques à la fois bien résolues et dynamiques. Nous avons vu qu’il était possible de coupler les informations de deux séquences pour en tiré le meilleur profit en résolution spatiale (séquence à effet angiographique) et temporelle (séquence SSFP dynamique). Finalement, les techniques d’extraction (automatiques ou manuelles) sont capables de passer d’un volume morphologique de bonne qualité à une reconstruction de surface des formes anatomiques.

Dans cette partie de mon travail je vais présenter la méthode qui permet d’animer le mouvement de ces formes anatomiques mais qui peut aussi bien animer tout autre objet qui soit placé dans un champ de déformation qui sera estimer à partir des images dynamiques à basse résolution spatiale. Nous sommes ici face à une technique hybride qui présente des avantages énormes pour la mécanique des fluides numériques. Ces avantages seront développés dans le chapitre suivant.

4.1 Principe de déformation

Une déformation est une opération qui déplace un objet dans son espace vers une nouvelle position de cet espace. Une déformation peut être linéaire (homothétie) ou bien non-linéaire (à chaque point lui correspond un vecteur de déplacement différent et indépendant). Nous avons décidé de ne pas utiliser les lois de comportement pour la paroi décrites dans la littérature car nous estimons qu’il existe un risque très grand de faire des erreurs sur son comportement biomécanique (rhéologie). Face aux structures très variées comme l’os, les parois malades (athérosclérose, anévrisme) et les organes environnants, nous avons opté pour imposer sa position à tout instant.

Pour partir d’une géométrie statique de la paroi et reconstruire l’évolution de cette paroi au cours d’un cycle cardiaque, nous devons alors imposer à chaque point de cette paroi et pour chaque phase cardiaque, un vecteur de déplacement indépendant. Les avantages de cette approche sont bien établis par rapport à l’extraction de plusieurs géométries, avec leurs reconstruction de surface indépendante :

• Une seule segmentation à faire, au lieu d’avoir à faire une par phase cardiaque (~ 20).
• Une cohérence topologique intrinsèque, par le principe de déformer la même structure, au lieu de travailler avec des structures non conformes les unes avec les autres. C’est l’avantage principal pour les calculs CFD à base de formulations ALE (Sec.??).

     

Le volume d’images dynamiques est un ensemble de volumes statiques acquis à des instants différents et qui donnent l’impression d’un effet ciné. Chaque volume est identifié sur le cycle cardiaque grâce au temps d’acquisition après l’onde R de l’ECG. Le cycle cardiaque est segmenté en 20 phases dont une seule correspond de très près au délai de lancement de l’acquisition statique (séquence injectée). Nous retrouvons alors une correspondance temporelle entre les deux séquences. Néanmoins la correspondance géométrique doit être satisfaite pour éviter de faire des erreurs pendant l’application du champ de déformation. Ce dernier est estimé de façon statistique entre l’information contenu dans une phase appelée "source" et une phase appelée "cible". La phase "source" sera toujours celle qui correspond au même délai après l’onde R que pour la séquence statique (séquence injectée). L’estimation d’un champ de déformation pour une phase se fait donc entre deux "volumes phase". L’alignement de tous les champs de déformation permet de reconstituer l’évolution du champ spatial dans le temps.

4.2 Méthode de déformations non-linéaires

Le programme de départ est une boîte à outils basée sur des algorithmes de probabilité (bayessiens). Spécifiquement conçu pour faire des tests statistiques sur des hypothèses d’activation cérébrale, SPM (Statistical Parametric Mapping, K. Friston et J. Ashburner) est développé par les membres et collaborateurs du Wellcome Department of Imaging Neuroscience (UK, http://www.fil.ion.ucl.ac.uk/spm/). Mais si nous utilisons quelques fonctionnalités de ce programme, nous n’allons pas opter pour l’approche statistique, nous allons plutôt voir les déformations comme un problème d’"optimisation".

4.2.1 Principe d’optimisation

Concrètement, l’approche "optimisation" est beaucoup plus naturelle. Nous cherchons à trouver le meilleur compromis entre transformation efficace par la fonction F₂ (eq.21b) et transformation simple par la fonction F₁ (Jacobien proche de 1)(eq.21a). L’optimisation est retrouvée en minimisant une combinaison linéaire de F₁ et F₂ (eq.21c).

(21a)

(21b)

(21c)

L’optimisation d’une transformation idéale est un travail itératif pendant lequel le domaine de calcul (le champ de déformation) peut aller dans les mauvaises solutions (minima locaux). Il faut le voir comme si on devait se rendre quelque part et on était pas sur de connaître le bon chemin. Il y a une infinité de chemins mais seul quelques uns sont vraiment optimum. Il peut aussi nous arriver de tomber dans une impasse et de se retrouver tout prêt de la destination voulu mais de ne pas y avoir parvenu. Les paramètres qui décrivent ce problème d’optimisation sont alors les suivants :

• La solution n’est pas unique
• La solution doit être bijective.
• La transformation est non linéaire, en général elle est représentée par un champ vectoriel de déformation.
• Un algorithme existe pour l’analyse statistique des images cerébrales. Il nous sera utile pour trouver la bonne solution.
• La contrainte sur le Jacobien n’est pas suffisante. Il faut éviter l’écrasement des mailles car elles peuvent entraîner des désaccord numériques lors des calculs en MFN. L’utilisation d’un filtrage dégressif (gaussien) limite ce problème.

Pour obtenir une transformation T idéale, le processus itératif cherche à maintenir le déterminant du Jacobien |J| de T le plus proche possible de l’unité. F₁ et F₂ sont obtenues numériquement dans un système de mailles aux éléments finis (champ de déformation). Leur gradients par rapport à T(x_i) sont obtenus par differences finies.

Un simple algorithme de descente est alors utilisé:

(22a)

(22b)

Si T est idéale, le volume source transformé s’ajuste parfaitement au volume cible. L’algorithme statistique analyse la variance obtenue entre les deux et donne une information pour ajuster au mieux les nœuds du champ de déformation. La combinaison de la transformation T et de la transformation inverse T^-1 est égale à l’identité si T est idéale. Ce qui veut dire qu’un volume transformé avec T puis avec T^-1 correspond "exactement" au volume de départ, sinon, T n’est pas la solution idéale. Théoriquement ceci est vrai sur les équations mais dans la pratique, les erreurs d’interpolation inhérents aux techniques discrètes, ne permettent pas de retrouver exactement le point de départ.

4.2.2 Algorithme

Le processus d’optimisation se déroule par une suite d’algorithmes que j’ai développé pendant mon travail de thèse. Je me suis inspiré des travaux de J.Ashburner [12, 13] et j’ai créer un code qui re-utilise certaines interfaces et routines mex-files (codées en c) disponibles dans la toolbox de SPM.

Le nouveau code a été dédié à estimer les déformations de la paroi vasculaire entre la phase native et les autres phases cardiaques (Fig.??). Ceci revenait à estimer les déformations de l’espace dans une région d’intérêt. L’algorithme a été organisé en modules pour avoir un accès facile aux résultats intermédiaires et pouvoir développer indépendamment chaque partie du code source.

Nous avons défini les tâches suivantes, qui sont en pratique chaque partie de l’algorithme général. Elles sont constituées de contrôles _ui (User Interface) pour faciliter l’échange entre l’utilisateur et le programme :

• makeroi_ui a deux fonctions très simples. Il permet de couper le volume brut pour en garder que la région d’intérêt et permet de fabriquer un masque sur les zones d’entré et sortie de flux. Ceci annule l’information du signal et ne provoque pas de déformations. (Utile pour les "patchs" d’entrée et sortie lors d’un calcul CFD). C’est une façon de fixer les structures des zones sélectionnées. Les deux opérations se font dans une fenêtre active et les coordonnées des points repères sont sauvegardés pour éviter de refaire la même manipulation pour chaque phase.
• multifiltre_ui permet d’enlever le bruit de fond; il se sert du filtre de flou selectif avec un SD = 14. Il se sert du filtre de diffusion anisotropique en 3D pour rendre homogènes les régions qui correspondent à la lumière vasculaire. Ces tâches sont répétées sur tous les volumes (les phases) sélectionnés.
• multi_nldf_ui lance le processus itératif qui détermine le champ de déformation non linéaire entre le volume natif et les autres phases cardiaques. C’est une boucle qui permet de combiner les paramètres de filtrage gaussien (LMH, Fig.2(b)) et de régularisation avec des valeurs dégressives, afin de trouver le minimum général pour la valeur résiduelle, signe d’une bonne transformation.
• multidef_ui calcule les nouvelles coordonnées des points qui ont été chargés pour la déformation. Elle utilise le champ de déformation pour re-écrire les nouvelles positions sans changer la topologie de l’objet.

Les modules avec le prefixe multi* ont été programmés avec l’intention de couvrir toutes les phases du cycle cardiaque. C’est alors qu’à partir de la phase de référence initiale appelée "native", l’algorithme va estimer toutes les autres phases cardiaques suivant la logique décrite à la figure 4.2.2. La construction de tous les maillages se fera donc à partir du premier maillage. Afin d’intégrer la référence temporelle, l’ensemble de tous les maillages sera organisé suivant le numéro de la phase cardiaque à laquelle chacun est associé. Nous verrons par la suite à la section ?? que l’intitulé "native" change de principe d’utilisation dans le déroulement du calcul CFD.

---

Figure 15: Logique de construction des maillages synchronisés sur le cycle cardiaque à partir du maillage natif et les transformations T.

---

Pour chaque paire de volumes analysés, notre algorithme donne à la sortie :

• Un fichier contenant le champ de déformation (Fig.20(d)).
• Un fichier *.txt contenant les nouvelles coordonnées des points demandés en entrée.
• Un graphique avec l’évolution de la variance résiduelle obtenue au cours du processus itératif (Fig.18(f)).

C’est important de retrouver une organisation homogène des mailles dans le nouveau champ car l’écrasement (mailles très petites à côté des grandes) est une cause d’erreur dans la construction du nouveau maillage pour les CFD. Il existe donc un intérêt à partir avec une valeur de LMH équivalente à 10 fois la résolution de l’image (en mm).

4.2.3 Probabilité de trouver la "bonne" déformation

Les images sont recalées en estimant un ensemble de paramètres qui maximise la probabilité d’avoir obtenu une transforation idéale. Ceci implique de commencer par un ensemble d’évaluations de départ, et de faire de petits ajustements de sorte à faire augmenter le potentiel de probabilité à la fin de la boucle itérative. Dans chaque itération, les positions des points de contrôle (noeuds) sont mis à jour in situ, en balayant séquentiellement le volume cible. Lors d’une itération de recalage tridimensionnel, la boucle peut fonctionner du bas vers le haut, postérieur vers le devant, et de la gauche vers la droite. Dans la prochaine itération, l’ordre de la mise à jour est renversé (haut vers le bas, antérieur vers le postérieur, et de la droite vers la gauche). Cette séquence alternative est continuée jusqu’à ce qu’il n’y ait plus une réduction significative de la pente du potentiel de probabilité, ou pour un nombre fini d’itérations. Chaque itération de l’optimisation implique de déterminer le taux de changement du potentiel de probabilité par rapport aux changements minuscules de chaque élément de T. Pour la nième itération, ε est choisi suffisamment petit pour satisfaire la condition 22b.

---

(a) Évaluation bidimensionel

(b) Division en tétraèdres

Figure 16: Évaluation du déplacement des nœud sur le maillage du champ de déformation. Évaluations pour un cas 2d et sous divisions en tétraèdres pour l’espace entre huit centres (sphères) de voxels dans le volume.

---

Le taux de changement est estimé à chaque boucle de programme. Lors de la mise à jour, chaque noeud est déplacé le long de la direction qui diminue le plus rapidement possible le potentiel de probabilité à travers la relation 22a (méthode de descente de gradient). Pour le recalage bidimensionnel, ∇F dépend des changements de la matrice du Jacobien des six triangles adjacents représentés sur le schéma de la figure 16(a). Puisque les mathématiques de calcul de ces dérivés partielles sont algébriquement denses, un sous-programme en C (mex file) est utilisé sur la base de la toolbox de SPM. Pour le cas tridimensionnel, le déplacement du noeud dans le maillage va influencer les matrices du Jacobien des tetrahèdres qui ont un sommet à ce noeud (Fig.16(b)). Si un noeud est déplacé de façon exagérée, alors le déterminant d’une ou plusieurs des matrices du Jacobien lié à un triangle ou à un tétraèdre voisin peut devenir négatif. Ceci impliquerait une violation de la contrainte linéaire dans le tracé (puisque le tetrahèdre voisin occuperait le même volume). La tentative initiale déplace le noeud d’une petite distance. Si l’un des déterminant devient négatif, alors la valeur de ε est réduite de moitié et une autre tentative est faite pour déplacer le noeud d’une distance plus petite de son emplacement original. Ceci continue pour le noeud jusqu’à ce que les contraintes soient satisfaites. Un procédé semblable est alors répété par lequel la valeur de ε est divisé en deux jusqu’à ce que le nouveau résidu soit moins que ou égal à la valeur précédente. En incorporant ce procédé, le résidu n’augmentera jamais lors du déplacement du noeud, assurant ainsi que le potentiel de probabilité d’avoir trouvé la "bonne" déformation augmentera avec chaque itération.

4.2.4 Inverser un champ de déformation

Á la fin des itérations du module nltf_ui, il est nécessaire de calculer l’inverse du champ de déformation. Deux raisons principales pour cette action :

• La qualité du champ de déformation obtenu est estimée par l’opération Id = T × T^-1. La régularité du maillage dans Id est une bonne indication du résultat obtenu.
• Pour échantilloner les nouvelles coordonnées des nœuds, il est nécessaire de passer par l’inverse de la transformation. Ceci parce que la représentation du nouveau volume doit être transférée sur le même espace que le point initial, dans l’espace euclidien de départ.

4.2.5 Exemple 3D avec un Cas test : Modèle "dalle Carrée"/ "dalle Ronde"

---

(a) Volume source filtré

(b) Volume cible filtré

(c) Jacobien

(d) Champ de T à

(e) Champ de T-1 à

(f) Champ Id = T × T-1 à

Figure 17: Cas test sur deux volumes images. T correspond à la transformation "idéale" entre le volume "source" et le volume "cible", quelques erreurs sont retrouvés dans la matrice identité Id.

---

---

(a) Champ de T à

(b) Champ de T-1 à

(c) Maillage natif (source)

(d) Maillage déformé

(e) Maillage avec T × T-1

(f) Evolution de la valeur résiduelle

Figure 18: Cas test, champ de T et T^-1 avec une meilleure résolution spatiale (). Déformation du maillage natif vers son nouveau champ à travers T. Le retour en arrière ne donne pas exactement le même maillage de départ par erreur d’interpolation et d’optimisation, la valeur résiduelle n’est pas nulle à la fin du processus itératif

---

Pour exposer la méthode avec un exemple académique, nous avons choisi d’obtenir la transformation entre deux volumes d’images représentant respectivement une "dalle carrée" et une "dalle ronde" (Fig.17(a) et Fig.17(b)). L’algorithme filtre les volumes avec un filtre gaussien pour retrouver progressivement la transformation T à travers l’optimisation de leur Jacobien (Fig.17(c)). La figure 17(e) montre l’allure des mailles du champ de déformation avec une résolution équivalente à d’un voxel dans le volume image de référence. Pour éviter les problèmes d’interpolation lors de l’écriture des nouvelles coordonnées, les mailles peuvent être réduites à de la résolution d’origine (Fig.18(b)). La dalle "native" (Fig.18(c)) se déforme vers son nouveau espace sous forme d’un nouveau maillage avec la même topologie de référence (connectivité)(Fig.18(d)), seul les coordonnées des nœuds ont été modifiées. Mais l’opération inverse partant du nouveau maillage vers la "source" ne permet pas d’obtenir exactement la même "dalle" d’origine (Fig.18(d)). Des erreurs résiduelles sont encore présentes et la transformation (Fig.18(f)). Hormis ces problèmes d’optimisation, ce cas test nous a permis de démonter que la méthode était capable d’appliquer un champ de déformation vectoriel (non-linéaire) aux coordonnées d’un maillage "natif".

4.2.6 Exemple 3D avec un modèle de fantôme d’anévrisme

Une expérience IRM sur un fantôme d’anévrisme à paroi déformable nous a permis d’utiliser les informations d’une séquence dynamique b-SSFP multi-shot pour animer un maillage "natif" vers ses homologues, au cours d’un cycle de pulsation. Nous allons décrire la partie qui concerne uniquement les déformations du maillage natif à partir des volumes images. La description du montage et ses applications CFD seront présentées au chapitre suivant.

La géométrie est obtenue avec la méthode level set en 3D (Fig.20(c)) sur un volume d’imagerie dynamique centré sur l’anévrisme du fantôme correspondant à la phase diastolique (Fig.20(a)). Le volume est discrétisé avec AMIRA 4 vers un maillage de surface contenant 4360 faces triangulaires. Pour corriger les défauts aux bords des "patchs" d’entrée et sortie nous avons utiliser une combinaison entre un algorithme MATLAB et le programme AMIRA 4. Après avoir identifié les différentes zones (entrée, sortie et mur) sur le maillage de surface, un maillage volumique est réalisé avec AMIRA 4. Ce maillage appelé "natif" est constitué de 7455 nœuds et 36473 tétraèdres (Fig.20(e)).

Afin d’illustrer les résultats d’une transformation entre phase "diastolique" et phase "systolique", nous avons sélectionné un deuxième volume à l’opposé du cycle de pulsation (Fig.20(b)). La transformation nécessaire entre les deux volumes est optimisée après 200 itérations. Le nouveau maillage est écrit avec la même connectivité et sa correspondance spatiale est validée visuellement sur des coupes d’imagerie à travers l’interface graphique d’AMIRA (Fig.20(f)). De la même façon nous pouvons faire la même opération pour chaque phase du cycle de pulsation.

Nous avons eu des problèmes après avoir estimer la transformation pour une phase en particulier, car, à cause des artéfacts de flux, le filtrage n’arrive pas a enlever les irrégularités du signal à l’intérieur de la lumière vasculaire. Ce phénomène entraîne des conséquences sur le champ de déformation et donc sur les déplacements des nœuds dans un maillage. Pour cette raison, il est vivement conseillé de lancer des segmentation automatiques dans les phases à modéliser pour en extraire des références volumiques homogènes (Fig.2(e)). Le cas contraire, un champ contenant les conséquence du bruit est alors obtenu (Fig.4.2.6).

---

Figure 19: Champ de transformation avec l’influence d’un fort bruit dans l’image. Bruit d’origine artéfact de flux

---

---

(a) Volume "diastolique"

(b) Volume "systolique"

(c) Diastole par Levelset

(d) Champ de déformation

(e) Maillage "natif"

(f) Validation géométrique

Figure 20: Exemple 3D avec un modèle de fantôme d’anévrisme fusiforme à paroi élastique. L’extraction de la lumière vasculaire réalisée par la méthode Level Set est discrétisée dans un maillage volumique non-structuré. Les images "diastoliques" et "systolyques" permettent d’estimer une transformation T entre les deux volumes. Son champ de déformation est utilisé pour déplacer les nœuds du maillage "natif" vers ses nouvelles coordonnées. Une vérification visuelle de la nouvelle géométrie est disponible sur un interface en 3D.

---

5 Flux artériels, ajustement et correction

Les images de vélocimétrie par contraste de phase nous donnent l’information sur les conditions d’écoulement aux frontières de notre domaine numérique (maillage). Le repérage spatial est disponible sur les reconstructions MIP (Fig.??) à travers l’interface de navigation de l’appareil d’imagerie. Nous allons voir comment une fonction de Fourier (eq.23) ajuste les échantillons de vitesse (~ 30) à une fonction périodique dans le temps. Cette fonction sera utilisée comme CL au moment d’imposer les flux en entrée et en sortie du domaine numérique. Les informations sont accessibles dans un fichier qui se présente par un tableau de données, elles sont disponibles après le contourage des ROI sur l’image de vélocimétrie. Ces conditions représentent la dynamique d’un fluide qui est non-compressible. Pour rester cohérent avec cette hypothèse, nous allons considérer que la somme des flux à travers tous les "patchs" doit être nulle dans l’ensemble du cycle et dans le domaine considéré.

L’ajustement avec une série de fourier à 7 harmoniques nous a parue convenable pour représenter tout type de profil de vitesse sur la base des 30 échantillons disponibles.

Rythme card.           :       67 bpm  
Intervalle-RR          :      896 ms       (à partir du r.c.)  
 
coupe 1                 | AoAscd  AoDescd  
                        |  
---------------------------------------  
Débit (ml)              |   63.8   44.1  
Vol. flux avant (ml)    |   64.1    0.4  
Vol. flux arrière (ml)  |    0.3   44.5  
Fract. régurgitante (%) |    0.5    0.8  
Débit abs. (ml)         |   64.4   44.9  
Flux moyen (ml/s)       |   71.3   49.3  
Distance systolique (cm)|    6.8    5.9  
Vitesse moyenne (cm/s)  |    7.6    6.6  
 
AoAscd, coupe 1  
 
   Num   Délai   Flux   Aire  Aire   Vit.   Vit.   Vit.   Vit. Ecart  
           de                       moy.   max.   min. optima  type  
      déclenc                                                   vit  
           ms   ml/s    cm2 pixel   cm/s   cm/s   cm/s   cm/s  cm/s  
-------------------------------------------------------------------  
    1       0    9.1    8.5   545    1.1    6.3   -3.3    6.3   1.6  
    2      30   37.9    8.7   555    4.4   10.4   -3.9   10.4   2.4  
    3      60  171.0    8.9   567   19.3   33.4   -9.6   33.4   8.2  
.....  
.....  
   29     842   15.8    9.1   580    1.7    7.4   -4.4    7.4   1.9  
   30     872    9.9    8.9   570    1.1    7.4   -4.8    7.4   1.8  
===================================================================

La relation entre ces mesures et la fonction de flux est la suivante :

(23)

Où T_cardiaque est la période cardiaque et a₀ le décalage de cette courbe par rapport à zero.

---

textures.def (a) Fantôme 1

textures.def (b) Fantôme 2

textures.def (c) Fonction ajustée sur 30 échantillons de vitesse en entrée du fantôme 1

Figure 18: Mesures expérimentales sur fantôme. Une seule entrée pour un écoulement pulsé. Hypothèse de fluide nul absolu non retrouvée sur les mesures de vélocimétrie par contraste de phase. Décalage équivalent au premier harmonique de la relation.

---

Des mesures expérimentales sur fantôme (segment de chambre à air Fig.?? et un gant en latex (Fig.??) avec une seule entrée et un débit pulsé périodique, nous ont permis de faire des corrections sur les mesures de vélocimétrie par contraste de phase. Il s’agissait d’obtenir un flux nul en entrée du compartiment. Mais plusieurs mesures par IRM ont donné des valeurs non nulles pour le flux total pendant la durée de la pulsation. Ce qui voulait dire que notre fantôme aurait dû se collapser (flux négatif) ou bien exploser (flux positif). Pourtant la réalité expérimentale montrait que le fantôme avait un mouvement de paroi périodique et stable. L’analyse par harmoniques de Fourier (eq. 23) de la fonction de débit Q nous a donné une valeur non nulle pour la composante constante (a₀) de ce modèle (Fig.??) Nous avons repris l’expérience en plaçant une bouteille remplie d’eau à côté de chaque fantôme. La mesure de la vitesse de l’eau (immobile) dans la bouteille était non nulle, elle était très proche de a₀, l’offset retrouvé donné par le modèle de Fourier ajusté aux échantillons de mesure. Basé sur ces résultats nous avons fait des corrections sur les mesures de vélocimétrie sur la base de la composante constante (offset) d’une ROI placée sur des structures immobiles présentes dans la zone d’exploration (Fig.??). Il est évident que cette hypothèse ne peut pas être confirmée avant de faire une analyse statistique basée sur un très grande multitude de mesures, mais nous avons opté par cette solution qui nous a parue la plus conforme aux réalités observées.

L’autre type de correction que nous avons effectué est basée sur l’hypothèse de la somme des flux nulle dans tout le domaine exploré. Logiquement, l’addition ne donne jamais une valeur nulle à cause des pertes de flux non mesurées (petit) à travers les collatérales qui partent du domaine étudié. Pour remédier à ce problème d’incohérence nous avons corrigé le flux des sorties avec un coefficient de pondération.

Conclusion

Chaque section de la chaîne de post-traitement a été analysée et plusieurs techniques utilisées ont été développées dans le courant de ce chapitre. Les exemples cités sont la base d’un travail d’optimisation qui couvre l’étape entre les images brutes et l’ensemble d’informations nécessaires à la mécanique des fluides dans le domaine discret.

Nous avons pu donner des réponses aux problèmes qui se sont presenté lors de l’analyse des images médicales concernant le filtrage et le ré-haussement du contraste. Les solutions proposées, par combinaison des différents filtres développés pendant ce travail, nous ont permis d’engager le processus d’extraction de la géométrie native par deux techniques différentes. La méthode d’extraction automatique, basée sur la méthode Level Set en 3D, qui est puissante par la simplicité de son utilisation, mais qui reste d’application lente et parfois limitée par la qualité des images qui lui sont proposées. La méthode manuelle, qui est accessible aux utilisateurs cliniciens à partir d’une interface graphique et qui autorise les corrections et les modifications virtuelles de la géométrie obtenue. Les résultats des deux techniques sont validés visuellement sur la même interface par la superposition des données brutes avec le résultat de la segmentation, ce qui limite les erreurs dans le processus d’extraction de la géométrie "native". La dynamique de la paroi artérielle, qui correspond à celle d’un maillage mobile, a été associée à la géométrie "native" par l’intermédiaire d’une méthode de déformations non linéaires de l’espace. Cette technique de déformation estime ses transformations sur la base des variations du signal d’une séquence d’IRM dynamique. La qualité du processus de déformation est d’une précision équivalente à la résolution spatiale des images dynamiques mais garde les détails de surface avec une précision équivalente à la résolution spatiale de la séquence injectée. La validation des géométries composant le maillage mobile est effectuée visuellement sur une les volumes dynamiques à travers l’interface citée plus haut. Finalement, l’analyse des images de vélocimétrie, avec une technique d’ajustement basée sur une analyse par harmoniques de Fourier, complètent l’ensemble de données nécessaires pour résoudre les équations de mécanique des fluides numériques. Les conditions d’acquisition de transformation et de préparation des résultats, fournissent des paramètres réalistes et spécifiques pour chaque cas étudié.

Cette méthode générale ne peut être applicable actuellement que sur la grande circulation et collatérales de calibre centimétrique. L’application future de ces méthodes pour l’exploration des coronaires sera possible dès que l’IRM pourra fournir des images haute résolution infra-milimétriques de bonne qualité.

References

[1]   H Soltanian Zadeh Y. A., Windham JP. A multidimensional nonlinear edge-preserving filter for magneticresonance image restoration. IEEE Transactions on Image Processing, 4(2), 147 – 161, Feb 1995.

[2]   Peters R. A. A new algorithm for image noise reduction using mathematical morphology. IEEE Transactions on Image Processing, 4(3), 554–568, May 1995.

[3]   Gensanne D., Josse G., Lagarde J. M., and Vincensini D. A post-processing method for multiexponential spin-spin relaxation analysis of mri signals. Phys Med Biol, 50(16), 3755–3772, Aug 2005.

[4]   Perona P. and Malik J. Scale-space and edge detection using anisotropic diffusion. IEEE Transactions on Pattern Analysis and Machine Intelligence, 12(7), 629–639, July 1990.

[5]   Osher S. and Sethian J. A. Fronts propagating with curvature dependent speed: Algorithms based on hamilton-jacobi formulations. Journal of Computational Physics, 79, 12–49, 1988.

[6]   Sethian J. A. Level Set and Fast Marching Methods, Evolving Interaces in Computational Geometry, Fluid Mechanics, Computer Vision and Material Science. Dept. of Mathematics, University of California, Berkeley, 1999.

[7]   Sethian J. A. A fast marching level set method for monotonically advancing fronts. Proc Natl Acad Sci U S A, 93(4), 1591–1595, Feb 1996.

[8]   Bertalmio M., Sapiro G., and Randall G. Region tracking on level-sets methods. IEEE Trans Med Imaging, 18(5), 448–451, May 1999.

[9]   Watanabe M., Kikinis R., and Westin C.-F. Level set-based integration of segmentation and computational fluid dynamics for flow correction in phase contrast angiography. Acad Radiol, 10(12), 1416–1423, Dec 2003.

[10]   Deschamps T. and Schwartz P. Vessel segmentation and blood flow simulation using level-sets and embedded boundary methods. International Congress Series, 1268, 75–80, 2004.

[11]   Adalsteinsson D. and Sethian J. A. A fast level set method for propagating interfaces. Technical report, Laurence Berkeley Laboratory and Department of Mathematics, University of California, Berkeley, CA 94720, 1994.

[12]   Ashburner J., Andersson J. L., and Friston K. J. Image registration using a symmetric prior–in three dimensions. Hum Brain Mapp, 9(4), 212–225, Apr 2000.

[13]   Ashburner J., Andersson J. L., and Friston K. J. High-dimensional image registration using symmetric priors. Neuroimage, 9(6 Pt 1), 619–628, Jun 1999.

ref:

http://www.ocfia.org/pages_iframe/traitement/traitement.html

 http://www.ocfia.org/?page_id=76