FR2999325A1 - Procede iteratif de determination d'une image en deux dimensions ou trois dimensions a partir de signaux issus de tomographie par rayons x - Google Patents

Abstract

Procédé itératif de détermination d'une image en deux dimensions ou trois dimensions à partir de signaux issus de tomographie par rayons X, comprenant : - une étape de segmentation (1) de l'image par séparation de matériaux distincts respectivement par des lignes de niveaux pour une image en 2D ou par des surfaces de niveaux pour une image en 3D, et - une étape d'adaptation du maillage (2) respectivement triangulaire pour une image 2D ou tétraédrique pour une image 3D, dans laquelle on répartit des points de façon homogène respectivement sur une ligne de niveau pour une image 2D ou sur une surface de niveau pour une image 3D, servant à élaborer le maillage progressivement grossier respectivement triangulaire pour une image 2D ou tétraédrique pour une image à 3D, de sorte que le nombre de points répartis dépend respectivement, pour une ligne, du rapport entre la longueur de ladite ligne et la longueur de la ligne minimale parmi les lignes de l'image segmentée, ou, pour une surface, du rapport entre l'aire de ladite surface et l'aire de la surface minimale parmi les surfaces de l'image segmentée.

Description

Procédé itératif de détermination d'une image en deux dimensions ou trois dimensions à partir de signaux issus de tomographie par rayons X L'invention porte sur un procédé itératif de détermination d'images 5 en deux dimensions ou trois dimensions à partir de signaux issus de tomographie par rayons X. La reconstruction d'images par tomographie à partir de données de projection est très répandue dans le domaine de l'imagerie médicale. En utilisant un nombre suffisamment élevé de projection, les algorithmes à 10 rétroprojections filtrées, ou algorithmes FBP pour acronyme de "filtered backprojections" en langue anglaise, permettent d'obtenir des reconstructions d'images rapides et précises. Toutefois, dans le cas d'un faible nombre de vues (imagerie faible dose) et/ou d'angle limité (contraintes spécifiques liées à l'installation), les 15 données disponibles pour l'inversion ne sont pas complètes, le mauvais conditionnement du problème s'accentue, et les résultats montrent des artefacts importants. Pour aborder ces situations, une approche alternative consiste à discrétiser le problème de reconstruction, et à utiliser des algorithmes itératifs (ART pour acronyme de "Algebraic Reconstruction 20 Technique" en langue anglaise, SIRT pour acronyme de "Simultaneous Iterative Reconstruction Technique" en langue anglaise, et leurs variantes) ou une formulation statistique du problème afin de calculer une estimation de l'objet inconnu. Cette dernière permet d'incorporer une information a priori concernant l'objet et calculer le MAP (pour maximum a posteriori) ou des 25 solutions de vraisemblance régularisées. L'introduction d'a priori évite les divergences parfois observées dans les méthodes non régularisées (par exemple MLEM pour acronyme de "Maximum Likelihood ExpectationMinimization" en langue anglaise et ses variantes dont le nombre d'itérations doit être défini par l'utilisateur pour stopper le processus d'estimation), mais 30 l'utilisation de critères quadratiques montre l'inconvénient majeur de lisser les contours (gradients) de l'objet reconstruit. Des solutions utilisant des fonctions non quadratiques ont fait l'objet de nombreuses études ("Reconstruction d'images: régularisation avec prise en compte des discontinuités" par P. Charbonnier, Thesis, University of Nice-Sophia Antipolis, France 1994). Il a également été proposé d'introduire, dans un contexte, des variables binaires de Markov ("Stochastic relaxation, Gibbs 5 distributions and the Bayesian restoration of images" de Geman S, Geman D, IEEE Trans. lm. Froc. 1984;6(6):721-741) ou continues ("Constrained restoration and the recovery of discontinuities" de Geman D et Reynolds G, IEEE Trans, on Pattern Analysis and Machine Intelligence 1992, 14(3):367383 ; "Nonlinear image recovery with half-quadratic regularization and FFT's" 10 de Geman D et Yang C, IEEE Trans, lm Froc 1995, 4(7):932-946). Toutes ces méthodes sont classiquement basées sur une discrétisation du volume en un ensemble de voxels, et fournissent des cartes en trois dimensions de la densité de l'objet étudié. Les temps de calcul et la ressource de mémoire de ces méthodes itératives sont les principaux points 15 faibles qui rendent peu usuelles. Par ailleurs, quelle que soit l'application (médicale ou industrielle), l'approche commune consiste à appliquer des outils de traitement d'image sur les volumes reconstruits afin de segmenter les régions d'intérêt pour une analyse quantitative. Cependant, de nombreuses études ont été menées au cours des deux dernières décennies 20 conduisant à un large éventail d'outils de segmentation, basés sur différentes interprétations des contours et de fonctionnelles à minimiser. Les choix sont donc multiples et les résultats peuvent en dépendre. Dans ce contexte, des méthodes ont été proposées afin de mener la segmentation et la reconstruction de l'objet étudié de façon simultanée. 25 Des approches basées sur la forme, modélisent la forme et reconstruisent des objets en optimisant alternativement la position des frontières et les paramètres décrivant la distribution d'intensité. Plus récemment, des modèles non paramétriques, comme les contours actifs, développés et utilisés actuellement dans le domaine de la segmentation d'images, ont été 30 mis en oeuvre. 2 99932 5 3 Divers documents ("Polygonal and polyhedral contour reconstruction in comuted tomography" de Soussen C, Mohammad-Djafari A, IEEE Trans, Image Proc, 2004 November, 13(11)) ont mis l'accent sur la reconstruction d'une image binaire en trois dimensions ou 3D composée d'un objet uniforme compact totalement inclus dans un fond uniforme. La scène est reconstituée à l'aide de modèles de surfaces paramétriques (splines) dont les paramètres sont estimés à partir des données. Les résultats obtenus avec cette approche polyédrique restent fortement liés à l'initialisation, pour laquelle les auteurs proposent l'utilisation d'une méthode globale de surfaces harmoniques. D'autres documents ("Reconstruction 3D de vaisseaux à partir d'un faible nombre de projections à l'aide de contours déformables" de Senasli M, Garnero L, Herment A, Mousseaux E, GRETSI, 1997; "3D reconstruction of vessel lumens from very few angiograms by dynamic contours using a stochastic approach" de M. Senasli LGAHEM, Graphical Models, 2000, 62:105-127) ont proposé d'utiliser un contour actif basé sur un processus stochastique pour la reconstruction d'images tomographiques binaires dans le cadre de données réduites. Des fonctions splines sont utilisées pour décrire les contours des vaisseaux. Dans ce contexte des méthodes basées sur la forme de l'objet pour la reconstruction tomographique, il a été proposé ("Introducing shape priors in object-based tomography reconstruction" de Gaullier G, Charbonnier P, Heitz F, ICIP, 2009) d'introduire un a priori sur la forme, dont le descripteur basé sur les moments de Legendre est compact et hiérarchique. Cet a priori conduit à rendre la contrainte sur la longueur de la courbe appropriée à l'objet à reconstruire. D'autres auteurs proposent ("3D reconstruction of a patientspecific surface model of the proximal femur from calibrated x-ray radiographs: A validation study" de Zheng G, Schumann S, Med Phys, 2009, 36(4):1155-1166) d'utiliser un modèle de surfaces statistiques de l'objet et un modèle de recalage pour reconstruire un modèle de surfaces 3D spécifique au patient, à partir d'une paire de radiographies à rayons X 2D calibrées. La technique est basée sur une méthode itérative de recalage non rigide des formes caractéristiques extraites du modèle statistique de surfaces 3D avec les données radiographies. Les méthodes de segmentation par ensemble de niveaux (ou 35 "Ievel set" en langue anglaise) ont d'abord été introduites dans le domaine du 2 99932 5 4 traitement d'image pour capturer les fronts mobiles de surfaces et pour reconstruction de formes dans le cadre de problèmes inverses. Elles ont ensuite été appliquées pour estimer simultanément les valeurs dans les formes reconstruites, et ont montré des résultats prometteurs pour leur 5 utilisation dans les problèmes de reconstruction tomographique ("Reconstruction absorption and diffusion shape profiles in optical tomography by a level set technique" de Schweiger M, Arridge S, Dom n 0, Zacharopoulos A, kolehmainen V, Optics letters, 2006, 31(4):471-473; "3-D shape and contrast reconstruction in optical tomography with level sets" de 10 Schweiger M, Dom n 0, Arridge S, In: IPIA, 2008; "A curve evolution approach to object-based tomographic reconstruction" de Feng H, Karl W, Castafion D, IEEE tras lm Froc 2003, 12(1):44-57; "Tomographic reconstruction of piecewise smooth images" de Alvino C, Yezzi JA, In: CVPR, 2004; "Level set reconstruction for sparse angularity sampled data" de Yoon S, Pineda A, 15 Fahrig R, In: IEEE NSS Conf Rec, 2006; "level-set approach for the inversion and segmentation of x-ray tomography data" de Ramlau R, Ring W. A Mumford-Shah, J Comput Phys, 2007, 221(2):539-557; "Dynamic conebeam reconstruction using a variationnal level set formulation" de Keil A, Vogel J, Lauritsch G, Navab N, In: MICCAI, 2009). 20 Certains de ces documents ont appliqué la méthode des ensembles de niveaux pour reconstruire l'objet à partir de vues limitées bruitées par tomographie électromagnétique. Le document "Tomographic reconstruction using curve evolution" de Feng H, Castafion D, Karl W, IEEE, 2000, utilise un procédé d'ensembles de niveaux avec un petit nombre de coefficients de 25 texture pour reconstituer les images dans un contexte de vues limitées en tomographie. D'autres documents ("Level set reconstruction for sparse angularity sampled data" de Yoon S, Pineda A, Fahrig R, In: IEEE NSS Conf Rec, 2006; "Simultaneous segmentation and reconstruction: A level set 30 method approach for limited view computed tomography" de Yoon S, Pineda A, Fahrig R, Med. Phys. 2010 May, 37(5)) ont proposé un algorithme itératif pour la reconstruction tomographique qui segmente et reconstruit simultanément le domaine de l'image. La segmentation et la reconstruction sont réalisées en alternant une étape de segmentation qui fait évoluer une fonction à lignes de niveaux à deux phases et une étape de reconstruction 2 99932 5 5 qui calcule une estimation des valeurs d'intensité de l'objet inconnu à l'aide d'un algorithme de gradient conjugué. Le modèle à deux phases par morceaux a ses limites: il peut fournir des informations suffisantes lorsqu'on s'intéresse à un seul matériau (objet mono-matériau ou que les autres 5 structures de l'objet sont de moindre intérêt). Une extension de la méthode est nécessaire lorsque l'objet est composé de plusieurs matériaux (par exemple en ajoutant une deuxième fonction de lignes de niveaux). Une autre fonctionnelle d'énergie, telle l'énergie élastique en relation avec un détecteur de bord, et différents modèles géométriques ont 10 été envisagés à travers les courbes paramétriques comme les "snakes" en langue anglaise (littéralement traduit pas serpents paramétrés) ou les techniques d'ensembles de niveaux ("level set" en langue anglaise). En 2007, il a été proposé ("An Active Contour Approach for a Mumford-Shah Model in X-Ray Tomography" de Hoetzl E, Ring W, Advances in Visual 15 Computing, 2009, LNCS 5876:1021-1030) l'idée de segmenter et reconstruire simultanément à partir de données mesurées par tomographie X. Le travail est basé sur l'observation que les méthodes de segmentation échouent souvent lorsque l'image est très bruitée. Dans la méthode proposée, le contour segmenté et la densité correspondante sont estimés 20 par minimisation de la fonctionnelle de Mumford-Shah directement à partir des données mesurées. Toutes les méthodes précédemment citées utilisent une représentation standard de l'image qui consiste à diviser le volume de reconstruction en une grille régulière de pixels en deux dimensions (ou 25 voxels en trois dimensions) et supposer que le milieu reconstruit est constant à l'intérieur de chaque élément de la grille. Il a été proposé d'autres types de représentation afin d'améliorer la représentation commune. En 1992, il a été introduit ("Alternatives to voxels for image representation in iterative reconstruction algorithms" de Lewitt R, Phys Med 30 Biol 1992, 37:705-716) les fonctions de Bessel-Kaiser, également appelées blob, comme une représentation alternative d'une image. Il est décrit une famille d'éléments de volumes sphériques dont deux paramètres contrôlent le lissage de l'image (au contraire des pixels conventionnels qui sont discontinus). Par ailleurs, la symétrie de ces éléments sphériques conduit à 35 un calcul efficace des projections. Cette approche augmente la précision par 2 99932 5 6 rapport à l'approche standard de la représentation d'image basée sur des voxels. Dans le cadre de la Tomographie d'Emission Mono Photonique d'acronyme TEMP (ou SPECT pour acronyme de "Single-Photon Emission 5 Computed Tomography" en langue anglaise), il a été proposé ("Tomographic image reconstruction based on a content-adaptative mesh model" de Brankov J, Yang Y, Wernick M, IEEE Trans. Med. Imaging. 2004 February, 23(2)) une méthode de reconstruction d'images en 2D basée sur une discrétisation irrégulière. Le modèle de maillage est basé sur un 10 échantillonnage non uniforme dont la densité de noeuds s'appuie directement sur les détails de l'objet. La génération du maillage est basée sur une image de référence obtenue en utilisant une reconstruction FDK (pour FeldkampKress-Davis) à partir de laquelle les sommets ou noeuds du maillage sont placés à l'aide d'une étape d'extraction de carte de caractéristiques qui est 15 réalisée sur la grille régulière de pixels de l'image de référence. Le nombre de noeuds du maillage est déterminé par un critère de longueur de description minimale, ou MDL pour acronyme de "minimum description length" en langue anglaise, qui code une image de référence avec le nombre minimum de bits. Les algorithmes itératifs classiques comme la maximisation 20 de vraisemblance par espérance-minimisation d'acronyme MLEM pour "Maximum Likelihood Expectation-Minimization" en langue anglaise et le maximum a posteriori d'acronyme MAP pour maximum a posteriori en langue anglaise, ont été adaptés aux maillages triangulaires et utilisés pour la reconstruction d'image. Ces auteurs démontrent qu'un modèle de maillage 25 adaptatif de contenu de l'image fournit une discrétisation précise du problème avec moins d'inconnues qu'en utilisant des grilles de pixels réguliers, et peut donc améliorer le temps de calcul et la qualité d'image car la densité du maillage est adaptée aux détails de l'image. Il a été divulgué ("Tomographic reconstruction using 3D 30 deformable models" de Battle X, Cunnigham G, Hanson K, Phys. Med. Biol. 1998, 43:983-990) un procédé qui explore la déformation d'un maillage triangulaire par estimation directe des déplacements des sommets des triangles dans le cadre de l'approche Bayésienne. La méthode est appliquée sur des données SPECT. Les objets sont représentés par une surface 35 fermée et finement divisée en triangles renfermant une région uniforme 2 99932 5 7 (activité uniforme). Cette approche s'avère simple et efficace, mais soulève des questions complexes pour la reconstruction de plus d'une région. Il a également été proposé ("Three-Dimensional Attenuation Map Reconstruction Using Geometrical Models and Free-Form Deformations" de 5 Battle X, Le Rest C, Turzo A, Bizais Y, IEEE Trans. Med. Imaging, 2000, 19(5)) un autre modèle déformable qui rend indépendant les surfaces triangulées et la façon dont elle est déformée. Ils ont proposé une méthode à base de modèle itératif où la carte inconnue est modélisée comme un ensemble de régions géométriques ayant un coefficient d'atténuation 10 uniforme. La reconstruction consiste à déformer l'espace pour correspondre au mieux aux mesures. L'influence du choix du modèle de la matrice de projection et le temps très long de convergence de l'algorithme sont les deux principales faiblesses de l'approche. Il a été proposé ("Tomographic reconstruction using an adaptive 15 tetrahedral mesh defined by a point cloud" de Sitek A, Huesman R, Gullberg G, IEEE Trans. Med. Imaging., 2006 September, 25(9)) une méthode de reconstruction tomographique 3D utilisant une représentation du volume par un nuage de points et représenté par des tétraèdres sans chevauchement. La densité des noeuds définit la résolution locale de l'image, et l'intensité d'un 20 point quelconque à l'intérieur du volume est estimée par une interpolation linéaire entre les valeurs des quatre noeuds qui définissent chaque tétraèdre. Bien que le mode de calcul du système de matrice soit lourd et complexe, cette représentation de l'image est mieux adaptée pour la visualisation de volume à l'aide de processeurs graphiques connus, ou GPU pour acronyme 25 de "graphics processing units" en langue anglaise. Les principaux problèmes soulevés par les méthodes citées ci-dessus sont les suivants : - Le modèle de douceur par morceaux à deux phases est correct dans le cas d'un objet composé d'un seul matériau (ou dans le cas d'un objet 30 pour lequel certaines structures ne sont pas essentielles), mais la méthode atteint ses limites dès lors que l'objet étudié est composé de plusieurs matériaux. Les hypothèses ne sont en effet alors plus vérifiées, et les images reconstruites comportent des artéfacts. - La représentation de l'image usuellement utilisée (grilles régulières de 35 pixels) induit l'estimation d'un grand nombre d'inconnues ce qui engendre un coût de calcul et de mémoire de stockage conséquents, surtout dans le cas de reconstructions de hautes résolutions. Un but de l'invention est de palier aux problèmes précédemment 5 cités, en réduisant le nombre d'inconnues à estimer et en permettant l'étude d'objet composé de plusieurs matériaux ou régions d'intérêt. Aussi, il est proposé, selon un aspect de l'invention, un procédé itératif de détermination d'une image en deux dimensions ou trois dimensions à partir de signaux issus de tomographie par rayons X, comprenant trois 10 étapes qui alternent de façon itérative jusqu'à convergence: - une étape d'estimation de l'objet par un algorithme itératif de reconstruction adapté à une représentation de l'image sur une grille irrégulière de triangle en deux dimensions, de tétraèdres en trois dimensions, 15 - une étape de segmentation de l'image par séparation de matériaux distincts respectivement par des lignes de niveaux pour une image en deux dimensions ou par des surfaces de niveaux pour une image en trois dimensions, et - une étape d'adaptation du maillage respectivement triangulaire pour une 20 image à deux dimensions ou tétraédrique pour une image à trois dimensions, dans laquelle on répartit des points respectivement sur une ligne de niveau pour une image en deux dimensions ou sur une surface de niveau pour une image en trois dimensions, servant à élaborer le maillage respectivement triangulaire pour une image à deux dimensions 25 ou tétraédrique pour une image à trois dimensions, de sorte que le nombre de points répartis dépend respectivement, pour une ligne, du rapport entre la longueur de ladite ligne et la longueur de la ligne minimale parmi les lignes de l'image segmentée, ou, pour une surface, du rapport entre l'aire de ladite surface et l'aire de la surface minimale parmi les 30 surfaces de l'image segmentée. Un tel procédé permet de limiter drastiquement la quantité de calculs et donc le temps de calcul nécessaire à une détermination ou reconstruction d'images par tomographie par rayons X, et offre en résultat 35 une image classique de reconstruction tomographique en niveaux de gris (coefficients d'atténuation) associée à un maillage contenant la segmentation des différents matériaux ou régions d'intérêt de l'objet étudié. En outre, l'utilisation d'un maillage adapté au contenu de l'objet étudié permet de réduire le nombre d'inconnues, et de limiter la place 5 mémoire nécessaire. Dans un mode de réalisation, le nombre de points contraignant l'adaptation dudit maillage est respectivement proportionnel, pour une ligne (contour segmenté en deux dimensions), au rapport entre la longueur de ladite ligne et la longueur minimale desdites lignes (contours segmentés), ou, lo pour une surface, au rapport entre l'aire de ladite surface et l'aire minimale desdites surfaces. Ainsi, chaque contour, respectivement surface, est représenté par un nombre adéquat de points répartis de façon homogène sur chacun des dits contours, respectivement desdites surfaces, permettant de générer un 15 maillage fin contraint par ces dits points aux niveaux des contours respectivement desdites surfaces, et progressivement plus grossier en s'en éloignant. Selon un mode de réalisation, ladite étape de segmentation utilise une approximation de l'image à partir d'une fonction constante par morceaux. 20 Ainsi, chaque matériau ou région d'intérêt représenté par un coefficient d'atténuation donné peut être segmenté par une méthode d'ensemble de niveaux (level set). Dans un mode de réalisation, ladite approximation utilise la minimisation suivante : 1 C1C2,c min Lc + pt f (cl - 02 + II f (c2 - n+ 0 2 } n- 25 dans laquelle : - C représente une courbe pour une image en deux dimensions ou une surface pour une image en trois dimensions, qui représente la frontière entre deux régions constantes (C = afl) ; - c1 et c2 représentent respectivement les valeurs moyennes des régions 30 constantes par morceaux à l'intérieur fr et à l'extérieur fr de C ; - Lc représente un terme de régularisation permettant de contrôler la longueur de la courbe C ; et - pt. représente un terme de régularisation qui contrôle le détail de segmentation. Ainsi, par union de ces courbes en deux dimensions ou surfaces 5 en trois dimensions, on obtient l'ensemble des frontières sur lesquelles lesdits points de maillage vont être répartis. Selon un mode de réalisation, le procédé comprend, en outre, une étape d'initialisation du maillage, par un maillage de maille élémentaire 10 constante respectant le critère de Delaunay. L'invention sera mieux comprise à l'étude de quelques modes de réalisation décrits à titre d'exemples nullement limitatifs et illustrés par les dessins annexés sur lesquels : - les figures 1 et 2 illustrent schématiquement le procédé selon un 15 aspect de l'invention; - les figures 4a et 4b illustrent un exemple de mise en oeuvre du procédé en 2D mono-matériau, selon un aspect de l'invention; et - les figures 5a et 5b illustrent un exemple de mise en oeuvre du procédé en 2D multi-matériau, selon un aspect de l'invention 20 Sur l'ensemble des figures, les éléments ayants les mêmes références sont similaires. Sur la figure 1 est illustré schématiquement un procédé itératif de détermination d'une image en deux dimensions ou trois dimensions à partir 25 de signaux issus de tomographie par rayons X selon un aspect de l'invention, comprenant : - une étape 1 de segmentation de l'image par séparation de matériaux distincts respectivement par des lignes de niveaux pour une image en deux dimensions ou par des surfaces de niveaux pour une image en trois 30 dimensions, et - une étape 2 d'adaptation du maillage respectivement triangulaire pour une image à deux dimensions ou tétraédrique pour une image à trois dimensions, dans laquelle on répartit des points de façon équidistante respectivement sur une ligne de niveau pour une image en deux 35 dimensions ou sur une surface de niveau pour une image en trois dimensions, servant à élaborer le maillage contraint respectivement triangulaire pour une image à deux dimensions ou tétraédrique pour une image à trois dimensions, de sorte que le nombre de points répartis dépend respectivement, pour une ligne, du rapport entre la longueur de ladite ligne et la longueur de la ligne minimale parmi les lignes de l'image segmentée, ou, pour une surface, du rapport entre l'aire de ladite surface et l'aire de la surface minimale parmi les surfaces de l'image segmentée. L'utilisation d'un maillage adapté à l'objet étudié permet de réduire 10 le nombre d'inconnues, de mieux représenter l'objet tout en limitant la place mémoire de stockage nécessaire. La figure 2 illustre un exemple plus détaillé des étapes du procédé selon un aspect de l'invention. 15 Une étape d'initialisation 3 consiste à mailler le domaine de détermination ou reconstruction respectivement par des éléments de type triangles en deux dimensions ou 2D et par des éléments de type tétraèdres en trois dimensions ou 3D, en respectant les conditions de Delaunay contraint. 20 Le maillage initial peut être un maillage simple comportant Néléments et vérifiant la relation suivante : Ndet/100 Ndet / 1 < N < io , dans laquelle Ndet est le nombre de pixels détecteurs suivant la plus grande dimension du détecteur. Le maillage initial peut être, en variante, un objet maillé qui peut 25 être soit le résultat d'une reconstruction antérieure soit un fichier de type CAO. Une étape 4 effectue une reconstruction itérative de tomographie par rayons X utilisant un algorithme itératif usuel (algébrique ou probabiliste) mais adapté à la représentation par maillage non régulier composé de 30 triangles en 2D ou respectivement de tétraèdres en 3D. On considère le volume de reconstruction composé d'éléments triangulaires en reconstruction 2D ou d'éléments tétraédriques en reconstruction 3D, dont les valeurs inconnues associées correspondent à l'objet que l'on cherche à reconstruire, et des projections connues de l'objet 35 sur le détecteur correspondent aux mesures de chaque rayon de projection, et qui sont décrites sous la forme d'un vecteur de dimension M = Nproj x Ndet, Nproj étant le nombre d'angle de projections sous lesquelles l'objet est vu par le détecteur et Ndet le nombre de pixels détecteurs imageant chaque projection.

Le but est de déterminer l'objet à partir des projections acquises, chacune des valeurs détectées dans un pixel du détecteur étant une combinaison linéaire des valeurs des éléments de volume à reconstruire. En regroupant l'ensemble des mesures dans un vecteur p et en représentant l'objet à reconstruire par le vecteur f, on peut modéliser la 10 transformation permettant de passer de l'objet représenté par le vecteur f aux données observées p, par un modèle linéaire de la forme : p = Hf H représentant la matrice de projection (ou matrice système) dont les coefficients hii représentent la contribution de chaque élément du maillage à chaque rayon de projection. Dans le cas de la représentation proposée, ces 15 contributions correspondent à la longueur d'intersection entre un rayon de projection i et l'élément du maillage Ti (dans le cas 2D un élément de maillage correspond à un triangle, en 3D à un tétraèdre). La figure 3 est une illustration de cette contribution dans le cas 2D d'un maillage composé de triangles sur laquelle ith projection représente le 20 rayon de projection i, 1 représente ladite longueur d'intersection, et Ti représente l'élément de maillage. Pour résoudre ce problème, c'est à dire estimer f à partir de ses projections, il existe de nombreux algorithmes. Dans la présente invention, on exploite un algorithme itératif usuellement utilisé dans le cadre 25 d'une représentation conventionnelle respectivement pixélisée en 2D et voxélisée en 3D, dont l'opérateur de projection H est calculé en considérant dorénavant des éléments de maillage respectivement triangulaires en 2D et tétraédriques en 3D comme indiqué précédemment et illustré sur la figure 3 pour le cas 2D. 30 Dès lors que les éléments du maillage s'adaptent au contenu de l'image, les éléments de la matrice H sont recalculés. En 2D, l'élément de structure triangle est composé de trois sommets et d'une face à laquelle est affectée la valeur du niveau de gris ou coefficient d'atténuation estimé par l'algorithme de reconstruction de 35 tomographie.

En 3D, l'élément de structure tétraèdre est composé de quatre sommets et de trois faces auxquelles sont affectées la valeur de niveau de gris ou coefficient d'atténuation estimée par l'algorithme de reconstruction de tomographie.

L'étape 1 de segmentation est basée sur les résultats obtenus à l'étape précédente de reconstruction itérative 4. Elle consiste à délimiter les contours des différents matériaux composant l'objet étudié. Le résultat de cette étape est ensuite utilisé comme contrainte pour l'adaptation des éléments du maillage au contenu de l'image.

Il est pris pour hypothèse que chaque matériau composant l'objet étudié est homogène, et on utilise l'approche de segmentation basée sur le modèle de Mumford-Shah qui segmente une image en trouvant la meilleure approximation de l'image par une fonction constante par morceaux. Pour cela, on utilise les travaux de Chan et Vese [IEEE Trans. 15 lm. Proc. (2001)10(2) :266-277] qui proposent une formulation simple et très efficace de ce modèle à travers un modèle de contours actifs sans bords ou "active contours without edges" en langue anglaise. Il s'agit alors d'approximer l'image par une fonction prenant en compte deux valeurs au moyen du problème d'optimisation suivant : min F(ci, c2) = min Lc + ji f (cl_ _2 + 1.1 f (2 -D2 c1c2,C c1c2,C n- n+ 20 dans lequel : - C représente une courbe pour une image en deux dimensions ou une surface pour une image en trois dimensions, qui représente la frontière entre deux régions constantes (C = 25 - c1 et c2 représentent respectivement les valeurs moyennes des régions constantes par morceaux à l'intérieur fr et à l'extérieur fr de C, - Lc représente un terme de régularisation permettant de contrôler la longueur de la courbe C ; et 30 - pi. représente un terme de régularisation qui contrôle le détail de segmentation.

Pour résoudre le problème, on choisit une représentation par ensembles de niveaux ou "Ievel set" pour ensemble de niveau : C est remplacée par une fonction levet set (1)(x), fr et .(2+ sont définis par la fonction de Heaviside (g-C) Le critère à minimiser devient alors : E(c1c2, (p) v 1 1 VJ-C ( (x)) 1 dx + il f J-C ( (x)) (c 1 - f (x))2 + (1 - J-C ((x))) (C2 - f(x))2 dx 12 12 dans lequel .O. représente le domaine de segmentation, cp la fonction level set et y représente le terme de régularisation contrôlant les détails de segmentation. L'énergie E est alors minimisée en faisant alterner estimation de la courbe C en 2D et de la surface C en 3D, et mise à jour des valeurs c1 et c2. Ainsi, à (p(x) fixée, on peut écrire la minimisation de l'énergie E par rapport aux constantes cl et c2 comme : f ci(q)(x)) = f(x)g-r(q)(x))dx n fn 1 - g-r(q)(x))dx A cl et c2 fixées, la courbe optimale C en 2D et de la surface optimale C en 3D peut être estimée en utilisant le flux de gradient suivant qui décrit l'équation d'Euler-Lagrange associée à la fonction (p(x) et qui évolue 20 selon la direction normale à la courbe C: chp V'w 7= g-c'( ((p) y. ii((ci - 02 + (C2 - 02)) Afin d'éliminer toute difficulté liée aux modèles non convexes et de s'affranchir des minima locaux, on effectue une segmentation 25 convexe globale ou SCG basée sur l'observation qu'une solution stable du flux de gradient précédent coïncide avec l'état stable du flux simplifié : chp( v(p _ 7 - 17.1,1 ii((ci - 02 + (C2 - 02)) 10 15 Ce flux simplifié représente la descente de gradient pour minimiser l'énergie : E(q)) = IV(Pii + ii(q),r), Avec r = (c1 - n2 - (c2 - n2 Cette énergie n'a pas de minimum global unique (comme les courbes en 2D et surfaces en 3D n'ont pas de représentation unique par ensemble de niveau ou "Ievel set"). Afin de bien définir le minimum global, la solution doit être contrainte dans l'intervalle [0; 1] ce qui aboutit au problème d'optimisation suivant : min IV(Pli -Fii(P,r) o(.pi.

La minimisation est effectuée en utilisant l'approche de Split Brigman définie par Goldstein et al. [J. Sci. Comput. (2010)45 :272-293]. Après minimisation, le contour de la région segmentée est trouvé en seuillant la fonction d'ensemble de niveau de sorte que : C = {xlcp(x) = a} Avec x E et a E ] 0,1 [.

II: Dans le cas d'un objet mono-matériau, l'image à segmenter est composée de deux régions : le matériau et le fond. Le contour est alors max fcp(x)} obtenu en fixant a = . 2 Dans le cas d'un objet multi matériau, il est nécessaire d'affecter 25 différentes valeurs à a afin de segmenter les régions deux à deux. La segmentation finale S correspond à la fusion de l'ensemble des courbes obtenues avec les différentes valeurs de a : s = U ci iEI 30 Les valeurs de a sont fixées de telle sorte qu'elles permettent de segmenter les régions de faibles et de forts coefficients d'atténuation. Ces valeurs peuvent être fixées par l'utilisateur par sa connaissance a priori de l'objet à reconstruire, à l'aide d'histogramme de l'image en prenant les valeurs correspondant aux limites des régions constantes.

L'étape de segmentation faisant intervenir des calculs de gradient, il est nécessaire de passer à une représentation usuelle de l'image c'est-à-dire une grille régulière de pixels en 2D et de voxels en 3D. Le nombre d'éléments composant ce volume est : M = -8 . N lm Avec N = Ne représente le nombre de pixels d'une grille régulière usuelle G, de pixels en 2D ou respectivement de voxels en 3D dans le cadre de la reconstruction tomographique avec D E {2,3} la dimension du volume à reconstruire et Np = max (Aix, Ny ), N, et Ny étant le nombre de pixels détecteur selon les dimensions du détecteur ô représente la résolution spatiale de la grille régulière G, lm désigne la longueur moyenne des arêtes des éléments du maillage On associe alors à chaque élément de volume de la grille régulière la valeur de niveau de gris de la face de l'élément de volume de la grille irrégulière dont les coordonnées barycentriques minimisent la distance euclidienne entre les deux éléments de volume de chacune des deux grilles.

Durant l'étape 2, le maillage est adapté au contenu de l'image. L'adaptation du maillage au contenu de l'image utilise la segmentation obtenue à l'étape précédente. Cette segmentation peut être composée d'un contour ou de plusieurs contours. Afin que le nouveau maillage présente une maille fine (i.e. petits éléments de maillage en forte densité) au niveau des contours, devenant plus grossière en s'en éloignant (i.e. gros éléments de maillage en faible densité), il est nécessaire de répartir de façon homogène un nombre suffisant de points sur chaque contour segmenté. Ces points constitueront l'ensemble initial de points P à partir duquel le nouveau maillage est généré. L'ensemble P constitue une contrainte à la génération du nouveau maillage, en ce sens que les points le constituant sont fixes. Le critère proposé pour calculer le nombre optimal de points Ki pour chaque contour Ci est : - dans le cas 2D : L. K. = k , Lmin avec : minimale des longueurs de est fixée par l'utilisateur, Lci représente la longueur du contour Ci Lmin = mine.} représente la longueur tous les contours Ci k est une constante dont la valeur généralement fixée à 100 - dans le cas 3D: A. K. = k Amin avec : Aci représente la longueur du contour Ci Amin = minitAcil représente l'aire minimale des longueurs de tous les contours Ci k est une constante dont la valeur est fixée par l'utilisateur, généralement fixée à 100 Les Ki points de chaque contour sont ensuite répartis de façon homogène de telle sorte qu'ils soient équidistants les uns des autres. A partir de l'ensemble de points P = Uci{Pk}1,...,Ki ,on crée, par 25 exemple à l'aide l'outil "mesh generation" de la bibliothèque CGAL, les éléments du nouveau maillage en considérant deux critères: leur forme et leur taille afin de générer une triangulation de Delaunay contrainte. Le critère de forme est pris en compte par le critère de qualité B= - , r représentant le rayon du cercle circonscrit à l'élément de maillage 30 considéré et / la longueur de la plus courte arrête de ce même élément de maillage considéré.

Il existe une relation entre B et le plus petit angle lemin de l'élément de maillage considéré de l'élément de maillage : 1 B = = 1 2. Sin(emin) de laquelle on déduit la relation suivante : 1\ flmin = arcs in (-2B La convergence de l'algorithme est assurée pour B -\12, ce qui correspond à flmin = 20,7°.

Le critère de taille est contrôlé par le paramètre lm' qui définit la longueur maximale des arrêtes de l'élément de maillage. Le choix de lm' joue un rôle fondamental dans la création du maillage. Ainsi, pour que les mailles soit générées de façon à respecter les contraintes et permettre un maillage progressivement plus grossier en s'éloignant des contours, la valeur de lm' doit satisfaire la condition suivante : Lc lm'>>-Kc où Lc est la longueur du contour C, Kc le nombre des points décrivant le contour C. La longueur du contour est données par : Kc Lc où pi e P,i = {0, ...,Kc -1} i=0 Dans le cas où plusieurs contours sont décrits par l'ensemble P, la valeur optimale lopt est : Lc . lopt>> 'mn Kcmm O ù L. est la longueur du contour le plus court Cmin, K. mtn e nombre de points décrivant ce plus petit contour Cmin.

Les figures 4a et 4b illustrent une mise en oeuvre du procédé pour une image à deux dimensions à partir de données expérimentales avec la méthode proposée. Il s'agit ici d'une coupe transversale d'un objet mono matériau qui est un foret en acier de 150 pm de diamètre. La figure 4a représente l'image en niveaux de gris, et la figure 4b illustre la grille irrégulière correspondante composée d'une maille triangulaire adaptée au contenu de l'image faisant apparaître le contour final de l'objet segmenté. Les figures 5a et 5b illustrent une mise en oeuvre du procédé pour une image à deux dimensions à partir de données numériques avec la méthode proposée. Il s'agit ici d'une coupe transversale d'un fantôme numérique, ou en d'autres termes d'un objet défini numériquement, appelé Shepp-Logan, multi-matériaux. La figure 5a représente l'image reconstruite en niveaux de gris, et la figure 5b illustre la grille correspondante composée d'une maille triangulaire adaptée au contenu de l'image (i.e. aux matériaux la composant). Les étapes du procédé décrit ci-dessus peuvent être réalisées par un ou plusieurs processeurs programmables exécutant un programme informatique pour réaliser les fonctions de l'invention en agissant sur des 20 données d'entrée et en générant des données de sortie. Un programme informatique peut être écrit dans n'importe quel langage de programmation, tels les langages compilés ou interprétés, et le programme informatique peut être déployé sous n'importe quelle forme, y compris en tant que programme autonome ou comme un sous-programme 25 ou fonction, ou tout autre forme appropriée pour une utilisation dans un environnement informatique. Un programme d'ordinateur peut être déployé pour être exécuté sur un ordinateur ou sur plusieurs ordinateurs sur un seul site ou sur plusieurs sites répartis et reliés entre eux par un réseau de communication. 30

Claims (5)

1. REVENDICATIONS1. Procédé itératif de détermination d'une image en deux dimensions ou trois dimensions à partir de signaux issus de tomographie par rayons X, comprenant : - une étape de segmentation (1) de l'image par séparation de matériaux distincts respectivement par des lignes de niveaux pour une image en deux dimensions ou par des surfaces de niveaux pour une image en trois dimensions, et - une étape d'adaptation du maillage (2) respectivement triangulaire pour une image à deux dimensions ou tétraédrique pour une image à trois dimensions, dans laquelle on répartit des points respectivement sur une ligne de niveau pour une image en deux dimensions ou sur une surface de niveau pour une image en trois dimensions, servant à élaborer le maillage respectivement triangulaire pour une image à deux dimensions ou tétraédrique pour une image à trois dimensions, de sorte que le nombre de points répartis dépend respectivement, pour une ligne, du rapport entre la longueur de ladite ligne et la longueur de la ligne minimale parmi les lignes de l'image segmentée, ou, pour une surface, du rapport entre l'aire de ladite surface et l'aire de la surface minimale parmi les surfaces de l'image segmentée.
2. 2. Procédé selon la revendication 1, dans lequel ledit nombre de points répartis pour générer ledit maillage adaptatif est respectivement proportionnel, pour une ligne, au rapport entre la longueur de ladite ligne et la longueur de ladite ligne minimale, ou, pour une surface, au rapport entre l'aire de ladite surface et l'aire de ladite surface minimale.
3. 3. Procédé selon la revendication 1 ou 2, dans lequel ladite étape de segmentation utilise une approximation de l'image à partir d'une fonction constante par morceaux.
4. 4. Procédé selon la revendication 3, dans lequel ladite approximation utilise la minimisation suivante :min fl,c + pi f (ci_ - 02 + I-1- f (c2 - Di cic2,c n- n+ dans laquelle : - C représente une courbe pour une image en deux dimensions ou une surface pour une image en trois dimensions, qui représente la frontière entre deux régions constantes (C = - c1 et c2 représentent respectivement les valeurs moyennes des régions constantes par morceaux à l'intérieur fr et à l'extérieur fr de C, - Lc représente un terme de régularisation permettant de contrôler la longueur de la courbe C ; et - pt. représente un terme de régularisation qui contrôle le détail de segmentation.
5. 5. Procédé selon l'une des revendications précédentes, comprenant, en outre, une étape d'initialisation du maillage, par un maillage de maille élémentaire constante respectant le critère de Delaunay.20