FR2888105A1

FR2888105A1 - Dispositif d'imagerie x ou infrarouge

Info

Publication number: FR2888105A1
Application number: FR0507362A
Authority: FR
Inventors: Norbert Beyrard
Original assignee: Individual
Current assignee: Individual
Priority date: 2005-07-11
Filing date: 2005-07-11
Publication date: 2007-01-12
Anticipated expiration: 2025-07-11
Also published as: FR2888105B1

Abstract

Dispositif d'imagerie X d'un corps comprenant, un support (1,3,5,7) pour recevoir un corps à examiner (9), une source (11) émettant un faisceau de rayons X, un détecteur (13) irradié par le faisceau, un convertisseur (15) des intensités détectées en données, un moyen (17) faisant tourner d'un angle de rotation le support (1-7) monté mobile autour d'un axe de rotation (19) par rapport à la source (11 ) et au détecteur (13) et un ordinateur (27) dûment programmé pour moyenner les données en n et m valeurs moyennes, construire une matrice initiale (n,m) avec les n et m valeurs moyennes, ajuster le coefficient d'atténuation en chaque zone élémentaire par une méthode des moindres carrés en tenant compte des valeurs moyennes, répéter les étapes précédentes pour des données acquises avec différentes paires d'angles de rotation, et moyenner terme à terme les matrices ajustées pour aboutir à une matrice de synthèse exprimant une image des coefficients d'atténuation du corps examiné (9).

Description

DISPOSITIF D'IMAGERIE X OU INFRAROUGE

L'invention se rapporte à un dispositif d'imagerie X ou infrarouge d'un corps comprenant plus particulièrement, un support pour recevoir un corps à examiner, une source émettant un faisceau de rayons X ou lumineux suivant une direction de propagation pour irradier ou illuminer le corps à examiner, un détecteur irradié ou illuminé par le faisceau pour détecter une intensité atténuée en considération d'une traversée des rayons X ou lumineux à travers le corps à examiner, et un convertisseur analogique numérique pour convertir les intensités détectées en données permettant de déterminer une atténuation par le corps à examiner des rayons X ou lumineux.

Un tel dispositif est connu par exemple du brevet US 3 924 131 ou US 3 919 552. Rappelons que la scannographie (ou tomodensitométrie), a été découverte par un ingénieur de la firme EMI, G. N. Hounsfield, en 1968. Le brevet de 1972 s'intitule: "A method and apparatus for examination of a body by radiation such as X or gamma-radiation". Cette invention valut en 1979 le prix NOBEL à son inventeur.

Le principe est le suivant: Un faisceau de rayons X, balaye un plan défini, il traverse un organe de manière linéaire et frappe une plaque ou un détecteur radiographique. La traversée de l'organe provoque une atténuation du faisceau, dont la mesure peut être effectuée grâce au détecteur. Le balayage de manière croisée, dans le plan de coupe, produit une série d'informations traitée par des logiciels appropriés sur un ordinateur associé.

En effet suivant chaque axe de balayage dans un milieu hétérogène l'atténuation peut s'exprimer par une loi exponentielle, tenant compte de l'absorption photoélectrique et de la diffusion par effet Compton.

Soit 10, la valeur de référence x la valeur en un point X, On peut écrire la relation suivante: 1=loe-F In=IoE r Ln = O = f o A (x) dx Dont on déduit par discrétisation: t Io _ + + +A Ln In = 0 A (x) dx - A,X, A2X2 X Les valeurs successives, AI, A2, An correspondent aux valeurs de chaque segment définis par X,, X2, Xn.

On peut alors exprimer par une série d'équations, les profils de chaque balayage associé à un angle (ou à une position) précis.

On peut définir une échelle particulière par la valeur, relative à une valeur de référence du coefficient d'atténuation, celle de l'eau par exemple, ou toute autre molécule convenablement choisie.

L'échelle, utilisée le plus souvent est celle relative à une molécule, abondante dans tous les organismes vivants, l'eau.

Si l'on appelle A (h2o), le coefficient d'atténuation de l'eau, on peut utiliser une échelle relative telle que: Bn = [An - A(h2o)] * 1000 / A(h2o) La valeur du coefficient de l'eau peut être définie comme égale à 1 ou à 0, créant ainsi un système de notation facile à utiliser puisque l'eau est un 20 composant essentiel du corps humain.

D'autres systèmes peuvent d'ailleurs être utilisées, en fonction de la manière dont l'information obtenue sera exprimée (de manière visuelle). On choisit souvent la valeur de 1000 pour l'os et de -1000 pour l'air.

Le traitement informatique d'un nombre suffisant de balayages croisés, 25 définissant en fait des petites cellules ou zones élémentaires, permet la résolution d'un ensemble d'équations linéaires à condition que le nombre de balayages soit égal aux nombre de cellules.

L'édition et l'utilisation de l'information sont effectuées par un ordinateur associé.

L'ordinateur collecte l'ensemble des données et calcule ainsi la valeur du coefficient d'atténuation de chaque zone élémentaire.

L'information résultant des calculs est traduite par une carte du plan de coupe tomographique.

L'ensemble des cartes constitue l'image scanner en trois dimensions de l'analyse, ce qui autorise des coupes longitudinales ou transversales.

L'interprétation médicale s'appuie ainsi sur une véritable image intérieure des tissus.

De telles images permettent de vérifier l'état de certains os, aussi bien que l'état du cerveau, pour détecter tumeur ou autre anomalie.

Les explorations sont précédées ou complétées par d'autres explorations par exemple échographie ultrasonore ou imagerie par résonance magnétique.

Le scanner et les méthodes qu'il a initiées restent un outil essentiel de l'exploration médicale.

A l'origine on procédait à une série de déplacements angulaires de l'ordre de 3 degrés, répétés une centaine de fois.

Les perfectionnements apportés depuis permettent d'associer plusieurs faisceaux à des barrettes de détection d'une longueur suffisante pour multiplier le nombre de mesures effectuées au même instant grâce à des détecteurs multiples.

Dans les scanners de cinquième génération on utilise des barrettes de détecteur perpendiculaires aux plans de coupe pour éviter tout déplacement.

L'image obtenue résulte d'un processus par étapes: - Obtention des valeurs des atténuations pour chaque projection, - Calculs des valeurs d'un profil, - Représentation matricielle de chaque plan de coupe, - Traduction de chaque représentation par une carte particulière, - Etablissement d'un système cartographique spatialisé.

On arrive aujourd'hui à des volumes de chaque zone élémentaire de l'ordre du mm 3.

C'est cependant loin de l'échelle microscopique puisque le nombre de cellules vivantes est de l'ordre du milliard dans un mm 3.

La détection précoce du cancer suppose un gain considérable de la définition. Mais le temps d'utilisation du système pour un patient déterminé ne saurait dépasser un seuil économique évident. Mais surtout l'accroissement du nombre de profils augmente la dose d'irradiation.

On admet en effet qu'un nodule cancéreux accélère son développement lorsqu'il provoque une vascularisation associée, ce phénomène apparaît lorsqu'une taille critique est atteinte soit par exemple 50 microns. Dans les procédés classiques la dose d'irradiation et le temps de calcul sont multipliés par 8000, pour atteindre ce niveau de finesse.

Le but de l'invention est de modifier un dispositif connu dans l'état de la technique rappelé ci-dessus pour réduire à la fois la dose d'irradiation et le temps de calcul de l'ordinateur dans le traitement des données pour et ainsi permettre d'augmenter la définition de l'exploration et des images en résultant.

A cet effet, l'invention a pour objet un dispositif conforme à celui rappelé en introduction, caractérisé en ce qu'il comprend - un moyen pour faire tourner d'un angle de rotation le support monté mobile autour d'un axe de rotation par rapport à la source et au détecteur montés sur un bâti ou pour faire tourner d'un angle de rotation la source et le détecteur montés sur un bâti mobile autour d'un axe de rotation par rapport au support et, - un ordinateur dûment programmé pour effectuer les étapes suivantes: (1) moyenner les données provenant de la conversion des intensités détectées, en n valeurs moyennes à l'intérieur de n intervalles résultant d'une première partition des données correspondant à un premier découpage de l'objet à examiner en n couches élémentaires parallèlement à la direction de propagation du faisceau pour un premier angle de rotation, et moyenner les données provenant de la conversion des intensités détectées, en m valeurs moyennes à l'intérieur de m intervalles résultant d'une deuxième partition des données correspondant à un deuxième découpage en m couches élémentaires parallèlement à la direction de propagation du faisceau pour un deuxième angle de rotation, de préférence différent de 90 degrés du premier angle de rotation, les découpages en couches élémentaires réalisant un quadrillage en n x m zones élémentaires d'un plan de coupe de l'objet à examiner défini par la première et la deuxième direction de propagation du faisceau pour respectivement le premier et le deuxième angle de rotation, (2) construire une matrice initiale (n,m) avec les n et m valeurs moyennes en affectant à chaque zone élémentaire un terme de ligne et de colonne (Bij) représentant un coefficient d'atténuation défini par la somme de la valeur moyenne sur l'intervalle (i) de même ligne que celle du terme, divisée par le nombre (m) de colonnes de la matrice initiale et de la valeur moyenne sur l'intervalle (j) de même colonne que celle du terme, divisée par le nombre (n) de lignes de la matrice initiale, (3) ajuster le coefficient d'atténuation en chaque zone élémentaire par une méthode des moindres carrés en tenant compte de contraintes imposées par les valeurs de bordures que constituent les valeurs moyennes en utilisant la formule suivante: Cij = Bij + n* (pj EBij) + m * (ci E; 1Bij) - nmJ ( E1,Pj IijBij) où, dans cette formule, Cij = la valeur recherchée 25 Bij = la valeur estimée initialement (n) = le nombre de lignes de la matrice initiale (m) = le nombre de colonnes de la matrice initiale n Ei -1Cij = pj pour toutes les valeurs de i, la contrainte de la colonne j E m Cij = ci pour toutes les valeurs de j, la contrainte de la ligne i, J =1 pour aboutir à une matrice redressée, (4) répéter les étapes (1) à (3) pour des données acquises avec différentes paires d'angles de rotation, et (5) moyenner terme à terme les matrices redressées pour les différentes paires d'angles de rotation pour aboutir à une matrice de synthèse exprimant une image des coefficients d'atténuation du corps examiné (9) sous une définition déterminée par le quadrillage.

De préférence, la source possède une embout d'émission large, par exemple de plusieurs centimètres de diamètre pour un embout cylindrique, émettant un faisceau également large pour irradier ou illuminer en une seule impulsion de commande de la source, une large étendue du corps à examiner et du détecteur, - l'ordinateur étant dûment programmé pour effectuer les étapes supplémentaires suivantes: (6) enregistrer les données provenant de la conversion des intensités détectées dans toute l'étendue irradiée ou illuminée du détecteur, (7) appeler les données correspondant à un plan de coupe particulier en sélectionnant, parmi les données enregistrées, celles qui proviennent de la conversion des intensités détectées dans une tranche de l'étendue irradiée ou illuminée des détecteurs, et (8) effectuer les étapes (1) à (5) à partir de ces données sélectionnées.

La première amélioration concerne l'acquisition des données, relatives aux coefficients d'atténuation ponctuels. L'auteur a choisi d'utiliser de préférence la projection de l'image sur un détecteur. Cette méthode est destinée normalement à obtenir une radiographie de l'objet à examiner, conduisant à produire sur un écran d'ordinateur, puis sur une imprimante une image de type radiographique.

L'acquisition de l'information peut se faire de deux manières: en produisant l'image radiographique puis en traitant celle-ci, ou en prenant l'information à la sortie du détecteur, en la traitant par conversion de l'information analogique en information digitale puis en traitant cette information.

On dispose de cette manière d'une masse considérable d'informations, puisque comme nous le verrons plus loin, un détecteur existant sur le marché permet de recueillir plusieurs millions d'informations élémentaires correspondant à des pixels dont la dimension est de l'ordre de 25 microns. La multiplication de l'information obtenue instantanément par cette technique conduit à des besoins de traitement eux mêmes considérablement amplifiés, ce qui a conduit à revoir le processus de traitement.

La deuxième amélioration concerne donc le traitement de l'information en substituant au traitement classique par l'algèbre linéaire, une autre méthode se résumant ainsi: a) on relève, par exemple, dans une première image obtenue une coupe ou tranche correspondant à une épaisseur de par exemple 1 mm et peut descendre jusqu'à 25 p m, cette tranche étant divisée en couches élémentaires réalisant un quadrillage de zones élémentaires de 1 mm2 de 1 mm. On peut grâce au logiciel de lecture de l'image mesurer la moyenne des intensités élémentaires, en déduire la moyenne des coefficients d'absorption par le corps à examiner, et ce pour chaque carré élémentaire et tirer le coefficient d'atténuation associé à chaque carré. De chaque coupe on déduit ainsi un premier vecteur.

Puis on prend une deuxième image obtenue sous un angle différent décalé de par exemple 90 degrés par rotation de l'objet ou de la prise de vue. On découpe de la même manière une bande située dans le même plan de coupe dont on tire un deuxième vecteur. On dispose ainsi de deux vecteurs orthogonaux.

b) Ces deux vecteurs orthogonaux (ou en tous cas sécants), peuvent permettre d'engendrer une première matrice initiale de la manière suivante: chaque vecteur est considéré comme une première série des valeurs des bordures de la matrice. Chaque ligne ou colonne de la matrice est divisée par un nombre égal au nombre de termes de l'autre vecteur, ce qui donne deux séries de termes élémentaires associés aux différentes lignes ou colonnes. On peut alors évaluer chaque terme de la matrice en prenant la moyenne arithmétique (ou géométrique) du terme correspondant à la ligne et du terme correspondant à la colonne sur lequel le terme est situé. On a ainsi créé une matrice initiale correspondant à un découpage millimétrique si les vecteurs sont eux mêmes obtenus ainsi.

c) Cette matrice initiale est ajustée en considérant les termes de chaque vecteur comme des contraintes de ligne ou de colonnes, par la méthode d'ajustement qui est décrite ci-dessous.

On peut traduire cette matrice de résultats par une image sur un écran donnant déjà une première image. Pour parfaire la recherche, on peut réaliser ne série de couples d'images décalées de 90 , d'où l'on tire autant de matrices élémentaires que de couples. On peut alors calculer une matrice de termes égaux à la moyenne des termes homologues des matrices individuelles. On peut enfin calculer un écart type pour chaque terme et si l'ensemble est satisfaisant produire l'image correspondant à la matrice des moyennes.

L'avantage du dispositif selon l'invention est double: -l'acquisition des données est beaucoup plus rapide et moins contraignante que par le système des balayages point à point. Si le détecteur a une définition des pixels de 25 microns approximativement (par exemple le détecteur de la société ATMEL dont la taille est de 23x6 cm offre sur sa surface 22 millions de zones de 25 microns de côté), il permet d'obtenir plusieurs millions de points même si l'on tient compte des espaces libres.

Par la haute résolution, un seul flash durant quelques millisecondes remplace plusieurs millions de balayages il en résulte que l'analyse par balayage physique est remplacée par une série de balayages informatiques.

Ainsi la quantité de radiations mise en jeu est particulièrement réduite.

- le traitement des données peut être extrêmement simplifié par le calcul de l'ajustement utilisé ainsi d'emblée, puis réitéré pour obtenir un nombre de matrices ajustées égales à p/2 pour p prises de vue. Ces prises de vue sous des angles régulièrement décalés au cours d'un tour complet sont appariées par couple décalés de 90 . L'image de synthèse est obtenue par superposition dans un même plan après rotation, et, par calcul des moyennes des valeurs homologues obtenues.

Pour une matrice de 20 000 termes, le temps de calcul est inférieur à la seconde au lieu d'être de l'ordre de l'heure sur un ordinateur PC disponible sur le marché.

Mais on objectera qu'une image obtenue par un faisceau large peut être de moindre qualité, en effet chaque rayon du faisceau légèrement conique produit une image influencée par les images des rayons voisins ou par les échos dus aux singularités fortes (par exemple aux inclusions métalliques) .

L'expérience montre que cette influence n'est pas considérable. De plus, un traitement mathématique peut faciliter l'analyse, en effet les valeurs obtenues sur une ligne ou une colonne de la matrice résultats peuvent être traitées comme valeurs redressables par un ajustement par exemple polynomial, ce qui a deux conséquences, les erreurs aléatoires sont lissées, mais on peut aussi améliorer la définition par interpolation entre les points mesurés et utilisation de la fonction de redressement en tenant compte de deux ou plusieurs interpolations croisées.

Cette recherche aboutit ainsi à une très grande économie de moyens techniques ou informatiques dans un domaine où le coût reste aujourd'hui trop élevé pour envisager dans beaucoup de pays des explorations systématiques dans la recherche de maladies comme le cancer.

La description qui va suivre expose donc la méthode de calcul d'ajustement proprement dite selon l'invention Cette méthode joue un rôle important dans le traitement des signaux résultant de la mesure, par les détecteurs radiographiques, de l'intensité ou de la valeur résiduelle du faisceau élémentaire produit par le dispositif à rayons X, après le parcours dans l'organisme à étudier.

Si l'on veut traiter une matrice de dimensions n lignes et m colonnes si on appelle BU la valeur estimée à la ligne i et à la colonne j, si on appelle Cij la valeur la plus probable du terme correspondant de la matrice, si on appelle pj la somme des termes de la colonne j, si on appelle cj la somme des termes de la ligne i.

L'estimation de Bij résulte soit du processus d'amplification matricielle, soit de toute autre méthode permettant une telle estimation, notamment à partir des techniques de l'ajustement linéaire ou polynomial.

Dans le cas présent, on recherchera la solution des valeurs de Cij en tenant compte des contraintes de lignes et de colonnes, c'est à dire le minimum de la fonction: E (Cij Bij)' Pour toutes valeurs de i et de j sous les contraintes: E Cij = pj pour toutes les valeurs de j Cij = ci pour toutes les valeurs de i La recherche d'un minimum de la fonction sous contraintes sera faite en utilisant la méthode des multiplicateurs de Lagrange, le Lagrangien s'écrira: L=- Eij (Cij Bij)2+E,, (Ln =, (Cijpj))+ gi( (Cij Cette fonction est composée de deux parties, une première n'a pas un 20 caractère gauche et la deuxième est un ensemble de relations linéaires.

Le Lagrangien est donc dérivable pour les variables Cij et j et fci, multiplicateurs de Lagrange associés aux contraintes de lignes et colonnes (nous disposons en effet de deux groupes de contraintes, les contraintes de ligne et les contraintes de colonne).

Dans ces conditions nous sommes en mesure d'obtenir un ensemble de relations linéaires concernant les Cij par dérivation du Lagrangien et un ensemble de valeur de relations relatives aux valeurs de contraintes, ce qui s'écrit: En précisant que le dL/dCij désigne une dérivée partielle, de la fonction L pour la variable Cij.

dL=-2(Cij Bij)+j+ i=0 1 dCij et les contraintes Cij = pj, pour tout j Cij = ci, pour tout i 1 a Cij = Bij+ X,j + i 10 2 L'ensemble des n*m relations correspondant aux dérivées partielles plus les n+m relations de contraintes est linéaire et n'admet qu'une solution correspondant aux nm+ n+m variables.

Si par exemple on veut traiter une matrice où n, le nombre de lignes, est 15 égal à 25 et m, le nombre de colonnes, est égal à 30, la solution par l'algèbre linéaire consiste à traiter: 750 variables Cij variables correspondant aux multiplicateurs de lignes, les p i variables correspondant aux multiplicateurs de colonnes, les a j Un premier objectif est déjà atteint puisque seulement 55 profils doivent être établis au lieu de 750.

Nous disposons au total de 750 relations correspondant aux dérivés partielles et de 55 relations correspondant aux contraintes, pour 805 variables. La résolution de ce problème en faisant appel au calcul matriciel est la solution la plus évidente mais elle implique des calculs très lourds, légèrement plus lourds que ceux impliqués par les méthodes classiques. L'auteur a d'abord tablé sur l'amélioration n r=1 m Ei=i rapide des processus de calculs, mais au delà de ce qui était son objectif essentiel, la limitation de taux d'irradiation en cours d'examen, il a poursuivi sa recherche pour tenter d'améliorer le temps de calcul.

De la combinaison des relations 1 et 2 on déduit: 5 En1Bij + 2 *+(En1,u i/2)=p1 lm Bij + *pi+(>m j/2)=ci 2 i=l On peut déduire de ces relations: xj= iJ *(2*(pj- n1Bj)- ni) n 2(ci-)-j) ( m) pi = Dans ces conditions en reportant par exemple la valeur de Xj dans pi, on obtient: Pour tout j = J*((pj - En,Bij) - n,u i) n Pour tout i pi _ W(Ein.=1,ui) + 21 (ci - m fil

E *

Si l'on précise que p- = (1/n) E (i=1à n) pi soit la moyenne des multiplicateurs liée à la contrainte des lignes, on est conduit aux deux relations suivantes: Pour tout j = 2 * (pj - En, Bi) - p - n Pour tout i pi = p- + ? * ( ci - En' Bij - mJ '-' En effet: m En' * En pi est égal à p - Dans ces conditions en reportant dans la relation: Cij = Bij + ('/2) * ( aj+pi) on est conduit à la relation algébrique Cette formule d'ajustement permet de déduire la matrice des Cij de la matrice des Bij, par un calcul terme à terme C i j = B i j + \n *(pj ElBij)+( * ( c i E m BU)- fl *( 1;0=1 p - leu) L'auteur a ainsi réussi à aboutir de manière tout à fait surprenante à un calcul de nature algébrique ne nécessitant pas le recours au calcul matriciel.

La méthode algébrique autorise le traitement partiel de la matrice de référence qui dans beaucoup de cas peut suffire.

La validation numérique de cette formule de traitement de signaux et 15 d'établissement des valeurs de définition de l'image recherchée sur le plan médical, apparaît ci après.

Exemple d'application de la méthode, sur un modèle réduit Soit donc une matrice de n lignes et m colonnes dans laquelle n=3, m=4

MATRICE INITIALE

1 2 3 4 lignes c 1 22 24 18 16 80 78 2 24 22 18 20 84 85 3 26 20 22 24 92 93 colonnes 72 66 58 60 256 P 70 67 59 60 256 Dans cette matrice, les valeurs estimées sont inscrites, sur les trois lignes et les quatre colonnes, les contraintes de lignes sont inscrites dans la colonne C. Les contraintes de colonnes sont inscrites dans la dernière ligne P. L'application de la formule ci dessus se simplifie puisque le total des contraintes de colonnes (ou de lignes) est égal à la somme des termes et conduit à:

EQUILIBRE APRES CALCULS

1 2 3 4 E2 C A 1 20.83333 23.83333 17.8333 15.5 78 78 0 2 23.5833 22. 58333 18.5833 20.25 85 85 0 3 25.5833 20.5333 22.583 24.25 93 93 0 2 70 67 59 60 256 P 70 67 59 60 0 0 0 0 On peut vérifier, en tenant compte de calculs réalisés par une simple calculatrice que la valeur des termes sommés verticalement ou horizontalement satisfont les contraintes, mais aussi qu'il conduisent bien aux résultats recherchés.

Equilibrage en utilisant la méthode de l'algèbre linéaire Cette méthode exprime directement les relations linéaires entre les variables Cij et Bij, et les variables X j et pi.

La méthode classique de résolution d'un système linéaire implique l'inversion de la matrice des coefficients des relations entre les variables et la multiplication par cette matrice inverse, du vecteur exprimant les seconds membres des relations.

a) II existe 12 relations entre les variables résultant de l'expression des dérivées partielles du LAGRANGIEN de la forme: Cij Xj / 2 - pi / 2 = Bij b) II existe 4 relations de contraintes relatives aux colonnes et 3 relations de contraintes relatives aux lignes, de la forme: E: Ci = nCa =ci i=1 J La méthode classique nécessite alors l'inversion d'une matrice d'une taille égale à n * m + n + m soit dans notre cas 19x19, dont le temps de calcul est 20 évidemment beaucoup plus élevé.

Pour permettre de vérifier le fonctionnement du dispositif, l'auteur a réalisé un prototype comportant les équipements suivants: - un faisceau à rayons X, de 70 kv, issu d'un tube en plomb d'une source 11 comprenant dont le diamètre est de 6 cm et constitue l'embout 29, - un détecteur 13 de petite dimension, de 45 x 30 mm, - un moyen 17 permettant de faire tourner un objet à examiner 9 placé entre le faisceau et le détecteur 13, avec un pas de 1, (voir figures 1 et 2), - un micro ordinateur 27 sur lequel est logé un logiciel d'imagerie médicale (de la société KODAK), qui permet d'obtenir une image radiographique à différents pas de rotation de l'objet à examiner 9, a) Dans un premier temps, 36 images ont été réalisées avec des pas de 10 , c'est à dire: 0, 10, 20, 90, 180, 270, 360 Une image de contrôle à 360 a donc été réalisée qui se superpose exactement à l'image à 0 , permettant de vérifier que la rotation a été effectuée convenablement. A titre d'exemple, nous joignons en annexe les images suivantes: figures 3a à 3d, les images à 0 , 90 , 180 , 270 des intensités détectées.

b) Dans un deuxième temps, on a procédé à une analyse par couples rectangulaires: - 0-90 ; 90- 180 ; 180-270 ; 270-360 (0) c) En effet on découpe à une certaine distance du plan de référence (que l'on voit apparaître dans les images ci-jointes en annexe, une tache noire qui est l'extrémité de la tige de rotation), une tranche étroite correspond à une épaisseur de 1 mm, soit un trentième de la largeur de chaque image. On obtient ainsi quatre bandes correspondant à 0 , 90 , 180 et 270 .

Chaque bande est à son tour découpée en carrés élémentaires correspondant dans l'espace à une zone élémentaire de 1 mm2, dans la partie utile de l'image qui est de 30 mm (approximativement). On lit grâce à un logiciel approprié la valeur moyenne de l'intensité des pixels dans chacune des zones de 1 mm2. On obtient alors un vecteur de trente (ou trente et un) termes, correspondant aux moyennes de chacune des zones élémentaires.

Comme on le constate à la lecture du tableau 4, les valeurs varient entre 0 et 255. Pour 255, le carré élémentaire est complètement blanc; dans ces conditions on admet que le coefficient d'atténuation le long de la ligne aboutissant au terme correspondant du vecteur, est égal à 255 moins la valeur de l'intensité lumineuse de la tache homologue. Dans ces conditions, s'il n'y a pas d'objet interposé, la valeur du coefficient d'atténuation est égale à 0. Si le carré est complètement noir, la valeur du coefficient d'atténuation est égale à 255, la totalité du rayonnement du faisceau a été absorbée.

On renouvelle le même traitement pour les bandes à 90 , 180 et 270 ce qui apparaît dans le tableau 4, dans les colonnes titrées VAL CORR (valeur corrigée). Les valeurs corrigées apparaissent pour les quatre coupes, 0 , 90 , 180 et 270 . Sur ce tableau, les moyennes sont estimées par un logiciel approprié ainsi que les écarts types qui restent dans des limites acceptables.

Dans la mesure où le nombre de pixels reconnu est élevé, on peut poursuivre l'analyse jusqu'à une définition largement inférieure au millimètre. On redresse les valeurs obtenues pour que la somme des valeurs de chacune des colonnes soit égale à la valeur moyenne, puisque dans les mêmes conditions de distance, l'absorption d'un corps donné est constante.

d) On divise par 31 le terme précédent pour en déduire la valeur de contrainte horizontale (90 ) ou verticale (0 ), ces valeurs apparaissent dans le tableau 5. Le tableau 5 (matrice initiale) fait donc apparaître: les sommes des lignes et des colonnes, - les contraintes tirées du tableau 4, - les écarts entre sommes et contraintes, - les écarts réduits par la division par 31 des écarts ci-dessus.

Cette matrice initiale apparaît dans le tableau 5, après calcul de chaque terme.

f) Si l'on veut obtenir rapidement une image à la définition de 1 mm, oneffectue le calcul d'ajustement, ce qui conduit au tableau de résultat 6.

On peut reprendre le processus que nous avons utilisé, qui consiste à créer 18 images à partir de 36 coupes effectuées dans 36 images radiophotographiques élémentaires et en calculant couple par couple les données estimées comme indiqué précédemment pour le découpage pour le premier angle de rotation égal à zéro associé au découpage pour le deuxième angle de rotation égal à 90 .

On superpose alors les 18 images en faisant effectuer une rotation de chacune des images (sauf la première), d'un angle égal à l'opposé de l'angle de rotation à partir du point 0 . Nous montrons à titre d'exemple 4 images, avant rotation, obtenues pour les couples: - 0 - 90 (figure 7), 90 - 180 (figure 8), 180 - 270 (figure 9), 270 - 360 (figure 10).

Si l'on agrége les tableaux après rotation, on est conduit aux tableaux 11, 12, et à l'image 13 (synthèse des images des figures 7, 8, 9 et 10).

Pour obtenir une image à haute définition, il a été nécessaire de créer un logiciel ayant les fonctions suivantes: - acquérir et enregistrer les images élémentaires et les données numériques associées comme celles apparaissant sur les figures 3a à 3d.

- définir graphiquement une tranche de l'étendue du corps à examiner 9 irradié par le faisceau de la source 11 à rayons X (voir figure 14).

- tirer dans cette tranche les valeurs associées aux images situées dans ce même plan, par exemple 36 coupes effectuées dans 36 images radiophotographiques élémentaires prisées à 36 angles de rotation, de 0 à 360 incrémentées de 10 .

- effectuer les calculs pour générer, dans ce cas 18 images, en associant les valeurs deux à deux pour des couples d'angle de rotation dont le décalage est de 90 .

- faire tourner les images en fonction du décalage d'un angle égal à l'angle de rotation de sorte à superposer de manière homogène les 18 images.

- produire une image de synthèse telle qu'apparaissant sur les figures 15 et 16 qui montrent clairement la coupe d'un os de poulet examiné dans le plan de coupe choisi. L'image peut être blanche sur fond noir ou inversement. Un léger halo apparaît qui peut être éliminé notamment si l'on tient compte du fait que l'axe de rotation et l'axe de symétrie de l'image sont légèrement décalés.

On peut rappeler que le calcul terme à terme résulte de la formule: H* Cij=Bij+(pj EBij)+ m *(ci 1Bij)-n l

J

(E7,Pj ijBIJ) Pour chaque couple ou paire d'angles de rotation décalés de 90 , nous avons obtenu une image différente que l'on peut associer en superposant les différentes images après rotation conformément à un procédé qui sera expliqué ici.

En effet il est possible grâce à un calcul simple d'obtenir un opérateur 5 simple.

Cet opérateur est constitué par la matrice inverse de la matrice suivante: COSINUS A - SINUS A 0 SINUS A COSINUS A 0 0 0 1 Donnons un exemple, exprimé par les trois tableaux ci dessous: Prenons par exemple un point de coordonnées égales respectivement à X=-2; Y= 1; Z= 0.

Les tableaux suivants montrent le processus de calcul: ANGLE EN DEGRE 270 ANGLE EN RADIAN 4.7124 SINUS -1 COSINUS 0 Matrice initiale 0 1 0 -1 0 0 0 0 1 Matrice inverse 0 -1 0 1 0 0 0 0 1 Cette matrice est multipliée par le vecteur X, Y, Z, ce qui conduit aux coordonnées suivantes: X'=-1;Y'=-2, Z'=0 Dans le cas exposé plus haut pour les planches 7, 8, 9, 10, nous avons été conduits à constater un décalage de l'axe de rotation et fait tourner cet axe, voir le tableau 17 pour retrouver ses coordonnées dans chaque tableau décalé de 90 , 180 , 270 de sorte à recentrer les tableaux et produire une image de synthèse ad hoc.

En tenant compte de ces décalages, nous avons généré une matrice initiale 29 x 29 (tableau 11) et une matrice résultat (tableau 12) qui se traduit par le graphe multicolore de synthèse (figure 13).

Le tableau 11 est obtenu, après rotation des matrices obtenues pour créer et recentrer les tableaux de résultats avec des angles de rotation.

/ 0 pour le graphe 90 180, / 270 pour le graphe 90 180, / 180 pour le graphe 180 270, / 90 pour le graphe 270 O. Ces rotations ont été simplifiées par l'utilisation d'un opérateur simple qui permet de calculer un tableau après une rotation de 90 , en prenant la matrice transposée du tableau initial en multipliant cette matrice transposée par la matrice dont tous les termes sont égaux à 0, à l'exception de la deuxième diagonale, comme par exemple la matrice ci-dessous.

0 0 0 1 0 0 1 0 0 1 0 0 1 0 0 0 Une rotation de 270 est ainsi le résultat de trois rotations successives de 90 .

II résulte de cette expérience que l'on ne peut envisager un nouveau type de dispositif d'imagerie X ou infrarouge tel que décrit à la figure 18 et constitué comme suit: - un plateau tournant 19, sur lequel une personne (ou un objet) à examiner est placé. Les personnes peuvent être normalement debout, mais si elles sont souffrantes, on peut les placer sur un siège articulé pour leur permettre de s'asseoir tout en restant inclinées pour que la région à examiner puisse l'être commodément grâce aux articulations du siège; - un faisceau émis par une source à rayons X 11 à axe horizontal déplaçable le long d'une potence verticale 23; - un détecteur 13 pouvant être placé verticalement ou horizontalement le long d'une potence verticale 25 selon la taille de la région à examiner; - un ordinateur 27 qui reçoit, via un convertisseur 15 analogique-numérique, instantanément les informations provenant du détecteur 13 pour chacun des flashs du faisceau au moment où le détecteur 13 et faisceau sont alignés sur un axe passant par l'axe de rotation 19 du plateau tournant 31.

Ce dispositif peut être complété par un deuxième couple faisceau détecteur placé à 90 du premier et fonctionnant de manière synchrone; Le fonctionnement de ce système peut être le suivant: - le plateau tournant 31 tourne pas à pas de sorte à ce que l'angle de rotation permette de prendre des images juste accolées les unes aux autres. Si par exemple la distance de l'axe de rotation 19 au détecteur 13 est de 75 cm, celui-ci effectuera un tour complet sur un cercle de l'ordre de 4.70 m; si le détecteur 13 est une plaque de 23 cm, une vingtaine de pas suffiront au cours d'un seul tour pour obtenir 10 d'images élémentaires à partir de couples ou paires pris à 90 l'un de l'autre; si le pas est de seulement 6 cm, le mouvement circulaire sera découpé de manière plus précise en environ 80 pas, permettant de produire 40 images à la définition voulue.

Dans le cas où l'on veut examiner des régions peu étendues verticalement on choisira de placer le détecteur 13 en position horizontale, si l'on peut se contenter de zones plus étendues verticalement on placera le détecteur 13 en position verticale, et dans certains cas on déplacera le patient pour bien cibler la zone à analyser.

Mais on peut aussi utiliser d'autres configurations de détecteurs pour toutes applications particulières sans que le système soit réellement alourdi.

On peut ajuster la taille du détecteur 13 en utilisant des plaques de détection différentes dont le coût est aujourd'hui raisonnable pour obtenir toute une gamme de possibilités.

Rappelons qu'un tour peut être effectué sans dommages en quelques secondes et dans le cas examiné on recueille 20 ou 80 fois de l'ordre de 12 millions de pixels d'information permettant de produire un nombre considérable de coupes densitométriques à la définition que l'opérateur pourra définir à son gré.

Bien entendu plus la définition est élevée, plus le temps de calcul le sera.

L'opérateur pourra choisir son protocole une fois les images stockées pour obtenir dans un premier temps des images à définition basse ou moyenne et affiner pour telle ou telle partie de l'organisme sa recherche en conservant les images créées au cours du balayage, le temps de poursuivre son investigation après avoir libéré son patient.

De toute manière, pour des images à basse définition elles sont susceptibles d'apparaître sur l'écran quelques secondes après l'investigation, et le cas échéant, l'investigation peut se poursuivre pour d'autres régions de l'organisme étudié.

La possibilité de disposer de deux ensembles faisceau détecteur décalés de 90 , permet d'obtenir en quelques fractions de secondes des images constituées d'informations prises au même moment, et en environ 2 à 3 secondes, le temps d'une rotation, on peut obtenir 10 ou 40 images de telle ou telle partie de l'organisme et dans ce cas disposer d'une possibilité de produire une série d'images décalées de quelques fractions de seconde avec des vitesses de défilement élevées.

On peut ainsi réaliser une vue dynamique de l'intérieur d'un organisme en percevoir les mouvements, en déduire des informations supplémentaires, par exemple sur la respiration, ou les mouvements cardiaques.

Le même dispositif est naturellement applicable à l'observation d'objets pour permettre: - un examen non destructif d'un objet quelconque; - et même une série d'images lorsque l'on veut observer des mouvements à l'intérieur d'objets clos, - en effet un objet physique peut être mis en rotation rapidement par exemple à une vitesse de 10 ou 25 tours par seconde.

- en utilisant un détecteur de 6 cm de large on réalisera 40 images par tours associables pour créer une image de synthèse parfaite et pour 25 tours par seconde obtenir une image dynamique équivalente à celle produite sur un écran de télévision.

Mais même pour des personnes on peut améliorer le caractère dynamique de la prise d'images en association plusieurs couples travaillant de manière synchrone deux à deux pour obtenir un nombre suffisant d'images de synthèse en un tour de une ou deux secondes de durée.

De tels dispositifs plus complexes seront dans un premier temps utilisés par la recherche.

Les faisceau X peuvent t être remplacés par un faisceaux lumineux par exemple infrarouge à fort pouvoir de pénétration.

Pour vérifier le fonctionnement du dispositif décrit, une routine a été ajoutée au logiciel de base. Cette routine permet de calculer point par point l'écart type relatif à ce point, et en divisant cet écart type par la valeur de référence, d'obtenir le niveau d'erreur relatif au point, puis la moyenne de ces niveaux d'erreur donc aussi la précision de l'image dans son ensemble ou par des zones particulières.

La planche 19 donne la carte des écarts types.

La planche 20 obtenue en tenant compte de cette carte pour toute la zone centrale donne une image dont la qualité peut être évaluée. Les résultats sont les suivants: l'erreur moyenne pour l'ensemble de l'image est de 1, 07 %, l'erreur moyenne pour la zone centrale où se situe la partie utile est de 0,08 %, soit une précision jamais atteinte. Si l'on réduit le nombre d'images de base de 18 à 9, l'écart type est multiplié par 2.

Claims

REVENDICATIONS

1. Procédé d'imagerie X ou infrarouge d'un corps dans lequel, un corps à examiner (9) étant reçu par un support (1,3,5,7;31), - on irradie ou illumine le corps à examiner (9) à l'aide d'une source (11) émettant un faisceau de rayons X ou lumineux suivant une direction de propagation, on détecte une intensité atténuée en considération d'une traversée des rayons X ou lumineux à travers le corps à examiner (9), à l'aide d'un détecteur (13) irradié ou illuminé par le faisceau, - on convertit les intensités détectées en données permettant de déterminer une atténuation par le corps à examiner des rayons X ou lumineux, à l'aide d'un convertisseur (15) analogique numérique, - on fait tourner d'un angle de rotation le support (1-7;31) monté mobile autour d'un axe de rotation (19) par rapport à la source (11) et au détecteur (13) montés sur un bâti (21, 23,25) ou tourner d'un angle de rotation la source et le détecteur montés sur un bâti mobile autour d'un axe de rotation par rapport au support, et - à l'aide d'un ordinateur (27) dûment programmé, on effectue les étapes suivantes: (1) moyenner les données provenant de la conversion des intensités détectées, en n valeurs moyennes à l'intérieur de n intervalles résultant d'une première partition des données correspondant à un premier découpage de l'objet à examiner en n couches élémentaires parallèlement à la direction de propagation du faisceau pour un premier angle de rotation, et moyenner les données provenant de la conversion des intensités détectées, en m valeurs moyennes à l'intérieur de m intervalles résultant d'une deuxième partition des données correspondant à un deuxième découpage en m couches élémentaires parallèlement à la direction de propagation du faisceau pour un deuxième angle de rotation, de préférence différent de 90 degrés du premier angle de rotation, les découpages en couches élémentaires réalisant un quadrillage en n x m zones élémentaires d'un plan de coupe de l'objet à examiner (9) défini par la première et la deuxième direction de propagation du faisceau pour respectivement le premier et le deuxième angle de rotation, (2) construire une matrice initiale (n,m) avec les n et m valeurs moyennes en affectant à chaque zone élémentaire un terme de ligne et de colonne (Bij) représentant un coefficient d'atténuation défini par la somme de la valeur moyenne sur l'intervalle (i) de même ligne que celle du terme, divisée par le nombre (m) de colonnes de la matrice initiale et de la valeur moyenne sur l'intervalle (j) de même colonne que celle du terme, divisée par le nombre (n) de lignes de la matrice initiale, (3) ajuster le coefficient d'atténuation en chaque zone élémentaire par une méthode des moindres carrés en tenant compte de contraintes imposées par les valeurs de bordures que constituent les valeurs moyennes en utilisant la formule suivante: \n) ( iJBI> ) où, dans cette formule, Cij = la valeur recherchée Bij = la valeur estimée initialement (n) = le nombre de lignes de la matrice initiale (m) = le nombre de colonnes de la matrice initiale Ein -1Cij = pj pour toutes les valeurs de i, la contrainte de la colonne j m Cij = ci pour toutes les valeurs de j, la contrainte de la ligne i, j =1 Cij Bij + (i lt *(pj ,Bi)+ *(ci 77! ;_iBij)- \ \MI f- 1 l' nm pour aboutir à une matrice redressée, (4) répéter les étapes (1) à (3) pour des données acquises avec différentes paires d'angles de rotation, et (5) moyenner terme à terme les matrices redressées pour les différentes paires d'angles de rotation pour aboutir à une matrice de synthèse exprimant une image des coefficients d'atténuation du corps examiné (9) sous une définition déterminée par le quadrillage.

2. Procédé selon la revendication 1, caractérisé en ce qu'on irradie ou illumine en une seule impulsion de commande de la source, une large étendue du corps à examiner (9) et du détecteur (13), à l'aide d'un embout d'émission (29) de la source (11) large, par exemple de plusieurs centimètres de diamètre pour un embout cylindrique, émettant un faisceau également large et, - à l'aide de l'ordinateur (27) dûment programmé, on effectue les étapes supplémentaires suivantes (6) enregistrer les données provenant de la conversion des intensités détectées dans toute l'étendue irradiée ou illuminée du détecteur (13), (7) appeler les données correspondant à un plan de coupe particulier en sélectionnant, parmi les données enregistrées, celles qui proviennent de la conversion des intensités détectées dans une tranche de l'étendue irradiée ou illuminée des détecteurs, et (8) effectuer les étapes (1) à (5) à partir de ces données provenant de la conversion des intensités détectées dans la tranche de l'étendue irradiée ou illuminée.

3. Programme d'ordinateur destiné à être chargé dans un ordinateur pour la mise en oeuvre d'un procédé selon la revendication 1 ou 2 au cours duquel, un corps à examiner (9) étant reçu sur un support (1,3,5,7;31), on irradie ou illumine le corps à examiner (9) à l'aide d'une source (11) émettant un faisceau de rayons X ou lumineux suivant une direction de propagation, on détecte une intensité atténuée en considération d'une traversée des rayons X ou lumineux à travers le corps à examiner (9) à l'aide d'un détecteur (13) irradié ou illuminé par le faisceau, on convertit les intensités détectées en données permettant de déterminer une atténuation par le corps à examiner des rayons X ou lumineux à l'aide d'un convertisseur (15) analogique numérique, et on fait tourner d'un angle de rotation le support (1-7;31) monté mobile autour d'un axe de rotation (19) par rapport à la source (11) et au détecteur (13) montés sur un bâti (21,23,25) ou tourner d'un angle de rotation la source et le détecteur montés sur un bâti mobile autour d'un axe de rotation par rapport au support, le programme d'ordinateur effectuant les étapes uivantes: (1) moyenner les données provenant de la conversion des intensités détectées, en n valeurs moyennes à l'intérieur de n intervalles résultant d'une première partition des données correspondant à un premier découpage de l'objet à examiner en n couches élémentaires parallèlement à la direction de propagation du faisceau pour un premier angle de rotation, et moyenner les données provenant de la conversion des intensités détectées, en m valeurs moyennes à l'intérieur de m intervalles résultant d'une deuxième partition des données correspondant à un deuxième découpage en m couches élémentaires parallèlement à la direction de propagation du faisceau pour un deuxième angle de rotation, de préférence différent de 90 degrés du premier angle de rotation, les découpages en couches élémentaires réalisant un quadrillage en n x m zones élémentaires d'un plan de coupe de l'objet à examiner (9) défini par la première et la deuxième direction de propagation du faisceau pour respectivement le premier et le deuxième angle de rotation, (2) construire une matrice initiale (n,m) avec les n et m valeurs moyennes en affectant à chaque zone élémentaire un terme de ligne et de colonne (Bij) représentant un coefficient d'atténuation défini par la somme de la valeur moyenne sur l'intervalle (i) de même ligne que celle du terme, divisée par le nombre (m) de colonnes de la matrice initiale et de la valeur moyenne sur l'intervalle (j) de même colonne que celle du terme, divisée par le nombre (n) de lignes de la matrice initiale, (3) ajuster le coefficient d'atténuation en chaque zone élémentaire par une méthode des moindres carrés en tenant compte de contraintes imposées par les valeurs de bordures que constituent les valeurs moyennes en utilisant la formule suivante: X11 /I*(PJ ,Bij)+ \n l r *(ci Bij) - -i nm Cij=Bij+ ( \m, * ( Ej, p j E JBiJ) où, dans cette formule, Cij = la valeur recherchée Bij = la valeur estimée initialement (n) = le nombre de lignes de la matrice initiale (m) = le nombre de colonnes de la matrice initiale n Ei -1= pj pour toutes les valeurs de i, la contrainte de la colonne j m Cij = ci pour toutes les valeurs de j, la contrainte de la ligne i, =1 pour aboutir à une matrice redressée, (4) répéter les étapes (1) à (3) pour des données acquises avec différentes paires d'angles de rotation, et (5) moyenner terme à terme les matrices redressées pour les différentes paires d'angles de rotation pour aboutir à une matrice de synthèse exprimant une image des coefficients d'atténuation du corps examiné (9) sous une définition déterminée par le quadrillage.