Jean-Baptiste Keck : Reconstruction efficace d'un volume 3D isotrope à partir d'images ultrason 2D localisées

De Ensiwiki
Aller à : navigation, rechercher


Reconstruction efficace d'un volume 3D isotrope à partir d'images ultrason 2D localisées

Labo TIMC-IMAG [1]
Equipe GMCAO [2]
Encadrants Matthieu.Chabanas@imag.fr


Etudiant

KECK Jean-Baptiste (2A - Modélisation Mathématique, Image et Simulation)

Contexte de recherche

Une sonde échographique miniaturisée a été développée pour acquérir des images intra-articulaires du cartilage du genou, sous conditions de chirurgie arthroscopique. L'objectif est de détecter des lésions cartilagineuses (arthrose) et de quantifier leur taille et localisation précise, pour aider le clinicien dans son diagnostic et sa prise de décision thérapeutique.

La sonde ultrasons est traquée par un localisateur optique, ce qui permet de localiser (position + orientation) chaque image 2D dans un référentiel 3D commun. De 2000 à 5000 images sont acquises lors du balayage d'une surface cartilagineuse. Un outil intéressant est de reconstruire un volume 3D isotrope (voxels cubiques) de l'espace étudié à partir de l'ensemble des images.

Mon but était de permettre une reconstruction rapide ce volume, dans l'optique d'une application clinique, à l'aide du GPGPU (General Purpose Computing on Graphic Processing Units). Le GPGPU est le fait d'effectuer un calcul générique sur un processeur graphique afin de bénéficier de leur capacité de traitement massivement parallèle. Le GPGPU suit par essence le paradigme SIMD (Single Instruction Multiple Data), la même instruction est appliquée simultanément à plusieurs données pour produire plusieurs résultats. Ici le choix a été fait d'utiliser CUDA mais tout est transposable en OpenCL.

Images US
Exemple d'images ultrasons localisées par rapport à la surface cartilagineuse
Setup
Interfaces eau/cartilage et cartilage/os sur une image

Reconstruction du volume

Traitements préliminaires

Il faut tout d'abord filtrer les images, ainsi que leur positions successives. L'on peut également recadrer toutes les images sur une région d’intérêt ou effectuer des conversions de données. Dans cet IRL, je n'ai qu'effectué la partie conversion de données afin de réduire l'empreinte mémoire des images dans la mémoire des GPUs, qui est souvent plus limitée que la mémoire dédiée du CPU (RAM).

Détermination de la grille

Le volume a reconstruire est un ensemble discret de voxels cubiques, il faut donc tout d'abord construire une grille de voxel. Pour déterminer cette grille, il faut tout d'abord trouver la boîte englobante de toutes nos images ultrasons. Une fois cela effectué, il faut choisir une taille de voxel, typiquement de 0.1 à 0.2mm.

Mais cela ne suffit pas. En effet la grille de voxel est l'élément qui va consommer le plus de mémoire dans l'algorithme de reconstruction, d'autant plus qu'il faut en stocker plusieurs en mémoire. La solution est de découper la grille principale en sous-grilles qui seront remplies indépendamment les unes des autres. Pour simplifier le code, on prend la plus petite grille de taille alignée sur des puissance de deux qui contient la grille initiale, puis on la découpe par une puissance de deux sur chacun des axes. Toutes les sous-grilles résultantes ont ainsi la même taille et satisfont des contraintes d'alignement intéressantes pour une application GPGPU.


Détermination de la boîte englobante, image extraite et modifiée de T. Wen et al. [1]
Setup
Choix de la taille des voxels cubiques
Setup
Sur-grille de taille alignée sur une puissance de deux
Setup
Découpe de la grille en sous-grilles

La découpe en sous-grilles ne permet pas seulement de s'affranchir des contraintes de consommation mémoire, mais permet aussi la mise en place aisée d'une solution multi-GPU, d’équilibrage de la charge et des transferts de mémoire entre CPU et GPU.

Les transferts de mémoire entre CPU et GPU sont les transferts de mémoire les plus lents qui interviennent en GPGPU. Heureusement, CUDA fournit un mécanisme pour faire travailler le processeur graphique tout en effectuant un transfert de mémoire, en passant par des "streams".

Setup
Exemple de streams sur deux GPUs différents. En GPGPU, le CPU maître est appelé le host, et les carte graphiques esclaves sont appelées les devices. Les fonctions exécutées sur les processeurs graphiques sont appelés kernels.

On voit que plus l'on va découper les données à transférer, plus l'on va pouvoir se faire chevaucher les kernels (en bleu) et les transferts de mémoire (en vert et jaune).

Algorithmes de reconstruction

Il existe plusieurs algorithmes de reconstruction qui partent d'images 2D localisées, les plus connus étant le PNN (Pixel Nearest Neighbor) et le VNN (Voxel Nearest Neighbor).

  • Dans la méthode du PNN, on va itérer sur chaque voxel de la grille, et pour chaque voxel, on va déterminer le pixel qui lui est le plus proche. Cet algorithme est de coût élevé mais il à l'avantage de ne pas laisser de voxels vides, et donc ne nécessite pas l'application d'un algorithme d'interpolation.
  • La méthode VNN quant-à elle fonctionne de la manière inverse : Chaque pixel est parcouru est est affecté au voxel qui lui est le plus proche. La complexité est grandement abaissée, mais cet algorithme laisse des trous ce qui implique une deuxième étape d'interpolation des voxels vides grâce à leurs voxels voisins.
  • Il existe d'autres méthodes plus complexes comme les algorithmes de type Distance Weighted (DW) qui prennent en compte un voisinage de pixels proches avec une pondération de distance, et d'autres qui prennent en compte la direction des contours lors de la phase d'interpolation comme la Fast Marching Methode (FMM).

Pour plus d'informations et des comparatifs sur ces méthodes, l'on pourra se référer à R. Rohling et al. [2] et T. Wen et al. [1] pour la méthode FMM.

Implémentation

Comme l'algorithme PNN présentait une complexité moindre et se prêtait bien au paradigme SIMD, c'est par lui que j'ai choisi de commencer. Comme dit précédemment, cet algorithme se décompose en deux partie : Une partie dite de 'hole filling' où l'on utilise l'information présente dans les images pour remplir notre grille de voxel et une partie dite de 'bin filling' où l'on va combler les trous laissés par la première étape.

Tout a été implémenté en C++ et CUDA. Une interface graphique a été développée avec Qt4 afin de mieux visualiser le volume et les coupes générées. L'afficheur de voxel a été réalisé en OpenGL avec le QGLviewer et CUDA.

Bin-filling

Tous les pixels sont traversés et affectés au voxel le plus proche. On effectue un moyennage si il y a redondance d'information (plusieurs pixels qui vont dans le même voxel).

Voici quelques résultats que j'ai obtenu sur le jeu de données d'un demi plateau tibial (environ 1600 images) :

Exemple d'une image seule reconstruite par la première étape de l'algorithme. Le découpage de la grille a été exagéré.
Exemple de volume reconstruit. Les trous ne sont pas visibles mais ils sont bien présents dans le volume.
Coupe dans ce même volume, on voit bien les voxels vides laissés par cette première passe (en noir).

Au niveau de la carte graphique, on alloue 4 grilles de voxels (une pour sommer les contributions, une pour compter les contributions, et deux pour stocker le résultat de la moyenne). Il y a deux grilles de résultats pour appliquer la stratégie des "streams", en effectuant des calculs tout en transférant l'autre sous-grille précédemment calculée vers le CPU. On garde également de la place pour y placer les données des images et des transformations.

Organisation mémoire pour le CPU et un GPU quelconque. Les échanges mémoires sont représentés par les flèches. La grille représentée en jaune est la grille où l'on stocke la moyenne résultante, elle est allouée deux fois.


Hole-filling

L'étape précédente laisse des voxels vides. On effectue une interpolation de ceux-ci avec leur voisins. La région de voisinage choisie est une simple sphère, et l'on pondère les contributions avec la distance au centre de cette sphère.

Toujours pour le même jeu de données (demi plateau tibial, environ 1600 images), j'ai obtenu les résultats suivant :

Vue 2D du pattern d'interpolation sphérique avec R=2. Les carrés noirs représentent les voxels vides. Le carré rouge situe le voxel vide qui est entrain d'être interpolé.
Exemple de volume reconstruit. L'interpolation donne une épaisseur plus importante aux images qui sont distantes, ce qui est visible notamment sur le haut du volume.
Coupe dans ce même volume, il y a beaucoup moins de voxels vides (cela dépend du rayon d'interpolation). L’interpolation des bords n'a pas été implémenté par manque de temps d’où la non reconstruction bien visible des bords des sous-grilles.

Au niveau du GPU, l'organisation de la mémoire est plus simple car il ne reste plus que 3 grilles, la première étant la grille source, et les deux autres des grilles pour stocker le résultats, toujours pour appliquer la stratégie des "streams".

Organisation mémoire pour le CPU et un GPU quelconque. Les échanges mémoires sont représentés par les flèches. La grille source est la grille de l'étape précédente, la grille destination est la grille sans les trous, elle est allouée deux fois.

Pour tirer au maximum partie du matériel graphique, on utilise la mémoire partagée pour accélérer grandement cette étape. De plus cela permet d'atténuer l'impact du rayon d'interpolation sur les performances (qui dépendraient du cube du rayon en temps normal).

Résultats

Ces résultats ont été réalisés dans un environnement multi-GPU (deux GTX 770 comprenant 4Go de mémoire embarquée et 1536 cœurs de calculs chacune). Le paramètre d représente la taille des voxels.

Nombre de sous-grilles Rayon d'interpolation d = 0.5mm (128x128x256) d = 0.2mm (256x512x512) d = 0.1mm (512x1024x1024)
4 0 1.57s 1.62s 2.64s
4 1 1.89s 10.8s 88.4s
4 2 1.96s 12.5s 109.5s
4 3 2.36s 16.5s 146.2s
16 3 6.24s 15.1s 129.9s
64 3 2.98s 14.8s 105.9s
256 3 6.22s 15.4s 87.9s
256 4 6.4s 16.4s 112.3s
256 5 6.52s 19.2s 153.4s
256 6 6.59s 22.9s 205.3s

Le programme a également été testé sur une carte graphique qui possédait moins de mémoire (une GTX 590, 2x1.5Go) et la reconstruction était pleinement fonctionnelle.

On voit que les temps de reconstruction restent globalement acceptable pour une application clinique (moins de 5 minutes dans tous les cas). L'étape de bin-filling (R=0) est clairement l'étape la plus rapide (moins de 3s dans tous les cas) car il y a beaucoup moins de pixels que de voxels. Le gros du travail est fait dans le hole-filling (interpolation) qui représente plus de 90% du temps d’exécution dans quasiment tous les cas. A chaque fois que l'on double la précision de la grille, on multiplie par 8 le temps d’exécution, ce qui n'est pas surprenant car à taille de grille constante, le nombre de voxels est multiplié par 8. L'on voit également que l'utilisation de la mémoire partagée limite l'effet du rayon d'interpolation sur le coût comme prédit, ce qui est une bonne nouvelle car la dépendance serait cubique ! Enfin, on constate qu'un découpage plus important de la grille permet d'obtenir des meilleures performances car la charge est mieux répartie entre le différentes cartes graphiques, et les transferts de mémoire sont mieux recouvert par les kernels.

Bilan et conclusion

  • L'objectif de reconstruction du volume en un temps raisonnable a été atteint : Une application clinique pourra être mis en œuvre, même sur des cartes graphiques non professionnelles présentes dans les ordinateurs personnels.
  • L'algorithme scale très bien sur plusieurs cartes graphiques.
  • J'aurais aimé avoir plus le temps pour m’intéresser un peu plus au filtrage des images et des positions.
  • Il aurait également fallu coder une implémentation de référence en C++ à des fins comparatives.
  • Ceci ouvre des porte à l'adaptation d'algorithmes de reconstruction plus robustes (comme le FMM qui lui aussi semble bien s'adapter au paradigme SIMD)

Références

  • [1] T. Wen, Q. Zhu, W. Qin, L. Li, F. Yang, Y. Xie, and J. Gu, “An accurate and effective fmm-based approach for freehand 3d ultrasound reconstruction,” Biomedical Signal Processing and Control, vol. 8, no. 6, pp. 645–656, 2013
  • [2] R. Rohling, A. Geel, and L. Berman, “A comparison of freehand three-dimensional ultrasound reconstruction techniques,” Medical Image Analysis, vol. 3, no. 4, pp. 339–359, dec 1999.

Fichiers

Rapport : Fichier:Rapport irl keckj.pdf

Slides de la soutenance : Fichier:Slides irl keckj.pdf