Brevet FR2700909A1 — publié le 19940729

DISPOSITIF ET PROCEDE AUTOMATIQUE DE CALIBRATION GEOMETRIQUE D'UN SYSTEME D'IMAGERIE PAR RAYONS X La présente invention concerne un dispositif et un procédé automatique de calibration géométrique d'un système d'imagerie à rayons X. Des systèmes de ce type sont utilisés principalement dans le domaine médical dans le but d'analyser des structures internes de corps humains notamment de structures vasculaires. Le but de l'invention est de permettre, par une calibration aisée du système, son utilisation pour évaluer quantitativement, et non plus uniquement qualitativement, les dimensions des structures internes analysées. L'invention se rapporte essentiellement aux systèmes dont le détecteur de rayonnement est bidimensionnel, par exemple, sensiblement plan.Elle se rapporte, pour le procédé, aux systèmes tournants d'imagerie par rayons X. Dans le cas d'un système d'imagerie à rayons X composé d'une source de rayons X et d'un détecteur bidimensionnel, tournant ou non, l'opération géométrique intervenant dans la production de l'image est une projection conique d'un objet analysé, déployé dans un espace à trois dimensions (3D), sur un espace à deux dimensions (2D) qui est celui du plan de projection correspondant au plan de détection. Les paramètres géométriques décrivant complètement la projection conique sont au nombre de neuf. Ces neuf paramètres peuvent être présentés selon différentes paramétrisations. Une des paramétrisations possibles est la suivante: - trois premiers paramètres concernent les trois coordonnées de la source S de rayons X dans le repère R de l'espace 3D dans lequel est placé l'objet analysé; - trois autres paramètres concernent les trois angles d'Euler du changement de repère, entre le repère R de l'espace 3D de l'objet et un repère R' lié à la projection conique, c'est à dire formé d'un repère 2D du plan de projection P et dd'un troisième axe orthogonal au plan de projection; - trois derniers paramètres concernent les trois coordonnées d'un point arbitraire I du plan de projection P (par exemple le centre de l'image) dans le repère R. Pour simplifier les calculs dans l'invention on verra par la suite qu'il y a lieu de préférence de retenir des coordonnées cylindriques pour les trois premières et trois dernières coordonnées. Le repère R est alors un repère cylindrique. La connaissance de tout ou partie de ces paramètres est très souvent utile en radiologie, notamment lorsqu'une donnée quantitative sur un objet 3D est estimée à partir d'une mesure faite dans la projection 2D, ou dans plusieurs projections 2D, de cet objet. Des paramètres de chaque point de vue doivent alors être connus. Un point de vue concerne l'orientation du système par rapport à l'objet. A titre d'exemple, on retiendra le problème classique de l'estimation de la section d'un vaisseau à t partir d'une ou plusieurs projections 2D. Or il est souvent impossible, ou trop imprécis, d'accéder à ces paramètres de façon directe, c'est à dire, par exemple, en mesurant directement sur le système d'acquisition la distance entre la source de rayons X et le détecteur. On appelle calibration géométrique d'un système d'imagerie, l'opération qui aboutit à la connaissance indirecte précise des paramètres géométriques qui interviennent dans la production d'une image. Le principe, classique, est basé sur l'utilisation d'un fantôme géométrique connu dans l'espace 3D, et dont on acquiert la projection 2D. La suite des opérations effectuées dans ce but comporte les étapes suivantes: - on dispose d'un objet connu, un fantôme de calibration présentant un certain nombre de points caractéristiques dont la position dans l'espace est connue par des coordonnées mesurées par rapport à un repère propre à cet objet; - on acquiert l'image de ce fantôme dans les conditions géométriques d'un point de vue (ou incidence) que l'on veut calibrer; - on reconnaît les projections des points caractéristiques dans l'image.Pour cela on associe chaque point caractéristique de l'objet à sa trace dans l'image acquise, la projection; - on inverse au sens mathématique le système d'équations décrivant la projection; - et on obtient finalement l'ensemble des paramètres de la projection pour le point de vue donné. Dans l'état de l'art, la reconnaissance des points caractéristiques dans la projection du fantôme de calibration, est faite "manuellement" par un opérateur humain. Par exemple, une forme de fantôme de calibration géométrique souvent utilisée est celle d'un cube, aux huit coins duquel sont disposées des billes métalliques opaques aux rayons X. Parfois des billes supplémentaires sont rajoutées pour augmenter la précision géométrique de#la calibration. Or, la projection au sens des rayons X produit une image "en transparence" de l'objet, dans laquelle, selon l'orientation du cube, il peut être très difficile d'associer sans se tromper les traces 2D des billes avec leurs homologues 3D. On peut se référer, pour l'état de l'art de la calibration géométrique en radiologie aux articles suivants: (1) D.L. Parker, J.Wu, D.L Pope, R. Van Bree, G.R. Caputo and H.W. Marshall, "Three-dimensional reconstruction and flow measurements of coronary arteries using multi view digital angiography", in New Developments in Quantitative Coronary Arteriography, J.C. Reiber and P.W. Serruys Eds., pp. 225-247, Kluwer Academic Publishers, 1988; (2) D.J. Hawkes, A.C.F. Colchester and C.R. Mol; "The accurate 3-D reconstruction of the geometric configuration of the vascular trees from X-ray recordings," in Physics and Engineering of Medical Imaging, R. Guzzardi Ed., Nijhoff, 1987. (3) M. Garreau, JTL Coatrieux, R. Collorec and C. Chardenon, "A knowledge-based approach for 3-D reconstruction and labeling of vascular networks from biplane angiographic projections", IEEE Medical Imaging, vol. 10, N02, pp. 122-131, June 1991. Cette nécessité d'intervention humaine est un inconvénient majeur, et peut même devenir prohibitive dans certains cas, lorsque le nombre de points de vue à calibrer est important C'est le cas en particulier pour tout système d'acquisition de projections bidimensionnelles par rayons X dans le but d'une reconstruction 3D. Un tel système de reconstruction est par exemple décrit dans la demande de brevet français NO 2 644 590 déposée le 20 mars 1989. Quelle que soit la méthode de reconstruction 3D, il est nécessaire au préalable de connaître parfaitement les paramètres géométriques qui caractérisent chacune des projections. Dans le cas probable où une méthode de calibration indirecte est nécessaire, il semble indispensable d'effectuer cette calibration de façon automatique. La présente invention propose donc un dispositif permettant une calibration automatique, qui s'applique à tout système d'imagerie bidimensionnelle à rayons X, et en particulier à tout système utilisé à des fins de reconstruction 3D. Le dispositif de l'invention est tel que la reconnaissance, sans équivoque, des correspondances entre les billes et leurs traces dans les projections est automatique. Ceci est obtenu en choisissant un fantôme dans lequel les billes sont distribuées, de proche en proche, en une succession telle que des altitudes de billes, mesurées le long de l'axe de rotation du système d'imagerie, et surtout un axe du fantôme, soient monotones croissantes (ou décroissantes) avec un numéro dtordre des billes dans la succession. Dans ce but l'invention a donc pour objet un système d'imagerie par rayons X avec un détecteur bidimensionnel situé en vis à vis d'un tube à rayons X, ce détecteur et ce tube pouvant dans certains cas tourner ensemble autour d'un axe, ce système d'imagerie par rayons X comportant un dispositif de calibration géométrique muni d'un fantôme de dimensions connues et de moyens de mesure pour déduire, à partir d'images en projection du fantôme sur le détecteur bidimensionnel, des coefficients de calibration de ce système d'imagerie par rayons X, caractérisé en ce que le fantôme comporte une succession de proche en proche de structures cellulaires à absorption radiologique différente, plus élevée ou plus faible, par rapport à leur environnement, chaque structure cellulaire étant identifiable automatiquement par une caractéristique ordonnée d'identification de cette structure cellulaire. Un autre problème à résoudre pour l'invention est lié au nombre des calibrations à entreprendre. Un système d'imagerie tournant de ce type est en effet susceptible d'occuper 1000 positions angulaires, ou plus, sur un tour. Des calibrations pour toutes les positions angulaires, même menées automatiquement, sont trop longues à effectuer pour pouvoir être répétées périodiquement, toutes les semaines ou tous les mois. Il faut en fait prévoir une journée de travail. Le procédé de l'invention vise à résoudre cet autre problème en faisant l'hypothèse que, si les dimensions théoriques de construction du système d'imagerie ne peuvent pas être considérées comme respectées dans le système réel, on peut néanmoins admettre que le système d'imagerie est à rotation circulaire parfaite ou presque parfaite. L'idée du procédé est alors de rechercher les paramètres intrinsèques du système: en définitive ceux qui donnent en plus la direction de son axe de rotation. A partir de ces paramètres intrinsèques, qu'on peut connaître très précisément, on peut ensuite calculer pour chaque position angulaire des coefficients de calibration. On montre que, par exemple, une centaine de calibrations conduisent à une connaissance très précise de ces paramètres intrinsèques. Avec cette connaissance, on calcule 900 autres "calibrations". On en connaît alors 1000 au total en n'en ayant mesurées que 100. Le gain de temps est considérable. L'invention a donc également pour objet un procédé de calibration géométrique d'un système d'imagerie par rayons X à détecteur bidimensionnel tournant autour d'un axe, en vis à vis d'un tube à rayon X et comportant les étapes suivantes: - on place un fantôme de dimensions connues entre le tube et le détecteur, - on mesure, pour une position en rotation de l'ensemble tube détecteur par rapport au fantôme, dans l'image radiologique du fantôme projetée sur le détecteur bidimensionnel, des coordonnées de lieux images de points caractéristiques du fantôme, - on en déduit, pour cette position en rotation de l'ensemble tube détecteur, et dans un repère associé à la position du fantôme, des coefficients de calibration relatifs aux positions respectives d'un foyer de rayonnement du tube à rayons X et du détecteur bidimensionnel, - et on réitère ces deux dernières étapes pour des positions désirées différentes de l'ensemble tube détecteur, caractérisé en ce que - on traite les coefficients de calibration relatifs à un certain nombre de positions en rotation de l'ensemble tube détecteur pour en extraire des paramètres intrinsèques du système d'imagerie par rayons X qui sont indépendants de la position en rotation de l'ensemble tube détecteur, - et on déduit de ces paramètres intrinsèques des coefficients de calibration relatifs à une position quelconque en rotation de l'ensemble tube détecteur en évaluant les valeurs de fonctions analytiques de ces paramètres intrinsèques pour une valeur angulaire de position en rotation de l'ensemble tube détecteur. L'invention sera mieux comprise à la lecture de la description qui suit et à l'examen des figures qui l'accompagnent. Celles-ci ne sont données qu'à titre indicatif et nullement limitatif de l'invention. Les figures montrent: - Figure 1: un système d'imagerie utilisable pour mettre en oeuvre l'invention avec la présentation des différents repères d'acquisition; - Figures 2a à 2d: respectivement des vues de profil et de face du fantôme selon l'invention, la trace des billes dans les projections, et un détail de mise en place des billes par incrustation; - Figure 3: la mise en place du fantôme dans le système d'imagerie; - Figures 4a et 4b: des images en projection d'une forme particulière du fantôme après traitement; - Figures Sa à 5d: des représentations géométriques permettant une meilleure compréhension du procédé; La figure 1 montre un système d'imagerie utilisable pour mettre en oeuvre l'invention. Ce système comporte un tube 1 à rayons X avec un foyer S de rayonnement. Un détecteur 2 de rayonnement, par exemple du type amplificateur de luminance utilisé en radiologie, est placé en vis à vis du tube 1. Une normale au plan P de détection du détecteur 2 passe sensiblement par le foyer S et le centre I de ce détecteur. Ceci peut être une nécessité pour la formation d'images classiques, mais n'est pas directement lié au problème de la calibration. En particulier, on peut calibrer, avec l'invention, un tomographe du type à translation homothétique.L'ensemble peut tourner avec un mécanisme (non représenté) de type connu pour les tomodensitomètres, ou avec un arceau, autour d'un axe 3 montré ici verticalement mais qui, dans la pratique, est normalement horizontal. Dans les cas des tomographes, il n'y a que des translations et un centre d'homothétie, sans même nécessairement un axe. Il en est de même bien sûr des appareils statiques. Le système comporte également un processeur de traitement 4 qui fait subir aux signaux électriques mesurés par le détecteur 2 un traitement contenu dans un programme 5 stocké dans une mémoire programme 6. Des images, ou des résultats statistiques résultant de ces traitements, peuvent être visualisés sur un moniteur 7. L'invention permet de résoudre le problème de la reconnaissance automatique des points homologues 2D/3D, en utilisant comme fantôme de calibration un fantôme particulièrement adapté montré sur les figures 2a à 2d et 3. Ce fantôme est dans un exemple, un cylindre 8 de plexiglass, qui sert de support à une série de petites billes métalliques incrustées dans la paroi du cylindre, et disposées régulièrement selon une hélice. Les billes ont toutes un diamètre identique, par exemple 1 mm, sauf la bille du milieu, qui a un diamètre plus élevé que les autres, par exemple 3 mm. Les paramètres de l'hélice sont grossièrement adaptés au système d'acquisition des images, de façon à ce que sa projection soit la plus grande possible sur le détecteur tout en restant contenue dans l'image. De même, le pas de l'hélice et le pas d'échantillonnage des billes sur l'hélice sont choisis de telle sorte qu'un nombre suffisant de billes apparaisse dans la projection, sans qu'elles puissent se superposer. Les figures la à ld sont cotées en millimètres.Pour éviter toute déformation du fantôme 8, notamment par vieillissement du plexiglass le cylindre 8 est obturé des deux côtés par un couvercle 9 qui en assure la rigidité. Le fantôme 8 est adapté par exemple à un système d'acquisition constitué d'une source de rayons X et d'un détecteur de champ rond de diamètre 37 cm, situé à une distance d'environ 135 cm, alors que l'objet est à environ 100 cm du foyer S. Ces chiffres sont des exemples qui ne sont en rien restrictifs, de même que les dimensions du fantôme indiquées en annexe. Dans l'exemple il y a 61 billes. Une des coordonnées des billes, l'altitude mesurée sur un axe longitudinal Z, évolue de proche en proche de 5 mm en 5 mm avec une précision de 5/100ème. Il y a 20 billes réparties régulièrement par tour de l'hélice. Les trois particularités principales du fantôme 8 sont les suivantes. D'abord le fantôme 8 dispose d'un nombre suffisant de billes pour la calibration. Elles sont bien réparties dans l'espace 3D. Leurs projections sont également bien réparties dans le plan P (condition sine qua non pour une calibration 3D). La régularité de la distribution n'est cependant pas indispensable. Dans ce but les billes peuvent se numéroter facilement en fonction de leur cote selon l'axe Z. Ceci est obtenu ici selon une formule très simple: z = h.n+zO, z = coordonnées selon Z, z0 = coordonnée de la bille numéro 0, h = pas d'échantillonnage des billes selon Z, n = numéro de la bille. C'est en fait la relation d'ordre qui nous importe ici: si zl > z2, alors nl > n2 zl, z2 = coordonnées des billes 1 et 2 selon Z nl, n2 = numéro des billes. Si le fantôme est positionné de façon judicieuse devant le détecteur de rayons X, c'est à dire essentiellement de telle façon que l'axe Z soit à peu près parallèle au plan P du détecteur 2, et si la conicité du faisceau n'est pas trop importante, cette relation d'ordre se retrouve également dans les projections 2-D. En effet, figure 3, soit V l'axe dans le plan P de projection qui correspond à la projection de l'axe Z. On aura la relation suivante: si vl > v2, alors nl > n2 vl, v2 = coordonnées des projections des billes 1 et2 selon V nl, n2 = numéros des billes Dans la suite, et par souci de simplicité, on supposera que l'axe V coïncide sensiblement avec l'axe vertical dans les images de projection, c'est à dire est parallèle aux colonnes de l'image. Si ce n'est pas le cas, on peut introduire une rotation dans le plan de projection pour s'y ramener. Toujours par souci de simplicité pour la suite, on définit le repère dans l'espace 3D comme étant lié à l'hélice: c'est à dire que l'axe Z de l'hélice définit un des axes du repère R, les deux autres axes étant perpendiculaires à Z. Dans ce repère, les coordonnées des billes s'expriment très facilement. Ceci signifie que, une fois la calibration faite, le système possède un axe vertical qui est celui qu'occupait l'axe du fantôme au cours de la calibration. Des reconstructions 3D relatives à des calibrations différentes ne sont donc pas directement comparables. Ceci sera par contre possible avec le procédé de l'invention vu plus loin. On a représenté dans la figure 3 le fantôme 8 positionné entre la source S de rayons X, et le détecteur D de rayons X.Le fantôme est positionné de telle façon que son axe principal Z est à peu près parallèle au plan P, et tel qu'il se projette selon l'axe V du détecteur. Enfin les billes peuvent se numéroter de façon absolue. Même si le haut et le bas de la projection de l'hélice sont tronquées (c'est à dire que l'on ne voit pas la première bille ni la dernière bille du fantôme), on a une référence absolue qui est la grosse bille particulière du milieu. La bille particulière, ou deux billes particulières, ou autres, sont donc de préférence placées dans le fantôme à un endroit qui a de grandes chances d'être dans le champ: au milieu du fantôme. Pour qu'on comprenne le principe mis en oeuvre dans le procédé de calibration, on a décrit dans les paragraphes précédent un cas particulier de fantôme qui répond aux critères exposés. En fait, tout fantôme ayant les particularités suivantes permet une calibration automatique. Ces particularités sont: - Le fantôme doit disposer de structures cellulaires, pour simplifier des billes, visibles en rayons X, et bien réparties dans l'espace 3D. Leur contraste radiologique par rapport à leur environnement doit être élevé. Les projections des billes doivent être contenues dans l'image 2D au moins pour les billes centrales. Il n'est pas du tout nécessaires que le support des billes soit un cylindre creux de plexiglass. - Il doit exister une relation d'ordre entre un numéro affectable aux billes et leurs coordonnées, aussi bien dans l'espace 3D selon l'axe Z, que dans la projection selon un axe, par exemple l'axe V, et quelle que soit l'orientation du fantôme dans l'appareil. L'axe du fantôme doit être sensiblement parallèle au plan du détecteur. La solution des billes réparties régulièrement sur une hélice comme dans le cas particulier du fantôme décrit, n'est pas du tout la seule. - On doit pouvoir reconnaître au moins une des billes parmi les autres, de façon à permettre une numérotation absolue même si on ne voit pas la première et la dernière bille dans la projection. Cette bille de référence peut se distinguer par autre chose que par son diamètre, par exemple sa forme, son opacité aux rayons X, son absence, ou encore une irrégularité du pas de distribution, ou autre. A titre d'exemple, on citera une autre réalisation possible pour le fantôme: les billes sont réparties sur une hélice "dégénérée", qui est portée par une sphère de plexiglass au lieu d'être portée par un cylindre. Dans toute la suite, on raisonne avec le fantôme préféré "hélice" tel que décrit ci-dessus. Toute la méthode exposée ci-après s'applique également avec un autre fantôme répondant aux critères. On indique ci-après la suite des opérations à effectuer pour calibrer le système. 1. Le fantôme 8 décrit ci-dessus est positionné de telle façon que son axe principal (axe Z) soit à peu près parallèle au plan du détecteur. C'est le cas par exemple lorsqu'on le positionne couché sur une table d'examen d'un système angiographique, et que l'arceau qui porte le tube à rayons X et le détecteur, est positionné dans un plan perpendiculaire à la table. Dans cette position, les projections des billes vérifient la relation d'ordre ci-dessus. 2. Une projection numérisée du fantôme est acquise dans les conditions d'incidence que l'on veut calibrer. 3. L'image est transférée au processeur 4 de traitement d'image. 4. L'image est corrigée des distorsions géométriques éventuelles, par exemple avec une méthode telle que celle décrite dans la demande de brevet français NO 2 633 793 déposée le ler juillet 1988. 5. L'image corrigée I est ensuite soumise à un algorithme de détection, localisation, et labellisation automatique des billes décrit ci-après. 6. Les coordonnées des billes dans le repère 3D, R, et de leurs projections dans le repère 2D du plan P, sont exploitées pour estimer les paramètres géométriques de la projection: les coefficients de calibration. Un algorithme basé sur une minimisation d'erreur par la technique du gradient conjugué est proposé et également décrit ci-après dans ce but. On décrit ici un exemple possible de l'algorithme de calibration, qui permet de passer de l'image en projection du fantôme de calibration, figure 4a - image I, à des paramètres géométriques intrinsèques de la machine. On montre tout au long de cette description les images résultant des différents traitements. 1. D'abord on segmente l'image I. L'objectif est ici de supprimer dans l'image I le signal de la projection du cylindre de plexiglass qui supporte les billes. On utilise pour cela des outils de morphologie mathématique: transformation chapeau haut de forme, avec un élément structurant de forme ronde et de taille supérieure à la taille des billes dans l'image I. On obtient une image S1 qui ne contient plus que le signal correspondant à l'absorption des billes. Les outils de morphologie mathématique sont de type connu et publiés dans de nombreux ouvrages. 2. L'image S1 doit être purifiée des artefacts résiduels, provenant des bords du cylindre. On utilise pour cela de nouveau une transformation de morphologie mathématique visant à supprimer les structures verticales dans l'image. L'image résultat est une image S2 visible figure 4b en négatif. L'image S2 contient donc toutes les billes (signal positif) sur un fond à valeur nulle, au bruit près. 3. On segmente ensuite l'image S2. L'objectif est ici de ne garder que le signal provenant de la grosse bille centrale du fantôme. On utilise pour cela une érosion de l'image S2, d'une taille supérieure à la taille des petites billes mais inférieure à la taille de la grosse bille. On obtient une image S3, qui ne contient que la grosse bille. 4. On détecte ensuite toutes les billes dans l'image S2. On veut ici localiser au pixel près les projections de toutes les billes dans l'image. On commence à la ligne 0 de l'image (en haut du détecteur). On parcourt toute l'image du haut vers le bas. Un signal positif rencontré dans l'image est déclaré être une bille si ce signal est supérieur à un certain seuil 1. On donne le numéro 0 à la première bille rencontrée, puis le numéro est incrémenté à chaque bille rencontrée, jusqu'en bas de l'image. On est sûr que l'ordre des numéros est le bon, grâce aux propriétés de l'hélice décrite ci-dessus. 5. On calcule ensuite les coordonnées exactes des billes. L'objectif ici est de connaître les coordonnées de chacune des billes dans l'image, et ceci avec une précision sub-pixel. Pour chaque bille dans l'image S2, on effectue un calcul de centre de gravité en niveaux de gris de l'ensemble des pixels autour de la bille détectée, et tels que leur niveau de gris soit supérieur à un deuxième seuil, seuil T2 (seuil T2 6. On détecte ensuite la grosse bille dans l'image S3; on connaît alors ses coordonnées au pixel près. On compare ses coordonnées avec celles de toutes les billes détectées dans l'image S2, et on reconnaît ainsi dans la liste de toutes le bille celle qui est la grosse bille. 7. La grosse bille étant reconnue dans la projection, on peut maintenant renuméroter toutes les billes de façon absolue, et ainsi être assuré d'une correspondance exacte, pour chaque bille, entre son numéro et ses coordonnées 3D, connues par ailleurs, et les coordonnées 2D de sa projection. Connaissant, pour chaque bille, ses coordonnées 3D et les coordonnées 2D de sa projection, on peut calculer l'erreur de reprojection associée à un choix de valeurs arbitraires pour les 9 paramètres de la projection conique. Cette erreur est obtenue par la différence entre la position calculée avec ces valeurs arbitraires et la position trouvée dans l'image du fantôme de calibration. On peut donc estimer les paramètres ou coefficients de calibration associés au point de vue, à la position en rotation, que l'on veut calibrer en choisissant les paramètres minimisant un critère d'erreur de reprojection pour l'ensemble des billes détectées. L'algorithme proposé d'estimation des paramètres utilise la technique de minimisation du gradient conjugué appliqué au critère d'erreur quadratique moyenne. Cette technique classique nécessite le calcul du critère et de ses dérivées premières en fonction des paramètres. Cet algorithme est publié aux noms de W.H. Press, B.P. Flannery, S.A. Teukolsky and W.T. Vetterling," Numerical Recipes in C", Cambridge University Press, 1988, chapitre 10.6, page 317 et suivantes. D'autres algorithmes, décrits au même endroit dans cet ouvrage, et dénommés Simplex et Newton sont aussi utilisables. Le premier, Simplex, à l'avantage de s'adapter à n'importe quel système cohérent de paramétrisation. Il ne converge cependant pas rapidement. L'algorithme du gradient conjugué converge plus vite mais nécessite un choix de paramétrisation tel que les expressions analytiques des paramètres soient dérivables. C'est le cas de la paramétrisation présentée ci-dessus. Les calculs relatifs à cet algorithme on été présentés au congrès "3D Advanced Imaging Processing in Medicine" à Rennes en France en 1992 au cours de la 14ème conférence annuelle IEEE EMBS, du 2 au 4 Novembre 1992, par Anne Rougee, Catherine Picard, Cyril Ponchut et Yves Trousset. Les étapes en sont les suivantes: 1. Initialisation: les paramètres doivent être initialisés à des valeurs réalistes, ce qui permet de lancer la minimisation dans une région proche de la solution. 2. Normalisation: autour de cette valeur initiale on peut calculer la valeur d'un critère correspondant à un déplacement d'une unité sur chaque axe de l'espace des paramètres, puis effectuer un changement de variable dans cet espace par changement d'échelle sur chaque axe afin que tout déplacement unitaire sur chaque axe provoque le même effet sur le critère. 3. Minimisation: partant de la valeur initiale, et tenant compte d'un changement d'échelle trouvé au cours de l'étape de normalisation, on applique la technique itérative du gradient conjugué pour minimiser l'erreur quadratique moyenne sur l'ensemble des billes détectées. A chaque itération on calcule une nouvelle estimation des paramètres en fonction des dérivés du critère. Le processus est stoppé lorsque l'écart relatif entre l'estimation courante et la suivante est inférieure à un seuil donné. On appelle machine de reconstruction 3D une machine qui effectue des radiographies d'un objet sous différentes incidences obtenues en tournant autour d'un axe. Ces radiographies sont ensuite utilisées pour reconstruire une information 3D. Cette machine peut par exemple être construite à base d'un anneau rigide sur lequel tourne l'ensemble tube détecteur auquel cas l'axe de rotation est fixe et ne change jamais au cours du temps. Ou bien elle peut être construite à base d'un arceau inclinable supportant l'ensemble tube détecteur, auquel cas l'axe de rotation peut être positionné de différentes façons dans l'espace. Plutôt que de calculer pour chaque point de vue, pour chaque position en rotation de l'ensemble tube détecteur, les paramètres dits coefficients de calibration, on peut selon un procédé de l'invention, calculer des paramètres intrinsèques du système tournant. Ces paramètres intrinsèques sont alors des fonctions analytiques de la position en rotation de l'ensemble tournant. Il suffit ensuite, pour obtenir les paramètres de calibration relatifs à une position de calculer la valeur des paramètres intrinsèques pour une valeur donnée d'orientation du système tournant. Les systèmes de l'art connu délivrent déjà normalement une mesure, le plus souvent sous forme d'un signal électrique, représentant leur orientation dans l'espace. Il suffit d'utiliser cette mesure pour obtenir les coefficients de calibration relatifs à cette orientation en la portant dans les expressions analytiques des paramètres intrinsèques. Pour estimer les paramètres intrinsèques décrivant sous forme condensée la géométrie globale d'une machine de reconstruction 3D, on peut procéder de la façon suivante: estimer les paramètres de la projection conique, vue par vue, lors d'une première étape, puis estimer les paramètres intrinsèques à partir de l'ensemble de ces paramètres lors d'une deuxième étape. Première étape Dans le cas d'une machine de reconstruction 3D, la méthode de calibration automatique est particulièrement bien adaptée pour estimer les coefficients de calibration vue par vue au cours d'une rotation. On positionne l'hélice sur la table d'examen, son axe Z approximativement centré sur l'axe de rotation de la machine. Pour chacun des points de vue d'acquisition de la machine au cours d'une rotation, on va acquérir la projection de l'hélice, et traiter ainsi toutes les images acquises, avec la méthode décrite précédemment. La paramétrisation proposée - avec des coordonnées 3D cylindriques - est adaptée à l'estimation des paramètres le long d'une trajectoire circulaire autour de l'axe Z. Les étapes d'initialisation et de normalisation sont effectuées pour la première vue seulement. Les valeurs des paramètres trouvés pour la première vue sont ensuite utilisées pour initialiser l'estimation des paramètres de la vue suivante et ainsi de suite jusqu'à la dernière vue. Ceci permet d'accélérer les calculs. On obtient donc les paramètres géométriques pour chaque point de vue, qui sont stockés sur le disque numérique de la machine, et serviront dans toute étape de reconstruction 3D. Deuxième étape Il est possible selon l'invention de passer de la connaissance des coefficients de calibration, vue par vue, de la machine, à la connaissance de paramètres intrinsèques décrivant sous forme condensée la géométrie globale de la machine. Soit N le nombre de points de vues calibrés. Par exemple N vaut 100. Dans le cas où chaque point de vue est calibré indépendamment de son voisin, on dispose de 9N paramètres pour décrire la géométrie des acquisitions. Or dans le cas d'un système tube détecteur liés, qui effectue une rotation autour d'un axe fixe, on voit bien que ces 9N paramètres sont redondants. La géométrie d'acquisition de la machine de reconstruction 3D peut se paramètrer avec un nombre beaucoup plus limité de paramètres, qui sont; de préférence: - les coordonnées d'un isocentre O du système (3 paramètres); - la direction de l'axe Z' de rotation passant par l'isocentre 0.Elle est donnée par 2 angles (2 paramètres); - les rayons rI et rs des trajectoires de I et de la source S (2 paramètres); - La différence de cote (suivant l'axe de rotation Z') entre la trajectoire de la source et celle de I (1 paramètre); - l'écart angulaire E entre l'angle polaire de S par rapport à Z' et l'angle polaire de I par rapport à Z' (1 paramètre); - les 2 tilts a1 et a2 du détecteur, les angles d'orientation du détecteur par rapport à l'axe Z' (2 paramètres); - l'angle a3 entre les colonnes de l'image et la projection de l'axe Z' dans l'image (1 paramètre) On obtient donc un total de 12 paramètres. D'autres choix pour ces douze paramètres sont possibles. Celui-ci est préféré pour des considérations analogues aux précédentes: simplification des calculs.On notera que dans ce cas le référentiel dans lequel sont exprimées les valeurs de ces paramètres intrinsèques est indépendant de la position du fantôme. Il est donc transposable d'une calibration à l'autre. Une méthode pour estimer ces 12 paramètres en fonction des 9N paramètres de départ peut être par exemple la suivante: 1. Trouver la direction de l'axe de rotation Z': On considère les N points In (n=l,N) formant la trajectoire de I. En théorie, cette trajectoire est plane et orthogonale à Z'. On cherche donc le plan le plus proche, au sens des moindres carrés, de l'ensemble des In et on en déduit la direction de Z', qui est la normale à ce plan. Soient (Xn,Yn, Zn) les coordonnées d'un point In. On cherche le plan d'équation z = ax + by + c qui soit le "plus proche" de l'ensemble des In. On appelle d(a,b,c) la distance entre le plan z = ax + by + c et l'ensemble des points In. On recherche les coefficients (a,b,c) qui minimisent la distance d(a,b,c) . On définit: d(a,b,c) = z (axn + byn + c - zn)2 et on trouve alors très facilement le triplet (a,b,c) qui minimise d(a,b,c) par régression linéaire. 2. Trouver les coordonnées de l'isocentre, et le rayon de la trajectoire de I: L'isocentre est le centre à la fois de la trajectoire de I et de celle de S. On va donc rechercher ce centre en analysant la trajectoire de I. Remarque: dans cette étape, on cherche les paramètres (centre et rayon) de la trajectoire de I, puis dans l'étape suivante on cherchera le rayon de la trajectoire de S. On aurait pu, en toute logique, faire l'inverse, c'est à dire commencer par chercher le centre et le rayon de la trajectoire de S, puis en déduire le rayon de la trajectoire de I. Cette solution est toutefois moins robuste que celle exposée ici, car le point S étant plus éloigné que le point I de l'axe de rotation, les informations sur la position de S sont plus entachées d'erreur que celles sur la position de I. On utilise donc les informations "sûres" (celles sur I) comme base de la méthode et non pas celles sur S, qui sont plus incertaines. Une autre solution consisterait également à trouver de manière simultanée les rayons et centres des trajectoires de I et de S. Pour les mêmes raisons de robustesse, cette méthode n'a pas été retenue. Comment trouver les coordonnées de l'isocentre, c'est-à-dire du centre de la trajectoire de I? On définit un plan ff comme étant orthogonal à Z' et ayant une côte suivant Z' égale à la moyenne des côtes des points In. On projette orthogonalement sur ce plan, en des points In tous les points In . La trajectoire de I étant supposée circulaire, les In sont à priori situés sur un cercle. On cherche alors le cercle le plus proche, au sens des moindres carrés, de l'ensemble des I'n ce qui nous fournit: - le rayon rl de la trajectoire de I (c'est le rayon du cercle que l'on a trouvé), - les coordonnées de l'isocentre, qui est ici défini comme l'intersection de 7r et de l'axe de rotation (c'est le centre du cercle que l'on a trouvé). Comment définir et trouver le cercle le plus proche? Soit C(A,r) un cercle de centre A et de rayon r. Ce qu'il faut, c'est définir une distance d(A,r) entre le cercle C(A,r) et l'ensemble des 1,n Plusieurs méthodes sont possibles. Méthode 1 - figure Sa Soient (x'n,y'n} les coordonnées de Inv On définit: d(A,r) = z 1,non le point J'n étant défini comme l'intersection du cercle C(A,r) et de la droite AI'n. Méthode 2 - figure 5b Quand la machine de reconstruction 3D émet des tirs X suivant des positions angulaires régulières, c'est à dire tous les se degrés (Ne constant), on peut utiliser cette connaissance pour être plus robuste dans la détermination des paramètres du cercle. On va définir d(A,r) par: d(A,r) = z I'n J'n, le point J'n étant défini cette fois-ci comme le point du cercle C(A,r) qui possède, par rapport à une direction de référence, un angle polaire: e = e0 + (n-l) * se. Que l'on adopte la méthode 1 ou la méthode 2 pour définir la distance entre le cercle et le nuage de points, la suite est la même : on recherche, par une méthode quelconque de minimisation (gradient conjugué par exemple) les paramètres A et r du cercle (oui pour la méthode 2, les paramètres a, r et e0) qui minimisent la distance au cercle. Une fois trouvé ces paramètres, c'est fini, puisque - Le rayon de la trajectoire de J est connu : c'est r - L'axe Z' de rotation est entièrement connu, puisque l'on a déterminé son orientation à l'étape précédent et que l'on sait maintenant que cet axe passe par A. 3. Trouver le rayon de la trajectoire de S On connaît le centre de la trajectoire de S, puisque c'est par définition le même que celui de la trajectoire de I, que l'on vient de trouver à l'étape précédente. Cherchons maintenant le rayon. On considère les points S'n, projections des points S sur le plan v. Le rayon de la trajectoire de S est défini comme la moyenne des distances entre llisocentre O et les points Sln rayon de I = 1/N* S S'n 4. Trouver la différence de côte suivant l'axe Z' entre la trajectoire de la source et celle de I. C'est la moyenne des côtes suivant l'axe Z' des points Sn: voir figure 5c. On en déduit la différence de cote par la distance au plan n. 5. Trouver l'angle E: On travaille à partir des projections I' et S'n dans le plan n. Pour chaque n, on calcule les angles polaires (par rapport à O) de 1,n et S'n et on appelle En leur différence (voir figure 5d). On estime E comme étant la moyenne des en. 6. Trouver les tilts a1 et a2 du détecteur Pour chaque vue, on calcule dans le repère lié à l'axe de l'hélice, les coordonnées de la normale au détecteur. Ensuite, par application d'une formule de changement de repère, entre le repère lié à l'axe de l'hélice et celui lié à l'isocentre et à l'axe de rotation Z', on en déduit les coordonnées de la normale au détecteur dans le repère lié à l'isocentre et à l'axe Z' de rotation. Cette formule de changement de repère est déduite de la connaissance du repère lié à l'hélice. Par ailleurs on connaît l'isocentre puisque c'est le centre du cercle trouvé précédemment. On connaît également l'orientation de l'axe Z'. Les deux autres axes de ce repère sont figés à une rotation près. De préférence on fait en sorte, pour les figer, que la première vue soit affectée dans le nouveau repère d'un angle e0 lié à une position privilégiée de la machine. Dans ce nouveau repère, ces coordonnées ne dépendent que des 2 angles a1 et a2. Pour chaque n, on calcule ces 2 angles notés aîn et a2n. Ensuite a1 est défini comme la moyenne des aln et idem pour a2n. 7. Trouver l'angle a3 Pour chaque vue, on calcule dans le repère lié à l'axe de l'hélice, les coordonnées du vecteur unitaire parallèle aux colonnes du détecteur. Ensuite, par application de la formule de changement de repère, on en déduit ses coordonnées dans le repère lié à l'isocentre et à l'axe Z' de rotation. Or, dans ce nouveau repère, ces coordonnées ne dépendent que de l'angle a3. Pour chaque n on calcule cet angle, noté a3n. Ensuite a3 est défini comme la moyenne des a3n. Troisième étape On connaît maintenant 12 paramètres intrinsèques qui définissent la géométrie du système. A partir de la connaissance de ces 12 paramètres intrinsèques, il est très facile de recalculer, pour chaque vue n, les 9 paramètres locaux de calibration définissant cette vue. Ce calcul est d'autant plus simple qu'il a lieu dans le repère intrinsèque de la machine, c'est-à-dire dans le repère lié à l'isocentre et à l'axe de rotation. On appelle (rsn, asn, zsn, thetan, phin, psin, rin, ain, zin) les 9 paramètres locaux relatifs à la vue n, tels qu'ils sont définis dans l'article publié au Congrès de Rennes.Il vient alors rsn = rayon de la trajectoire de la source (= constante) asn = e0 + 180 + (n-l) * se zsn = différence de cote calculée précédemment (= constante) thetan = e 0 + (n-l) * se + al phi, = a2 psi, = a3 rin = rayon de la trajectoire de S (= constante) ain = e0 + (n-l) * se zin = 0 (= constante) Ces neuf paramètres ou coefficients de calibration sont donc déterminés à partir de l'incidence déduite de e0 + (n-l) se. Ils sont connus avec beaucoup de précision compte tenu des évaluations statistiques dont ils résultent. Avec des mesures précises de l'incidence on arrive alors à des coefficients de calibration calculés très précis. L'incidence elle même est déduite de l'indice n. Cet indice est connu par la division de la circonférence de la machine en un nombre donné de positions préétablies. On notera que, de préférence, on recalculera, à partir des expressions analytiques selon les coefficients intrinsèques, les coefficients de calibration relatifs aux N vues calibrées. De ce fait, les coefficients de calibration pour ces N (100) vues sont évalués avec autant de précision que pour les autres vues (les 900 autres). Le résultat d'ensemble est encore meilleur. Une autre méthode possible consisterait à estimer directement les paramètres globaux de la machine à partir de l'ensemble des coordonnées des billes détectées dans toutes les projections de l'hélice. Mais cette méthode est longue. REVENDICATIONS 1 - Système d'imagerie par rayons X avec un détecteur bidimensionnel situé en vis à vis d'un tube à rayons X, ce système d'imagerie par rayons X comportant un dispositif de calibration géométrique muni d'un fantôme de dimensions connues et de moyens de mesure pour déduire, à partir d'images en projection du fantôme sur le détecteur bidimensionnel, des coefficients de calibration de ce système d'imagerie par rayons X, caractérisé en ce que le fantôme comporte une succession, de proche en proche, de structures cellulaires à absorption radiologique contrastée par rapport à leur environnement, chaque structure cellulaire étant identifiable automatiquement par une caractéristique ordonnée d'identification de cette structure cellulaire. 2 - Système d'imagerie par rayons X selon la revendication 1, caractérisé en ce que, ce détecteur et ce tube tournant ensemble autour d'un axe, le fantôme comporte également un axe, un support, et une distribution des structures cellulaires sur ce support, de telle façon que des altitudes de ces structures cellulaires, mesurées selon cet axe du fantôme, soient monotones et croissantes avec un numéro d'ordre de ces structures cellulaires, l'axe du fantôme étant destiné à être orienté dans le système d'imagerie par rayons X, en phase de calibration, sensiblement parallèlement à l'axe de ce système d'imagerie par rayons X, et de telle façon que cette relation d'ordre se retrouve dans la projection du fantôme. 3 - Système d'imagerie par rayons X selon la revendication 2, caractérisé en ce que le support est un cylindre, de préférence circulaire, dont les génératrices sont parallèles à l'axe de ce cylindre et en ce que la succession ordonnée des structures cellulaires sur ce support a sensiblement la forme d'une hélice. 4 - Système d'imagerie par rayons X selon la revendication 3, caractérisé en ce que le support cylindrique est muni à ses extrémités de plaques de calage pour empêcher sa déformation. 5 - Système d'imagerie par rayons X selon l'une des revendications 1 à 4, caractérisé en ce que la caractéristique ordonnée d'identification de chaque structure cellulaire est constituée par un repérage de cette cellule par rapport à au moins une structure cellulaire particulière du fantôme. 6 - Système d'imagerie par rayons X selon la revendication 5, caractérisé en ce qu'une structure cellulaire particulière du fantôme possède une taille ou une absorption radiologique différente de celle des autres structures cellulaires et est placée sensiblement au milieu de la distribution des structures cellulaires. 7 - Système d'imagerie par rayons X selon l'une des revendications 1 à 6, caractérisé en ce que les structures cellulaires possèdent des absorptions radiologiques insuffisantes pour empêcher totalement le passage du rayonnement X à travers elles. 8 - Système d'imagerie par rayons X selon l'une des revendications 1 à 7, caractérisé en ce que les structures cellulaires comportent des billes métalliques placées au fond de puits réalisés dans un support. 9 - Procédé de calibration géométrique d'un système d'imagerie par rayons X à détecteur bidimensionnel tournant autour d'un axe, en vis à vis d'un tube à rayon X et comportant les étapes suivantes: - on place un fantôme de dimensions connues entre le tube et le détecteur, - on mesure, pour une position en rotation de l'ensemble tube détecteur par rapport au fantôme, dans l'image radiologique du fantôme projetée sur le détecteur bidimensionnel, des coordonnées de lieux images de points caractéristiques du fantôme, - on en déduit, pour cette position en rotation de l'ensemble tube détecteur, et dans un repère associé à la position du fantôme, des coefficients de calibration relatifs aux positions respectives d'un foyer de rayonnement du tube à rayons X et du détecteur bidimensionnel, - et on réitère ces deux dernières étapes pour des positions désirées différentes de l'ensemble tube détecteur, caractérisé en ce que - on traite les coefficients de calibration relatifs à un certain nombre de positions en rotation de l'ensemble tube détecteur pour en extraire des paramètres intrinsèques (~1, a2, a3...) du système d'imagerie par rayons X qui sont indépendants de la position en rotation de l'ensemble tube détecteur, - et on déduit de ces paramètres intrinsèques des coefficients de calibration (rsn, asn...) relatifs à une position quelconque en rotation de l'ensemble tube détecteur en évaluant les valeurs de fonctions analytiques de ces paramètres intrinsèques pour une valeur angulaire de position en rotation de l'ensemble tube détecteur. 10 - Procédé selon la revendication 9, caractérisé en ce que pour déduire des coefficients de calibration relatifs à une position désirée en rotation de l'ensemble tube détecteur, - a) on simule, avec des coefficients de calibration arbitraires, des coordonnées de lieux théoriques qui soient des images, sur le détecteur bidimensionnel, de points caractéristiques du fantôme, - b) on compare les coordonnées de ces lieux théoriques avec celles des lieux mesurés dans l'image du fantôme, - c) on modifie en correspondance les coefficients de calibration arbitraires, - d) et on réitère les deux étapes a) et b) jusqu'à ce que la comparaison montre un différence négligeable, auquel cas on retient comme coefficient de calibration pour cette position en rotation de l'ensemble tube détecteur les derniers coefficients de calibration modifiés. 11 - Procédé selon la revendication 10, caractérisé en ce que - on remplace les coefficients de calibration arbitraires par des coefficients de calibration déduits de la construction mécanique du système d'imagerie par rayons X. 12 - Procédé selon la revendication 10 ou la revendication 11, caractérisé en ce que - pour déduire des coefficients de calibration relatifs à une nouvelle position désirée en rotation de l'ensemble tube détecteur, décalée par rapport à une position en rotation précédente, on utilise à la place des coefficients de calibration arbitraires relatifs à cette nouvelle position en rotation, les derniers coefficients modifiés obtenus pour cette position en rotation précédente. 13 - Procédé selon l'une des revendications 9 à 12, caractérisé en ce que les coefficients de calibration sont exprimés par des fonctions analytiques des paramètres intrinsèques, ces fonctions analytiques étant décrites dans un repère cylindrique. 14 - Procédé selon l'une quelconque des revendications 9 à 13 caractérisé en ce que les paramètres intrinsèques sont obtenus par évaluation statistiques à partir d'un nombre restreint d'images du fantôme pour diverses orientation du système d'imagerie. 15 - Procédé selon l'une quelconque des revendications 9 à 14 caractérisé en ce que on déduit, à partir des paramètres intrinsèques, des coefficients de calibration relatifs à des positions en rotation de l'ensemble tube détecteur pour lesquelles on a effectué des mesures des coordonnées de lieux images de points caractéristiques du fantôme, et pour lesquels on en avait déduit des coefficients de calibration.