FR2830962A1

FR2830962A1 - Dispositif et methode de traitement d'image pour detection de lesions evolutives

Info

Publication number: FR2830962A1
Application number: FR0115780A
Authority: FR
Inventors: Nicolas Ayache; Gregoire Malandain; David Rey; Jonathan Stoeckel
Original assignee: Institut National de Recherche en Informatique et en Automatique INRIA
Current assignee: Institut National de Recherche en Informatique et en Automatique INRIA
Priority date: 2001-10-12
Filing date: 2001-12-06
Publication date: 2003-04-18
Anticipated expiration: 2021-12-06
Also published as: US7822240B2; EP1440413A1; FR2830962B1; US20050141757A1; WO2003044719A1

Abstract

L'invention concerne un dispositif de traitement d'images, ayant une entrée pour recevoir une suite temporelle de jeux de données représentant des images numériques volumiques comparables, chaque donnée ayant une composante de position et une composante d'intensité, des moyens de pré-traitement (18) pour modifier les jeux de données afin d'obtenir des images recalées géométriquement et en intensité, et des moyens de traitement comparatif (20) pour examiner des jeux de suites temporelles d'éléments d'image et d'y détecter des signes de variations. Les moyens de traitement (20) comprennent une fonction de modélisation (12) pour ajuster un modèle paramétrique séparément sur certains des jeux de suites temporelles d'éléments d'image, pour obtenir des couples d'informations (élément d'image, temps) et une fonction d'analyse statistique (14) de ces couples d'informations pour isoler les couples d'informations représentatifs d'une variation significative.

Description

Dispositif et méthode de traitement d'image pour détection de lésions évolutives L'invention concerne le traitement d'images, notamment médicales.

Le traitement d'images permet notamment de comparer des images entre-elles. Selon certaines techniques actuelles, les images sont comparées deux à deux. D'autres techniques permettent de comparer les images d'une série. Selon le contexte, différentes catégories de problèmes se posent pour rendre les images comparables, comme on le verra.

Ces comparaisons ont des applications très diverses, par exemple dans le domaine médical pour l'aide au diagnostic de pathologies. Dans le cas de la comparaison des images d'une série temporelle prises sur une période donnée, il est possible de détecter des zones à évolution pathologique.

Dans une application particulière du traitement d'images, par exemple dans le cas de lésions de sclérose en plaques, des images médicales peuvent être analysées temporellement et automatiquement de façon à établir une détection de zones à évolution pathologique. Des techniques d'analyses temporelles automatiques sont particulièrement utilisées pour fournir une aide au diagnostic d'un patient. La performance des techniques d'analyses temporelles et automatiques actuellement utilisées est limitée, notamment dans le cas de la détection des zones à évolution pathologique. Dans le domaine médical, ces techniques sont mal adaptées par exemple pour quantifier a posteriori l'effet de médicaments administrés dans le cas d'essais thérapeutiques. Plus généralement, il s'agit de détecter d'une image à une autre des évolutions fines dans des zones à évolution pathologique. Il s'agit de pouvoir le faire sur des prises de vues espacées dans le temps, et également

sur un phénomène en cours, c'est-à-dire dont le début n'a pas été observé.

De façon générale, la présente invention vient améliorer la situation.

L'invention concerne un dispositif de traitement d'images, comprenant : - une entrée pour recevoir une suite temporelle de jeux de données représentant des images numériques volumiques comparables, chaque donnée comprenant une composante de position et une composante d'intensité, - des moyens de pré-traitement pour modifier les jeux de données de façon à obtenir des images recalées géométriquement et en intensité, et - des moyens de traitement comparatif, propres à examiner des jeux de suites temporelles d'éléments d'image, afin d'y détecter des signes de variations.

Selon une caractéristique principale de l'invention, les moyens de traitement comprennent : * une fonction de modélisation, capable d'ajuster un modèle paramétrique séparément sur certains au moins des jeux de suites temporelles d'éléments d'image, ce qui fournit des couples d'informations du genre : élément d'image, temps ; et * une fonction d'analyse statistique de ces couples d'informations, propre à isoler ceux des couples d'informations qui sont représentatifs d'une variation significative.

L'invention concerne un procédé de traitement d'images comprenant les étapes suivantes : a-recevoir une suite temporelle de jeux de données représentant des images numériques volumiques comparables, chaque

donnée comprenant une composante de position et une composante d'intensité, b-modifier les jeux de données de façon à obtenir des images recalées géométriquement et en intensité, Selon une caractéristique principale de l'invention, le procédé comprend les étapes complémentaires c-ajuster un modèle paramétrique séparément sur certains au moins des jeux de suites temporelles d'éléments d'image, ce qui fournit des couples d'informations du genre : élément d'image, temps, d-isoler ceux des couples d'informations qui sont représentatifs d'une variation significative par analyse statistique.

D'autres caractéristiques et avantages de l'invention apparaîtront à l'examen de la description détaillée ci-après, ainsi que des dessins annexés sur lesquels : - la figure la illustre un dispositif informatique comprenant un dispositif de traitement d'image selon la technique antérieure, - la figure 1b illustre un dispositif informatique comprenant un dispositif de traitement d'image selon l'invention,

- la figure 2a représente un ordinogramme du procédé automatique d'analyse d'images, - la figure 2b représente une partie détaillée de l'ordinogramme de la figure 2a, - la figure 2c représente une partie détaillée de l'ordinogramme de la figure 2a selon un mode de réalisation particulier,

- la figure 3a est une image de cerveau représentant une alternance de cubes d'une première image 3D et d'une deuxième image 3D ayant un biais temporel avec la première image, - la figure 3b est une image de cerveau représentant une alternance de cubes de la première image 3D de la figure 3a et de la deuxième image 3D corrigée, - la figure 3c représente l'histogramme lié à la figure 3a, cet histogramme représentant les intensités de la première image par rapport aux intensités de la deuxième image, - la figure 3d représente l'histogramme lié à la figure 3b, cet histogramme représentant les intensités de la première image par rapport aux intensités de la deuxième image corrigée, - la figure 4a est une image de cerveau représentant une première image 3D, - la figure 4b est une image de cerveau représentant une deuxième image 3D ayant un décalage spatial et un biais temporel avec la première image, - la figure 4c est une image de cerveau représentant la deuxième image 3D de la figure 4b dont le décalage spatial a été corrigé, - la figure 4d est une image de cerveau représentant la deuxième image 3D de la figure 4c dont le biais temporel a été corrigé, - les figures 5a, 5b, 5c, 5d, 5e, 5f sont des profils d'intensité I évoluant au cours du temps T de six points de cerveau sélectionnés,

- la figure 5g est un exemple du modèle paramétrique avec certains paramètres, - la figure 5h est un modèle moyen de profil d'intensité calculé selon entre autres les profils d'intensité des figures 5a à 5f, - la figure 5i est le modèle moyen de profil d'intensité de la figure 5h superposé aux mesures discrètes normalisées, - la figure 5j est une discrétisation du modèle moyen de profil d'intensité de la figure 5i, - les figures 6a à 6f représentent des profils d'intensité I sur lesquels les temps de montée sont calculés, - la figure 7 est une représentation schématique du modèle général linéaire, - la figure 8 représente un graphe de distribution des tailles de groupes de voxels selon le mode de réalisation particulier.

De plus, la description détaillée est accompagnée de l'annexe I, dans laquelle est présenté un exemple de modèle paramétrique, un exemple du modèle linéaire et le calcul de la variable de Student.

Les dessins contiennent, pour l'essentiel, des éléments de caractère certain. Ils pourront donc non seulement servir à mieux faire comprendre la description, mais aussi contribuer à la définition de l'invention, le cas échéant.

Le demandeur a déjà formulé des propositions dans EP 0927 405, EP 0 977 148 et EP 0 564 318.

Différentes publications ont également été produites par le demandeur : "New feature Points Based on Geometric Invariants for 3D Image Registration", Jean-Philippe Thirion, International Journal of Computer Vision, 18 (2) : 121-137.

"Recalage d'images médicales par inférence statistique", Alexis Roche, Thèse de Doctorat, Université de Nice-Antipolis, janvier 2001.

"Etude de la symétrie bilatérale en imagerie cérébrale volumique", Sylvain Prima, Thèse de Doctorat, Université de Paris XI, mars 2001.

"Automatic Quantification of changes in the Volume of Brain Structures", Guillaume Calmon, Neil Roberts, Paul Eldridge and Jean-Philippe Thirion, First International Conference on Medical Image Computing And Computer-Assisted InterventionMICCAI'98, Lecture Notes in Computer Science 1496 (1998), 964- 973.

"Automatic Detection and Segmentation of Evolving Processes in 3D Medical Images : Application to Multiple Sclerosis", David Rey, Gérard Subsol, Hervé Delingette et Nicholas Ayache, Information Processing in Medical Imaging-IPMI'99, Lectures Notes in Computer Science 1613 (1999) 154-167.

"Using Continuum Mechanics Operators for Detection and Quantification of Evolving Processes in 3D Medicals Images", David Rey, Hervé Delingette, Gérard Subsol, Nicholas Ayache, EUROMECH'99,1999, 185-188.

La performance des techniques d'analyses temporelles et automatiques actuellement utilisées est limitée, notamment dans le cas de la détection des zones à évolution

pathologique. Dans le domaine médical, ces techniques ne permettent pas, a posteriori, de quantifier rigoureusement l'effet de médicaments administrés dans le cas d'essais thérapeutiques.

De façon générale, les traitements d'images en question nécessitent un pré-traitement. Pour permettre une analyse temporelle et automatique d'images, ces techniques de prétraitement d'images sont au préalable appliquées aux images à analyser. Ces techniques, présentées ci-après, ont pour objectif de rendre comparables les images entre-elles et de permettre l'analyse de ces images.

Le rapprochement d'images nécessite en général deux opérations de pré-traitement. Ces opérations permettent une normalisation (ou une calibration) des images afin de les comparer. Dans l'imagerie médicale, ces deux opérations de pré-traitement s'appliquent dans le cas d'un rapprochement d'images provenant d'acquisitions sur patients.

Une première opération est un recalage géométrique des images par rapport à une image de référence. Suivant les cas, ce recalage peut être rigide, c'est-à-dire qu'il ne change pas les formes, ou bien affine, où il modifiera les formes soit selon un certain nombre de degrés de liberté ou encore une transformation plus générale. Il convient d'observer au passage qu'il s'agit en général d'images tridimensionnelles, dont l'élément de volume est également appelé"voxel". Toutefois, il peut également s'agir d'images bidimensionnelles.

Ainsi, on comprendra qu'il est très difficile d'obtenir la mise en correspondance exacte en volume d'une partie du corps d'un patient, par exemple son cerveau, d'un examen à un autre, sans un recalage préalable des images. De plus, la partie du corps examinée peut avoir légèrement évoluée dans sa forme.

Les techniques de recalage rigide ou affine permettent de ramener toutes les images à une même géométrie de base. Un exemple d'une telle technique est décrit dans EP 0 977 148, déjà cité.

A côté de cela, les images peuvent présenter des"biais spatial"et"biais temporel". Ainsi, l'image peut présenter un"biais spatial", c'est-à-dire qu'il existera dans certaines des images des zones dont l'intensité sera plus important, c'est à dire des zones plus lumineuses (plus claires) et des zones dont l'intensité sera moins importante, c'est à dire des zones moins lumineuses (plus sombres).

Ce"biais spatial"peut être traité, par exemple de la manière décrite dans : -"Etude de la symétrie bilatérale en imagerie cérébrale volumique", Sylvain Prima, Thèse de Doctorat, Université de Paris XI, mars 2001.

-"Automated Model-Based Bias Field Correction for MR Images of the Brain", K. Van Leemput, F. Maes, D. Vandermeulen, et P. Suetens, IEEE Transactions on Medical Imaging 18 (10) : 885- 896,1999.

Reste le problème du biais temporel", moins souvent abordé.

Dans notre cas, on peut le définir comme le fait que l'échelle de gris varie d'une image à l'autre.

Certaines techniques d'imagerie permettent une calibration que l'on peut qualifier d'absolue, en ce sens qu'il est possible de retrouver les mêmes niveaux de luminosité, de contraste dans des images provenant d'examens espacés de plusieurs semaines, voire de plusieurs mois. C'est le cas dans certaines techniques d'imagerie médicale par exemple.

D'autres techniques, qui ne permettent pas cette calibration absolue, vont donc fournir des images entachées d'un"biais temporel". C'est le cas, notamment, des machines d'imagerie par résonance magnétique nucléaire (IRM).

Une proposition pour tenter de corriger le"biais temporel" a été faite dans : -"New variants of a method of MRI Scale Normalization", Laslo G. Nyul, Jayaram K. Udupa, Information Processing in Medical Imaging-IPMI'99, Lecture Notes in Computer Science 1613 (1999), 490-495.

Mais, quelle que soit la qualité de ces pré-traitements, l'existence d'un"biais temporel"rend particulièrement difficile le rapprochement automatique d'images prises de façon échelonnée sur plusieurs semaines ou plusieurs mois.

C'est particulièrement gênant, pour suivre différentes affections, par exemple la sclérose en plaques. L'invention vise notamment à améliorer le pré-traitement du biais temporel, cette sous caractéristique de l'invention a un intérêt en tant que telle et est susceptible d'être revendiquée séparément. Cette amélioration n'est pas la principale caractéristique de l'invention, cette dernière se rapportant plutôt à l'analyse temporelle, quantitative et automatique des images.

D'autre part, l'analyse temporelle automatique des images médicales permet d'établir a posteriori une détection de zones à évolution pathologique, par exemple dans le cas de lésions de sclérose en plaques. Cette détection est particulièrement utile pour suivre a posteriori l'effet de certains médicaments sur l'organisme dans le cas d'essais thérapeutiques. Des méthodes d'analyse temporelle des images médicales ont été proposées dans les ouvrages suivants :

- brevet EP 0 927 405"Dispositif électronique de traitement d'images pour la détection de variations dimensionnelles", Guillaume Calmon, Jean-Philippe Thirion, -"Automatic Quantification of changes in the Volume of Brain Structures", Guillaume Calmon, Neil Roberts, Paul Eldridge and Jean-Philippe Thirion, First International Conference on Medical Image Computing And Computer-Assisted InterventionMICCAI'98, Lecture Notes in Computer Science 1496 (1998), 964- 973.

-"A Computerized Approach for Morphological Analysis of the Corpus Callosum", Chris Davatzikos, Marc Vaillant, Susan M.

Resnick, Jerry L. Prince, Stanley Letovsky and R. Nick Bryan, Journal of Computer Assisted Tomography-JCAT, vol. 20, pp 88- 97,1996.

-"Automatic Detection and Segmentation of Evolving Processes in 3D Medical Images : Application to Multiple Sclerosis", David Rey, Gérard Subsol, Hervé Delingette et Nicholas Ayache, Information Processing in Medical Imaging-IPMI'99, Lectures Notes in Computer Science 1613 (1999) 154-167.

-"Using Continuum Mechanics Operators for Detection and Quantification of Evolving Processes in 3D Medicals Images", David Rey, Hervé Delingette, Gérard Subsol, Nicholas Ayache, EUROMECH'99,1999, 185-188.

Ces techniques permettent notamment de détecter et de quanti- fier les zones évolutives dans des séries d'images. Dans la plupart des cas, l'analyse est effectuée sur des couples d'images successives. Ainsi, ces techniques permettent d'effectuer une analyse à court terme de l'évolution des lésions pour aider au diagnostic immédiat du patient.

Une méthode d'analyse temporelle automatique des images médicales a été proposée dans l'ouvrages suivant : -"Exploring the discrimination power of the time domain for segmentation and characterization of active lesions in serial MR", Guido Gerig, Daniel Welti, Charles R. C. Guttmann, Alan C. F. Colchester, Gabor Székely, Medical Image Analysis-MedIA, 4 (1) 31-42,2000.

Cette technique d'analyse temporelle automatique permet sur une série d'images d'un patient de détecter des zones à évolution pathologique. Toutefois, cette technique est limitée à la prise en considération d'informations temporelles pour chaque unité de volume considérée (appelée voxel), par exemple le profil temporel d'intensité pour chaque unité de volume.

Une statistique descriptive (application de moyenne, d'écarttype) ainsi qu'une logique floue sont appliquées à ces informations temporelles. Ces applications ne permettent pas de prendre en compte la dimension spatiale des zones à évolution pathologique.

On se réfère à la figure la décrivant l'état de l'art et aux figure 4a, 4b, 4c et 4d. Un écran de visualisation 1, par exemple un écran d'ordinateur, est relié à un système d'exploitation 3, lui-même relié à un dispositif de traitement d'image 5. Ce dispositif 5 comprend une mémoire de données 6, un module de pré-traitement 8 et un module de traitement comparatif 10. La mémoire de données 6 est adaptée pour stocker des données de séries d'image en 3D, par exemple des données de voxels représentés par une position spatiale de coordonnées spatiales (x, y, z) pour chaque voxel et une intensité I (x, y, z, t1) à un instant donné tj pour chaque voxel. j est un entier pouvant varier de 1 à n, n étant un entier représentant le nombre d'instants différents auxquels des images 3D d'un même patient sont acquises. Ces séries d'images 3D peuvent être des séries d'images d'un même patient

prises selon la technique de l'imagerie IRM. Comme les images sont acquises à des instants éloignés dans le temps (une semaine, voire plus), le patient ne prend pas la même position. Ceci est illustré par les images 4a et 4b représentant respectivement une première image acquise à un instant donné et une deuxième image acquise à un instant différent présentant un décalage spatial par rapport à la première image 4a.

Ainsi, comme indiqué sur les figures 4a et 4b, une coordonnée (x, y, z) d'un voxel d'une image ne correspond pas au même point anatomique sur toutes les images de la série. De plus, les techniques d'acquisition d'images médicales peuvent introduire des inhomogénéités dans l'image elle-même (biais spatial) ou entre les images acquises à des instants différents (biais temporel) dans chaque image et/ou entre les images d'une série. Ceci est illustré par les images 4a et 4c représentant respectivement une première image acquise à un instant donné et une deuxième image, acquise à un instant différent, alignée sur la première image mais présentant un biais temporel par rapport à la première image 4a.

Aussi, comme vu précédemment, le module de pré-traitement 8 est adapté pour effectuer différents pré-traitements : - aligner les images d'une série sur l'une des images prise en référence par un recalage (rigide, affine ou autre) comme présenté sur la figure 4c qui représente l'image 4b alignée sur la figure 4a, - corriger éventuellement le biais spatial et le biais temporel selon des techniques exposées dans les ouvrages précédemment cités. La correction du biais temporel est présentée sur la figure 4d qui représente l'image 4c corrigée par rapport à la figure 4a.

D'autre part, le module de traitement comparatif 10 est adapté, après pré-traitement d'image, pour une analyse temporelle automatique des images, notamment des images de cerveau dans le cas de la sclérose en plaques. Ce module de

traitement comparatif 10 est apte à fonctionner selon l'une des techniques des ouvrages précédemment cités.

Toutefois, ces techniques ne permettent pas une analyse statistique quantitative d'une série temporelle d'images. Une telle analyse permettrait de retrouver les zones correspondant à une évolution pathologique grâce à un calcul de probabilités. En particulier, ceci est indispensable pour mener une étude statistique rigoureuse de l'effet de médicaments administrés dans le cas d'essais thérapeutiques.

La figure 1b est une réalisation du dispositif selon l'invention. Associé à l'écran de visualisation 1 et au système d'exploitation 3, un dispositif de traitement d'image 11 comprend une mémoire de données 16, un module de pré-traitement 18 et un module de traitement comparatif 20. La mémoire de données 16 est apte à stocker par exemple en mémoire, pour chaque image en 3D formée de voxels de dimensions 1x1x3 mm, 54 coupes d'images 2D ayant comme format 256x256. Le module de pré-traitement 18 correspond au module de pré-traitement 8 de la technique antérieure. Avantageusement, le module de pré-traitement 18 comprend une fonction de correction de biais temporel des intensités de voxels IRM basée sur un algorithme permettant la correction de l'histogramme conjoint entre deux images. Cet histogramme représente les intensités d'une image par rapport aux intensités de l'autre. L'algorithme est basé sur la recherche d'une droite de régression par la méthode des moindres carrés orthogonale robuste, méthode développée dans la publication "CRC Concise Encyclopedia of Mathematics", E. Weisstein. CRC Press LLC, Novembre 1998.

La droite de régression J=a. I+b est obtenue, J représentant l'intensité d'une première image, I l'intensité d'une deuxième image. La correction est obtenue en appliquant la transformation ainsi trouvée à la deuxième image. Ainsi, l'image 3a est formée d'une alternance de cubes d'une première image 3D de

référence et d'une deuxième image présentant un biais temporel par rapport à la première image de référence. Quant à l'image 3c, elle représente l'histogramme conjoint correspondant à l'image 3a. Par l'application de l'algorithme de correction de biais temporel, l'image 3b ne présente qu'une différence entre la première image de référence et la deuxième image corrigée, différence représentée par le cube B. L'histogramme 3d correspond à cette image 3b représentant une comparaison entre l'image de référence et la deuxième image corrigée.

Le module de traitement comparatif 20 comprend une fonction de modélisation 12 et une fonction d'analyse statistique 14.

La fonction de modélisation 12 est apte à modéliser une courbe représentant le profil d'intensité de voxels de zones évolutives au cours du temps. Cette fonction s'appuie sur un modèle mathématique paramétrique pouvant représenter la forme des différents profils d'intensité.

A partir de cette modélisation décrite plus loin, la fonction d'analyse statistique 14 procède à une inférence statistique. Selon cette inférence statistique, il est possible de détermi-

ner des amas de points ou de voxels suivant"significativement", au sens statistique, un modèle d'évolution patholo- gique.

Ainsi, la fonction de modélisation 12 et la fonction d'analyse statistique 14 du module de traitement comparatif 20 sont présentés ci-dessous selon un premier mode de réalisation préférentiel illustré par les figures 5a à 5j, 6a à 6f et 7.

Un second mode de réalisation préférentiel est également décrit plus loin.

Par exemple, pour la sclérose en plaques, une telle courbe comprend typiquement une partie montante et une partie descendante au cours du temps. Ainsi, il a été choisi de modéliser une courbe de Gauss asymétrique avec cinq paramètres

comme présentée à l'annexe I-a.. Les paramètres de cette courbe représentent - pl, un paramètre d'amplitude maximale, - p2, un paramètre d'amplitude minimale, - p3, un paramètre de temps au pic, - p4, un paramètre de largeur de référence et - p5, un paramètre de dissymétrie.

Dans ce mode de réalisation, une méthode semi-manuelle est utilisée pour calculer un modèle moyen de courbe de type Gauss asymétrique, comme représenté en figure 5g, en fixant les cinq paramètres. Ce modèle moyen sera indifféremment désigné comme modèle moyen de profil d'intensité ou modèle moyen paramétrique. Par ce modèle moyen continu, la suite temporelle discrète de valeurs d'intensité I (x, y, z, tl) mesurées en certains points (x, y, z) de l'espace est approximée au mieux. Ainsi, un opérateur peut sélectionner, sur les images 3D préalablement corrigées, un jeu de points correspondant à des parties de lésions évolutives. Ainsi, les figures 5a à 5f illustrent six profils d'intensité évoluant au cours du temps de six points sélectionnés. Le"temps de montée", défini comme le temps que met l'intensité à passer de 5% de son maximum à son maximum, c'est-à-dire le temps entre la ligne verticale DM et le maximum d'intensité M sur les figures 6a à 6d, est calculé pour chaque profil. Les profils correspondant à un "temps de montée"de moins d'une semaine ou de plus de dix semaines sont rejetés, ne correspondant pas à un profil d'intensité considéré comme pathologique. Les profils 6e et 6f sont aberrants et sont donc rejetés. Les profils retenus sont centrés et normalisés selon leur maximum respectif. Le modèle moyen est alors calculé selon les profils normalisés et selon une estimation classique des moindres carrés. Les cinq paramètres sont calculés selon une estimation des moindres carrés utilisant une méthode de convergence quadratique de Powell, présentée dans l'ouvrage suivant :

[1]- W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery. Numerical Recipes. The Art of Computing. Cambridge University Press, 2nd edition, 1997.

Par exemple, dans une étude réalisée sur 54 profils extraits, 43 profils ont été retenus et le modèle moyen de profil d'intensité représenté sur la figure 5h a été calculé selon les profils d'intensité de ce jeu de points sélectionnés et selon la technique proposée ci-dessus. La figure 5i représente le modèle moyen superposé aux mesures discrètes normalisées. La figure 5j représente une discrétisation du modèle moyen obtenu.

La fonction d'analyse statistique est ensuite utilisée avec ce modèle moyen obtenu. Elle est répétée aux instants où on cherche les zones évolutives. Pour le cas de la sclérose en plaques, l'analyse statistique est appliquée à chaque instant discret t, avec j=l.. n, n représente le nombre d'images prises à des instants différents pour chaque voxel. Le modèle général linéaire permet de rechercher la présence du modèle moyen de profil d'intensité en tout voxel de l'image. Le modèle général linéaire est décrit dans l'ouvrage suivant : [2] R. S. J. Frackowiak, K. J. Friston, C. D. Frith, R. J. Dolan, and J. C. Mazziotta. Human Brain Function. Academic Press, 1997.

Chaque voxel est ajusté suivant un modèle général linéaire comme présenté en annexe I-b. La figure 7 illustre ce modèle général linéaire. Chaque colonne d'un vecteur de dimension nxx, x étant un entier, est composée de n composantes faisant référence chacune à un instant discret sur l'échelle des temps ti. Chaque colonne d'un vecteur forme ainsi un profil de valeurs discrètes en fonction du temps. Ainsi, le vecteur Y (dimension nxl) représente les données discrètes d'origine du profil d'intensité d'un voxel donné. La matrice X (dimension nx2) est appelée la matrice de conception. X a pour compo- sante, comme illustré sur la figure 7, les valeurs discrètes du modèle moyen de profil d'intensité, représenté par le

vecteur XI (dimension nxl), et une valeur constante au cours du temps, représentée par le vecteur X2 (dimension nxl), pour prendre en compte l'intensité moyenne du voxel au cours du temps ti. Les vecteurs X et Y sont connus. Les composantes ssl et ss2 du vecteur ss (2 xl) représentent respectivement la variable d'adéquation du modèle moyen de profil d'intensité et du niveau de référence du profil d'intensité d'un voxel donné. Le vecteur e représente le vecteur des termes de résidus. Les vecteurs ss représente les paramètres estimés du modèle général linéaire.

L'inférence statistique permet de détecter les voxels voisins ayant un même comportement pathologique évolutif. Pour cela, on calcule pour chaque voxel une variable t, appelée variable de Student telle que définit en annexe I selon l'ouvrage [2] et l'ouvrage suivant : - [3] D. J. Sheskin. Handbook of Parametric and Non-parametric Statistical Procedures. CRC Press, 1997.

Une carte de t est dressée. Dans la formule du calcul de t,

02 représente la variance résiduelle estimée par le carré moyen résiduel. Le vecteur c, défini en annexe I-c, permet de choisir la variable ss1 du vecteur ss. Dans la formule de l'annexe I, les valeurs de la variable ss1 sont divisées par la variance résiduelle. En estimant que les résidus sont distribués selon une distribution gaussienne, les valeurs de t suivent une distribution de Student (appelée aussi une "distribution en t").

Des voxels présentant des valeurs significativement grandes de t sur la carte de t obtenue sont sélectionnés lorsque leur valeur de t est supérieure à une valeur seuil Th choisie. Parmi ces voxels sélectionnés apparaissent des groupes de voxels voisins.

L'inférence statistique utilise la théorie des champs aléatoires développé dans l'ouvrage [2] et dans l'ouvrage suivant :

- [4] K. J. Worsley. Local maxima and the expected euler characteristic of excursion sets of x f and t fields.

Advances in Applied Probability, 26 : 13-42,1994.

Cette théorie permet de déterminer la probabilité d'obtenir des groupes de voxels d'une certaine taille. Les groupes de voxels ayant une probabilité très faible d'occurrence sont considérés comme "significatifs", c'est-à-dire, dans le cas d'une recherche de zones pathologique, ces groupes de voxels sont susceptibles de représenter une zone pathologique évolutive. Grâce à cette analyse statistique, on peut associer une probabilité à la taille d'un ensemble de voxels voisins.

Seuls les ensembles de voxels voisins pour lesquels cette probabilité d'occurrence est suffisamment faible sont considérés comme significatifs. Une probabilité très faible d'occurrence peut être une probabilité de 0,01 par exemple.

Afin de pouvoir appliquer la théorie des champs aléatoires, certaines hypothèses doivent être vérifiées : - l'homogénéité du bruit dans les images ; - une bonne représentation en treillis.

- l'indépendance temporelle du bruit ; - le bruit gaussien ; - l'homogénéité de la fonction de répartition de points effectifs.

Un second mode de réalisation est décrit ci-après. Ainsi, les courbes représentant le profil d'intensité de voxels de zones évolutives au cours du temps sont modélisées par un modèle paramétrique choisi. Dans cet exemple de réalisation, le modèle paramétrique de l'annexe I-a présentant une courbe de Gauss asymétrique avec cinq paramètres est choisi.

Ainsi, un modèle paramétrique est ajusté en chaque voxel sur un profil d'intensité temporel tel que ceux montrés sur les figures 5a à 5f. En d'autres termes, les paramètres d'un

modèle sont déterminés pour chaque voxel. Il peut être difficile de déterminer en même temps l'ensemble des paramètres d'un modèle lié à un profil d'intensité. Ainsi, avantageusement, il est possible de fixer par exemple le paramètre de largeur de référence p4 et le paramètre de dissymétrie p5. Ces paramètres sont en effet estimés sur un ensemble de profils normalisés en temps et en amplitude. Une fois ces paramètres fixés, les paramètres restant sont évalués pour chaque voxel, comme par exemple l'amplitude maximale pl, le paramètre d'amplitude minimale p2, et le paramètre de temps au pic p3. Les paramètres sont évalués au sens des moindres carré en utilisant la méthode de convergence quadratique de Powell présentée dans le référence suivante : - W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P.

Flannery. Numerical Recipes. The art of Scientific Computing.

Cambridge University Press, 2nd edition, 1997.

Ainsi, on obtient différentes courbes d'intensités selon les voxels. Ce résultat amène à une analyse non linéaire dans ce cas.

Afin de déterminer la probabilité d'obtenir des groupes de voxels d'une certaine taille, une analyse statistique est appliquée.

Cette analyse statistique est décrite plus en détail ci-après en référence à la figure 2c.

Les figures 2a et 2b illustrent le procédé automatique d'analyse d'images selon l'invention.

A l'étape 110, les images 3D d'un patient au cours du temps sont alignées sur une image de référence par recalage rigide comme précédemment expliqué. A l'étape 120, une correction temporelle est effectuée sur les images 3D, cette correction concerne le biais d'intensité. Une correction spatiale est

éventuellement effectuée. Les étapes 130 à 150 présentent des étapes générales particularisées par l'ordinogramme de la figure 2b. Ainsi, un modèle paramétrique de type Gaussien asymétrique comportant 5 paramètres est préparé à l'étape 130 pour modéliser les profils d'intensité des voxels des images 3D. De cette façon, le modèle paramétrique est ajuster sur les profils d'intensité des voxels des images 3D à l'étape 140. Ainsi, une (voire plusieurs) carte permettant de déterminer les groupes de voxels significativement pathologiques est élaborée à l'étape 150. Les zones pathologiques peuvent être détectées par inférence statistique à l'étape 160.

La figure 2b illustre plus particulièrement les étapes 130 à 150 du procédé.

A l'étape 132, un modèle moyen paramétrique est au préalable élaboré suivant un jeu de profils d'intensité de voxels. Selon une réalisation de l'invention, ce jeu de profils d'intensité de voxels est choisi par un opérateur comme étant un jeu de profils représentatifs de voxels présentant des lésions évolutives. Comme précédemment décrit, les cinq paramètres du modèle moyen paramétrique sont ainsi fixés. A l'étape 142, les paramètres du modèle général linéaire sont estimés en chaque voxel du cerveau. Ainsi, chaque profil d'intensité d'un voxel présenté sous forme d'un vecteur de données est représenté par une combinaison linéaire du modèle moyen paramétrique, d'une constante et des résidus.

A l'étape 152, un calcul de valeur de t est effectué pour chaque voxel, c'est-à-dire qu'on obtient une carte en t. A l'étape 154, la carte de t permet de déterminer les groupes de voxels présentant une valeur supérieure à une valeur seuil Th.

Par application de la théorie des champs aléatoires décrite dans l'ouvrage [4] cité ci-dessus, des groupes de voxels parmi

ceux déterminés présentent une probabilité d'occurrence très petite et forment des groupes de voxels significatifs à l'étape 156. De cette manière sont détectées et modélisées les zones pathologiques probables. Cette inférence statistique permet de détecter des zones pathologiques significatives statistiquement en prenant en compte la structure spatiale et temporelle des séries d'images 3D.

La figure 2c illustre plus particulièrement les étapes 130 à 150 du procédé selon un second mode de réalisation.

A l'étape 232, un modèle paramétrique est choisi, par exemple celui de l'annexe I-a. Dans ce mode de réalisation, à l'étape 242, ce modèle est déterminé indépendamment pour chaque voxel d'une série d'images de patient. Ainsi, les paramètres du modèle sont ajustés pour chaque voxel. Ceci permet d'obtenir un modèle avec des paramètres particuliers pour chaque voxel.

Avantageusement, avant l'étape 242, au moins un paramètre peut être déterminé selon une moyenne réalisée sur un échantillon de voxels afin de diminuer le temps de calcul des autres paramètres pour chaque voxel à l'étape 242.

Une analyse statistique est ensuite effectuée des étapes 252 à 258. A l'étape 252, une permutation des images est effectuée sur la série d'images considérée pour laquelle un modèle paramétrique particulier a été déterminé pour chaque voxel.

La permutation de ces images est avantageusement aléatoire. De cette façon, chaque image est affectée au hasard à un instant de la série d'images. Plus précisément, une permutation identique est appliquée à chaque voxel de chaque image pour préserver la structure de corrélation spatiale des images.

A l'étape 254, un modèle paramétrique est alors ajusté pour chaque voxel comme précédemment.

A l'étape 256, un seuil t ayant une valeur appropriée est appliqué à au moins un des paramètres déterminés précédemment pour les voxels de la série d'images permutées. De manière particulière, cette valeur t est appliquée par exemple à l'amplitude maximale des modèles obtenus. Les voxels dépassant cette valeur t sont détectés comme étant des voxels"hors seuil". Certains voxels hors seuil voisins forment des groupes de voxels dont la taille est alors déterminée. La valeur choisie de t n'influence pas la validité des statistiques mais la sensibilité des résultats.

Les étapes 252 à 256 sont répétées un nombre donné de fois, par exemple N fois à l'étape 257. Cet entier est déterminé de manière à correspondre à des exigences de temps de calcul corrects et à avoir une bonne approximation de la distribution des probabilités de densités. Par exemple, N peut être égal à 10.

Après un nombre donné d'itérations, les tailles significatives des groupes de voxels hors seuil sont déterminées par la méthode décrite ci-après. Ainsi, on cherche à obtenir une distribution des tailles de groupes de voxels sous une hypothèse dite nulle. L'hypothèse nulle est vérifiée quand il n'existe aucun processus évolutif pathologique cohérent. La distribution des tailles de groupes de voxels permet de détecter une lésion évolutive pour les tailles de groupes de voxels ayant une probabilité d'apparaître plus petite qu'une valeur critique a sous l'hypothèse nulle. Ainsi, à l'étape
258, les groupes de voxels significatifs sont déterminés par inférence statistique.

La figure 8 illustre une distribution de taille de groupes de voxels pour des séries d'images permutées. Ainsi, les barres représentent la densité obtenue pour chaque taille de groupes de voxels calculée en nombre de voxels par groupe. La densité représente le nombre de groupes d'une taille donnée divisée

par le nombre total de groupes considérés. Dans l'exemple de la figure 8, t a pour valeur 0,2 et a a pour valeur 0, 01. Dans cet exemple, les groupes de voxels significatifs ont une taille plus grande que 10 voxels.

L'hypothèse selon laquelle les images de la série sont échangeables permet de réaliser rigoureusement cette méthode statistique par permutations.

L'amélioration du pré-traitement du biais temporel est une sous caractéristique de l'invention qui a un intérêt en tant que telle et qui est susceptible d'être revendiquée séparément.

Bien entendu, l'invention n'est pas limitée à la forme de réalisation décrite précédemment à titre d'exemple, elle s'étend à d'autres variantes. Ainsi, on pourra utiliser d'autres modèles paramétriques. D'autres études statistiques peuvent être menées avec des tables de distribution différentes. Dans le second mode de réalisation, la méthode peut être appliquée sans permutation des images. L'invention s'étend également à d'autres pathologies que la sclérose en plaques.

Annexe I

Claims

Revendications 1. Dispositif de traitement d'images, comprenant : - une entrée pour recevoir une suite temporelle de jeux de données représentant des images numériques volumiques comparables, chaque donnée comprenant une composante de position et une composante d'intensité, - des moyens de pré-traitement (8 ; 18) pour modifier les jeux de données de façon à obtenir des images recalées géométriquement et en intensité, et - des moyens de traitement comparatif (10 ; 20), propres à examiner des jeux de suites temporelles d'éléments d'image, afin d'y détecter des signes de variations, caractérisé en ce que - les moyens de traitement (20) comprennent : * une fonction de modélisation (12), capable d'ajuster un modèle paramétrique séparément sur certains au moins des jeux de suites temporelles d'éléments d'image, ce qui fournit des couples d'informations du genre : élément d'image, temps ; et * une fonction d'analyse statistique (14) de ces couples d'informations, propre à isoler ceux des couples d'informations qui sont représentatifs d'une variation significative.
2. Dispositif selon la revendication 1, caractérisé en ce que le modèle paramétrique est du type gaussienne asymétrique, comprenant un paramètre d'amplitude maximum, un paramètre d'amplitude minimum, un paramètre de temps au pic, un para- mètre de largeur de référence et un paramètre de dissymétrie.
3. Dispositif selon la revendication 2, caractérisé en ce que le modèle paramétrique comprend un modèle moyen dont les paramètres sont établis à partir de jeux de suites temporelles d'éléments d'image mesurés sur des zones pathologiques évolutives.

<Desc/Clms Page number 26>
4. Dispositif selon la revendication 3, caractérisé en ce que lesdits paramètres du modèle moyen sont obtenus à partir de jeux de suites temporelles normalisées en temps et en amplitude.
5. Dispositif selon l'une des revendications 1 à 2, caractérisé en ce que la fonction de modélisation (12) comprend en outre un modèle linéaire général appliqué à chaque jeu de suites temporelles d'éléments d'image de façon à obtenir, pour chacun desdits jeux de suites, une combinaison linéaire d'au

moins le modèle moyen (Xl) et une constante (X2) suivant un vecteur (ss) comprenant au moins un premier et un deuxième coefficient.
6. Dispositif selon la revendication 5, caractérisé en ce que la fonction d'analyse statistique (14) comprend une opération du type T-test révélant des zones à évolution pathologique connexes où ledit premier coefficient est significativement différent de zéro.
7. Dispositif selon l'une des revendications 1 et 2, caractérisé en ce que la fonction de modélisation (12) comprend en outre d'ajuster au moins un modèle paramétrique pour au moins une suite temporelle d'élément d'image.
8. Dispositif selon la revendication 7, caractérisé en ce qu'au moins un paramètre du modèle est fixé à partir de jeux de suites temporelles d'éléments d'images normalisées en temps et en amplitude.
9. Dispositif selon la revendication 7, caractérisé en ce qu'au moins un paramètre du modèle est obtenu à partir du jeu de suite temporelle dudit élément d'image.

<Desc/Clms Page number 27>
10. Dispositif selon la revendication 1 en combinaison avec l'une des revendications 7 à 9, caractérisé en ce que la fonction d'analyse statistique (14) comprend de répéter, pour au moins une suite dans un ordre temporel différent dudit élément d'image, la fonction de modélisation (12) comprenant d'ajuster un modèle paramétrique pour une suite temporelle d'élément d'image.
11. Dispositif selon la revendication 10, caractérisé en ce que la fonction d'analyse statistique (14) comprend la détermination de groupes d'éléments d'image voisins correspondant à des suites d'éléments d'image ayant au moins un des paramètres de leur modèle supérieur à un seuil donné, ces groupes d'éléments d'image voisins formant des zones pathologiques évolutives d'une taille déterminée par la fonction d'analyse statistique.
12. Dispositif selon la revendication 11, caractérisé en ce que la fonction d'analyse statistique (14) comprend en outre une distribution desdites zones pathologiques évolutives suivant leur taille et une détermination de zones pathologiques évolutives significatives.
13. Dispositif selon l'une des revendications 3,5 et 7, caractérisé en ce que la fonction de modélisation (12) comprend une fonction d'optimisation appliquée aux paramètres du modèle paramétrique.
14. Dispositif selon l'une des revendications précédentes, dans lequel le prétraitement comprend un recalage rigide, suivi d'une correction du biais temporel, caractérisé en ce que la correction du biais temporel comprend la construction d'un histogramme conjoint 3D des intensités d'images prises deux à deux.

<Desc/Clms Page number 28>
15. Dispositif selon la revendication 14, caractérisé en ce que la correction du biais temporel comprend la construction d'une droite de régression sur ledit histogramme par une technique de moindres carrés orthogonale robuste.
16. Dispositif selon l'une des revendications 14 et 15, caractérisé en ce que la correction du biais temporel comprend la correction en conséquence de l'une des images par rapport à l'autre.
17. Procédé de traitement d'images, comprenant les étapes suivantes : a-recevoir une suite temporelle de jeux de données représentant des images numériques volumiques comparables, chaque donnée comprenant une composante de position et une composante d'intensité, b-modifier les jeux de données de façon à obtenir des images recalées géométriquement et en intensité, caractérisé en ce qu'il comprend les étapes complémentaires c-ajuster un modèle paramétrique (140) séparément sur certains au moins des jeux de suites temporelles d'éléments d'image, ce qui fournit des couples d'informations du genre : élément d'image, temps, d-isoler ceux des couples d'informations qui sont représentatifs d'une variation significative par analyse statistique (156).
18. Procédé selon la revendication 17, caractérisé en ce que le modèle paramétrique de l'étape c-est du type gaussienne asymétrique, comprenant un paramètre d'amplitude maximum, un paramètre d'amplitude minimum, un paramètre de temps au pic, un paramètre de largeur de référence et un paramètre de dissymétrie.
19. Procédé selon la revendication 18, caractérisé en ce que le modèle paramétrique de l'étape c-comprend un modèle moyen

<Desc/Clms Page number 29>

dont les paramètres sont établis à partir de jeux de suites temporelles d'éléments d'image mesurés sur des zones pathologiques évolutives.
20. Procédé selon la revendication 19, caractérisé en ce que lesdits paramètres du modèle moyen de l'étape c-sont obtenus à partir de jeux de suites temporelles normalisées en temps et en amplitude.
21. Procédé selon l'une des revendications 17 à 18, caractérisé en ce que l'étape c-comprend cl-appliquer un modèle linéaire général à chaque jeu de suites temporelles d'éléments d'image de façon à obtenir, pour chacun desdits jeux de suites, une combinaison linéaire d'au moins le modèle moyen et une constante suivant un vecteur comprenant au moins un premier et un deuxième coefficient (142).
22. Procédé selon la revendication 21, caractérisé en ce que l'étape d-comprend dl-appliquer aux couples d'information une opération du type T-test révélant des zones à évolution pathologique connexes où ledit premier coefficient est significativement différent de zéro (152).
23. Procédé selon l'une des revendications 17 et 18, caractérisé en ce que l'étape c-comprend cl-ajuster au moins un modèle paramétrique pour au moins une suite temporelle d'élément d'image.
24. Procédé selon la revendication 23, caractérisé en ce que l'étape cl-comprend en outre de fixer au moins un paramètre de modèles à partir de jeux de suites temporelles normalisées en temps et en amplitude.

<Desc/Clms Page number 30>
25. Procédé selon la revendication 23, caractérisé en ce que l'étape cl-comprend d'obtenir au moins un paramètre du modèle à partir du jeu de suite temporelle dudit élément d'image.
26. Procédé selon la revendication 17 en combinaison avec l'une des revendications 23 à 25, caractérisé en ce que l'étape d-comprend dl-répéter l'étape c-pour au moins une suite dans un ordre temporel différent dudit élément d'image.
27. Procédé selon la revendication 26, caractérisé en ce que l'étape d-comprend d2-déterminer des groupes d'éléments d'image voisins correspondant à des suites d'éléments d'image ayant au moins un des paramètres de leur modèle supérieur à un seuil donné, d3-déterminer la taille de zones pathologiques évolutives formées par ces groupes d'éléments d'image voisins.
28. Procédé selon la revendication 23, caractérisé en ce que l'étape d-comprend d4-distribuer lesdites zones pathologiques évolutives suivant leur taille, d5-déterminer les zones pathologiques évolutives significatives.
29. Procédé selon l'une des revendications 19,21 et 23 caractérisé en ce que l'étape c-comprend d'optimiser les paramètres du modèle paramétrique.
30. Procédé selon l'une des revendications précédentes, dans lequel l'étape b-comprend d'effectuer un recalage rigide, suivi d'une correction du biais temporel, caractérisé en ce que l'étape b-comprend d'effectuer la correction du biais temporel par la construction d'un histogramme conjoint 3D des intensités d'images prises deux à deux.

<Desc/Clms Page number 31>
31. Procédé selon la revendication 30, caractérisé en ce que l'étape b-comprend d'effectuer la correction du biais temporel par la construction d'une droite de régression sur ledit histogramme par une technique de moindres carrés orthogonale robuste.
32. Procédé selon l'une des revendications 30 et 31, caractérisé en ce que l'étape b-comprend d'effectuer la correction du biais temporel par la correction en conséquence de l'une des images par rapport à l'autre.