FR2700909A1 - Dispositif et procédé automatique de calibration géométrique d'un système d'imagerie par rayons X. - Google Patents

Abstract

On utilise un fantôme constitué d'une hélice. on montre que cette hélice permet une calibration géométrique automatique de n'importe quel système d'imagerie par rayons X à détecteur plan.

Description

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 < seuil T1). Un calcul de ce type est par exemple entrepris dans la deuxième demande de brevet français citée ci-dessus. Son principe repose sur le fait que les billes sont suffisamment fines pour ne pas absorber tout le rayonnement X. Alors, en l'absence de saturation d'absorption une précision inférieure au pixel d'image est obtenue.

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.

Claims (15)

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.