Projet image 2015 : Segmentation et recalage d'images 3D cardiaques

De Ensiwiki
Aller à : navigation, rechercher


Projet Image 2015 : Segmentation et recalage d'images 3D cardiaques
Projet Projet de spécialité 2A 2015
Thème Imagerie médicale
Outil MatLab
Étudiants
Krewcun Camille (MMIS)
Guéry Bastien (MMIS)
Némausat Tara (MMIS)
Encadrant Desvignes Michel


Le contexte

Un contexte d'actualité

Schéma du coeur

Les maladies coronaires chroniques constituent un problème de santé majeure à l’heure actuelle, et peuvent entraîner des crises cardiaques causées par un manque d’oxygénation et la mort d’une partie du muscle cardiaque, le myocarde.

Pour prévenir ce risque, les médecins souhaitent étudier l’alimentation en sang du myocarde (muscle cardiaque) en comparant son afflux sanguin en situation de repos et en situation de stress.

On utilise pour celà un marqueur radioactif, injecté dans le sang du patient, et on suit son avancée et sa concentration au cours du temps dans les différentes parties du coeur.




On dispose d’un jeu de données de 4 individus réalisé par le CHU de La Tronche, en dimension 4 (3D et le temps).

Le marqueur radioactif est repérable sur ces images en nuances de gris par son fort rayonnement (cf image de gauche ci-dessous).

L’objectif principal est alors de faire de la classification non supervisée (clustering) afin de délimiter les différentes régions du cœur (cf image de droite ci-dessous) et de pouvoir en déduire de manière automatisée dans laquelle se trouve le marqueur à chaque instant. Les médecins pourront alors estimer si, pour un patient donné, la quantité du marqueur qui atteint le myocarde est suffisante ou non selon ses états de repos et de stress.

Av.png Ap.png
Image 2D du coeur avant / après segmentation


Objectif

L’algorithme de clustering qui nous a été fourni n'est pas satisfaisant car il regroupe des régions du cœur éloignées qui ne constituent pas un même tissu. On obtient alors un résultat biaisé de la quantité de marqueur radioactif présent dans le myocarde à un instant donné.

L’objectif de notre projet est donc l'amélioration de l'algorithme de découpe des régions de voxels (pixels 3D) pour définir clairement les trois parties du cœur qui nous intéressent : ventricule droit, ventricule gauche, myocarde.

Les signaux temporels correspondant à la quantité de marqueur radioactif dans les différentes régions du cœur devront également correspondre au circuit sanguin: entrée du marqueur dans le ventricule droit (en bleu), puis passage dans le gauche (en vert) et arrivée dans le myocarde petit à petit (en rouge).

Profil idéal des classes et courbes pour le ventricule droit, gauche et le myocarde


Démarche

Les bases théoriques sur lesquelles nous nous sommes appuyés regroupent à la fois des notions de Fouille de Données (avec l’algorithme de K-means pour le clustering), de Traitement d'Image (avec le repérage de zones connexes et le lissage 3D), de la Théorie de l'Information (avec l'information mutuelle et la divergence de Kullback) et des documents scientifiques plus développés pour la deuxième partie du projet (algorithmes de segmentation).

Au niveau de la démarche en elle même nous avons procédé en deux étapes : le clustering des données, puis leur segmentation.

Dans un premier temps nous avons utilisé un algorithme de clustering (basé sur la méthode K-means) qui nous était fourni et nous avons tenté d’améliorer cet algorithme en modifiant séparément les différents paramètres qui le caractérisent (le nombre de classes initiales, la façon dont elles sont choisies et surtout la notion de distance entre les voxels).

Une fois l'algorithme de clustering opérationnel, nous nous sommes intéressés à la segmentation des données, en vue de pouvoir déterminer par la suite de façon automatique la localisation du myocarde sur des clichés SPECT.



Le travail réalisé

Traitement des données pré-clustering

Filtrage temporel avec filtre médian

Avant de lancer le clustering, nous avons implémenté la possibilité de choisir un filtrage temporel à appliquer sur les valeurs prises par un voxel au cours du temps. Nous avons appliqué ce filtrage à l'ensemble des voxels qui composaient notre image de cœur 3D à l'aide de la fonction cellfun.

Le premier filtrage étudié pour cela a été le filtrage médian. Il nous a permis de lisser les signaux temporels et donc de réduire le bruit présent sur nos images 3D.

Résultats avant / après filtrage médian

On peut observer sur ces courbes la façon dont la quantité de liquide radioactif évolue dans chaque classe au cour du temps (avec ou sans lissage des valeurs). On observe également grâce aux coupes 2D du cœur que ce lissage médian enlève une bonne partie du bruit présent avant le lissage.

Filtrage temporel avec ondelettes

De la même façon que pour le filtre médian, nous avons implémenté la possibilité de pré-traiter les signaux temporels avec un filtrage basé sur les ondelettes. Nous avons pour cela utilisé une décomposition en ondelettes (avec wavedec) puis une recomposition (avec waverec) en ayant au préalable mis à zéros les coefficients haute fréquence (correspondant au bruit).

La première version des ondelettes que nous avons implémentée utilisait une notion de filtrage naïve (et simpliste) qui consistait à mettre tous les coefficients correspondant aux hautes fréquences dans la décomposition en ondelettes à 0.

Afin d'affiner le filtrage, nous utilisons désormais la notion de seuil, en ne supprimant que les coefficients inférieurs au seuil fixé. Nous testons deux méthodes différentes, basées sur les fonctions wnoisest (pour estimer la déviation standard du bruit à partir des coefficients de détails) et wbmpen (pour calculer le seuil de filtrage). Pour dé-bruiter le signal, nous appliquons soit la même méthode que précédemment (en éliminant seulement les coefficients de haute fréquence inférieurs au seuil), soit la fonction wdencmp.

On remarque que la méthode manuelle (qui filtre tous les coefficients de détail de la même manière) a tendance a réduire l'intensité du signal de façon beaucoup plus importante qu'avec la fonction wdencmp. Si ceci est très efficace pour réduire le bruit (on passe par exemple d'un pic d'intensité 8000 à un pic d'intensité 800 pour un bruit de haute fréquence, contre un pic de 6000 avec wdencmp), on remarque cependant que les pics d'intensité qui nous intéressent au niveau des ventricules voient leur importance diminuer fortement, jusqu'à être comparable à celle du bruit environnant. Par conséquent on préférera l'utilisation de wdencmp qui a une action moins efficace sur le bruit, mais qui garantit la conservation du corps du signal.

Le résultat de ce filtrage est présenté ci-dessous.

Résultats avant / après filtrage par ondelettes


Le clustering de données

Nous avons implémenté différentes méthodes pour l'utilisation de l'algorithme du K-means. Chaque méthode utilise une nouvelle notion de distance entre les classes. Nous testons ces méthodes sur les différents patients afin de les valider ou bien de les rejeter.

Les résultats sont présentés dans le pdf ci-dessous.

Fichier:Methodes-K-means.pdf

On peut remarquer que seules la distance euclidienne et la divergence de Kullback sur les données non normalisées fournissent des résultats exploitables. Nous avons donc combiné ces résultats avec des méthodes de traitement pré et post-clustering afin d'améliorer le rendu.

Traitement des données post-clustering

Nous allons ici traiter les données obtenues après le clustering. Les trois étapes suivantes sont appliquées dans l'ordre. Elles permettent le regroupement de classes appartenant au même tissu ainsi qu'une élimination du bruit.

Élimination des petites zones connexes

Dans un premier temps, les signaux sont filtrés avec un filtre médian de fenêtre 5, puis on applique un k-means sur 6 classes avec la distance euclidienne (distance la plus efficace pour l'algorithme de clustering). On obtient alors 50 images 2D qui représentent le résultat du clustering (ce sont les 50 'couches' du cœur).

Pour chacune de ces 50 images 2D, on regarde les zones connexes (voxels connexes de même classe). Si une zone connexe contient moins de 20 voxels, on fusionne la zone à une autre, peu importe laquelle, qui lui est connexe.

Image post-clustering avant / après élimination des petites zones

Élimination des zones connexes proches du bord

On applique les modifications précédentes. On refait ensuite le même type d’itérations pour éliminer les objets connexes de moins de 150 pixels qui touchent les bords car ceux-ci sont des tissus différents des ventricules ou du myocarde. On reclasse alors ces régions de pixels en la classe de cardinalité maximale.

Image post-clustering avant / après élimination des zones proches des bords


Fusion de classes

Dorénavant, on ne considère les signaux que sur des temps faibles (du temps 1 au temps ou la courbe du ventricule gauche atteint son maximum + 2). On normalise les signaux puis on les divisent par leur maximum respectif, on compare plus aisément l'allure des courbes comme cela. On prend ensuite les N signaux dont les maxima sont supérieurs à 80% du maximum moyen des signaux (80% étant une proportion obtenue empiriquement) et qui ne possèdent pas plus de 6000 voxels. On réalise alors N-2 fusions de classes dans le but de n'obtenir plus que 2 classes pour ces N signaux (correspondant dans l'idéal aux 2 ventricules). Ces fusions ont lieu entre les 2 régions dont la distance euclidienne est la plus faible, ainsi cet algorithme fusionnera des signaux moyens homologues au niveau des premiers pics et regroupera ainsi les classes appartenant au même ventricule.

Image post-clustering avant / après fusion des classes dont les courbes sont semblables


Nous supprimons également les très petites classes qui ne correspondent à aucune classe réelle et qui sont dues au bruit. Nous fusionnons donc les classes qui contiennent moins de 800 pixels (cf. courbe bleue foncée par exemple sur le graphe ci-dessous) avec la classe la plus proche.

Le seuil 800 pixels, empirique, a été pris de manière à éviter de fusionner les classes correspondant aux ventricules.



La validation

Atlas

Nous avons créé notre propre image 3D des trois régions du cœur qui nous intéressent avec des connaissances médicales. Cette image 3D a les même dimensions (70x70x50) et la même résolution (taille de pixel) que les images fournies par le CHU.

Nous en présentons 6 couches ci-dessous.

6 couches de l'atlas

Nous avons ensuite attribué aux voxels de chaque zone une fonction correspondant à l'évolution de la concentration du marqueur radioactif dans la zone à laquelle il appartient (tel que nous l'espérons).

Fonctions associées aux parties de l'atlas

Images "fantômes"

Afin de tester la performance et l'efficacité de nos algorithmes nous avons détérioré notre atlas avec différentes fonctions Matlab. Voici les résultats pour un bruit additif blanc gaussien, un bruit multiplicatif de type 'speckle' et un lissage gaussien.

Différentes détériorations d'une image 2D

Nous avons choisi d'appliquer un bruit de speckle puis un lissage gaussien afin d'avoir des images 3D proche de celles du CHU. Ces images 3D plus ou moins bruitées sont appelées 'fantômes'.

Résultats

Nous avons enfin appliqué notre programme (lissage pré-clustering, K-means avec distance euclidienne, élimination des petites zones et des zones proches des bords puis fusion des courbes des ventricules) aux images 3D ainsi créées. Le but est de détecter le seuil de détérioration à partir duquel notre méthode ne fournit plus des résultats acceptables.

On utilise ici le coefficient de Dice : D = \frac{|X \cap Y|}{|X|+|Y|}

X représente une région de l'atlas (ex: tous les voxels du myocarde de l'atlas) et Y la même région obtenue après segmentation du fantôme (tous les voxels du myocarde de l'image 3D en sortie de l'algorithme).

Voici la variation du coefficient de Dice en fonction de la détérioration appliquée à l'atlas (pour la distance euclidienne à gauche et la divergence de Kullback à droite).

Distance euclidienne Divergence de Kullback
Variation du coefficient de Dice en fonction de la détérioration appliquée

On remarque que la détection des deux ventricules est peu affectée par le bruit ou le lissage, cela est dû au nombre élevé de voxels dans ces région ainsi qu'aux forts pics en début de courbes qui sont peu perturbés par la détérioration.

Le myocarde quant à lui est plus sensible aux perturbations. On peut voir que pour σ supérieur à 1.5 le coefficient de Dice devient très mauvais. L'algorithme n'est plus efficace.

Conclusion

Ce projet s’est révélé très enrichissant dans la mesure où il a consisté en une approche concrète du métier de chercheur, une première dans notre parcours à l’Ensimag. En effet, la prise d’initiative, l’application de ses connaissances et le travail en équipe seront des aspects essentiels de notre futur métier.

De plus, il nous a permis d’appliquer nos connaissances en fouille de données, en théorie de l’information et en traitement d’image à un domaine pratique et d’actualité. En effet, le secteur de la médecine est aujourd’hui en quête de méthodes rapides et peu onéreuses pour détecter des pathologies d’ischémie chez les patients.

Les principaux problèmes, que nous avons rencontrés concernaient la collecte de données, accessibles par l’intermédiaire du CHU de Grenoble. L’interne du CHU responsable de la prise de mesure était en effet indisponible durant toute la période du projet. Nous nous sommes donc contenté de peu de données pour mener nos recherches.

Ainsi, nous avons touché du doigt la difficulté d’extrapoler le jeu de données disponibles pour obtenir des valeurs cohérentes sur l’ensemble du domaine étudié, ce que pourtant nous serons vraisemblablement amenés à faire dans notre futur métier.