Brevet FR2886433A1 — publié le 20061201

La présente invention concerne une méthode de segmentation d'une image ou d'une séquence d'images tridimensionnelles en une image tridimensionnelle partitionnée en régions d'intérêt. L'invention s'applique d'une manière générale à tous types d'images, en particulier dans les domaines de la pharmaco-imagerie et de l'imagerie satellite. Elle s'applique également aux images bidimensionnelles, qui sont des cas particuliers d'images tridimensionnelles dans lesquelles la valeur d'une des dimensions est égale à l'unité. Une évolution récente et remarquable dans le développement des médicaments est l'utilisation de la pharmaco-imagerie pour mesurer les paramètres pharmacocinétiques. De manière connue, la pharmacocinétique est la composante de la pharmacologie qui décrit la disposition des médicaments dans l'organisme, et précise de façon qualitative et quantitative les processus d'absorption, de distribution, de métabolisation et d'élimination d'un principe actif. Le grand avantage de la pharmaco-imagerie par rapport aux techniques classiques, qui nécessitent le prélèvement répété d'échantillons biologiques, est sa capacité à mesurer la concentration du principe actif dans les organes du patient humain ou animal vivant. On évite ainsi de sacrifier des lots d'animaux, et l'on diminue de façon substantielle la dispersion des valeurs - toutes les mesures de la cinétique étant faites chez le même animal. De plus, on peut analyser la cinétique dans les organes profonds chez l'homme. En contrepartie, la pharmaco-imagerie requiert une mise en oeuvre qui peut être technologiquement lourde: d'abord le marquage, radioactif ou autre, du principe actif à suivre; ensuite, la preuve que la détection par imagerie donne bien une mesure quantitative de la concentration du principe actif; enfin, le repérage de la localisation anatomique du signal lié au principe actif marqué. D'une façon générale, mais non exclusive, une technique préférée de la pharmaco-imagerie est la tomographie d'émission de positons (TEP) qui, du fait du principe de détection sur lequel elle s'appuie, est parmi les techniques d'imagerie les mieux adaptées aux mesures quantitatives de concentrations moléculaires dans les organes de l'homme ou de l'animal. En tout état de cause, la mesure de la concentration d'un principe actif marqué (appelé ci-après traceur) en un point de l'organisme étudié n'est interprétable sur le plan pharmacocinétique que si l'on peut assigner ce point à un organe défini, soit sur une base anatomique, soit sur une base physiologique: le coeur, le foie, le rein, une tumeur, une région du cerveau, etc. Or, l'image est uniquement représentative de la localisation de l'indicateur (émetteur de positons radioactif dans le cas de la TEP) lié au principe actif étudié, et ne contient a priori aucune information anatomique ou physiologique, mais uniquement une information pharmacocinétique. Il est donc nécessaire de superposer les images successives de localisation du traceur à une ou des images renseignant sur l'anatomie et/ou la physiologie. Cette superposition ou recalage peut se faire de plusieurs manières: É Dans le cas le plus favorable et le moins fréquent, la distribution du traceur dévoile une anatomie reconnaissable. C'est le cas des traceurs se concentrant dans l'ensemble d'un système anatomique remarquable comme le squelette, qui ne nécessitent pas de recalage supplémentaire, toute l'information étant dans l'image. Inversement, c'est aussi le cas de certains traceurs qui diffusent très bien dans tous les organes (par exemple le fluorodésoxyglucose). Dans ce dernier cas, les organes sont délimités sur les images par un contraste dépendant du niveau de rétention du traceur, ce qui n'est pas forcément suffisant pour leur identification, par exemple dans le cas de deux organes accolés ayant le même niveau de rétention du traceur. D'une façon générale, ce cas où le recalage est peu ou pas utile ne concerne que les traceurs généralistes (ions, précurseurs métaboliques, etc.) et pas l'imagerie des principes actifs. É Dans de nombreux autres cas, les organes sont repérés par plaquage d'une image anatomique issue des connaissances de l'examinateur sur les pharmaco-images. Un biologiste peut souvent deviner la localisation des organes anatomiques, tels que le foie ou le coeur, sur les images dont l'échelle de contraste suit la dynamique de concentration du traceur. La précision de la localisation obtenue dépend alors de la capacité à identifier l'organe, donc de la distribution du traceur qui délimite plus ou moins bien les contours d'un organe, et de l'habileté ou de l'habitude de l'examinateur. On comprend que cette méthode sans image anatomique est inapplicable au cas d'un traceur qui ne serait présent qu'en un point. É De plus en plus souvent, l'examinateur fait appel à une deuxième image qu'il superpose à la pharmaco-image, via un deuxième mode d'imagerie dans lequel le contraste n'est pas fondé sur la localisation du traceur mais sur l'anatomie (par exemple la tomodensitométrie, ou l'imagerie par résonance magnétique ou IRM) qui est appliquée au même individu, de préférence lors de la même session d'imagerie. En superposant les deux images, on recale la pharmaco-image sur une image anatomique qui identifie les organes. Cette double modalité est de plus en plus largement utilisée pour l'imagerie clinique avec l'avènement récent du système PET-CT associant caméra TEP et tomodensitomètre X. Toutefois, la double modalité n'est pas totalement satisfaisante pour la pharmaco- imagerie pour des raisons à la fois fondamentales et pratiques: É Sur le plan fondamental, l'image anatomique ne fait que plaquer une information statique sur l'information pharmaco-cinétique donnée par la pharmaco- image. Le lien, autre que le sujet, entre l'image anatomique et la pharmaco-image, ne peut donc être établi que si la localisation pharmacocinétique du traceur suit exactement les contours des organes décrits par l'image anatomique. Par conséquent, l'information supplémentaire apportée par l'image anatomique est limitée à la résolution de la méthode d'imagerie anatomique utilisée, qui est loin d'être toujours suffisante pour identifier la nature de tous les pharmaco- organes. Par exemple, la tomodensitométrie apporte peu de contraste dans les tissus mous et le cerveau, I'IRM distingue mal les masses pulmonaires, etc. É Sur le plan pratique: 0 - Les caméras couplées de type PET-CT ne sont actuellement disponibles que pour l'homme, donc avec une grande ouverture et une résolution moindre que les caméras TEP dédiées aux petits animaux. Ces caméras sont nécessairement plus chères que les caméras TEP seules. - Dans le cas d'une imagerie anatomique par tomodensitométrie, la dose d'irradiation nécessaire pour obtenir une image chez le petit rongeur est loin d'être négligeable. Ceci pose des problèmes de toxicité etlou d'interprétation de résultats dans le cas de l'imagerie des tumeurs. - La qualité du recalage obtenue est tributaire de l'immobilité totale du sujet pendant toute la durée de la pharmacocinétique explorée, ce 10 qui n'est pas le cas lors du transport d'un animal d'une caméra à une autre. Par ailleurs, l'interprétation pharmacocinétique des images TEP nécessite généralement la délimitation de régions d'intérêt ( ROI ) représentatives des organes. Kimura et coll. [Kimura, 2002] et Zhou [Zhou, 2000] ont montré qu'un regroupement adéquat des éléments de volume (voxels) améliore la qualité de la quantification pharmacocinétique. En effet, le tracé des régions d'intérêt suppose que chaque organe ou partie d'organe présente un comportement homogène pour un traceur donné, qui peut être caractérisé par une cinétique physiologique. Par conséquent, quelle que soit la méthode de segmentation, les limites des organes tracés doivent être visibles pour éviter que les régions d'intérêt n'englobent deux organes ou parties d'organes aux fonctions différentes (appelés pharmaco-organes ci-après). La qualité de la segmentation de l'image en régions d'intérêt cohérentes sur le plan pharmacocinétique est donc cruciale. Or, les régions d'intérêt sont en général tracées manuellement, ce qui est un travail long, fastidieux, dépendant de l'opérateur et nécessitant un degré d'expertise certain. Un certain nombre de travaux, dont le but n'était pas limité à l'analyse pharmacocinétique d'images, ont proposé des méthodes de segmentation automatique permettant de s'affranchir de l'opérateur pour la segmentation des régions d'intérêt. Les images TEP souffrent d'un faible rapport signal sur bruit. Comme elles souffrent également d'une faible résolution spatiale, l'activité mesurée en un élément de volume donné de l'image est polluée par l'activité contenue dans les éléments de volumes voisins. En particulier, la mesure de la concentration du traceur dans les organes de petite taille peut être soit sous-estimée, soit surestimée selon que la concentration radioactive dans les structures environnantes est moins élevée ou plus élevée. Cet effet, dit de volume partiel peut en effet être modélisé comme un lissage de l'image [Frouin, 2002] et rend donc les contours des organes intrinsèquement flous. En conséquence, les parties des organes proches des frontières de ces derniers ne contiendront pas la cinétique d'un unique pharmaco-organe, mais une combinaison linéaire des cinétiques de tous les pharmaco-organes proches. Ainsi, la cinétique de pharmaco-organe la plus représentée au sein d'un élément de volume n'est pas forcément celle du ou des pharmaco-organe(s) réellement contenu(s) dans cet élément de volume. En outre, les images TEP dynamiques représentent de manière connue une grande masse de donnée (un demi-million de cinétiques), dont le traitement est prohibitif en termes de temps de calcul et de mémoire informatique requise [Guo, 2003]. Dans les images, l'organisme étudié représente entre 20% et 40% des données. La partie restante, extérieure à l'organisme, contient principalement du bruit. On sait par ailleurs que la cinétique au sein d'un élément de volume (ce terme désigne l'unité volumique élémentaire de l'image) n'a de sens qu'à condition que le pharmaco-organe étudié soit parfaitement immobile au cours du temps d'acquisition de la séquence d'images. Tout mouvement du pharmaco-organe imagé brise le lien entre la cinétique d'un élément de volume et les cinétiques de ce pharmaco-organe. Les mouvements physiologiques, extrêmement difficiles à corriger, sont de deux sortes. Certains mouvements sont périodiques, de période en général inférieure à la durée d'acquisition d'une image de la séquence. Les battements cardiaques et la respiration imposent un déplacement des organes voisins, générant un effet de flou non négligeable, mais à peu près constant d'une acquisition à la suivante. D'autres mouvements non périodiques sont imprévisibles, comme le déplacement des viscères lors de la digestion et sous l'effet de la respiration, ou le remplissage de la vessie au cours de l'examen: cette dernière peut voir son volume apparent décupler entre le début et la fin de l'examen. Pour les traceurs excrétés dans les urines, la concentration du traceur et des métabolites dans les urines devient très élevée. Un voxel correspondant à une région viscérale ou musculaire en début d'examen peut donc contenir une cinétique caractéristique de la vessie en fin d'examen. Comme le montre la figure 2 jointe, le remplissage de la vessie génère une famille de cinétiques composées pour les premiers temps d'acquisition de la cinétique d'un pharmaco-organe proche de la vessie, et pour les derniers temps de la cinétique de la vessie. L'hypothèse d'un nombre fixe de cinétiques contenues dans une image TEP, au bruit et à l'effet de volume partiel près, ne s'applique donc pas dans le cas d'un organisme imagé sujet à des mouvements physiologiques non périodiques. Enfin, il est connu que le bruit dans les images TEP diffère 20 selon les méthodes de reconstruction d'images utilisées. Au sein d'une image issue d'une reconstruction analytique par rétroprojection filtrée, le bruit est souvent considéré comme additif gaussien stationnaire dans l'image, non corrélé avec le signal. En revanche, le bruit d'une image de la séquence dépend de la durée d'acquisition de cette image. Cette dépendance peut être considérée comme linéaire [Guo, 2003] par rapport à l'inverse de la durée d'acquisition de l'image. Dans les images issues de méthodes statistiques itératives, comme OSEM ( Ordered Subset Expectation Maximization , méthode reconstruisant l'image en maximisant la vraisemblance de sa projection selon diverse incidences) ou AWOSEM ( Attenuation Weighted Ordered Subset Expectation Maximization , méthode opérant de façon similaire à OSEM en prenant en compte le phénomène d'atténuation des photons par l'organisme), on ne peut plus supposer un bruit stationnaire dans l'espace. Le bruit dépend du nombre d'itérations utilisées pour la reconstruction, phénomène stationnaire, mais est aussi corrélé au signal. On peut écrire: G2 = a2 x S, où o représente la variance locale du bruit à un instant donné, a une constante indépendante de la position et du temps, mais dépendant de la méthode de reconstruction utilisée ainsi que du nombre d'itérations choisi, et S le signal au point considéré. Dans ce cas de figure, le rapport signal sur bruit s'exprime comme S/G = a x Sv'. La séparation du bruit et du signal en un point quelconque de l'image permet le calcul de (X et permet de caractériser finement le bruit. Comme dans le cas d'une reconstruction par méthode analytique, il conviendra de corriger les données de l'influence de la durée d'acquisition de chaque image de la séquence. Plusieurs méthodes de segmentation des images TEP en agrégats regroupant des zones à cinétiques homogènes, sans connaissance anatomique a priori, ont récemment été proposées pour tenter de résoudre les problèmes précités. Ashburner [Ashburner, 1996] suppose qu'une image TEP contient un nombre connu de cinétiques une par agrégat , et décrit donc toute cinétique d'un voxel de l'image comme la cinétique d'un agrégat, multipliée par un facteur d'échelle, à laquelle s'ajoute un bruit gaussien multinormal de loi normale pour chaque image acquise identique pour chaque agrégat. Wong et coll. [Wong, 2002] supposent les cinétiques homogènes au sein d'un même agrégat, mais dissimilaires entre agrégats différents. Ils proposent une approche d'agrégation par la méthode des k-moyennes, qui maximise la variance des cinétiques entre les classes, tout en minimisant la variance des cinétiques au sein d'une même classe. La cinétique d'un agrégat est alors estimée comme la moyenne des cinétiques des individus le composant. Frouin F. et coll. [Frouin, 2001] utilisent aussi la méthode des kmoyennes afin de réaliser une segmentation du coeur sur des images de perfusion. Cependant, l'agrégation ne s'opère pas sur les cinétiques elles-mêmes, mais sur des facteurs extraits des cinétiques par analyse factorielle, assurant une meilleure robustesse de la segmentation. Cependant, l'interprétation des facteurs d'une analyse factorielle devient difficile au-delà de 4 ou 5 facteurs, ce qui exclut de l'utiliser directement sur un corps entier. On notera que la méthode décrite par Frouin F. et colt ne segmente pas un organisme en pha'rmaco-organes, mais un pharmaco-organe en zones de prééminence de facteurs tels que les cinétiques artérielle, veineuse et tissulaire. Acton [Acton, 1999] utilise la méthode des c-moyennes floues, proche des k-moyennes, mais permettant de mieux prendre en compte la nature intrinsèquement floue des données acquises en tomographie par simples photons, tout en assurant une meilleure robustesse. Kimura et coll. [Kimura, 2002] proposent une méthode d'agrégation sur les projections des cinétiques sur les vecteurs propres associés aux deux plus grandes valeurs propres d'une analyse en composantes principales opérée sur l'ensemble des cinétiques de l'organisme, afin d'extraire des paramètres de modélisation compartimentale. Brankov et coll. [Brankov, 2003] proposent d'utiliser une mesure de similarité définie comme le cosinus de l'angle formé entre deux vecteurs représentés dans l'espace des cinétiques, plutôt qu'une norme euclidienne ou une norme de variance totale. Cette mesure de similarité 2. 5 présente une forte sensibilité au bruit dans les zones de faible rapport signal sur bruit. Brankov et coll. présentent deux algorithmes de type attentemaximisation ( EM en anglais), l'un flou et l'autre binaire, qui déterminent des agrégats à l'intérieur desquels les individus cinétiques présentent une forte similarité. La méthodologie EM permet d'estimer de façon itérative une variable cachée dont l'image est une réalisation particulière conformément à un modèle choisi. Chaque itération comprend une première phase d'attente où l'on calcule la vraisemblance attendue des données complètes à partir de la distribution jointe des données observées et cachées, et une seconde phase de maximisation où l'on estime les paramètres du modèle qui maximise cette vraisemblance du modèle attendue. Le processus est répété jusqu'à convergence de l'algorithme. Brankov compare sa méthode notamment à une application des mélanges d'analyseurs en composantes principales probabilistes. Cette méthode, introduite par Tipping et Bishop [Tipping, 1999], modélise le signal au sein de la zone à segmenter par un mélange de projections sur des sousespaces de l'espace des cinétiques. Un inconvénient majeur de toutes ces méthodes est qu'elles sont initialisées aléatoirement. A chaque lancement, l'algorithme converge vers l'un de ses minima locaux. Or, la solution recherchée correspond a priori au minimum global de cette énergie. Plusieurs lancements du programme avec des initialisations à chaque fois différentes sont donc nécessaires pour approcher de la meilleure solution. Guo [Guo, 2003] propose une agrégation par classification ascendante hiérarchique permettant d'obtenir un nombre d'agrégats défini a posteriori, mais aussi de conserver les agrégats de petite taille. Pour le calcul de la distance entre cinétiques, la valeur aux éléments de volume considérés est pondérée à chaque instant de la cinétique par le rapport signal sur bruit attendu pour l'image correspondante. Ce rapport, pour une image donnée de la séquence, est supposé dépendre essentiellement de la durée d'acquisition de cette image. L'approche hiérarchique assure la convergence de l'algorithme vers un minimum optimal, mais au prix d'un temps de calcul ne permettant pas d'examiner tout le volume de données. De plus, toute fusion de voxels obtenue par cet algorithme est définitive. En conséquence, une affectation erronée d'un voxel lors des premières itérations, par exemple du fait du bruit, ne peut être corrigée lors des itérations successives. Guo et coll. [Guo, 2003] proposent une pré-segmentation par histogramme, afin d'obtenir des premiers regroupements accélérant une classification ascendante hiérarchique. Cette dernière opère par une succession de fusions d'individus dans un ordre optimal selon un critère choisi. Les deux individus typiquement des voxels de l'image les plus proches en terme d'une distance choisie sont agrégés, puis les agrégats ou individus les plus proches sont de nouveaux agrégés, et ainsi de suite jusqu'à ce qu'un critère d'arrêt soit vérifié ou bien qu'il n'existe plus qu'un unique agrégat regroupant tous les individus. L'histogramme utilisé par Guo et coll. peut être décrit comme une comptabilisation du nombre de voxels ayant une intensité donnée. Cette classification par histogramme se base sur les valeurs de la dernière image au sens temporel de la série acquise après administration du traceur. On suppose que les premières fusions ont peu 1 o d'incidence sur les agrégats finaux et on regroupe les voxels correspondant au même intervalle de l'histogramme de la dernière image de la séquence. La variation de la distribution spatiale est supposée minimale pour cette dernière image parmi toutes celles contenues dans l'intervalle de temps considéré, le traceur ayant eu le maximum de temps pour se distribuer dans les organes suivant son affinité. Toutefois, dans le cas des isotopes à courte période radioactive utilisés en TEP, du fait de la décroissance exponentielle de la radioactivité au cours du temps, ceci entraîne l'inconvénient d'un bruit accru du fait que la dernière série présente aussi le plus faible rapport signal sur bruit de toutes les images de la série. Parmi toutes ces méthodes, seule celle de F. Frouin est validée sur des organes en mouvement de type mouvements périodiques de faible période. Cependant, de par son principe, elle n'est applicable qu'à des zones de l'organisme comprenant très peu de pharmaco-organes. Toutes les autres méthodes précitées n'ont pas été validées dans le cas d'une problématique corps entier, et ne sont pas adaptées aux mouvements physiologiques propres à cette problématique. Le document de Brevet EP-A-1 365 356 présente une méthode de segmentation semi-automatique d'images acquises par TEP, qui requiert notamment la délinéation préalable d'une région d'intérêt et de voxels-modèles à extraire dans les images. On notera que la méthode présentée dans ce dernier document est limitée au domaine de l'oncologie et qu'elle ne permet pas de segmenter les images en autant de régions d'intérêt Ii que de pharmaco-organes, mais seulement en deux parties dont une seule contient une tumeur. Un but de la présente invention est de proposer une méthode de segmentation d'une image ou d'une séquence d'images tridimensionnelles de départ à base de voxels pour l'obtention d'une image tridimensionnelle de, segmentation comportant une partition en régions d'intérêt, ladite image ou séquence d'images comprenant des mesures, pour chaque voxel et au cours de n intervalles de temps (n 1), de l'évolution réelle d'un signal représentatif d'au moins une variable de type physique, chimique ou biologique de ladite image ou séquence, qui permette de remédier à ces inconvénients. A cet effet, la méthode de segmentation selon l'invention comprend essentiellement les étapes suivantes: a) une modélisation du signal comprenant la définition d'un modèle paramétrique d'évolution spatiotemporelle dudit signal, ce modèle comprenant des jeux de paramètres de sorte que lesdits jeux soient respectivement propres à des structures correspondant auxdites régions d'intérêts et que chaque jeu de paramètres soit indépendant des coordonnées 20 spatiales dans la structure correspondante (ces paramètres seront dits homogènes ci-après) ; b) une extraction d'échantillons de voxels de sorte que lesdits échantillons soient respectivement inclus dans lesdites structures, cette extraction comprenant: (i) un calcul, pour chaque voxel de ladite image ou séquence d'images de départ ou bien d'une zone d'intérêt de celle-ci, d'un critère de validité d'une hypothèse selon laquelle ledit modèle d'évolution de ladite variable au sein du voisinage de ce voxel est propre à une et une seule desdites structures, (ii) une extraction de voxels-modèles qui sont définis comme réalisant les maxima locaux dudit critère, 30 (iii) une définition, pour chaque voxel-modèle, de l'un desdits échantillons de voxels inclus dans ladite structure correspondante, de manière à ce que cet échantillon présente une évolution de ladite variable qui soit conforme à celle du modèle de la structure à laquelle appartient ledit voxel-modèle, puis (iv) une estimation, sur chaque échantillon ainsi défini, des paramètres dudit modèle d'évolution de la structure dans laquelle ledit échantillon est inclus, puis c) une fusion desdits échantillons regroupant les échantillons dont le modèle d'évolution est propre à la même structure, ladite fusion suivant, précédant ou incluant une classification de tous les voxels de ladite image ou séquence d'images ou d'une zone d'intérêt de celle-ci par agrégation à un groupe d'échantillons, de manière qu'une similarité maximale existe entre l'évolution de ladite variable pour ces voxels et l'évolution dudit modèle caractérisant ce groupe d'échantillons. Cette invention s'applique aussi bien aux images bidimensionnelles, qui peuvent être des cas particuliers d'images tridimensionnelles dans lesquelles une des trois dimensions adopte la valeur un. Comme indiqué cidessus, nous utilisons le terme tridimensionnel pour désigner les images bidimensionnelles, qui sont un cas particulier d'images tridimensionnelles, aussi bien que les images tridimensionnelles. II est essentiel de noter que la méthode selon l'invention permet une segmentation automatique en structures ou régions d'intérêt (e.g. des pharmaco-organes), qui est uniquement basée sur l'information (e.g. pharmacocinétique) présente dans une série séquentielle d'images, notamment grâce à l'initialisation non aléatoire qu'elle comprend à l'étape b) précitée où l'on extrait les paramètres des modèles des structures d'intérêt dans des zones positionnées automatiquement au coeur de ceux-ci. Cette initialisation est adaptée aux données à traiter. Cette méthode permet ainsi d'obtenir systématiquement le même résultat sur un même jeu de données avec les mêmes paramètres, du fait de son caractère déterministe. Par voxel , on entend dans la présente description un élément unitaire de l'image volumique 3D , qui est généralement choisi cubique. C'est le plus petit volume spatial dont on dispose au sein de l'image. Les pixels sont des éléments d'images bidimensionnelles et peuvent être considérés comme des cas particuliers de voxels. Par la suite, nous utiliserons le terme voxel pour englober les termes voxel et pixel . Par région , on entend une zone connexe de l'image, i.e. en un seul morceau. Deux régions seront dites connexes si elles se touchent. De préférence, ladite image ou séquence d'images de départ est obtenue par une technique d'imagerie choisie dans le groupe constitué par la tomographie par émission de positons (TEP), l'imagerie par résonance magnétique (IRM), la tomographie d'émission de photons (e.g. SPECT ), les imageries optiques, le scanner aux rayons X, l'imagerie histologique, l'imagerie autoradiographique, l'imagerie satellitaire et l'imagerie photographique. D'une manière générale, on notera que la technique d'imagerie utilisée peut être aussi bien de type 2D que 3D . A titre encore plus préférentiel, ladite séquence d'images est obtenue par la technique TEP. Dans ce cas, on met avantageusement en oeuvre à l'étape b) des analyses en composantes principales locales dans un espace à n dimensions, soit au voisinage d'un voxel, soit en un échantillon, pour, d'une part, calculer le critère de validité de l'hypothèse selon laquelle ledit modèle d'évolution au sein du voisinage dudit voxel est propre à une unique structure et, d'autre part, estimer à partir d'un échantillon les paramètres du modèle d'évolution de la structure dans laquelle cet échantillon est inclus. Une analyse en composantes principales (ACP) opérée sur un tel ensemble de mesures multidimensionnelles comprend un changement de repère, tel qu'un minimum d'axes rend compte d'un maximum de la variance du signal. On distingue les vecteurs propres ou vecteurs directeurs des axes dégagés par l'ACP, et les valeurs propres qui sont respectivement associées à ces axes et qui représentent chacune la variance du signal le long de l'axecorrespondant. En général, on trie les axes ou vecteurs propres par valeur propre décroissante. La reconstruction d'une mesure est représentée par les k premiers vecteurs propres de l'ACP, et on réalise la projection de cette mesure sur le sous-espace déterminé par les k vecteurs propres associés aux k plus grandes valeurs propres de l'ACP. Le résidu de la projection est égal à la mesure elle-même à laquelle est soustraite la projection de cette mesure sur l'ACP, et ce résidu, représentatif du bruit relatif aux mesures, correspond à une projection sur les axes de faible variance. La variance du signal non reconstruit au sein de l'ensemble de mesures est égale à la moyenne des plus petites valeurs propres de l'ACP non prises en compte pour la reconstruction du signal. Si le nombre de valeurs propres retenues est juste suffisant pour reconstruire le signal au sein de l'ensemble, la variance de signal non reconstruit représente la variance du bruit. Si ce nombre de valeurs propres est trop petit, la variance de signal non reconstruit contiendra à la fois le bruit et du signal. On définit en outre comme suit la variance de signal non reconstruit relative au signal. Si l'on divise l'écart type du signal non reconstruit par la norme de la moyenne du signal au sein de l'ensemble, on obtient le rapport signal sur bruit lorsque k est juste assez grand pour que les k premiers vecteurs propres de l'ACP reconstruisent le signal. Dans une zone sans signal, alors que la variance de signal non reconstruit sera identique à ce qu'elle serait dans une zone où le signal est bien reconstruit, la variance de signal non reconstruit relative au signal sera très forte. On définit également la proportion de signal non reconstruit comme étant égale à la norme du résidu de reconstruction de la mesure par 30 l'ACP, divisé par la norme de cette même mesure. Dans son application particulièrement avantageuse à la pharmaco-imagerie, on notera que la méthode selon l'invention répond aux exigences de l'imagerie corps entier. En effet, lesdites images de départ peuvent être avantageusement des images d'un organisme entier humain ou animal, l'image de segmentation segmentant ledit corps en une partition de pharmaco-organes dont les contours respectifs délimitent les régions d'intérêt. Dans cette application, on notera également que la méthode selon l'invention répond en outre aux exigences de l'imagerie corps entier en mouvement, dans la mesure où cette méthode est applicable à un organisme entier animé de mouvements physiologiques de type périodiques tels que la respiration et les battements du coeur, dont la période est réduite en comparaison de la durée d'acquisition de l'une des images de départ de ladite séquence, ou bien à des mouvements physiologiques non périodiques. Ainsi, une cinétique de concentrations pourra être constituée d'un mélange des pharmaco-organes qui auront traversé l'élément de volume où cette cinétique est lue. La cinétique lue à un instant donné contiendra le même ratio de chaque pharmaco-organe que la même cinétique à un autre instant. De plus, le modèle choisi à l'étape a) permet de prendre en compte ce type de mouvements périodiques lors de l'étape d'extraction des voxelsmodèles, excluant les zones de l'image traversées par plusieurs pharmacoorganes au cours du mouvement périodique. Quant aux mouvements non périodiques ou de période de l'ordre de la durée d'acquisition d'une image de la séquence (e.g. digestion, dilatation de la vessie, mouvement globaux du corps), aucun voxel-modèle ne pourra être extrait dans les zones affectées par de tels mouvements, du fait que les cinétiques d'aucune de ces zones ne pourront satisfaire localement le modèle à un seul organe. Ces zones seront donc exclues de la phase critique de détermination des modèles des pharmaco-organes (voir la figure 2 qui illustre des cinétiques relatives à des mouvements non périodiques de vessies de rats). 1 6 On notera que la présente invention permet de s'affranchir de la limitation qu'impose une définition purement anatomique des organes, tout en restant compatible avec les techniques d'imagerie bimodales. Par pharmaco-organe , on entend de manière connue dans la présente description une parcelle de l'organisme présentant une réponse identique à des traceurs, de sorte que les pharmaco-organes ne sont pas totalement identifiables aux organes. Par exemple, l'organe rein pourra être dissocié en au moins deux pharmaco-organes: le cortex du rein et le bassinet, dans lesquels les traceurs présentent des cinétiques du fait du temps nécessaire à la filtration du sang par le cortex du rein pour former l'urine dans le bassinet. Par traceur , on entend de manière connue dans la présente description une molécule entrant dans les mécanismes biologiques étudiés et marquée de façon à pouvoir être suivie dans l'organisme. Cette molécule est injectée par voie intraveineuse à dose traceuse , c'est-à-dire suffisamment importante pour être repérée, mais suffisamment faible pour ne pas perturber le fonctionnement de l'organisme. Un tel traceur, alors marqué par un isotope radioactif, est notamment utilisé dans les techniques d'imagerie moléculaire, notamment la tomographie d'émission de positions (TEP), la gammatomographie ( Single-Photon Emission Computed Tomography , SPECT en abrégé), l'imagerie optique, l'imagerie par IRM et l'imagerie par ultrasons. A titre d'exemple de variable de type physique, chimique ou biologique caractérisant ladite séquence, on utilise avantageusement la concentration radioactive en un instant donné to à tn d'au moins un principe actif marqué et injecté dans ledit organisme, l'ensemble des voxels à l'intérieur de chaque pharmaco-organe présentant des cinétiques biochimiques de répartition dudit principe actif qui sont similaires. Selon d'autres caractéristiques de la méthode de l'invention: - on admet qu'il existe dans chaque image de départ de ladite séquence une partition de l'espace en un nombre fini desdites structures, lesquelles sont chacune connexes et présentent chacune en leur sein un comportement homogène en réponse au phénomène étudié dont l'évolution de ladite variable est représentative; et - on détermine le nombre desdites structures a posteriori. Selon d'autres caractéristiques avantageuses de l'étape a) de modélisation selon l'invention: - ledit modèle comprend en outre des paramètres hétérogènes dépendant des coordonnées spatiales des voxels au sein d'une même structure, et l'on admet que lesdits paramètres homogènes et hétérogènes peuvent être estimés sur un ou plusieurs éléments de volume inclus dans cette même structure; - on admet en outre que lesdites structures présentent entre elles des réponses différentes audit phénomène étudié, à moins qu'elles ne soient pas connexes; - on admet en outre que le bruit, qui est inhérent auxdites mesures et qui est dû au mode d'acquisition de ladite séquence d'images de 20 départ et de segmentation, est additif; - on se fixe en outre les deux contraintes suivantes, pour prendre en compte des mélanges locaux de différentes évolutions temporelles dudit signal: (i) si la totalité des voxels qui sont voisins d'un voxel et qui contribuent à l'évolution de ladite variable relative à ce voxel est incluse dans la même structure, alors ledit signal correspondant est une réalisation du modèle de cette structure, et (ii) ledit jeu de paramètres homogènes pour ces voxels voisins est le même que celui dudit modèle propre à cette structure; et - dans le cas d'une technique d'imagerie avec injection de telle que la TEP ou la SPECT , l'imagerie optique, l'imagerie par t l'imagerie par ultrasons, ledit modèle peut être un modèle ls compartimentai de type à un ou plusieurs traceurs indépendants et à plusieurs compartiments, tels que des compartiments biologiques (on suppose que les traceurs se répartissent dans les compartiments, et ce modèle étudie la concentration du traceur dans chaque compartiment via des vitesses de transfert entre les compartiments). Selon un exemple préférentiel de réalisation de l'invention, le modèle choisi est un modèle compartimentai à quatre paramètres, mais la présente méthode peut être étendue sans modification à un modèle compartimentai à six paramètres. Cette modélisation suppose les organes sains. En cas d'organes présentant des zones pathologiques, ces dernières feront l'objet d'un modèle séparé, et seront donc segmentées séparément de la partie saine de l'organe. Les zones pathologiques des organes pourront donc être mises en relief par la présente méthode. On notera que la modélisation du signal qui est mise en oeuvre à cette étape a) passe par une étude des mécanismes essentiels du phénomène en jeu, par exemple en introduisant des paramètres physiologiques comme la concentration de traceur dans le plasma ou le tissu, ou encore la fraction de volume sanguin, dans le cas de la pharmaco-imagerie (par ( Ledit modèle utilisé à l'étape a) est notamment conçu pour permettre de calculer des expressions analytiques des mesures attendues, conformément au phénomène étudié. La confrontation de ces expressions analytiques avec les mesures des effets du phénomène permet d'en estimer les paramètres (e. g. pour le phénomène activation de neurone)), la consommation de glucose peut être mesurée indirectement et partiellement via la technique TEP, par injection puis acquisition de molécules de fluoro déoxy glucose marquées par un isotope émetteur de position). On notera que ces mélanges locaux de différentes évolutions temporelles du signal traduisent un effet dit de volume partiel , dû à la faible résolution intrinsèque de la caméra utilisée. Le signal mesuré en un voxel donné contient alors un mélange des cinétiques physiologiques des régions d'intérêt voisines (e.g. des pharmaco-organes voisins). On notera également que la méthode selon l'invention permet de déterminer de façon fiable, à l'effet de volume partiel près, les modèles d'évolution correspondant à chaque structure d'intérêt, tels que des modèles de cinétiques d'un pharmaco-organe. Ces modèles sont en effet estimés dans les zones de l'image où le modèle de cinétique à un seul organe est le plus valable. Il est donc inutile d'affiner les paramètres des modèles des pharmaco-organes au cours de multiples itérations: une simple agrégation de chaque élément de volume au modèle qui représente le mieux sa cinétique suffit. On choisit dans la présente méthode un modèle de bruit additif gaussien non stationnaire dans le temps. Le bruit est estimé sur l'ensemble des échantillons extraits lors de la segmentation, et il est corrigé de l'influence du signal dans le cas d'une séquence d'image reconstruite par une méthode itérative. La méthode selon l'invention prend en compte: - la non-stationnarité spatiale du bruit au sein des images reconstruites par une méthode itérative, en corrigeant les données de l'influence du signal lors de chaque étape de notre algorithme, et -la non-stationnarité temporelle du bruit, ce qui permet de donner aux différentes images de la séquence des pondérations différentes rendant compte de la variabilité réelle des données autour du modèle. Selon d'autres caractéristiques de l'étape b) d'extraction de 30 voxelsmodèles selon l'invention: - cette étape b) comprend une sélection préalable, dans un espace à n dimensions (n 1) correspondant à ladite séquence d'images de départ (par exemple n=38 pour une séquence d'images TEP données), de voxelsmodèles pour que chacune desdites structures d'intérêt contienne en son coeur au moins un voxel-modèle ainsi que le voisinage de celui-ci; - cette étape b) comprend, suite à ladite sélection préalable, la définition d'une métrique prévue pour définir les distances dans ledit espace et ledit critère d'extraction local; et - on admet à l'étape b) que pour toute structure à segmenter, il existe un voxel dont le voisinage est compris au sein de la région d'intérêt correspondante. Selon une autre caractéristique de l'invention, cette étape b) comprend: une détermination par exemple itérative, pour chaque structure d'intérêt, d'un échantillon de voxels (i.e. une liste de voxels) appartenant à cette structure, cette appartenance étant déterminée par une similitude de l'évolution de ladite variable avec le modèle relatif à cette structure, puis - une estimation, pour chaque échantillon, des paramètres homogènes du modèle qui correspondent à la structure. On entend par distance dans la présente description la distance entre deux mesures de l'évolution temporelle dudit signal (e.g. mesures de cinétiques dans le cas de mesures de concentrations radioactives par la technique TEP, notamment). En général on utilise une distance pondérée de Minkowski, dont les distances euclidiennes, de Manhattan et distance maximum sont des cas particuliers. Le facteur de pondération, qui est choisi pour compenser une éventuelle non-stationnarité du signal, est en général pris égal à 1 pour chaque valeur de mesure (e.g. un instant déterminé de la cinétique des concentrations). 3 0 On entend par seuillage d'une image l'opération qui consiste à sélectionner les voxels dont la mesure dans l'image est supérieure à un seuil et, par extension, l'opération où l'on sélectionne les voxels dont la mesure dans l'image est comprise entre deux seuils: un seuil bas et un seuil haut. Un seuillage peut donc consister à sélectionner les voxels dont la mesure dans l'image est égale à une valeur donnée. Selon une autre caractéristique de l'invention, on notera que la bordure de chaque structure est exclue de l'extraction desdits échantillons qui est mise en oeuvre à l'étape b) précitée selon l'invention. Selon une autre caractéristique de l'invention, l'étape c) de fusion et de classification comprend une fusion desdits échantillons correspondant à une même structure d'intérêt, avant ou après ladite classification des voxels. Cette fusion est mise en oeuvre par exemple par une classification ascendante hiérarchique. Quant à la classification des voxels qui est également comprise dans l'étape c), elle est par exemple mise en oeuvre par une analyse discriminante, étant entendu que toute méthode de classement pourrait convenir pour cette classification. A l'issue de cette étape c) de fusion et de classification, on obtient une image de label indiquant à quel groupe d'échantillons un voxel donné est apparenté. Bien que l'on puisse se satisfaire de cette image, il est nécessaire, pour vérifier la définition d'une image de segmentation, de séparer chaque composante connexe pour chaque label. Selon une autre caractéristique de l'invention, ladite méthode de segmentation comprend en outre, suite à l'étape c), une étape de séparation de composantes connexes qui est mise en oeuvre pour chacun desdits labels. On notera que cette étape de séparation des composantes connexes permet de séparer entre elles des structures non connexes qui présentent des évolutions similaires de ladite variable, et que cette étape de séparation est nécessaire pour passer d'une image de classification à une image de segmentation en régions d'intérêt. Selon une autre caractéristique avantageuse de l'invention, ladite méthode comprend en outre une étape optionnelle de pré-segmentation qui est mise en oeuvre avant ou après l'étape a) et qui consiste à séparer lesdites images de départ en une première partie contenant lesdites structures d'intérêt et en une seconde partie correspondant à un fond d'image exclu de la segmentation. On notera que cette étape de pré-segmentation permet non seulement d'accélérer la mise en oeuvre de la méthode selon l'invention, mais encore d'éviter la pollution de la segmentation de la zone d'intérêt de l'image par des cinétiques non significatives de structures. Selon une autre caractéristique avantageuse de l'invention et indépendante de la caractéristique optionnelle précédente (i.e. non conditionnée par celle-ci), ladite méthode comprend en outre une autre étape optionnelle d'optimisation des échantillons qui est mise en oeuvre suite à l'étape b) et avant l'étape c) et qui consiste à rechercher par un algorithme par exemple itératif, pour chaque échantillon, une forme par exemple géométrique, paramétrique ou libre qui minimise en son sein le signal d'autres structures que celle où le voxel-modèle a été extrait. On utilisera par exemple pour cette étape d'optimisation des échantillons une croissance de région , i.e. une agrégation des voxels voisins des agrégats à ces derniers, par ordre décroissant de similitude avec ledit modèle. On propose, pour exemple, l'algorithme suivant permettant 25 d'extraire un échantillon y/; minimisant en son sein le signal provenant d'autres structures que i. A la première itération, l'échantillon y, recherché est égal à l'ensemble des cinétiques contenues dans Vj. A chaque itération, le modèle de la structure i est estimé au 3( sein de la forme yr;. Ensuite, une croissance de région est opérée à partir de F ment de volume j, en maximisant à chaque nouvelle agrégation un critère c similarité entre la cinétique de l'élément de volume considéré et le modèle estimé de la structure i. La croissance de région cesse lorsque le nombre d'éléments de volume inclus dans la région réalisée devient suffisant pour estimer les paramètres {O (i)}, Le processus itératif cesse lorsque yr; cesse d'évoluer d'une itération à l'autre. On notera que cette étape optionnelle d'optimisation des échantillons améliore la qualité et la robustesse de la segmentation obtenue de chaque structure détectée, car elle assure que l'échantillon est bien inclus dans la structure correspondante. Selon une autre caractéristique avantageuse de l'invention et indépendante de chacune des deux caractéristiques optionnelles précédentes, ladite méthode comprend en outre une étape optionnelle de représentation hiérarchique desdites régions d'intérêt obtenues qui est mise en ceuvre suite à l'étape c) et avant ladite étape de séparation des structures non connexes. On notera que cette étape optionnelle de représentation hiérarchique des organes permet à l'utilisateur de choisir le niveau de fusion optimal pour une structure donnée. Un niveau de fusion trop élevé ne permettrait pas de détecter la structure, qui serait fusionné avec une autre, tandis qu'un niveau de fusion trop bas scinderait la partie correspondante de l'image en deux zones représentant la même structure. Cette représentation hiérarchique permet, de plus, de choisir a posteriori le nombre de structures au sein de l'image. La présente méthode propose de conserver toutes les images de segmentation pour chaque niveau de fusion et de laisser le choix à l'utilisateur, e.g. pour un pharmaco-organe donné, du niveau de fusion où il est le mieux segmenté. On notera que cette méthode fusionne non pas les éléments de volume, mais les échantillons via les modèles d'évolution (e.g. de cinétiques de pharmaco-organes), dont le nombre est avantageusement plus faible de quatre ordres de grandeur. La navigation entre les divers niveaux de fusion est donc possible et peut être interprétée comme une recherche du niveau de fusion qui est tel que tous les modèles correspondant au même pharmaco-organe soient fusionnés. L'identification entre agrégat et région d'intérêt (e.g. pharmaco-organe) devient donc totale. D'une manière générale, on notera que la méthode selon l'invention est applicable à la segmentation de tous types d'images ou de séquences d'images de départ, e.g. des images satellites. D'une manière générale, on notera que la méthode selon 15 l'invention permet de traiter une grande masse de données. En effet, le temps de calcul reste compatible avec les besoins des utilisateurs, même dans le cas d'acquisition de données par une caméra de type TEP à haute résolution spatiale. La complexité des étapes critiques de l'algorithme extraction des voxels-modèles et segmentation est en effet proportionnelle au nombre d'individus à classer. Plusieurs millions d'individus sont ainsi classés sans que l'augmentation de ce nombre d'individus ne nécessite d'effectuer d'approximations lors de la segmentation. De plus, la présente méthode ne nécessite pas d'hypothèses sur la prévalence de certaines images de la séquence sur d'autres en termes de stabilité des phénomènes physiologiques. On notera également que la méthode selon l'invention bénéficie d'un coût en temps de calcul qui est fixe et est notamment compatible avec le temps de traitement de données issues de la technique TEP. On notera en outre que la présente méthode permet de segmenter des régions de tailles très différentes, en ne faisant aucune hypothèse sur la taille des structures à segmenter, telles que des pharmaco-organes. Les caractéristiques précitées de la présente invention, ainsi que d'autres, seront mieux comprises à la lecture de la description suivante de plusieurs exemples de réalisation de l'invention, donnés à titre illustratif et non limitatif, ladite description étant réalisée en relation avec les dessins joints, parmi lesquels: la figure 1 est un diagramme à blocs illustrant les différentes étapes obligatoires ou optionnelles de la méthode de segmentation selon l'invention, la figure 2 est un graphique illustrant notamment les familles de cinétiques non physiologiques qui sont générées par le déploiement de 15 vessies de rats, dans une image acquise par la technique TEP, les figures 3a et 3b sont des schémas illustrant un exemple de modélisation compartimentale à un seul traceur d'un processus physiologique d'interaction entre le traceur et un organe, via trois compartiments et quatre paramètres, la figure 4 est un histogramme illustrant, dans le cas d'une pré-segmentation selon l'invention, le nombre réduit de valeurs négatives d'une cinétique obtenue par la technique TEP qui est relative à un élément de volume inclus dans un organisme de rat, la figure 5 est un histogramme illustrant, dans ce cas optionnel de pré-segmentation, le nombre plus élevé de valeurs négatives d'une cinétique relative à un élément de volume extérieur à cet organisme de rat, les figures 6a et 6b illustrent deux analyses d'histogrammes de pré-segmentation illustrant l'extraction de l'indice du minimum d'un histogramme situé entre les deux pics correspondant respectivement à l'organisme de rat de la figure 4 et au fond d'image de la figure 5, cet indice correspondant à une valeur de seuillage permettant de séparer les éléments de volume de l'organisme de rat de ceux du fond, les figures 7a et 7b sont deux images masquées obtenues par pré-segmentation à partir d'oligo-nucléotides injectés à ce rat, ces deux 5 images correspondant respectivement aux figures 6a et 6b, les figures 8a et 8b sont deux images de coupes coronales d'organismes de rats obtenues lors de l'étape d'extraction de voxels-modèles selon l'invention, l'image 8a illustrant une coupe coronale passant par le coeur et les reins et l'image 8b illustrant une coupe coronale passant par la vessie, les figures 9a et 9b sont deux images de coupes coronales qui dérivent respectivement des figures 8a et 8b et qui illustrent les minima locaux des variances de données non reconstruites ayant été obtenus au coeur des organes par analyse en composantes principales, les figures 10a et 10b sont deux images de coupes coronales qui dérivent respectivement des figures 9a et 9b et qui ont été obtenues suite à l'étape d'extraction par une étape d'optimisation des échantillons selon l'invention mise en oeuvre par croissance de région, les figures 11 a et 11 b sont deux images de coupes coronales qui dérivent respectivement des figures 9a et 9b et qui ont été obtenues suite à une étape de classement des cinétiques selon l'invention en cent classes, la figure 12 illustre un exemple d'arbre de fusion obtenu par une étape de fusion des échantillons selon l'invention qui fait suite à l'étape de classement des cinétiques et qui est mise en oeuvre par une classification ascendante hiérarchique, les figures 13a et 13b sont des images de coupes coronales qui dérivent respectivement des figures 11 a et 11 b et qui ont été obtenues par fusion des échantillons via une classification ascendante hiérarchique, où dix-huit classes ont été conservées, les figures 13c et 13d sont des images de coupes coronales qui dérivent respectivement des figures 13a et 13b et qui ont été obtenues par des labels de l'image de classement des cinétiques via une cation ascendante hiérarchique, où dix- huit classes ont été conservées, les figures 14a, 14b et 14c sont des images 3D de segmentation respectivement obtenues par cette étape de fusion selon des vues coronale, sagittale et oblique, et la figure 15 est un diagramme à blocs illustrant une 5 représentation hiérarchique des pharmaco-organes selon l'invention qui peut être mise en oeuvre suite à cette étape de fusion des échantillons. Dans l'exposé suivant de la méthode selon l'invention, on a utilisé les notations suivantes: É i: un organe, un label ou un échantillon et l: le nombre d'échantillons. É o: un organe. É t et t' : des instants ou des indices d'images de la séquence. É 4, : la durée d'acquisition de l'image t de la séquence. É j et k: des éléments de volume (voxels). É s: un traceur injecté et S le nombre de traceurs injectés. É r: un rayon. É C; : la variation de la concentration au sein de l'élément de volume j au cours du temps, que l'on nommera par la suite cinétique. É : le voisinage de l'élément de volume j. É (p; : un échantillon de cinétiques représentatives du modèle de la structure i. É Cp(t) : la cinétique plasmatique, constante dans l'organisme. É Ca(t) : la cinétique artérielle, constante dans l'organisme. É Cf,(t) : la cinétique du compartiment libre. É Cb,(t) : la cinétique du compartiment lié. É Ct,(t) : la cinétique au sein du tissu: Ct,(t) = Cf,,(t) + Cb,(t) . É M;.j(t) : la cinétique prédite par modèle à un seul organe, l'organe i au sein de l'élément de volume j. É (t) : le bruit de détection et de reconstruction à l'élément de volume j. É a: l'écart type du bruit. É K,, k2, k3 et k4: les paramètres d'un modèle à trois compartiments. É Vb,1: le ratio de volume plasmatique au sein de l'élément de volume j de l'organe i. É ,u,,(t) : la cinétique moyenne du modèle M,.j(t) au sein de l'organe i. É /Jj,k: le coefficient de contamination de la cinétique de l'élément de volume j par celle de l'élément de volume k, et inversement. É i(k) : la notation signifiant que i dépend de k. Il s'agit ici de l'organe i physiologiquement présent dans l'élément de volume k. É ni: le nombre de valeurs nulles ou négatives de la cinétique Cj(t) à l'élément de volume j. É Modefond: le pic de l'histogramme correspondant au fond. É Modeorg: le pic de l'histogramme correspondant à 20 l'organisme. É nmin: m'indice du minimum de l'histogramme situé entre Modeorg etModefofd, correspondant à une valeur de seuil permettant de séparer les éléments de volume de l'organisme imagé de ceux du fond. É Bj,r: la boule de rayon r centrée en j. É fr;(t) : la moyenne des cinétiques Ck au sein de la boule B;.,. É e;,/(t) : le terne vecteur propre et 2y,, : la ferre valeur propre de l'analyse en composantes principales calculée sur les cinétiques des éléments de volumes inclus dans la boule Bj,r. La méthode selon l'invention, dont les principales étapes obligatoires ou optionnelles sont illustrées à la figure 1, est préférentiellement mise en oeuvre via la technique d'imagerie TEP. On notera toutefois que cette méthode de l'invention n'est nullement limitée à cette technique. Le diagramme de la figure 1 illustre les huit étapes de cette méthode, cinq de ces étapes 10, 30, 50a, 50b et 70 étant obligatoires (symbolisées par des cases en traits pleins) et trois autres étapes 20, 40 et 60 étant facultatives (symbolisées par des cases en pointillés). Une première étape 10, essentielle dans la méthode selon l'invention, consiste en une modélisation du signal au sein de l'image, via la définition d'un modèle adapté au phénomène étudié. Ce modèle est défini une fois pour toutes pour un phénomène donné, et n'a donc pas à être redéfini pour une nouvelle image rendant compte de ce phénomène. La donnée de ce modèle permettra par la suite non seulement de déterminer si un volume donné contient les cinétiques provenant d'une unique structure, mais aussi de définir une distance entre une cinétique et un échantillon de cinétiques supposé représentatif d'une structure donnée. Si nécessaire, une seconde étape optionnelle 20 de pré-20 segmentation de l'image séparera cette dernière en deux parties, l'une de ces parties correspondant à la zone étudiée, telle qu'un organisme avec la technique TEP, et l'autre partie correspondant au fond de l'image, jugé sans intérêt. Une troisième étape 30, également essentielle dans la méthode selon l'invention, consiste en une extraction de voxels-modèles au coeur de chaque structure d'intérêt, rendue possible par la mesure de la présence d'un unique modèle de cinétique de structure au sein du voisinage de chaque voxel de l'image. Une même structure peut alors contenir plusieurs voxels-modèles. Il est possible de conserver les voisinages définis ci-dessus des voxelsmodèles extraits comme échantillons de points situés dans les organes. Cependant, il peut être judicieux, pour un voxel-modèle donné, de mettre en oeuvre une quatrième étape optionnelle 40 d'optimisation des échantillons dans laquelle on recherche un échantillon dont la forme est plus adaptée à celle de la structure au sein de laquelle le voxel-modèle a été extrait. Ce nouvel échantillon peut par exemple être extrait par une méthode connue de croissance de région. On met ensuite en oeuvre deux étapes 50a et 50b toutes deux obligatoires, selon l'un ou l'autre des deux modes de réalisation selon l'invention décrits ci-après. Selon un premier mode de réalisation, on met en oeuvre l'étape 50a de classement des cinétiques suite à l'étape 30 ou 40, dans laquelle chaque élément de volume de l'image est agrégé à l'échantillon dont le modèle est le plus proche de celui de cet élément. Puis on met en oeuvre l'étape 50b de fusion des échantillons, où les échantillons appartenant à un même organe sont fusionnés. Selon un second mode de réalisation, on met en oeuvre l'étape 50b de fusion des échantillons suite à l'étape 30 ou 40, puis on met en oeuvre l'étape 50a de classement des cinétiques. Si l'étape 50b de fusion utilise une méthode hiérarchique, cette hiérarchie de fusion peut être utilisée pour décrire les structures dans une étape optionnelle 60 de représentation hiérarchique des organes, une structure donnée apparaissant sur une plage de niveaux de fusion donnés. Dans un cas de figure idéal, l'intersection de toutes les plages d'apparition des structures est non vide et peut être extraite automatiquement. Dans un cas contraire, le choix peut être laissé à l'utilisateur, pour chacune des structures, de la plage de niveaux de fusion pour lequel elle apparaît le mieux. On obtient à cette étape 60 une image d'agrégation des éléments de volumes aux fusions d'échantillons dont ils sont les plus semblables. Cependant, comme deux structures peuvent présenter des réponses identiques au phénomène étudié tout en étant disjointes, on met en oeuvre préférentiellement une étape optionnelle 70 de séparation des composantes connexes, qui permet de séparer de telles structures. Au terme de la mise en oeuvre de cette méthode selon l'invention, on dispose d'une segmentation de l'image conforme au modèle de signal choisi. Plus précisément, la méthode de segmentation selon l'invention peut être décomposée comme suit. On étudie un phénomène dont les effets sont mesurables en chaque point de l'espace discrétisé image et dont la mesure comporte T valeurs (T > 1) ces valeurs pourront par exemple être les mesures du phénomène à des moments successifs. Dans la suite, on appellera cinétique la suite ordonnée des T valeurs de la mesure. Appelons S la mesure sans bruit de ce phénomène en chaque 15 point de l'espace. Appelons C la mesure réelle, soumise aux limitations du mode de détection utilisé. On va détailler ci-après les principales hypothèses et contraintes sur lesquelles se base cette méthode selon l'invention: Hypothèse H1: il existe une partition de l'espace en un nombre fini de structures présentant chacune en son sein un comportement homogène et telle que chacune de ces structures soit connexe. La présente méthode cherche donc à réaliser une telle partition au sein de l'image en se basant sur l'information contenue dans les mesures réelles C. Pour ce faire, nous commencerons par définir un modèle de signal basé sur la modélisation du phénomène lui-même. Ce modèle nous permettra d'effectuer une mesure de la présence de plusieurs structures dans un petit volume. Cette mesure est nécessaire à l'étape d'extraction des vc x:i:ls-modèles. D'autre part, la probabilité d'appartenance d'un élément de e à une structure sera déterminée par la similarité de la cinétique de cet t de volume avec le modèle de la structure. Etape 10 de modélisation du siqnal: Définissons un modèle paramétrique M(0,0), (p (i)) du 5 phénomène, dont certains paramètres {Ou(i)}1 u ne dépendent que de la structure étudiée i et sont par conséquent homogènes au sein de cette structure. Les autres paramètres {cp,,0)}, v du modèle M dépendent de la position j dans l'espace. Les paramètres N)v(j)}, La réponse d'une structure au phénomène étudié est caractérisée par l'ensemble des paramètres {O (i)}, u, dont la valeur est propre à la structure i et homogène au sein de celle-ci. M sera donc choisi de telle sorte que deux structures aux comportements différents présentent chacune un jeu de paramètres {0 (i)}, Hypothèse H2: les paramètres de M, peuvent être estimés sur un échantillon de mesure inclus dans i et de taille suffisamment grande. Hypothèse H3: les structures dont on cherchera à délimiter les contours présentent entre elles des réponses différentes au phénomène étudié, à moins qu'elles ne soient pas connexes. Deux structures connexes présentant la même réponse au 2.5 phénomène seront donc considérées comme une seule et même structure. Hypothèse H4: Le bruit de mesure dû au mode d'acquisition et de reconstruction de l'image est additif. On en conclut que la mesure C1 au sein de l'élément de 30 volume j peut s'exprimer comme: C; = S1 + si où le bruit E1 peut être non stationnaire tant d'un point de vue spatial que temporel. Par conséquent, si j est inclus dans i, C; s'écrit: C1(t) = M,1(t) + ì(t) . Une fois estimés les paramètres RN) , u du modèle de la structure i, il devient possible de calculer pour tout élément de volume j de 5 l'image l'erreur de reconstruction 0) = Cj(t) M;,j(t). Cette erreur de reconstruction est un indicateur de la présence de la structure i au voxel j. On va examiner à présent le problème des mélanges locaux de cinétiques. Dans le cas où Si n'est pas égal M;,j, mais à une fonction f connue des réalisations des modèles des structures i(vk) présents aux éléments de volume vk, les {Vk}, Dans le cas d'une fonction fj linéaire, M; doit dépendre linéairement de j. Dans le cas où f est un barycentre à coefficients positifs, M; doit être une fonction convexe de j. Etape 20 de pré-segmentation: Dans certains cas, l'image à segmenter contient un grand nombre d'éléments de volume dont la mesure n dimensionnelle correspond à une absence de signal. Ces éléments de volumes peuvent alors être regroupés en une même région. Le modèle de cinétique n'a pas de sens dans de telles régions, et les modèles estimés ne correspondront à aucun phénomène, si ce n'est le bruit. Toute analyse faite dans de telles régions est une perte de temps, qui peut être considérable lorsque leur taille devient non négligeable par rapport à la taille totale de l'image. Du fait de la nature aléatoire de la mesure C au sein de telles régions, on peut s'attendre à y extraire un grand nombre de voxels-modèles lors de l'étape 30 d'extraction. Les modèles extraits étant eux-mêmes aléatoires, ils seront fort différents les uns des autres et la probabilité pour qu'ils soient tous fusionnés dans une même région est faible. L'application de l'algorithme aux zones de l'image ne présentant pas de signal le ralentira et perturberont l'étape d'extraction 30 et l'étape 50b de fusion des échantillons. Il est donc avantageux d'isoler, lors de cette étape 20, les zones ne présentant pas de signal, que nous nommerons extérieur , des zones présentant un signal, que nous nommerons intérieur D. On notera que nombre d'opérations permettent d'effectuer une telle séparation, comme par exemple un seuillage. - Etape 30 d'extraction des voxels-modèles: Soit HOv l'hypothèse selon laquelle le voisinage V de l'élément de volume j ne contient que des cinétiques provenant du modèle d'une seule structure. Cette étape cruciale de la méthode selon l'invention va extraire des points dont le voisinage maximise un critère d'inclusion dans une structure. L'étape 10 de modélisation du signal présent dans l'image permet de définir ce critère comme la mesure de la validité de l'hypothèse HOvj. Si un modèle à une seule structure est incapable de capter le signal au sein du voisinage considéré, alors il est probable que ce dernier contient les modèles de cinétiques de plusieurs structures. Cette étape nécessite cependant une hypothèse supplémentaire, qui assure d'extraire des points à l'intérieur de chaque structure: Hypothèse H5: Pour toute structure, il existe un élément de 10 volume j dont le voisinage V; est compris au sein de l'organe. Selon cette hypothèse, toute structure i à segmenter contient une forme identique à celle de V connue et déterminée a priori. Le non-respect de cette condition engendrerait une incertitude sur la détection de cette structure. Dans le cas où le signal présent au sein d'un élément de volume j se calcule comme une fonction f des modèles des structures aux éléments de volume voisins, l'hypothèse H5 devient: Hypothèse H5' : Pour toute structure, il existe un élément de volume dont le voisinage est compris au sein de l'organe, et dont et tel que pour tout k appartenant à ce voisinage, le domaine d'application de fk est lui aussi inclus dans la structure. Indicateur de l'intérieur des organes: Soit une mesure r 1 de la validité de l'hypothèse HOvi au sein du voisinage V; de j comparable entre deux points de l'espace. L'éventuelle non-stationnarité du bruit -tant spatiale que temporelle - devra donc être prise en compte dans le calcul de cette mesure. Cette mesure Iv; est calculée pour chaque élément de volume de l'image ou de la zone d'intérêt extraite lors de l'étape 20. Extraction: Les fortes valeurs Tvi correspondent à des éléments de volume j dont chaque point du voisinage V; selon l'hypothèse H5 ou dont le domaine d'application de fk pour tout k appartenant à V; selon l'hypothèse H5' est inclus dans la même structure. En effet, l'hypothèse HOvj devient non valide au sein de BIr lorsque B r est à cheval entre deux structures i et o caractérisées par des réponses différentes au phénomène étudié. On extrait donc les maxima locaux de I vi, la mesure de validité de HOv, calculée pour chaque élément de volume j de la zone d'intérêt de l'image (un maximum local j d'une mesure T est tel que, pour tout k voisin de j, on aTvjjvk). Cette étape 30 de la méthode selon l'invention détectera des points dans le voisinage duquel seul un modèle de structure est réalisé. Dans le cas d'un mélange local de cinétique par une fonction f, les points pour lesquels l'hypothèse H5' n'est pas vérifiée, bien que le voisinage Vj soit inclus dans l'organe, n'ont donc pas la certitude d'être détectés. A fortiori, cela est vrai pour les structures dont la taille est de l'ordre de la résolution de l'image (et non du dispositif d'acquisition), ainsi que pour les structures fines et allongées ne présentant aucun renflement. Cependant, l'extraction de maxima locaux du critère permet d'extraire des points au sein de régions de taille ou d'épaisseur inférieure - mais non négligeable - à la taille du voisinage considéré. On notera toutefois que l'estimation du modèle dans de telles régions sera polluée par les modèles de cinétiques des structures voisines. Etape 40 d'optimisation des échantillons: Pour une géométrie de voisinage V, donnée et conformément au modèle M du phénomène étudié, les cinétiques contenues dans Vj constituent un échantillon localement optimal de cinétiques de la structure i contenant j. Les cinétiques aux éléments de volume contenus dans Vi peuvent donc être choisies comme échantillons de cinétique de la structure i. Cependant, certaines structures peuvent avoir une géométrie très différente de celle de V et ne pas disposer de renflements suffisants pour qu'il existe un voisinage Vi vérifiant H5 ou H5'. Le voisinage Vi entourant un voxel j extrait lors de l'étape 30 pourra donc servir d'initialisation à un algol 'thme itératif cherchant une forme qui minimise en son sein le signal d'autres structures que i. Cette forme pourra être géométrique (ellipse, cylindre), paramétrique (beta-splines) ou encore de forme libre. Selon la présente invention, un algorithme permettant d'extraire un échantillon yr; minimisant en son sein le signal provenant d'autres structures que i peut être le suivant. A la première itération, l'échantillon yr; recherché est égal à 15 l'ensemble des cinétiques contenues dans V. A chaque itération, le modèle de la structure i est estimé au sein de la forme y/. Ensuite, une croissance de région est opérée à partir de l'élément de volume j, en maximisant à chaque nouvelle agrégation un critère de similarité entre la cinétique de l'élément de volume considéré et le modèle estimé de la structure i. La croissance de région cesse lorsque le nombre d'éléments de volume inclus dans la région réalisée devient suffisant pour estimer les paramètres {04i)}, u du modèle M pour la structure i par exemple égal à celui contenu dans Vp yr; est alors calculée comme l'ensemble des cinétiques des 25 éléments de volume de la région réalisée. Le processus itératif cesse lorsque yr; cesse d'évoluer d'une itération à l'autre. On a ainsi déterminé un nouvel échantillon y; de cinétiques de la structure i, dont la forme optimise localement l'estimation du modèle M du 7 n5F: omène étudié pour la structure i. Etape 50a de classement des cinétiques: La définition des échantillons donne directement accès aux modèles de cinétiques des structures i qui composent l'image. La classification des cinétiques processus généralement itératif au cours duquel les paramètres du modèle sont estimés se réduit ici à un simple classement: chaque élément de volume de l'image est agrégé à la structure dont le modèle rend le mieux compte de sa cinétique. En cas de bruit non stationnaire d'un point de vue spatial, le bruit pourra être blanchi , à moins que la variance du bruit ne soit estimée au sein des échantillons. Elle peut encore être calculée analytiquement si la dépendance spatiale de la variance du bruit suit une loi connue a priori. Notons que le bruit pourra être estimé au sein des échantillons par la relation: 0-;'(1)_ C,(1) Etape 50b de fusion des échantillons: L'étape d'extraction 30 peut extraire plusieurs voxels-modèles au sein d'une même structure. Aussi, l'image d'agrégation résultante est peu exploitable en tant que telle, car un élément de volume appartenant à une structure est assigné de façon aléatoire à l'un des échantillons associés à cette même structure. On dispose de l échantillons, dont plusieurs peuvent se trouver au sein de la même structure. Les échantillons inclus dans la même structure doivent donc être fusionnés, afin que chaque classe fusionnée obtenue corresponde à une structure. Il convient de déterminer une mesure de similarité entre les modèles M; et g, estimés sur les échantillons i et i'. Deux échantillons i et o contenus dans la même structure seront caractérisés par des valeurs de paramètres {O (i)), u et {O (i')}, La fusion réunira prioritairement les échantillons dont la similarité entre leurs modèles respectifs est maximale. Une fusion 1 o d'échantillons sera considérée comme un nouvel échantillon, ce qui permettra de réitérer l'étape de fusion. Les paramètres {0,0)}I Il existe donc théoriquement une étape de fusion après laquelle à chaque structure détectée correspond exactement à un groupe d'échantillons fusionnés. Cette étape de fusion particulière peut par exemple coïncider avec une soudaine diminution de la mesure de similarité entre les modèles des deux échantillons fusionnés. Cette étape optimale de fusion peut donc être déterminée par une étude de la mesure de similarité ou encore de ses dérivées par rapport à l'étape de fusion. Etape 60 de représentation hiérarchique des orqanes: L'attribution d'un même label aux régions fusionnées par l'étape de fusion 50b permet théoriquement d'obtenir une image de classement dans laquelle chaque structure correspond exactement à une fusion d'échantillons, et ce pour une étape de fusion identique pour chaque structure étudiée. Toutefois, les réponses de certaines structures au phénomène peuvent être très proches sans être identiques, ce qui s 'curerait pas que les échantillons correspondant à des structures 30 différentes soient fusionnés après ceux qui correspondent à la même structure. Il n'existera alors pas de niveau de fusion tel que chaque structure ou chaque groupe de structures présentant exactement la même réponse au phénomène observé corresponde exactement à un échantillon. II est aussi envisageable que ce niveau de fusion optimal existe, mais qu'il soit impossible à déterminer à partir de l'étude de la mesure de similarité entre modèles de cinétiques d'échantillons à fusionner. Afin de pallier ces défauts, la présente méthode propose de représenter l'image sous la forme d'un arbre dont chaque feuille représente un agrégat obtenu lors de l'étape 50a de classement des cinétiques. Chaque noeud de l'arbre représente une fusion ou un agrégat issue d'une fusion, et les feuilles de fils représentent les agrégats qui composent cette dernière. Par une reconnaissance de structures utilisant un algorithme automatique ou bien encore par une reconnaissance manuelle, cette étape 60 se propose de déterminer le noeud de l'arbre correspondant exactement à une structure donnée. L'identification de chaque structure dans l'arbre permet d'obtenir une image de classement telle qu'à une structure corresponde une fusion d'échantillons, mais à une étape de fusion dépendant de la structure étudiée. Etape 70 de séparation des composantes connexes: D'après l'hypothèse H3 précitée, il est possible que deux structures aient la même réponse au phénomène étudié, mais soient non connexes. Les étapes 50a de classement et 50b de fusion, voire l'étape 60 de représentation hiérarchique, peuvent donc réunir sous un même label de telles régions. Une simple séparation des composantes connexes, classe par classe, permet de séparer les structures qui présentent la même réponse au phénomène mais qui sont non connexes, et permet d'obtenir une image de segmentation, résultat final de la méthode selon l'invention. Exemple de mise en oeuvre de la méthode selon l'invention pour la segmentation d'images issues d'une Tomographie par Emission de Positions (TEP) Bien que la problématique corps entier du petit animal, tel qu'un rat, soit abordée dans cet exemple de réalisation, ce dernier pourrait tout aussi bien être appliqué à des images de zones spécifiques, comme le cerveau, et chez d'autres animaux ou chez l'homme. En effet, la segmentation du corps entier est particulièrement problématique du fait des importants mouvements physiologiques affectant l'image à traiter. On a choisi une segmentation de type corps entier et chez un petit animal pour montrer le bon fonctionnement de la méthode malgré la petite taille des organes à segmenter et malgré l'occurrence de mouvements physiologiques. On montrera ici que le modèle de cinétique choisi vérifie les hypothèses de travail du principe général, et l'on décrira la méthode suivie comme une application directe du principe général énoncé ci-dessus. a) Etape 10 de modélisation du signal: Signal en TEP La Tomographie par Emission de Positons permet de suivre le devenir d'une ou plusieurs molécule(s) présentant des propriétés biologiques d'intérêt, que l'on aura marquée au préalable par un isotope radioactif émetteur de positons. Ce traceur est injecté par voie intra-veineuse au sujet imagé installé dans le champ de vue de la caméra. Lors de la désintégration dr, !' %tope, un positon est émis qui, après avoir perdu son énergie cinétique f tocs successifs avec des électrons, s'annihile avec un électron pour deux photons partant dans des directions diamétralement opposées. Un système de détection composé de couronnes de détecteurs permet de détecter les photons émis et de reconstituer la ligne sur laquelle s'est produite la désintégration. A partir de ces événements, un algorithme reconstruit une image rendant compte, en chaque endroit de l'espace discrétisé, du nombre de désintégrations ayant eu lieu dans cet élément de volume. Un pharmaco-organe peut être défini comme une structure de l'organisme dont tous les éléments présentent une réponse identique au traceur. Cela ne signifie pas pour autant que la cinétique mesurée au sein de tous les éléments de volumes d'un même organe soit forcément identique. On peut définir un modèle de cinétique M;,; permettant de rendre compte de la cinétique physiologique de chacun des éléments de volume du pharmacoorgane. Par cinétique physiologique, on entend le signal émis par les molécules marquées et non le signal mesuré. Ainsi, si l'on suppose le bruit additif et si l'on ignore l'effet de volume partiel, on peut décrire la cinétique Ci, pour un élément de volume j au sein d'un organe i comme: C;=M;,i+ì Si l'on prend en compte l'effet de volume partiel, qui mélange au voxel j les cinétiques physiologiques des pharmaco-organes i(k) au sein 20 des voxels k voisins de j, Ci s'écrit: = 1(/3,.k x Mrrti). ti)+ kerors) A,k est un coefficient de mélange déterminé par la fonction de réponse de I'imageur TEP. Le modèle de cinétique physiologique à un seul organe devra être choisi de façon à ce que, pour deux pharmaco-organes i et o présentant des réponses différentes au traceur, le nombre d'éléments de volume j tels que M,4(t) = Mo,;(t) soit faible. 3C Modélisation compartimentale à un seul traceur: Modélisons les processus physiologiques d'interaction entre le traceur et un organe i donné par trois compartiments et quatre paramètres. Ce traceur traverse la barrière cellulaire avec une constante de vitesse K,. Il peut de la même façon repasser la barrière pour retourner dans le sang avec une constante de vitesse k2. Il est transformé dans les cellules avec une constante de vitesse k3. Si cette transformation est réversible la transformation inverse est caractérisée par une constante de vitesse k4. Les trois compartiments étudiés sont Cp le compartiment plasma-sérum, Cf, le compartiment libre de l'organe i et Cb, le compartiment lié de l'organe i. Soit Cp(t), Cf,(t) et Cb,(t) les concentrations du traceur dans chacun de ces trois compartiments. Ce processus peut être schématisé comme illustré à la figure 3a. Dans ce modèle, la variation de concentration des composés dans les différents compartiments s'écrit: dÇy; (t) _ K,CP(t) (k.2 +k,)Cf (t)+ kaCb,(t) dt dCb (r) = k,C'f (r) k4Cb;(t) dt Dans un premier temps, plaçons-nous dans le cas où un seul traceur est injecté à l'organisme. Soit Ct,(t) la concentration d'un traceur dans le tissu dans l'organe i, Ct,(t) s'écrit: Ct,(t) = Cf,{t) + Cb,(t) Cp(t), la concentration du traceur dans le sérum plasmatique, n'est pas accessible directement. En effet, le compartiment du sérum plasmatique est inclus dans un autre compartiment, le compartiment artériel Ca(t). Soit donc Ca(t) la concentration du traceur dans le compartiment artériel. Comme illustré à la figure 3b, le compartiment comprend un compartiment Ccelt contenant les cellules du sang, ainsi que le compartiment plasmatique Cptasma, comprenant lui-même le compartiment du sérum plasmatique Cp entrant dansle modèle compartimentai choisi et les protéines du sang Cprot. 20 Le modèle M,;(t) cinétique du signal TEP au sein d'un élément de volume j donné compris dans l'organe i s'écrit: M (t)=((l Vb, )Ct,(t)+Vb, x Ca(t)) où Vb,,j est le rapport de volume artériel au sein de l'élément de volume j de l'organe i. Ca(t) est homogène dans l'organisme entier et Ct(t) est homogène au sein d'un pharmaco-organe sain donné. Dans le cas d'un pharmaco-organe au fonctionnement pathologique, celui-ci peut être séparé en une zone saine et une zone pathologique. Ces zones se comportant différemment par rapport au traceur, elles sont pharmaco-cinétiquement différentes et seront considérées comme deux pharmaco-organes différents. Ca(t) et Ct(t) étant donc homogènes au sein d'un pharmaco- organe donnés, M,,,(t) appartient au segment Ca(t)Ct,(t), et donc par extension à la droite Ca(t)Ct,(t), lorsque j varie tout en restant au sein de l'organe i. On peut donc modéliser la cinétique physiologique au sein d'un organe par une droite D, dans l'espace des cinétiques passant par Ca(t) et Ct,(t). La connaissance des modèles M,,/t) de cinétiques physiologiques pour quelques organes i devrait donc permettre de calculer Ca(t), comme l'intersection de toutes les droites D,. En revanche, Ct,(t) ne peut être explicitement calculée. Modélisation compartimentale à plusieurs traceurs On montre aussi (voir l'annexe A joint), que lorsque plusieurs traceurs s au nombre de S sont injectés à l'organisme, le signal résultant S,,s,;(t) dans l'élément de volume j décrit toujours une droite dans l'espace des cinétiques, cette fois la droite passant par,(t) et 1 ces différents traceurs n'interagissent pas entre eux. 0), à condition que Le modèle de cinétique reste donc valide dans le cadre du modèle compartimenta) à quatre paramètres, quel que soit le nombre de traceurs indépendants injectés. > Mélanges locaux de cinétiques: - Impact de l'effet de volume partiel sur le signal: L'effet de volume partiel peut être modélisé [Frouin, 2002] comme un lissage de l'image. Le signal S; mesurée dans l'image au sein d'un élément de volume j est donc une combinaison linéaire des cinétiques Sk des éléments de volume k voisins. Lorsque plusieurs traceurs sont injectés: S Si(r)=l s-11( AeI?rsmaeF(l) Mi(kL..k (1)) où /3 vérifie =1. AeI oisin aerl j) Sj(t) appartient toujours à une droite de l'espace des cinétiques, passant par,,(t) et Ca (() (voir Annexe A) lorsque le voisinage Vj de j dont les cinétiques physiologiques contribuent au signal Sj est inclus dans l'organe i. Lorsque ce voisinage contient plusieurs organes modélisés par des droites différentes, Sj(t) appartient à un sous espace de l'espace des cinétiques centré en Ca(t) et de dimension strictement supérieure à 1. - Impact des mouvements physiologiques sur le signal: Les mouvements physiologiques non périodiques ou dont la période ne peut être négligée par rapport à la durée d'acquisition d'une image de la série génèrent des cinétiques non interprétables. En effet, à chaque instant t, j ne contient pas les mêmes sites physiologiques. Les cinétiques des régions subissant un mouvement physiologique périodique de période négligeable par rapport à la durée d'acquisition d'une image de la série, comme la respiration ou les battements du coeur, redeviennent interprétables (voir annexe A). - Conclusion: On a un modèle de cinétique à un unique pharmaco-organe. Soit un volume A fixe dans le repère de la caméra. Soit un volume B fixe dans le repère d'un l'organe, placé au sein du pharmaco-organe i, de façon à ce que le voisinage de B (Vois(B)) sur lequel s'étend l'effet de volume partiel soit lui aussi inclus dans le pharmaco-organe i. Supposons A et B superposés au temps to: A ne contient en to que du signal issu de l'organe i. Soit Vois(A) le voisinage de A sur lequel s'étend l'effet de volume partiel. Lors d'un mouvement physiologique périodique de faible période, B se déplace par rapport à A, si bien que les deux volumes ne sont plus superposés. Cependant, le modèle M;,j(t) à un unique organe i est valide au sein de A, fixe dans le référentiel de la caméra, si et seulement si Vois(A) reste inclus dans l'organe i lors de la période du mouvement physiologique (voir Annexe A). Vérification des hypothèses du principe général: - Partition de l'espace en pharmaco-organes: 25 L'organisme peut être séparé en un nombre fini de pharmacoorganes, correspondant chacun à un organe, à un groupe d'organes connexes ou à une sous-partie d'un organe. Par définition du modèle compartimentai, chaque pharmaco-organe est caractérisé par une cinétique 2ç rielle et une cinétique tissulaire toutes deux homogènes dans le r 'jaco-organe. L'hypothèse H1 précitée est donc vérifiée. - Estimation des paramètres du modèle: Le modèle M,1 est une droite dans l'espace des cinétiques. Il peut être estimé au moyen d'une analyse en composantes principales opérée sur les cinétiques des éléments de volume de l'échantillon et centrée sur la moyenne du signal au sein de l'échantillon. L'hypothèse H2 précitée est donc vérifiée. - Réponses différentes des pharmaco-organes au traceur: Un pharmacoorgane est défini par des cinétiques tissulaire et artérielle homogènes en son sein. En revanche, rien n'empêche deux pharmaco-organes voisins i et o de montrer une réponse identique au traceur, ou bien telle que Cto appartienne à la droite CaCt;. Rien n'assure donc que l'hypothèse H3 précitée soit vérifiée. Considérons cependant comme un seul et unique pharmacoorgane deux pharmaco-organes voisins présentant des réponses identiques. L'hypothèse H3 est alors vérifiée. La méthode sera toutefois incapable de séparer les deux parties de ce pharmaco-organe. - Bruit de mesure et de reconstruction: En TEP, pour des méthodes de reconstruction itératives ou analytiques, le bruit peut être considéré comme additif gaussien, 25 conformément à l'hypothèse H4 précitée. - Mélange local de cinétiques: Comme nous venons de le voir, l'effet de volume partiel et les mouvements physiologiques génèrent au sein d'un élément de volume donné un mélange local de cinétiques de pharmacoorganes voisins. Dans le cas de l'effet de volume partiel, nous venons de voir que le mélange de cinétique, opéré au sein d'un élément de volume j donné, est un barycentre des modèles de cinétiques physiologiques aux éléments de volume voisin. De même, les mouvements physiologiques périodiques de faible période opèrent au sein du voxel j un barycentre des modèles de cinétiques des pharmaco- organes qui ont traversé j au cours du mouvement. Plaçons-nous dans le cas où tous les voisins de j appartiennent au pharmaco-organe i. M; étant une droite D; dans l'espace des cinétiques, tout barycentre de cinétiques appartenant à la droite D, appartient à la droite D,. La contrainte Cl précitée est donc vérifiée. b) Etape 20 de pré-segmentation: Comme indiqué ci-dessus dans l'exposé général de la méthode selon l'invention, les zones présentant une absence de signal nuisent à la rapidité de l'algorithme et à la qualité des résultats. Ces zones sans signal correspondent en TEP à l'extérieur de l'organisme: le fond. L'image TEP peut donc être découpée en deux zones, l'une représentant l'organisme et l'autre le fond. Les cinétiques du fond sont caractérisées par un faible signal, engendrant un grand nombre de valeurs nulles ou négatives, qui, étant non physiologiques, ne peuvent représenter une cinétique organique. Ces dernières, caractérisées par une activité en générale importante par rapport au bruit, s'annulent rarement et ont une moyenne plus forte que les cinétiques du fond. L'analyse d'histogramme est une méthode rapide permettant de déterminer des seuils pour séparer plusieurs populations d'éléments de volume. L'histogramme représente le nombre d'éléments de volume ayant une activité donnée, ou plus généralement dont l'activité appartenant à un intervalle donné. Les intervalles sont jointifs et souvent pris de même largeur. 3G > Méthodes usuelles: En TEP, un histogramme de l'image de transmission permettra de segmenter aisément l'organisme, mais dissociera difficilement ce dernier des objets présents dans le champ de vue: lit, dispositif pour rongeurs, éventuel respirateur, etc. L'image usuellement utilisée pour réaliser un histogramme à partir de l'émission est la séquence d'images elle-même, moyennée sur le temps ou bien dont une image particulière a été extraite. Sur un tel histogramme, le premier pic représente le fond ou une partie de celui-ci. Du fait de la grande disparité des valeurs au sein de l'image, une analyse par histogramme des valeurs de l'image, si elle permet de détecter un pic correspondant au fond, n'assure pas la détermination d'un seuil séparant ce fond du reste de l'image (voir figures 6a, 6b et 7a, 7b) . Histogramme du nombre de valeurs négatives des cinétiques En revanche, l'histogramme du nombre de valeurs négatives ou nulles de la cinétique présente deux pics appelés modes distincts, l'un pour l'organisme, qui présente une radioactivité non négligeable par rapport au bruit, et l'autre pour l'air entourant l'organisme, qui présente une radioactivité négligeable par rapport au bruit. La cinétique d'un élément de volume inclus dans l'organisme présentera un faible nombre de valeurs nulles ou négatives, tandis que la cinétique d'un élément de volume du fond présentera de nombreuses valeurs négatives, comme le montrent les graphiques des figures 4 et 5 (respectivement des exemples de cinétiques d'un élément de volume inclus dans l'organisme et d'un élément de volume extérieur à l'organisme). Soit ni le nombre de valeurs négatives prises par une cinétique C; donnée. On s'attend à ce que ni soit grand dans le fond et petit dans l'organisme. L'histogramme des valeurs prises par n; dans l'image présente deux modes, l'un correspondant au fond (Mode,o d)et l'autre à l'organisme (Modeorg). Celui-ci présentant généralement peu d'involutions (corps entier, crâne), sa surface représente un nombre d'éléments de volume négligeable par rapport à son volume et au volume du fond. Cet histogramme présente donc au moins un minimum d'indice nmin situé entre Modeorg et Modefo d. Une analyse multi-échelle de l'histogramme [Mangin, 1998] assure l'extraction du minimum d'indice nmin de l'histogramme entre les deux modes Modeorg et Modefo d. Tous les éléments de volume j de l'image qui sont tels que n; A l'intérieur de l'organisme, certains organes, comme les poumons, peuvent présenter une très faible activité, et être exclus de la zone d'intérêt. Afin de pallier cet inconvénient, on ajoute à la zone d'intérêt toute zone de fond non connexe avec les bords de l'image. c) Etape 30 d'extraction des voxels-modèles: Dans la suite nous prendrons comme voisinage t' d'un élément de volume j la boule Bj,r de rayon r. De plus, nous ne détaillerons l'exemple de réalisation que dans le cas d'une image reconstruite par une méthode analytique, pour laquelle le bruit au sein de l'image peut être supposé stationnaire dans l'espace. La variation de cet exemple de réalisation dans le cas d'une image reconstruite par une méthode itérative, pour laquelle la variance du bruit est liée au signal, est donné dans l'annexe C jointe. Selon le modèle M, JO) choisi précédemment pour rendre compte de la cinétique physiologique au sein du pharmaco-organe i et pour tous les éléments de volume j inclus dans i, M;,j(t) appartient à une droite dans l'espace des cinétiques. Plaçons-nous dans le cas d'une reconstruction analytique, pour laquelle le bruit est stationnaire dans l'espace. Supposons l'hypothèse H5 vérifiée et que le bruit au sein d'une image donnée de la séquence, non corrigée de la décroissance radioactive, dépende essentiellement de la durée d'acquisition de cette image. Supposons cette dépendance linéaire 6(t) = / At, où Z est une quantité stationnaire dans le temps et l'espace et Ar est la durée d'acquisition de l'image t de la séquence. Soit c,.(r) = c, (r) x A, la cinétique à l'élément de volume k, dont la valeur au temps t a été pondérée par la durée d'acquisition de l'image t de 5 la séquence. L'écart-type du bruit de l'ensemble des cinétiques c,'(1) ainsi pondérées, est égal à Z. L'étude des cinétiques ch (u) permet de s'affranchir de la non stationnarité temporelle du bruit. - Indicateur de l'intérieur des organes: L'hypothèse HOv1 précitée est équivalente, dans le cas présent, à une hypothèse selon laquelle le signal au sein de la boule B1,r centrée au voxel j et de rayon r peut être modélisé par une droite dans l'espace des cinétiques, plus du bruit additif. Considérons non plus les cinétiques, mais les cinétiques 15 c,'(1) = c, (r) x A,, le bruit pouvant être de plus supposé gaussien. Soit une analyse en composante principale centrée (ACP) réalisée sur les cinétiques c,(t) au sein de la boule B11,. Soit {Ad, , Discussion: La mesure I'vj de validité de l'hypothèse HOv, est une approximation de la moyenne des proportions de signal non reconstruit au sein de B;,,. Cette approximation présente l'avantage d'être peu bruitée, d'être un bon indicateur de l'inclusion de B;,, dans un pharmaco-organe, mais aussi de détecter les points dont le voisinage est, en l'absence de points complètement inclus dans un organe donné, au moins majoritairement inclus dans ce dernier. Ainsi, des pharmaco-organes ne présentant aucune cinétique non polluée par les signaux physiologiques des organes voisins peuvent néanmoins être détectés. Les figures 8a et 8b illustrent des cartes des IR'i dans le cas d'une image reconstruite par une méthode itérative (la figure 8a étant une coupe coronale passant par le coeur et les reins du rat, et la figure 8b une coupe coronale passant par la vessie du rat). - Extraction des voxels-modèles Le rayon r de la boule B;,, est choisi de façon à ce que le nombre de cinétiques qu'elle contient suffise au calcul de l'ACP. Il sera pris par exemple égal à n x 3. Nous calculons donc la carte des Iv;, dont les minima locaux correspondent à des points situés au coeur des organes. Le voisinage de ces points est peu affecté par l'effet de volume partiel. Ces minima locaux, qui peuvent être multiples au sein du même organe, sont extraits automatiquement. Afin d'éliminer les maxima locaux dus au bruit, l'image des est filtrée par un filtre gaussien dont la variance, choisie par un algorithme itératif, est juste suffisante pour que le nombre de minima locaux extrait soit inférieur à un seuil fixé a priori. Ce seuil sera fixé relativement haut pour éviter qu'un trop fort lissage ne déplace ou ne fasse disparaître des minima locaux correspondant à l'inclusion de la boule dans l'organe. Les figures 9a et 9b illustrent les minima locaux de la carte e:._ v;, par extraction des voxels-modèles (la figure 9a étant une coupe coronale passant par le coeur et les reins du rat, et la figure 9b une coupe coronale passant par la vessie du rat). d) Etape 40 d'optimisation des échantillons: 5 La géométrie choisie en forme de boule Bir peut ne pas être appropriée à certains pharmacoorganes aux formes allongées, plates ou en forme de pelure d'oignon . Nous allons donc définir un échantillon de points connexes à j minimisant un critère de non-reconstruction par le modèle à un organe calculé sur les cinétiques des éléments de volume inclus dans Bj,r (deux points connexes sont tels qu'ils appartiennent à la même région et qu'il existe un chemin inclus dans cette région que les relie) . La cinétique corrigée c; (t) = C,(t) x o, à l'élément de volume k s'écrit: c (t)=(s,(t) +%z_(t)) x0. Nous utiliserons ici la méthode par croissance de région proposée dans le principe général afin d'optimiser la forme des échantillons yr,. Cette croissance de région, dont l'échantillon initial est formé des cinétiques des éléments de volume inclus dans Bj,r, utilisera comme mesure de similarité entre une cinétique C'(t)et l'échantillon yr, l'inverse de la proportion de signal de C'^ (t) non reconstruit par le vecteur propre associé à la plus grande valeur propre de l'ACP centrée calculée sur yr,. Soit p'i,k cette proportion de signal non reconstruite: P j.k(t)= ( x, A ,_ (1)x max(A,.) / -- est de l'ordre de (t) C: (t) C (t) 1 au sein de l'organe, et de l'ordre de 1 en dehors de l'organe, où C (t)L E,k(t) contient encore du signal non reconstruit. Elle est aussi grande en bordure de l'organisme, car y décroît fortement alors que fjk(t) y reste c (r) de l'ordre de Elle permet donc non seulement d'écarter les cinétiques trop éloignées du modèle, mais encore d'éviter, dans le cas d'un organe i présentant une forte activité polluée par les activités d'organes voisins, de préférer une cinétique faible à une cinétique de l'organe i. A chaque étape de raffinement de l'échantillon, une ACP est donc calculée sur y, puis une croissance de région par propagation de front est opérée à partir de l'élément de volume j, seul constituant de la région au départ. A chaque étape, l'élément de volume k du front qui minimise pj, k(t) est agrégé à l'échantillon i. L'introduction d'un modèle d'Ising permet de régulariser les contours des échantillons obtenus. La taille finale de chaque échantillon est prise égale au nombre d'éléments de volume inclus dans la boule Les figures 10a et 10b illustrent le résultat de la croissance de région sur les échantillons (la figure 10a étant une coupe coronale passant par le coeur et les reins du rat, et la figure 10b une coupe coronale passant par la vessie du rat). e) Etape 50a de classement des cinétiques: Le bruit des cinétiques Ck est considéré comme stationnaire dans l'espace, et le bruit des cinétiques CA, est considéré comme stationnaire dans l'espace et le temps. La probabilité a posteriori de la donnée CA connaissant le modèle NI, R s'écrit: r2 5 p(C, / M; ,_) = y x e " , où y est un scalaire indépendant de l'organe i comme de l'élément de volume k considéré. Une fois les modèles M; calculés à partir des cinétiques C,' in de chaque échantillon par une analyse en composantes principales dont on ne conservera que le vecteur propre associé à la plus grande valeur propre, la maximisation de p(C; M A) revient à trouver l'échantillon i qui minimise CA ti1, A On attribue donc à l'élément de volume k le label o tel que: o= arg min C, M;ti Régularisation du classement des cinétiques: Afin de limiter les chances d'une agrégation erronée de k du fait du bruit, on définit le processus de segmentation de l'élément de volume k non plus comme: o= arg min C M;_A mais comme o=argmax où Vois(k) désigne Axe l'ensemble des éléments de volumes voisins de l'élément de volume k, et A.k un facteur de pondération dépendant de la distance entre les éléments de volume j et k. Les figures 11 a et 11 b illustrent le classement des cinétiques en 100 classes (la figure 11 a étant une coupe coronale passant par le coeur et les reins du rat, et la figure 11 b une coupe coronale passant par la vessie du rat). f) Etape 50b de fusion des échantillons: L'étape d'extraction des voxelsmodèles peut extraire plusieurs voxels-modèles au sein d'un organe donné. Aussi l'image de segmentation résultante est peu exploitable en tant que telle, car un élément de volume appartenant à un pharmaco-organe est assigné de façon aléatoire à l'un des échantillons associés à ce même pharmaco-organe. Nous disposons de I échantillons, dont plusieurs peuvent se trouver au sein du même organe. Une classification ascendante hiérarchique 5 (HL) [Jain, 1988] permet de fusionner les échantillons. Deux échantillons ont de bonne chance de représenter le même organe si les cinétiques de l'un sont bien reconstruites par le modèle à un seul organe estimé sur l'autre, et si les cinétiques moyennes de chacun d'entre eux présentent des pics d'activité similaires (en position et en intensité). Dans le cas d'une reconstruction analytique, la variance du bruit associé aux cinétiques cF (t) = cyt) x n, est égal à 2, stationnaire dans l'espace et le temps. Soit ï; = t la moyenne des T-1 plus petites Tt,_2 valeurs propres d'une ACP calculée sur les cinétiques (-Jr) au sein de 15 l'échantillon r. 42 peut être estimée comme la moyenne des A2 sur les échantillons: La fonction de coût de fusion de deux échantillons a et b peut s'écrire ainsi: Cam /jxmax N 1 min +tSx ninl >harmaco-organes aux cinétiques proches, mais dont les modèles de tes présentent un léger déphasage. Les poids fi, 8 et y sont fixés de façon heuristique à fi = 0,5, 6 = 0,25 et y= 0,5. Entre deux groupes d'échantillons Ga et Gb pouvant contenir chacun plusieurs échantillons, on définit la fonction de coût suivante CG,,(b = max (C h), dont le calcul est extrêmement rapide, une fois Cam aega_b eb calculé pour tout couple (a, b). É Initialisation Initialement, chaque échantillon forme un groupe au sein duquel il se trouve seul. Le coût de fusion de chaque couple d'échantillon est calculé. É Etape de fusion A chaque étape de la classification ascendante hiérarchique, 15 les deux groupes d'échantillons Ga et Gb pour lequel le coût de fusion CGa,cb est minimum sont fusionnés en un nouveau groupe Gc. É Fin de l'algorithme Le processus cesse lorsque tous les échantillons ont été fusionnés. La séquence des fusions entre échantillons est conservée sous la 20 forme d'un arbre illustré à la figure 12. Le niveau E de fusion correspond à une image de classement des cinétiques entre E échantillons. Afin d'obtenir une image de segmentation associée, il est nécessaire de dissocier les zones rattachées au même échantillon, mais non connexes entre elles. Cette opération est réalisée par une extraction des composantes connexes pour chaque label. On ne conservera que les composantes connexes présentant un volume supérieur au plus petit volume attendu pour un pharmaco-organe. L'affectation d'un label différent pour chacune de ces composantes connexes permet d'obtenir nn partition de l'image en régions, c'est-à-dire une segmentation. Les figures 13a et 13b résultent d'une fusion des échantillons par classification ascendante hiérarchique, et les figures 13c et 13d résultent d'une fusion des labels de l'image de classement des cinétiques par classification ascendante hiérarchique, 18 classes ayant été conservées pour chacune de ces figures (les figures 13a et 13c étant des coupes coronales passant par le coeur et les reins du rat, et les figures 13b et 13d des coupes coronales passant par la vessie du rat). Les figures 14a, 14b et 14c sont des vues 3D de segmentation après fusion, respectivement coronale, sagittale et oblique. 10 g) Etape 60 de représentation hiérarchique des orqanes: L'arbre de fusion des régions permet de représenter l'image segmentée non pas sous forme d'une image de label, mais d'un arbre d'images de label, dont chaque noeud ou feuille correspond à un groupe de pharmaco-organes à une étape de la fusion. Il est donc possible de suivre, moyennant un logiciel de visualisation, la fusion des diverses zones de label d'un même organe. Nous laissons le choix à l'opérateur du niveau de fusion qu'il estime optimal pour la reconnaissance d'un pharmaco-organe donné. En effet, ce niveau de fusion optimal peut être différent pour deux pharmacoorganes différents. La figure 15 illustre un exemple de représentation hiérarchique des organes des rats examinés. h) Etape 70 de séparation des composantes connexes: Les composantes connexes de chaque agrégat obtenu par les étapes précédentes sont séparées par une méthode de fusion de labels. A chaque élément de volume est attribué un numéro. Les numéros de deux éléments de volumes voisins appartenant au même agrégat sont fusionnés en seul numéro, jusqu'à ce que cette opération ne soit plus possible. On obtient alors un numéro identifiant chaque composante connexe de chaque classe, ce qui constitue l'image de segmentation finale. Applications de la méthode selon l'invention 1) Applications réalisées: a) Fantôme de rat: Fantôme numérique MCAT : Afin de valider la méthode, un fantôme numérique de rat a été réalisé à partir du fantôme MCAT [Segars, 2004] et simulé sur une caméra Siemens de type HR+ . Les organes de ce fantôme, segmentés sur une image de scanner X, sont modélisés par des surfaces NURBS . Outre un lissage des contours des organes, cette modélisation permet un calcul rapide de déformation. Le fantôme de rat prend donc en compte les mouvements dus à la respiration et aux battements du coeur. Ces mouvements périodiques à haute fréquence génèrent un mélange de cinétiques des organes en mouvements. Un élément de volume de l'image, fixe dans le repère de la caméra simulée, contient donc un mélange des cinétiques des organes qui l'ont traversé lors de ces deux mouvements physiologiques. Le remplissage de la vessie lors de l'examen, mouvement physiologique handicapant la segmentation bien plus que les mouvements périodiques, comme nous l'avons vu précédemment, n'est pas pris en compte par le fantôme MCAT . Nous l'avons simulé par une dilatation progressive de la zone de l'image correspondant à la vessie le long de l'examen simulé. 3u EL!, -c:i comprend 37 images successives acquises sur une durée d'une chacune. Détermination des cinétiques physiologiques: Un calcul préalable, pour chaque trame, de la probabilité d'appartenance de chaque élément de volume aux différents organes permet de recomposer rapidement une image TEP non bruitée, à partir de l'équation suivante [Kamasak, 2003] : Ci (i) =1(p x[Vb,j xCa(i)+(] Vb,)xCi;(t)), où p;" est la probabilité d'appartenance de l'élément de volume j à l'organe i, Ca(t) est la cinétique vasculaire, identique pour tous les éléments de volume j et tous les pharmaco-organes i, Ct;(t) est la cinétique tissulaire de l'organe i, identique pour tous les éléments de volume j et Vb;j est un scalaire représentant le ratio de volume sanguin dans l'élément de volume j de l'organe i. Une fois la cinétique artérielle Ca(t) déterminée, Ct;(t) est calculée comme une convolution de la cinétique artérielle par une somme 15 pondérée d'exponentielles décroissantes [Kamasak, 2003] . ( CL(t)=Ca(i) * 1[a,,,.xe Les paramètres a;,W et b;,W sont déterminés aléatoirement pour chaquefonction exponentielle w et chaque pharmaco-organe i. Projection et reconstruction: Cette image 4D est ensuite projetée analytiquement, bruitée et reconstruite par un algorithme de reconstruction soit analytique, soit itératif. Les différents phénomènes de détection ou de reconstruction tels que les coïncidences fortuites, la diffusion Compton et le spill-over sont aussi simulés. Segmentation et critère de qualité de segmentation: La segmentation de l'image et la comparaison des résultats avec les probabilités d'appartenance aux divers organes permet d'établir, une fois l'identification i(1) faite entre le label 1j et l'organe i correspondant, un critère global de qualité de segmentation: 1 x Ce critère 1=1 sera égal à 1 pour une segmentation parfaite et à 0 pour une segmentation catastrophique. On peut aussi définir un critère de détection des organes, correspondant au nombre d'organes pour lesquels une correspondance directe d'organe à organe est impossible, du fait d'une fusion ou d'une scission de deux organes dans l'image de segmentation, divisé par le nombre d'organes défini dans le fantôme. b) Projet PNA : Les ( Par conséquent, ce type d'oligonucléotide semblait idéal pour développer des méthodes permettant de moduler leurs paramètres pharmacocinétiques de façon prévisible. La fonctionnalisation (modification chimique) du squelette peptidique des PNA a permis de synthétiser 12 dérivés possédant toujours leur activité anti-sens. C'est dérivés ont été marqué au fluor 18 et leurs pharmacocinétiques ont été analysées par jerie TEP chez le rat. Pour chaque PNA , quatre animaux ont été par injection simultanément. En appliquant la méthode de entation automatique d'image TEP à cette série de 12 PNA plus un ( PNA de contrôle, il a été possible d'isoler de façon automatique les organes et ainsi de d'obtenir les paramètres pharmacocinétiques par organe pour chaque type de PNA . Les cinétiques obtenues automatiquement ont montré une bonne concordance aux cinétiques obtenues par une segmentation manuelle des organes discernables sur l'image TEP. Cette étude montre que les pharmacocinétiques et pharmacodistributions d'oligonucléotides peuvent être fortement modifiées. Ces résultats présentent une avancée forte dans le développement de ( 2) Applications envisagées de la méthode de l'invention: La méthode de segmentation automatique d'images selon l'invention utilisant l'information pharmacocinétique peut être appliquée avantageusement à toute étude basée sur une mesure localisée dans l'espace, au sein d'un organisme, et présentant une variation au cours du temps. Les modalités concernées sont notamment la médecine nucléaire (techniques TEP et SPECT , notamment), mais aussi l'IRM fonctionnelle et l'imagerie optique. On peut notamment citer les applications suivantes. a) Pharmaco imagerie: Imagerie TEP corps entier: La présente méthode de segmentation permet d'obtenir une bonne segmentation, malgré les mouvements physiologiques du sujet imagé. - Pharmaco-imagerie sans a priori: La présente méthode de segmentation ne nécessite aucune connaissance préalable sur la réponse des divers organes au traceur. Tri in vivo de molécules: Devant les capacités de la chimie combinatoire à générer de grandes banques de molécules, se pose le problème de la sélection de molécules d'intérêt. Cette sélection peut être opérée in vitro ou en fonction de caractéristiques chimiques, de polarité ou de caractéristiques géométriques de la molécule, mais rien n'assure alors le passage des barrières systèmes et l'efficacité in vivo de la molécule. La bio-distribution in vivo renseigne sur les capacités de ciblage de la molécule injectée au sujet imagé. Il devient donc possible de sélectionner une molécule parmi un groupe de molécules d'intérêt. Cependant, le traitement manuel des données d'un grand nombre d'examens est une tâche longue et fastidieuse, ce qui limite dans les faits la quantité de molécules testées. La présente méthode de segmentation permet, par son caractère automatique et sa rapidité, de traiter rapidement un grand nombre d'examens, permettant un tri in vivo dont le débit est désormais limité par l'acquisition des examens et non le traitement des données issues de ceuxci. > Etude de nouveaux traceurs: La présente méthode ne requérant pas de connaissance a priori du comportement des traceurs le modèle appliqué reste valable pour tout traceur , il permet l'obtention rapide de renseignements sur la 25 pharmacocinétique d'un nouveau traceur, aux propriétés inconnues. Ciblage d'organes: Certaines molécules présentent un intérêt in vitro du fait de leur principe actif, mais n'entrent pas dans un organe cible. Elles ne présentent donc aucun intérêt in vivo. Cependant, par ajout d'étiquettes chimiques ne modifiant pas le principe actif d'une telle drogue , ou en ha I'encapsulant dans un vecteur, il devient possible de modifier ses capacités de ciblage. Une telle molécule ainsi déclinée génère une famille de molécules au principe actif supposé identique, mais dont les pharmacocinétiques seront différentes. Les molécules ainsi générées peuvent subir un tri in vivo, comme décrit ci-dessus, pour sélectionner celles qui présentent les capacités de ciblage requises. b) imagerie en physiologie; c) Détection d'analoqie de réponse des orqanes à un traceur: La présente méthode de segmentation ne permet pas de distinguer deux pharmaco-organes aux fonctionnements identiques. Dans l'image résultant du procédé de segmentation, ces deux organes se verront attribuer le même label. Il est ainsi possible de repérer des analogies entre les fonctionnements de deux organes. De façon analogue, certains organes sont présents dans l'organisme par paire (yeux, reins, poumons). Si le procédé de segmentation permet de les dissocier, il est fort probable que leur réponse au traceur soit différente, par exemple du fait du dysfonctionnement de l'un d'entre eux. d) Détermination de modèles de cinétiques: L'étape 40 de définition d'échantillons voire même directement l'étape 30 qui extrait les voxelsmodèles permet d'estimer les modèles de cinétiques propres à chaque organe. Dans l'exemple de réalisation donné, les modèles des tous les organes, des droites dans l'espace des cinétiques, ont un unique point d'intersection: la cinétique artérielle. Celle-ci peut donc théoriquement être déterminée par un calcul de l'intersection des droites-modèles dans l'espace des cinétiques. L'estimation de ces modèles étant biaisée par l'effet de volume partiel et le nombre d'échantillons relativement faible, les droites n'auront pas probablement pas de point d'intersection, mais passeront près de la cinétique artérielle. L'estimation des modèles peut donc permettre une estimation directe de la cinétique artérielle ou bien servir d'initialisation à une méthode permettant de le faire. e) Seqmentation préalable à une reconstruction des cinétiques des structures: 10 Parmi les utilisations de la segmentation, on compte désormais la reconstruction des cinétiques dans l'organisme, non pas élément de volume par élément de volume, mais organe par organe. En effet, le rapport signal sur bruit dépend de la finesse de l'échantillonnage temporel et de la résolution spatiale. En dégradant la résolution spatiale, il est possible d'affiner l'échantillonnage temporel. Plutôt que de dégrader la résolution spatiale de manière arbitraire, il paraît naturel de regrouper les éléments de volume ayant une réponse identique aux traceurs. Plutôt que de segmenter des images anatomiques et de les recaler sur la modalité TEP afin de déterminer ces regroupements, le procédé de segmentation proposé permet de se passer de la modalité anatomique en utilisant comme regroupement les pharmaco-organes issus de la segmentation. En effet, le recalage entre modalité anatomique et fonctionnelle est rendu difficile dans le cas du corps entier, et la modalité anatomique n'est même pas toujours disponible. f) Utilisation des modèle de cinétique en reconstruction d'image: La méthode de segmentation selon l'invention permet d'accéder à un modèle de cinétique de chaque structure présente dans l'image. L'accès à ce modèle ne nécessite pas de segmenter entièrement l'image. Ce modèle de cinétique peut avoir de nombreuses applications, comme par exemple la reconstruction 4D en TEP [Nichols, 2002]. Plus généralement, l'information temporelle pourra être utilisée comme a priori, ou en information supplémentaire pour améliorer la qualité des images, notamment en terme de fluctuation statistique. g) Reconnaissance des structures en TEP: La méthode de segmentation selon l'invention permet de regrouper sous le même label les éléments de volume appartenant au même pharmaco-organe. Il reste alors à reconnaître les organes correspondant aux différents labels. La représentation hiérarchique proposée ici devrait faciliter cette étape de reconnaissance, en repérant le niveau de fusion permettant de faire apparaître un organe précis isolé des autres et en un seul morceau. h) Vérification de la cohérence fonctionnelle d'un organe anatomiquement segmenté : La corrélation entre les contours anatomiques d'un organe et l'identité des réponses de ses sous-parties à un traceur donné n'est pas assurée lorsqu'une tumeur est présente par exemple. Le procédé proposé permet de séparer automatiquement un organe en sous-parties aux réponses distinctes au traceur. i) Séquences temporelles en imagerie optique: L'imagerie optique permet d'obtenir des images anatomiques, mais aussi fonctionnelles. Ces images présentent une très haute résolution spatiale et temporelle. La méthode selon l'invention permet de segmenter automatiquement de telles images. j) Images statiques: La présente méthode peut être appliquée à des images purement 2D ou 3D , moyennant un modèle adéquat choisi lors de l'étape 10. Références bibliographiques citées - Ashburner et coll., A cluster analysis approach for the characterisation of dynamic PET data. Quantification of Brain Function using PET, pages 301-306, San Diego, Academic Press (1996). - Acton, P.D., Pilowsky, L.S., Kung, H.F., Ell, P.J. Automatic segmentation of dynamic neuroreceptor single-photon emission tomography images using fu77y clustering. European Journal of Nuclear Medecine and Molecular Imaging, Springer-Verlag, Heingelberg, 26(6):581-590 (1999). - Brankov, J.G., Galatsanos, N.P., Yongyi Yang Wernick, M.N., Segmentation of dynamic PET or fMRI images based on a similarity metric. IEEE Transactions on Nuclear Science, 50(5):1410-1414 (2003). - Guo, H., Renaut, R., Chen, K., Reiman, E. Clustering huge data sets for parametric PET imaging,. BioSystem 71, 81-92 (2003). - Frouin, F., Boubacar, P., Frouin, V., De Cesare, A., Todd-Pokropek, A., Merlet, P., Herment, A. 3D regularisation and segmentation of factor volumes to process PET H215O myocardial perfusion studies. Functional Imaging and Modelling of the Heart (FIMH'2001). Lecture Notes in Computer Sciences, Springer Verlag, Heidelberg, 2230: 91-96 (2001). - Frouin, V., Comtat, C., Reilhac, A., Grégoire, M.-C., Correction of Partial-Volume Effect for PET Striatal Imaging: Fast Implementation and Study of Robustness. Journal of Nuclear Medicine 43- 12:1715-1726 (2002). - Jain, A.K., Dubes, R.C., Algorithm for Clustering Data. Advanced Reference. Prentice-Hall, Englewood Cliffts, NJ (1988). - Kimura, Y., Senda, M., Alpert, N. Fast formation of statistically reliable FDG parametric images based on clustering and principal 3 0 components. Phys. Med. Biol. 47(3):455-468 (2002). - Mangin, J.-F., Coulon, O., Frouin V., Robust Brain Segmentation Using Histogram Scale-Space Analysis and Mathematical Mor r'hology. MICCAI, Springer Verlag, Heidelberg, LNCS 1496, pp. 1230 124 i (1998). - Minka, T.P., Automatic choice of dimensionality for PCA. M.I.T. Media Laboratory Perceptual Computing Section Technical Report No 514 (2000). - Nichols, T.E., Qi, J., Asma, E., Leahy, R.M., Spatiotemporal Reconstruction of List-Mode PET Data. IEEE Transactions on Medical Imaging, 21(4):396-404 (2002). - Segars, W.P., Tsui, B.M.W., Frey, E.G., Johnson, G.A., Berr, 10 S.S., Development of a 4D Digital Mouse Phantom for Molecular Imaging Research. Molecular Imaging and Biology (2004). - Tipping, M.E., Bishop, C.M., Mixture of probabilistic principal component analysers. Neural Computation, 11(2):443-482 (1999). - Wong, K.-P., Feng, D., Meikle, S.R., Fulham, M.J. 2001. 15 IEEE Transactions on Nuclear Science 49:200-207 (2002) - Zhou, Y., Huang, S.-C. , Bergsneider, M., Wong, D.F. A non linear regression with spatial constraint for generation of parametric images in dynamic PET studies. Journal of Nuclear Medecine 42 (5), 100. ANNEXE A: Allure des cinétiques dans un modèle compartimentai à quatre paramètres. > Modèle de cinétique à un seul traceur: Le modèle M;,1(t) du signal TEP au sein d'un élément de volume j donné compris dans l'organe i s'écrit, après correction de la décroissance radioactive: M (t) = ((l Vb, )x Ct;(t) + Vb, x Ca(t))x D où D est la dose injectée et Vb,,1 est le ratio de volume sanguin au sein de l'élément de volume j de l'organe i. Lorsque j se déplace dans l'organe i, M,.j(t) décrit un segment 15 de la droite Ca(t)Ct,(t) dans l'espace des cinétiques. Modèle de cinétique à plusieurs traceurs: Lorsque plusieurs traceurs { T r s} , ket'oisinnee( j) - Au sein d'un organe: Lorsque le voisinage de j sur lequel s'étend l'effet de volume 15 partiel ne contient que l'organe i, le signal S1(t) s'écrit: ) xl(D xCtiS(t))+ 1L/3i.k x VU k kal'oisinngi{ J) s x x Cas (1)) s=1 S,(î)= I(A.k kdoinnnoel j) S1(t) appartient toujours à la droite passant par 1[D x Ct s (t)) s et 1 (Ds x Cas (1)). - En bordure d'organe: Lorsque, au contraire, ce voisinage de j contient I (I > 1) organes dont les modèles sont différents, C1 se calcule comme une combinaison linéaire des signaux appartenant chacun à une des droites M;,1. L'ensemble des signaux contenus dans le voisinage de j ne peut donc pas être décrit par une droite. En effet: x 1(D, x CI (t))+ s=1 x lJlJ A AEI ns5nner I i l'+ )) 1(P,_A Si (0= bd I 1orsin ngai i la \ D xCn,(t)) x s=I \ / A)) X1(DS XCt(t)) + L[Nr. A xJ h,_A] _ s I Irl Ael'oignngft j P A S, (t) _ IO, A) Considérons une zone de l'image affectée par un mouvement périodique de période très inférieure à la durée d'acquisition. Le signal est accumulé pendant une durée au cours de laquelle un grand nombre de périodes sont achevées. Le temps d'acquisition At d'une image peut être découpé en: 4t = n x P + cSt, ou P est la période du mouvement physiologique. Si P At, on a aussi ()t At. Si l'on suppose continue la cinétique des traceurs, le signal enregistré pendant la durée vt est négligeable par rapport au signal enregistré pendant la durée n x P. Ce sont donc, à peu de chose près, les mêmes sites physiologiques qui sont observés au sein d'un volume donné pour chaque image de la série étudiée. La démonstration du caractère linéaire du modèle à un organe malgré l'effet de volume partiel est 2886433 73 aisément transposable au problème du mouvement physiologique périodique rapide. Soit k une unité de volume de centre x se déplaçant avec les structures le long d'une courbe P et j une unité de volume fixe dans le champ 5 de vue. - Mouvement au sein du même organe i: Le signal Sj(t) à l'élément de volume j s'écrit: S1(t)= f (%j;_A, x(1 Vb,A,,,)) 'x 2 s x 10, x C't;, (r)) + s x1 (D, xCas(t)) s Sj(t) appartient toujours à la droite passant par 1(DS x Cr (t)) s=1 s et 1 (Ds xCas(t))_ - Mouvement passant par plusieurs organes: ( s f (/3 0 x VIz 1Y x 1(D X Ct,_s (0) x1(DS xCa(t)) , f (Pj A1 X Vb A1,1)(./Y rz-, + Sj(t) ne décrit donc pas une droite, mais un sous espace de dimension d telle que 1 Mouvement physiologique et effet de volume partiel: On montre de manière connue que le modèle à un organe i n'est conservé que si le voisinage de k sur lequel s'étend le volume partiel, durant tout le trajet SZ de k de centre x dû au mouvement physiologique périodique, est inclus dans i. Toujours avec l'hypothèse qu'il n'existe pas deux organes différents qui sont tels que leurs modèles de cinétiques sont identiques, le signal S1(t) à l'élément de volume j fixe dans le champ de vue appartient à la droite Ca(t)Ct,(t) pour tout j, si et seulement si le voisinage de j sur lequel s'étend l'effet de volume partiel reste inclus dans l'organe i durant le mouvement physiologique périodique. - Conclusion: 5 Si l'on exclut les zones de mouvement physiologique non périodique ou à longue période, et si l'on suppose qu'aucun couple d'organes ne possède le même modèle de cinétique, le modèle à un organe i est donc valide au sein d'un volume B, fixe dans le référentiel de la caméra, si et seulement si le voisinage Vois(B) de ce volume sur lequel s'étend l'effet de volume partiel reste inclus dans l'organe i lors de la période du mouvement physiologique. ANNEXE B: Choix de la mesure de non-validité locale du modèle à un organe. 1) Vraisemblance du modèle de signal choisi: Si Bj,r est incluse dans un pharmaco-organe i, Ho_ q(t) peut s'écrire C,(t) = ,u, (t) + y, / x e, (t) + g) (t) La vraisemblance de Ho est un bon indicateur de la validité du modèle sur Bi.r. Lorsque la boule Bi,, est incluse dans un organe, la dimension du signal est égale à 1. La sélection de modèle dans le cadre Bayésien, si CL, est l'ensemble des cinétiques contenues dans la boule BI, et K est le nombre d'organes présent dans Bj,,, revient à maximiser: p(C, / k) = J p(Ci, t ) p(/ K)dO pour K = 1. Dans le cas d'une modélisation par ACP (analyse en composantes principales) avec un bruit additif gaussien, l'approximation de Laplace de l'évidence p(C, /k) pour le modèle à K vecteurs propres se révèle [Minka, 2000] être un excellent indicateur de la dimension réelle du 20 modèle. L'approximation de Laplace de la vraisemblance des données selon le modèle à K dimensions s'écrit: p(0) i(0)(2, )ro srA) 2 A avec 0=argrnax(/(9)), A= Ud' logf (0) d.dO. _0_net 0 _ (U,A,v). U est la matrice des vecteurs propres de la matrice de corrélation, A la matrice des valeurs propres et v une estimation de la variance du bruit et du signal non reconstruit par l'ACP au sein de la boule B;,,. Soit U; L' et v' les valeurs des paramètres U, A et v maximisant f. d N A /=K +I n(d k) 2 Les valeurs diagonales de A s'écrivent À'1= NA. + a pour 1 N 1 +a K et,i,=V sinon, avec a représentant l'importance de l'a priori. p(C) /K) p(U) ( K -N 2 2i y-.F(d-k) 2(2rl(m+k) 2 i2N-k2 Maxima locaux de la vraisemblance du modèle à un organe: Soit j tel que la vraisemblance des données selon le modèle 10 de l'organe i soit maximale sur Bi.r. j vérifie p(CI /K = l) > p(Ck /K =1) pour tout k appartenant au voisinage de j. Après un calcul de la carte spatiale des p(C /K =1), les maxima locaux sont extraits. Chaque maximum local situé dans un organe correspond à une vraisemblance maximale de la validité du modèle. Ces points sont donc très peu entachés de volume partiel, auquel cas la dimension du modèle nécessaire à la caractérisation du signal serait supérieure à 1. Discussion: Cette mesure s'avère en pratique extrêmement 20 bruitée et peu sélective de l'intérieur des organes. Carte du nombre d'organes présents dans le voisinage: D; = argmax(p(C /K)) est la dimension la plus vraisemblable données au sein de Bi,r. Les zones Di est égale à 1 sont des zones 6=gis à l'intérieur des organes. Discussion: Cet indice présente des points de bruit impossibles à dissocier des réels éléments de volume j tels que B;r soit incluse dans un organe. De plus, les organes de petite taille ou de forme allongée ne présentant pas de points non entachés par le signal des pharmaco-organes voisins du fait de l'effet de volume partiel ils ne sont pas détectés. 2) Analyses en composantes principales locales: Soit {e;,,(t)},_, =T les vecteurs propres d'une analyse en composante principale (ACP), centrée en 4), opérée sur les cinétiques des éléments de volume inclus dans la boule Bj,r. Soit o-2ji la variance au temps t du bruit s;(t) et {} ,,,T les valeurs propres de l'ACP. La variance moyenne de Ej(t) s'écrit: 6-2 = I x T Plaçons-nous pour commencer dans le cas d'une reconstruction par une méthode analytique. Le bruit étant stationnaire tant dans l'espace que dans le temps, sa variance est constante: 6;,,2 = a2. Au sein r-1 avec /' _ de l'organe, 0- 2 = T 1 i=z J. 112 semble donc un bon indicateur de la position de l'élément de volume j par rapport aux frontières des organes. En effet, lorsque l'élément de volume j contient, du fait de l'effet de volume partiel, les cinétiques de plusieurs pharmaco-organes, les valeurs propres de l'ACP non prises en compte pour modéliser le signal contiennent à la fois de l'information de pharmaco-cinétique et du bruit. On aura alors u2 Toutefois cette mesure ne permet pas de discriminer les points intérieurs à l'organisme des points extérieurs. En effet, lorsque j est situé en dehors de l'organisme, l'extérieur de l'organisme ne contenant pas de traceur, j ne contient que du bruit, d'où = 1. De même, entre deux organes i et o, i présentant un fort signal et o un faible signal de l'ordre de la variance du bruit, %i2/o-2 sera peu différent de 1. - Moyenne de la proportion de signal non reconstruit: Définissons la mesure suivante d'inhomogénéité locale des cinétiques relative au signal: 1 PNR = --- X N jA AEB Soit j tel que Bj,, c organe i, et m tel que #(Bm, rn Bi,r) #( Bm,r - Bm,rn Bi,r) et que Bm.r Bm.rn Bir sont constitué de points extérieurs à l'organe, dont les cinétiques ne peuvent être modélisées par les K premiers vecteurs de l'ACP calculée suri. Supposons 6, 1) Pour k tel que k E Bm,n i, 2, avec 2 m > ai, comme nous l'avons vu plus haut. ", alors que fnLA C 2) Pour k tel que k E Bm,, Bm,,n i, 1, car Ck est mal Fm.A C'A reconstruite par l'ACP caractérisant l'organe i. 3) Pour k tel que k B,,,,,ni, Pk appartient à l'organe i: 21. PNR,,, - PNR i = 1 x N {k ke/3,,) ik k B,,,,) {k,ke8f,r)B ) j x N 11E/ Ick B 2 2 A - N, x ( ' _ - 6 /dl > O. PNR - PNR. _ > 0. Pour tout élément de volume k en bordure de l'organe i, et pour tout j tel que Bj,r est incluse dans l'organe i, PNRk > PNR. Avec les hypothèses que ai et que les organes voisins de l'organe i sont /ki caractérisés par des cinétiques présentant une forte proportion de signal non reconstruite par le modèle défini sur i, il existe donc un minimum PNR de 15 PNR au sein de l'organe i, tel que Bj.r soit incluse dans l'organe i. Cette mesure d'inhomogénéité nécessite le calcul d'une ACP pour chaque Bj, r, mais aussi celui de la projection de la cinétique de chaque élément de volume inclus dans Bj,r sur le premier vecteur propre de cette ACP. Or, le critère étant calculé sur chaque élément de volume inclus dans 2 G l'organisme imagé, ce surcoût de calcul se révèle prohibitif. - Variance relative des résidus: Moyennant l'approximation qu'au sein et en bordure d'un organe, 2, définissons la mesure suivante d'inhomogénéité locale Ck des cinétiques relative au signal: IR= où fy est la moyenne du signal sur la boule B incluse dans l'organe, et /12, = est la variance du signal non reconstruit par une ACP opérée sur les cinétiques des éléments de volume 10 appartenant à la boule Bj,1 centrée en j et de rayon r, dont seules les K premiers vecteurs propres ont été conservés. Variance relative non reconstruite corrigée pour les reconstructions itératives: Dans le cas d'une reconstruction itérative avec un nombre fixe d'itérations pour chaque image de la séquence reconstruite, la variance du bruit 62(t) à l'élément de volume j et à l'instant t peut s'écrire: a 2(t)= a2 x S;(t) / 4,. o 2(t) est non stationnaire, tant d'un point de vue spatial que temporel, car proportionnel au signal. Afin de s'affranchir de la contribution du signal à la variance du bruit, interdisant toute comparaison de variance de bruit inter voxels, nous étudierons non pas 6 j(t), mais a2 = Gj(t)2 / ,S (t) x 4, . En effet, a2 est une quantité indépendante du signal de la durée d'acquisition des images de la séquence, donc stationnaire dans l'espace et le temps. Si l'on suppose une variation continue et faible du signal dans l'espace, S;,k(t), pour k E B,,,, peut être approché par la cinétique moyenne,cil, au sein de la boule Bi,,. De même, S, (r) = S; (/)x A, peut être approché par u;(1) , la moyenne des cinétiques c1 (r) au sein de la boule B,,. si Soit l'analyse en composantes principales K(t), {e,,,(t)}) calculée sur les cinétiques C,,( ) = C (r)/liu; (r) = CA (t) x o, /J/I (,) . C;A(t)=S. A(t)xA,Nf',(r)+N(0,a2)=S"(t)+N(0.a'2), avec aa, stationnaire dansl'espace. Quand B1,1 est incluse dans un organe i, 5 S A. (r) = M 'A (t) et décrit une droite dans l'espace des cinétiques lorsque j se déplace au sein de i. Le critère de variance relative non reconstruite devient Ti, = , qui mesure la validité du modèle à un organe au sein de la boule B1,1. ANNEXE C: Variations de l'exemple de réalisation dans le cas d'images reconstruites par un algorithme itératif. Dans le cas d'un algorithme itératif, la variance du bruit est corrélée au signal. Afin de pouvoir comparer deux indices basés sur une mesure de l'erreur de reconstruction par le modèle de cinétique de pharmacoorgane, il semble donc essentiel de corriger de cette influence du signal à chaque fois que c'est possible. Dans le cas d'une reconstruction itérative avec un nombre fixe d'itérations pour chaque image de la séquence reconstruite, la variance du bruit (et) à l'élément de volume j et à l'instant t peut s'écrire: a2(t)= a2 x Sj(t) / 4,. aj2(t) est non stationnaire, tant d'un point de vue spatial que temporel, car proportionnel au signal. Afin de s'affranchir de la contribution du signal à la variance du bruit, interdisant toute comparaison de variance de bruit inter voxels, nous étudierons non pas 6 j(t), mais a2 = 61(t)2 / S,(t) x 4t. En effet, a2 est une quantité indépendante du signal de la durée d'acquisition des images de la séquence, donc stationnaire dans l'espace et le temps. Indicateur de l'intérieur des organes: Si l'on suppose une variation continue et faible du signal dans 2.5 l'espace, S;,k(t), pour k E Bj,r, peut être approché par la cinétique moyenne jy,, au sein de la boule Bj,r. De même, s, A (t) = S; r (t) x A, peut être approché par u, (i) , la moyenne des cinétiques ch (t) au sein de la boule BI,. Soit l'analyse en composantes principales ('j(t), {e',)(t)}) r.ée sur les cinétiques ci A, (t) = C;_ (t) (t) = Ch (t) x o, /Vf,, (t) . C,, _(t) =S1A(t)x4,Np;(r)N(O.a)=5(t)+N(O.a'), avec a a, stationnaire dans l'espace. Quand B1,,- est incluse dans un organe i, Si 1. (t) _/1;h (t) et décrit une droite dans l'espace des cinétiques lorsque j se déplace au sein de i. Le critère de variance relative non reconstruite devient = , qui mesure la validité du modèle à un organe au sein de la "; boule B1,1. - Optimisation des échantillons: La proportion de signal non reconstruite calculée sur B11, doit être corrigée de l'influence du signal sur le bruit. Elle devient p''j.k Ei-ti(t)2 s,h(t)) = où S; (t) est le signal estimé en k comme la reconstruction de Ch (t) par le modèle à un organe. s, h (t) s'écrit comme S, = fp. (t) (e, . C, )x e (t) . p'j(t) est la moyenne des cinétiques Ch (t) au sein de la boule B1,r et e,., le vecteur propre associé à la première valeur propre d'une ACP calculée sur les cinétiques Ch (t) au sein de B1,1 et centrée en,u1(t). 2 0 > Classement des cinétiques: Soit le modèle M, définit sur les cinétiques C; (I) = C, (I)/-f,u, (I) = C; (t) x A. / V,1, (t) au sein de l'échantillon L M,.i décrit toujours une droite dans l'espace des cinétiques lorsque j varie au sein de i. L'erreur de reconstruction moyenne au sein de l'échantillon i de la cinétique 5 C;(I) par le modèle suit une loi stationnaire pour./ appartenant à i. Une fois les modèles M, calculés à partir des cinétiques C,", au sein de chaque échantillon, la segmentation de l'image revient à trouver, pour tout élément de volume k, l'échantillon i qui minimise C; -11L1, On attribue donc à l'élément de volume k le label o tel que: o=argmin -M,"., Fusion des échantillons: La variance du bruit associé aux 'cinétiques C, , ( I) = C; (I) x o /V/1:(I) , pour l'élément de volume j inclus dans l'échantillon i, et égal à a'2, stationnaire dans l'espace comme dans le temps. Soit _ 1 T la moyenne des T-1 plus petites valeurs propres d'une ACP calculée sur les cinétiques c, h (I) au sein de l'échantillon i. a peut être estimée comme la moyenne des 2,2 sur les échantillons: a''=1/? ,_, La fonction de coût de fusion de deux échantillons a et b devient:'max, 'max, t max, t maxr, Ca,b fi xmax +8x(1 minulf1oll_/II/tbll>_>1I bII,/If REVENDICATIONS 1. Méthode de segmentation d'une image ou séquence d'images tridimensionnelles de départ à base de voxels pour l'obtention d'une image tridimensionnelle de segmentation comportant une partition en régions d'intérêt, ladite image ou séquence d'images comprenant des mesures, pour chaque voxel et au cours de n intervalles de temps (n 1), de l'évolution réelle d'un signal représentatif d'au moins une variable de type physique, chimique ou biologique de ladite image ou séquence, caractérisée en ce qu'elle comprend les étapes suivantes: a) une modélisation (10) du signal comprenant la définition d'un modèle paramétrique d'évolution spatio-temporelle dudit signal, ce modèle comprenant des jeux de paramètres homogènes de sorte que lesdits jeux soient respectivement propres à des structures correspondant auxdites régions d'intérêts et que chaque jeu de paramètres soit indépendant des coordonnées spatiales dans la structure correspondante; b) une extraction (30) d'échantillons de voxels de sorte que lesdits échantillons soient respectivement inclus dans lesdites structures, cette extraction comprenant: (i) un calcul, pour chaque voxel de ladite image ou séquence d'images de départ ou bien d'une zone d'intérêt de celle-ci, d'un critère de validité d'une hypothèse selon laquelle ledit modèle d'évolution de ladite variable au sein du voisinage de ce voxel est propre à une et une seule desdites structures, (ii) une extraction de voxels-modèles qui sont définis comme réalisant les maxima locaux dudit critère, (iii) une définition, pour chaque voxel-modèle, de l'un desdits échantillons de voxels inclus dans ladite structure correspondante, de manière à ce que cet échantillon présente une évolution de ladite variable qui 30 soit conforme à celle du modèle de la structure à laquelle appartient ledit voxel-modèle, puis (iv) une estimation, sur chaque échantillon ainsi défini, des paramètres dudit modèle d'évolution de la structure dans laquelle ledit échantillon est inclus; puis c) une fusion desdits échantillons regroupant les échantillons dont le modèle d'évolution est propre à la même structure, ladite fusion suivant, précédant ou incluant une classification de tous les voxels de ladite image ou séquence d'images ou d'une zone d'intérêt de celle-ci par agrégation à un groupe d'échantillons, de manière qu'une similarité maximale existe entre l'évolution de ladite variable pour ces voxels et l'évolution dudit modèle caractérisant ce groupe d'échantillons. 2. Méthode de segmentation selon la revendication 1, caractérisée en ce que l'on admet qu'il existe dans chaque image de départ de ladite séquence une partition de l'espace en un nombre fini desdites structures, lesquelles sont chacune connexes et présentent chacune en leur sein un comportement homogène en réponse à un phénomène étudié dont l'évolution de ladite variable est représentative. 3. Méthode de segmentation selon la revendication 2, caractérisée en ce que l'on détermine le nombre desdites structures a posteriori. 4. Méthode de segmentation selon une des revendications précédentes, caractérisée en ce que ledit modèle comprend en outre des paramètres hétérogènes dépendant des coordonnées spatiales des voxels au sein d'une même structure, et en ce que l'on admet pour l'étape a) que lesdits paramètres homogènes et hétérogènes peuvent être estimés sur un ou plusieurs éléments de volume inclus dans cette même structure. 5. Méthode de segmentation selon la revendication 4, caractérisée en ce que l'on admet en outre pour l'étape a) que lesdites structures présentent entre elles des réponses différentes à un phénomène étudié dont l'évolution de ladite variable est représentative, à moins qu'elles ne soient pas connexes. 6. Méthode de segmentation selon la revendication 5, caractérisée en ce que l'on admet en outre pour l'étape a) que le bruit, qui est inhérent auxdites mesures et qui est dû au mode d'acquisition de ladite séquence d'images de départ et de segmentation, est additif. 7. Méthode de segmentation selon la revendication 6 caractérisée en ce que l'on se fixe en outre pour l'étape a) les deux contraintes suivantes, pour prendre en compte des mélanges locaux de différentes évolutions temporelles dudit signal: (i) si la totalité des voxels qui sont voisins d'un voxel et qui contribuent à l'évolution de ladite variable relative à ce voxel est incluse dans la même structure, alors ledit signal correspondant est une réalisation du modèle de cette structure, et (ii) ledit jeu de paramètres homogènes pour ces voxels voisins est le même que celui dudit modèle propre à cette structure. 8. Méthode de segmentation selon la revendication 2, caractérisée en ce que, dans le cas où ladite image ou séquence d'images de départ est obtenue par une technique d'imagerie par injection de traceur, ledit modèle utilisé à l'étape a) est un modèle compartimenta) de type à un ou plusieurs traceurs indépendants et à plusieurs compartiments, tels que des compartiments biologiques. 9. Méthode de segmentation selon une des revendications précédentes, caractérisée en ce que l'étape b) d'extraction des voxelsmodèles comprend une sélection préalable, dans un espace à n dimensions (n 1) correspondant à ladite image ou séquence d'images de départ, de voxels- modèles pour que chacune desdites structures d'intérêt contienne en son coeur au moins un voxel-modèle ainsi que le voisinage de celui-ci. 10. Méthode de segmentation selon la revendication 9, caractérisée en ce que l'étape b) comprend, suite à ladite sélection préalable, la définition d'une métrique prévue pour définir les distances dans ledit espace et ledit critère d'extraction local. 11. Méthode de segmentation selon la revendication 9, caractérisée en ce que l'on admet à l'étape b) que pour toute structure à segmenter, il existe un voxel dont le voisinage est compris au sein de la région d'intérêt correspondante. 12. Méthode de segmentation selon une des revendications précédentes, caractérisée en ce que l'étape b) comprend: - une détermination, pour chaque structure d'intérêt, d'un échantillon de voxels appartenant à cette structure, cette appartenance étant déterminée par une similitude de l'évolution de ladite variable avec le modèle relatif à cette structure, puis - une estimation, pour chaque échantillon, des paramètres homogènes du modèle qui correspondent à ladite structure. 13. Méthode de segmentation selon une des revendications précédentes, caractérisée en ce qu'il comprend en outre une étape de pré-segmentation (20) qui est mise en oeuvre avant ou après l'étape a) et qui consiste à séparer lesdites images de départ en une première partie contenant lesdites structures d'intérêt et en une seconde partie correspondant à un fond d'image exclu de la segmentation. 14. Méthode de segmentation selon une des revendications précédentes, caractérisée en ce que l'étape c) de fusion (50b) et de classification (50a) comprend une fusion desdits échantillons correspondant à une même structure d'intérêt, avant ou après ladite classification des voxels. 15. Méthode de segmentation selon la revendication 9, caractérisée en ce que la fusion desdits échantillons correspondant à une même structure d'intérêt est opérée par une classification ascendante hiérarchique. 16. Méthode de segmentation selon une des revendications précédentes, caractérisée en ce qu'elle comprend en outre une étape d'optimisation (40) des échantillons qui est mise en oeuvre suite à l'étape b) et avant l'étape c) et qui consiste à rechercher, pour chaque échantillon, une forme géométrique, paramétrique ou libre qui minimise en son sein le signal d'autres structures que celle dont le voxel-modèle a été extrait. 17. Méthode de segmentation selon une des revendications précédentes, caractérisée en ce qu'elle comprend en outre, suite à l'étape c), une étape de séparation (70) des structures non connexes qui présentent des évolutions similaires de ladite variable. 18. Méthode de segmentation selon la revendication 17, caractérisée en ce qu'elle comprend en outre une étape (60) de représentation hiérarchique desdites régions d'intérêt obtenues qui est mise en ceuvre suite à l'étape c) et avant ladite étape de séparation des structures non connexes. 19. Méthode de segmentation selon une des revendications précédentes, caractérisée en ce qu'elle exclue la bordure de chaque structure de l'extraction desdits échantillons mise en oeuvre à l'étape b) . 20. Méthode de segmentation selon une des revendications précédentes, caractérisée en ce que ladite image ou séquence d'images est obtenue par une technique d'imagerie choisie dans le groupe constitué par la tomographie par émission de positons, l'imagerie par résonance magnétique, la tomographie d'émission de photons, les imageries optiques, le scanner aux rayons X, l'imagerie histologique, l'imagerie autoradiographique, l'imagerie satellitaire et l'imagerie photographique. 21. Méthode de segmentation selon une des revendications précédentes, caractérisée en ce que ladite image ou séquence d'images est obtenue par la technique de tomographie par émission de positons. 22. Méthode de segmentation selon la revendication 21, caractérisée en ce que l'on met en oeuvre à l'étape b) des analyses en composantes principales locales dans un espace à n dimensions pour déterminer, en fonction des valeurs propres et des vecteurs propres obtenus, les voxels qui présentent une évolution temporelle similaire de ladite variable. 23. Méthode de segmentation selon une des revendications précédentes, caractérisée en ce que lesdites images de départ sont des images d'un organisme entier, ladite méthode de segmentation segmentant ledit organisme selon une partition en pharmaco-organes. 24. Méthode de segmentation selon la revendication 23, caractérisée en ce que ledit corps est animé de mouvements physiologiques soit de type périodique, dont la période est réduite en comparaison de la durée d'acquisition de l'une des images de départ de ladite séquence, soit de type non périodique. 25. Méthode de segmentation selon la revendication 23 ou 24, caractérisée en ce que ladite variable représente la concentration radioactive en un instant donné to à tn d'au moins un principe actif marqué et injecté dans ledit organisme, l'ensemble des voxels à l'intérieur de chaque pharmacoorgane présentant des cinétiques pharmacologiques de répartition dudit principe actif qui sont similaires.