FR 3117652 A1 20220617 FR 2013342 A 20201216 La présente invention concerne un procédé de recalage d’un modèle 3D d’une pièce, ledit modèle 3D étant un nuage de points 3D qui constituent un échantillonnage de la surface 3D de la pièce ; de préférence associé à une triangulation de Delaunay, par rapport à une image radiographique 2D de ladite pièce, ledit procédé comprenant les étapes suivantes, mises en œuvre dans une unité de traitement : a) acquisition (E0) d’une image radiographique d’une pièce (1) ; b) détermination (E11), à partir de l’image radiographique d’une image , dite image d’épaisseur, dont les valeurs des pixels correspondent à une fonction linéaire de l’épaisseur de matière traversée par les rayons X visant le pixel ; c) détermination (E12) d’une image 2D des bordures de l’image 2D radiographique de la pièce à partir de l’image d’épaisseur ; d) détermination (E4) d’une position 3D de la pièce pour laquelle l’image 2D a été acquise, ladite position étant celle pour laquelle l’image 2D des bordures et des points silhouette projetés 2D de la pièce pour cette position diverge le moins. Figure pour l’abrégé : Fig. 1 Procédé de recalage pour le contrôle non destructif d’une pièce aéronautique DOMAINE DE L'INVENTION La présente invention concerne le domaine du contrôle non destructif de pièces industrielles, notamment dans le domaine de l’aéronautique au moyen d’images 2D obtenues par radiographie. Et plus précisément, l’invention concerne le recalage à l’aide de transformations de modèles informatiques sur ces images. ETAT DE LA TECHNIQUE Dans un certain nombre de domaines, il est important de connaître précisément la position d’un objet. C’est par exemple le cas dans le domaine du contrôle non-destructif (CND) de pièces notamment aéronautiques, et notamment lors de CND par radiographie numérique. Le CND par radiographie d’un objet comprend classiquement l’acquisition d’une image radiographique de l’objet, au cours de laquelle l’objet à contrôler est positionné entre une source de rayons X et un plan comprenant des détecteurs. Ainsi, chaque pixel formant l’image radiographique comporte l’information de la quantité de matière de l’objet traversée par les rayons X entre la source et le pixel. En pratique, on observe un décalage entre l’objet sur l’image radiographique et l’objet réel. Un tel décalage est souvent dû à un positionnement approximatif de la pièce par rapport au système d’imagerie. Or, il peut être important d’avoir une estimation précise de la position tridimensionnelle de l’objet à contrôler au cours de l’acquisition de la radiographie. En effet, cette information est cruciale pour l’application de certains procédés d’analyse automatique et l’obtention de résultats de CND précis et fiables. Pour pallier ce problème, une solution connue est d’effectuer un recalage rigide 3D/2D d’un modèle informatique de la pièce tridimensionnel avec une image radiographique bidimensionnelle de la pièce. Dans ce cas, le recalage rigide peut être défini comme la recherche de la transformation géométrique, impliquant rotation et translation, permettant de passer du modèle informatique à l’image radiographique de la pièce bidimensionnelle. Autrement dit, il s’agit de trouver la transformation géométrique qui permet de maximiser l’adéquation entre la structure contenue dans le modèle informatique et la structure contenue dans l’image radiographique. La transformation géométrique obtenue permet ensuite de déterminer la position précise de l’objet réel. Les procédés de recalage peuvent être extrinsèques, intrinsèques ou par calibration. Les procédés extrinsèques consistent à placer dans le champ de vision des objets-marqueurs, en plus de l’objet d’intérêt, dans l’unique but de recaler avec précision. Les procédés par calibration s’appuient sur la mise en œuvre de méthodes en amont de l’acquisition pour maîtriser au mieux la géométrie d’acquisition. Toutefois, les procédés de recalage extrinsèques et par calibration sont très spécifiques à l’objet imagé et sont potentiellement inutilisable sur l’objet d’intérêt. On utilise donc généralement les procédés intrinsèques qui n’utilisent que le contenu de l’image radiographique et du modèle informatique pour effectuer le recalage. Ces procédés consistent à mettre en correspondance certains éléments de l’image radiographique bidimensionnelle avec certains éléments du modèle informatique tridimensionnel. Parmi les procédés intrinsèques, on connait les procédés qui opèrent avec des radiographies numériques reconstituées (en anglais, « Digitally Reconstructed Radiography », (DRR)), qui sont des simulations de l’image radiographique obtenues par projection en perspective du modèle 3D de la pièce sur un plan bidimensionnel, Les procédés de recalage intrinsèques avec DRR cherchent la transformation géométrique du modèle informatique qui minimise une mesure de corrélation d’image entre la radiographie mesurée et les DRR, et requièrent donc une production intensive itérative de DRR par un simulateur de DDR. Le problème est que, dans le cadre des procédés de recalage intrinsèques avec DRR, la production intensive itérative de DRR engendre un coût calculatoire important et par conséquent, généralement un temps de calcul important. De plus, les simulations DRR présentent généralement des approximations et des incertitudes qui se propagent ensuite dans l’estimation du recalage et ne permettent donc pas d’obtenir la position précise de l’objet. Il existe aussi des procédés de recalage intrinsèques qui opèrent sans DRR mais nécessitent l’ajout de produit de contraste dans le sujet (domaine médical en l’occurrence), ce qui n’est pas envisageable dans le cadre du CND. Un but de l’invention est de remédier aux inconvénients exposés ci-avant en proposant, un procédé de recalage entre un modèle informatique 3D d’une pièce et une image radiographique 2D de la pièce permettant de déterminer la position précise de la pièce radiographiée lors de l’acquisition. Il est à cet effet proposé, selon un premier aspect de l’invention, un procédé de recalage d’un modèle 3D d’une pièce, ledit modèle 3D étant un nuage de points 3D qui constituent un échantillonnage de la surface 3D de la pièce ; de préférence associé à une triangulation de Delaunay, par rapport à une image radiographique 2D de ladite pièce, ledit procédé comprenant les étapes suivantes, mises en œuvre dans une unité de traitement : a) acquisition d’une image radiographique d’une pièce ; b) détermination, à partir de l’image radiographique d’une image , dite image d’épaisseur, dont les valeurs des pixels correspondent à une fonction linéaire de l’épaisseur de matière traversée par les rayons X visant le pixel ; c) détermination d’une image 2D des bordures de l’image 2D radiographique de la pièce à partir de l’image d’épaisseur ; d) détermination d’une position 3D de la pièce pour laquelle l’image 2D a été acquise, ladite position étant celle pour laquelle l’image 2D des bordures et des points silhouette projetés 2D de la pièce pour cette position diverge le moins. L’invention est avantageusement complétée par les caractéristiques suivantes, prises seules ou en une quelconque de leur combinaison techniquement possible : - l’étape d) comprend une minimisation d’une fonction de mesure de la divergence entre l’image 2D des bordures et les points silhouettes 2D projetés pour une position du modèle 3D de la pièce de manière à déterminer une position * du modèle 3D minimisant ladite fonction , la position * correspondant à la position 3D de la pièce pour laquelle l’image 2D a été acquise. - l’image est obtenue au moyen de l’opération suivante : avec la densité moyenne de la pièce, une image transmittance définie par où est l’intensité reçue par le capteur et est l’intensité émise par la source peut être estimée, par exemple, selon l’une des deux techniques décrites plus haut. - l’image des bordures est obtenue en filtrant (E12) l’image I L au moyen d’un filtre de Canny appliqué sur l’image d’épaisseur . - la fonction est minimisée par un algorithme d’optimisation globale et/ou par un schéma numérique de descente. - la détermination de la fonction comprend le calcul d’une transformée en distance DT d’un point silhouette à l’image des bordures. - la transformée en distance DT est définie à partir d’une distance, par exemple, une distance euclidienne, distance de Tchebichev ou une distance de type cityblock. - le calcul de la transformée en distance est effectué au moyen d’un interpolateur calculé à cet effet avant l’étape d’optimisation, conduisant à des évaluations rapides de la mesure de divergence f durant l’étape d’optimisation. - la fonction pour une position est définie par avec - un ensemble des points silhouettes 2D projetés ; - s un point de la silhouette projetée, appartenant à l’ensemble ; - un point de l’image des bordures ; - un poids associé au point s ; - une transformée en distance du point s par rapport à l’image de bordure. Selon un deuxième aspect, l’invention concerne un procédé de contrôle non destructif d’une pièce, de préférence aéronautique, comprenant les étapes suivantes : - acquisition par imagerie en 2D d’une image d’une pièce, l’image 2D comportant plusieurs pixels ; - traitement de recalage d’un modèle 3D sur l’image 2D acquise de manière à maitriser finement la géométrie d’acquisition de la pièce imagée, ledit traitement de recalage étant mis en œuvre au moyen d’un procédé selon le premier aspect de l’invention afin d’obtenir une estimation fine de la position de la pièce au moment de l’acquisition de l’image radiographique ; - analyse de l’image 2D couplée à la géométrie d’acquisition affinée au moyen du traitement de recalage, de manière à détecter dans l’image et localiser finement dans le modèle des anomalies ou à dimensionner finement certains éléments de la pièce visibles dans les images L’invention permet de réaliser un recalage performant d’un modèle informatique 3D d’une pièce par rapport à une image radiographique 2D au moyen d’une « approche silhouette ». La position de la pièce lors de l’acquisition de l’image radiographique 3D peut être déterminée de manière fiable et rapide. En effet, elle ne nécessite pas de simuler le processus radiographique et donc ne nécessite pas la production de DRR, ce qui permet de diminuer le coût calculatoire. De plus, l’invention ne s’appuyant pas sur des points pré-détectés, elle n’est pas spécifique et est donc utilisable sur toute pièce d’intérêt. L’invention ne nécessite pas d’apprentissage statistique. Avantageusement, le procédé de recalage selon l’invention est compatible avec une optimisation locale par schéma numérique de descente ou avec une optimisation globale. DESCRIPTION DES FIGURES D’autres caractéristiques, buts et avantages de l’invention ressortiront de la description qui suit, qui est purement illustrative et non limitative, et qui doit être lue en regard des dessins annexés sur lesquels : La illustre de façon schématique le dispositif selon un premier mode de réalisation de l’invention. La illustre un organigramme d’étapes d’un procédé selon un mode de réalisation de l’invention. La illustre un organigramme d’étapes d’un procédé selon un mode de réalisation de l’invention ; La illustre en détail des étapes du procédé de la ; La illustre en détail des étapes du procédé de la ; La illustre la notion de points silhouettes ; La illustre la notion d’image de bordure par rapport à la notion de points silhouettes ; Sur l’ensemble des figures, les éléments similaires portent des références identiques. DESCRIPTION DETAILLEE DE L'INVENTION Système de contrôle non-destructif La figure 1 illustre un système de contrôle non destructif d’une pièce. Un tel système comprend un système d’imagerie 10 permettant d’acquérir des images 2D d’une pièce 1 à contrôler. Des rayons X, issus d’une source S sont émis vers la pièce 1. Après avoir traversé la pièce 1, les rayons sont détectés par un détecteur D et formant une image 2D. L’acquisition est donc mise en œuvre par le détecteur D et au-dessus duquel la pièce 1 est placée. Il s’agit notamment d’un détecteur de rayons X. Ensuite, l’image 2D et un modèle 3D sont communiqués à une unité de traitement 11 pour mettre en œuvre un procédé de contrôle non destructif et en particulier un procédé de recalage qui va être décrit ci-après. Un tel modèle est un modèle CAO (pour Conception Assistée par Ordinateur) standard qui est généralement fourni par le concepteur de la pièce (un bureau d’étude par exemple). On peut noter qu’un tel model est également requis pour toute méthode intrinsèque de recalage 3D/2D, en particulier pour la génération de DRR. Il s’agit de préférence d’un modèle d’une pièce associé à une triangulation de Delaunay. On considère que le modèle 3D de la pièce 1 est un nuage de points 3D qui constituent un échantillonnage de la surface 3D de la pièce. La source de rayons X est préférentiellement équipée d’un collimateur (non représenté) de sorte que le rayonnement est émis sous la forme d'un faisceau plan d'axe horizontal, définissant un plan de coupe X, dit aussi plan de tir, passant par le foyer F de la source X et le plan X, positionné à l'horizontal. La source de rayon X peut, par exemple, émettre les rayons X selon une géométrie de projection conique (en anglais, « conebeam »). De manière préférée, une pluralité d'images radiographiques est obtenue par radiographie à rayons X obtenues selon des géométries d’acquisition différentes. Ces images sont alors successivement contrôlées afin de détecter d'éventuels défauts de la pièce aéronautique. Présentation générale d’un procédé de contrôle non-destructif d’une pièce par radiographie numérique En relation avec la un procédé de contrôle non-destructif d’une pièce, de préférence aéronautique comprend les étapes générales suivantes. Une pièce à contrôler est imagée au moyen d’un système d’imagerie par rayons X (étape E0) pour obtenir une image 2D radiographique de la pièce, l’image 2D comportant des pixels. Ceci permet d'obtenir une représentation 2D de la pièce 1. Un modèle 3D de la pièce à contrôler est également obtenu (étape E0’) de manière connue de l’homme du métier. Ensuite, un traitement (étape REC3D2D) de recalage du modèle 3D avec l’image 2D est mis en œuvre de manière à s’assurer que l’image 2D acquise est une représentation fidèle de la pièce ainsi recalée. Comme on va le décrire ci-après, ce traitement de recalage permet d’estimer avec précision la position de la pièce ayant engendrée l’image radiographique. En particulier, le recalage met en œuvre un critère de divergence entre le modèle 3D et l’image 2D acquise de la pièce afin d’obtenir une transformation affinée minimisant ce critère de divergence. Le critère de divergence est donc minimal (nul) lorsque les objets comparés sont identiques. Ainsi, pour tout point de l’image radiographique 2D, une correspondance peut être établie avec une projection 2D du modèle informatique 3D auquel a été appliquée la transformation optimisée. De façon plus générale, le recalage permet d’identifier avec précision la colonne de matière traversé par le rayon X ayant atteint un pixel particulier de l’image 2D. Une analyse (étape EANA) de l’image radiographique couplée à l’information de positionnement affiné de la pièce obtenue par le recalage peut faire l’objet de traitements complémentaires utiles pour le contrôle non destructif en tant que tel, tels que la détection d’anomalies et leur localisation 3D dans la pièce ou le dimensionnement d’éléments visibles dans l’image radiographique. Procédé de recalage La illustre des étapes d’un procédé de recalage. Le procédé de recalage consiste à déterminer une transformation géométrique rigide (translation et rotation globales) du modèle telle que la projection 2D simulée du modèle transformé correspond à l’image 2D acquise. Pour ce faire, il consiste à recaler une image des bordures de l’image radiographique 2D et des points silhouettes du modèle 3D M 3D transformé par une transformation T de paramètre dont on cherche l’optimum. L’ensemble des points silhouette 3D est déterminé à partir du modèle 3D M 3D transformé et correspond à l’ensemble des points de la surface 3D de la pièce pour lesquels les rayons X sont tangents à la pièce. Ainsi, on entend par silhouette un ensemble de points 3D déterminés à partir du modèle 3D et qui correspond à l’ensemble des points de la surface de la pièce dont les rayons X associés sont localement tangents à la surface de la pièce. De manière complémentaire, les points silhouettes 3D peuvent être projetés dans le plan image (selon la même géométrie de projection conique de production d’image radiographique) et forment l’ensemble des points dits points silhouettes2D projetés. L’image 2D des bordures est, quant à elle, déterminée à partir de l’image radiographique 2D I 2D et correspond à un masque binaire où à chaque pixel de l’image 2D est attribué une première valeur s’il appartient à une bordure de la pièce et une deuxième valeur s’il n’appartient pas à une bordure de la pièce Au cours de ce procédé, on optimise donc un paramètre permettant de faire coïncider le plus possible l’image des bordures à l’image de silhouette. Il s’agit du critère de divergence déjà mentionné plus haut. On décrit ci-après les différentes étapes du procédé de recalage. A partir de l’image 2D obtenue (étape E0), une image 2D des bordures est déterminée (étape E1). L’image 2D des bordures est calculée (étape E12) à partir d’une image I L (étape E11) dont chaque pixel est proportionnel à la colonne de matière traversée par le rayon X pour atteindre ledit pixel. De manière complémentaire à l’étape E1, on peut déterminer et exploiter les vecteurs normaux de chaque pixel appartenant à une bordure, dits normales au bordures, (étape E2), ce qui permet d’améliorer le recalage. Nous y reviendrons. En outre, on peut mettre en œuvre, de manière complémentaire mais non obligatoire, un calcul d’un interpolateur 4D des transformées en distances (étape E3), nous y reviendrons également. Ensuite, une détermination (étape E4) d’une position 3D de la pièce pour laquelle l’image 2D a été acquise est mise en œuvre, ladite position étant celle pour laquelle l’image 2D des bordures et une image silhouette du modèle 2D de la pièce pour cette position diverge le moins. En particulier, il s’agit de minimiser une fonction de mesure de la divergence entre l’image 2D des bordures et les points silhouettes 2D projetés pour une position du modèle 3D de la pièce de manière à déterminer une position * du modèle 3D minimisant cette fonction , la position * correspondant à la position 3D de la pièce pour laquelle l’image 2D a été acquise. Il s’agit donc de la position optimisée. On détaille ci-après chaque étape du procédé de recalage. Calcul de l’image des bordures E1 Comme illustré sur la figure 3, le procédé de recalage comprend une étape E1 de détermination d’une image 2D des bordures de l’image radiographique 2D de la pièce à partir d’une image proportionnelle à chaque colonne de l’image 2D radiographique. De manière avantageuse, l’image I L correspond à une image de la quantité de matière traversée par les rayons X. Détermination de l’image I L E11 L’image est déterminée à partir de l’image radiographique 2D I 2D . L’image radiographique 2D I 2D comporte une pluralité de pixels. La position de chaque pixel est donnée soit par son couple d’indices ligne et colonne, , soit par son couple de coordonnées continues dans le plan image. Le capteur est modélisé géométriquement comme une sous-partie rectangulaire du plan image décomposée en n rows lignes et un nombre de n cols colonnes, et j correspondent aux indices de ligne et colonne, respectivement. A chaque pixel de l’image radiographique 2D I 2D est associée une valeur de luminance qui correspond à la luminance reçue par le capteur associé au pixel. La luminance associée au pixel est donnée par la loi de Beer-Lambert et est définie de la manière suivante : , avec est l’intensité émise de la source de rayons X et représente le champ 1D de densité traversée par les rayons X sur une colonne de matière d’épaisseur avant d’atteindre le pixel (i, j). On se place dans un cadre approximatif d’une densité binaire, c’est-à-dire dans un cas l’on considère que la densité du matériau de la pièce à contrôler est constante dans la pièce. L’équation de Beer-Lambert donnant la luminance associée à un pixel peut être simplifiée de la manière suivante : , avec si la densité de la pièce est homogène et l’épaisseur cumulée de matière traversée par le rayon X associé au pixel e= . L’image , est notamment fonction d’une image de transmittance . L’image de transmittance est définie telle que, à chaque pixel e est associée une transmittance, la transmittance étant définie comme le rapport entre la luminance associé au pixel e et l’intensité des rayons X émis par la source : . Ensuite, à partir de l’image de transmittance , une image de niveau de gris est déterminée qui est linéairement liée à l’image de luminance : avec g un gain et Off un offset (ces paramètres sont intrinsèques au capteur utilisé). De manière avantageuse, l’image de niveau de gris est normalisée de sorte à obtenir une image de transmittance normalisée . L’image de transmittance normalisée peut être estimée par la formule suivante : Il est possible que les valeurs du gain g, de l’offset Off et de ne soient pas connues. Ainsi, en variante, l’image de transmittance associée à chaque pixel peut être estimée par : où est un paramètre tel que et devant être fixé ; est la valeur maximale de pixel rencontrée dans l’image de niveau de gris ou une valeur positive fixée arbitrairement correspondant au niveau de gris maximum que l’on s’attend à trouver dans l’image; est la valeur minimale de pixel rencontrée dans l’image de niveau de gris ou une valeur fixée arbitrairement entre et max[ ] et qui correspond au niveau de gris minimum que l’on s’attend à trouver dans l’image. Ainsi, l’image de transmittance est à valeurs dans un intervalle J = et, par construction, les pixels sont à valeurs dans l’intervalle J. Le paramètre est avantageusement introduit strictement positif pour anticiper un passage au logarithme lors d’une étape de linéarisation. Par la suite, l'image de transmittance normalisée est linéarisée pour enfin obtenir l’image proportionnelle à la colonne de l’image radiographique 2D. La linéarisation de l’image de transmittance normalisée permet de corriger la géométrie de l’image radiographique, déformée à cause de la fonction exponentielle de la loi de Beer-Lambert, par l’utilisation de la fonction du logarithme. On se place dans le cas où les matériaux de la pièce sont connus, et le paramètre peut être fixé. Dans un mode de réalisation, la pièce à contrôler peut-être constituée d’un matériau de densité constante, ce qui est souvent le cas des pièces industrielles élémentaires (non-assemblées). En variante, la pièce peut comprendre un matériau de densité variable ou plusieurs matériaux. Dans ce cas, le paramètre peut, par exemple, être une moyenne des densités des matériaux constitutifs de la pièce ou une densité maximum. Dans ce mode de réalisation, peut être fixé arbitrairement à 1 car l’important est que l’image soit proportionnelle à la quantité de matière traversée. L’image peut être alors estimée par la formule suivante . Calcul de l’image des bordures (masque binaire) (étape E12) Par la suite, l’image des bordures associée à l’image radiographique 2D est calculée (étape E12) à partir de l’image proportionnelle à chaque colonne de l’image 2D radiographique. L’image des bordures calculée à partir de l’image est plus précise que l’image des bordures calculées à partir de l’image de niveau de gris En effet, l’image linéarisée au cours de l’étape de linéarisation ne présente plus de déformation de la géométrie due à la présence de la fonction exponentielle dans la loi de Beer-Lambert, contrairement à l’image de niveau de gris . Dans cette application, l’image des bordures est obtenue par un filtrage de Canny, bien connu dans l’état de l’art appliqué à l’image . L’image des bordures ainsi obtenue est de même dimension que l’image et comporte des pixels p tels que : Pour chaque pixel de l’image des bordures : si le pixel de l’image associé au pixel de l’image des bordures appartient à une des bordures, alors le pixel est associé à une première valeur, par exemple 1 : . si le pixel de l’image associé au pixel de l’image des bordures n’appartient pas à une des bordures, alors le point-pixel p est associé à une deuxième valeur, par exemple 0 : p . Dans un mode de réalisation, on considère des points quadri dimensionnelles. Ainsi, à chaque pixel bordure de l’image 2D des bordures, c’est-à-dire tel que , on peut associer, en plus des coordonnées positionnelles 2D, les coordonnées 2D du vecteur gradient de l’image calculé en ce point. Calcul de normales aux bordures (étape E2) Dans un mode de réalisation 4D, le procédé de recalage peut comprendre le calcul des vecteurs normaux aux bordures (étape E2). Dans ce mode de réalisation 4D, à partir de l’image 2D des gradients, on calcule donc pour chaque pixel de bordure, un vecteur normal. Puis on construit une image 4D à valeur binaire (0 ou 1) dont les deux premières dimensions sont les dimensions spatiales et les deux dernières dimensions décrivent les coordonnées (régulièrement échantillonnées) du vecteur normal associé. Ces coordonnées sont échantillonnées régulièrement en 10 niveaux par exemple de manière à obtenir une image 4D de dimensions nrows x ncols x 10 x 10. L’image 4D vaut donc 0 presque partout sauf aux points de bordures ayant un certain vecteur normal ainsi calculé. Ce mode de réalisation 4D permet à la fois d’augmenter la précision du recalage et le bassin de convergence pour l’optimisation. Alors que dans le mode 2D, un point de bordure ou de silhouette projetée est seulement caractérisé par sa position 2D (coordonnées ligne et colonne), dans le mode 4D il est caractérisé par sa position et sont vecteur normal unitaire. En effet, les normales extérieures des points homologues de bordure et de silhouette projetée doivent avoir la même position et le même vecteur normal extérieur. Cette signature enrichie est donc plus discriminante et chaque point de la silhouette projetée a une plus grande probabilité d’être attiré par son point homologue de la bordure de l’image radiographique. Pour calculer l’image des normales aux bordures, un gradient est calculé en tout pixel de l’image appartenant aux bordures, c’est-à-dire, à tout pixel de l’image auquel est associée la première valeur dans l’image des bordures. Le gradient au pixel de l’image peut être calculé suivant la formule suivante : . Pour chaque pixel appartenant à la bordure, on estime un vecteur normal extérieur : - ;Dans un mode de réalisation, notamment dans le cas où l’image est trop bruitée, l’étape de calcul de l’image gradient peut être précédée d’une étape de filtrage de l’image . Le filtrage peut, par exemple, être un filtrage gaussien basse fréquence. Calcul de l’interpolateur 4D des transformées en distance (étape E3) La valeur du champ de bordure est connue sur une grille régulière (2D ou 4D selon le mode de réalisation). Par la suite, il sera nécessaire d’évaluer sa transformée en distance en tout point (2D ou 4D) du même espace de manière intensive. La réalisation intensive du calcul exact de la transformée en distance est globalement coûteuse et peut être avantageusement décomposé en deux étapes : un calcul préliminaire exact sur une grille d’échantillonnage et, pour l’évaluation intensive en tout point, un calcul rapide approximatif exploitant le résultat du calcul préliminaire. Plus précisément, le calcul préliminaire consiste à évaluer la transformée en distance du champ de bordure 2D ou 4D sur la même grille régulière que celle du champ de bordure (sans perte de généralité), puis à construire une fonction d’interpolation de ce champ de transformée en distance, appelée interpolateur. Lors du calcul intensif, l’estimation de la transformée en distance aux points (2D ou 4D selon le cas) d’intérêt appellera l’interpolateur (qui estime la transformée en distance en n’importe quel point 2D ou 4D à partir des points de la grille les plus proches). Dès lors à chaque point de la fonction de transformée en distance une valeur est associée. La grille 4D d’échantillonnage du domaine peut, par exemple, avoir les dimensions suivantes : avec une valeur permettant de paramétrer l’échantillonnage du vecteur normal en échantillonnant régulièrement et de -1 à 1 sur valeurs régulièrement espacées. Le vecteur est donc échantillonné sur valeurs vectorielles. Il ne peut néanmoins prendre qu’une certaine sous partie des valeurs possibles car le vecteur normal est unitaire. Sans perte de généralité, on peut fixer . Tout famille d’interpolateurs peut être utilisée sans perte de généralité (linéaire, splines, polynomiale, ...). Position 3D de la pièce (étape E4) Une fois obtenus l’image des bordures et/ou l’ensemble des vecteurs normaux aux bordures ainsi que l’interpolateur 4D selon les cas, le procédé comprend la détermination (étape E4) d’une position 3D de la pièce pour laquelle l’image radiographique a été acquise. En particulier, et comme déjà indiqué, il s’agit de minimiser une fonction de mesure de la divergence entre l’image 2D des bordures et les points silhouettes 2D projetés pour une position du modèle 3D de la pièce de manière à déterminer une position * du modèle 3D minimisant ladite fonction , la position * correspondant à la position 3D de la pièce pour laquelle l’image 2D a été acquise. C’est notamment par une stratégie d’optimisation itérative que cette position * est obtenue. Cette optimisation nécessite d’évaluer fonction en tant que telle faisant appel à des points silhouettes 2D projetés comme illustré sur la . Cette évaluation consiste à déterminer (étape E410) pour une position du modèle 3D de la pièce, un ensemble de points 3D silhouette du modèle 3D de la pièce ; puis de projeter ces points dans le plan image (étape E411) de manière à obtenir l’ensemble des points silhouette 2D projetée, et ensuite de déterminer (étape E412) une fonction f de mesure de la divergence entre l’image 2D des bordures et l’ensemble des points silhouette 2D projetée pour la position . Ces étapes sont impliquées dans l’évaluation de la fonctionnelle pour une position . Ces étapes sont répétées dans le cadre d’une stratégie d’optimisation itérative permettant de d’explorer itérativement plusieurs valeurs de jusqu’à convergence vers une position * du modèle 3D de la pièce minimisant la fonction f, la position * correspondant à la position 3D de la pièce pour laquelle l’image 2D a été acquise comme cela est illustré » sur la . Points silhouettes On introduit ici les points silhouettes qui sont utilisés dans la fonction . Ces points sont obtenus à partir du modèle informatique 3D M 3D , On se place dans le cas où le modèle informatique 3D M 3D est dans une position paramétrée par un vecteur . Le paramètre décrit la position du modèle informatique de la pièce dans un certain référentiel. Typiquement peut présenter six composantes ; trois pour le décalage 3D, trois pour les rotations 3D. Les points silhouettes de la pièce pour une position dans un système d’acquisition radiographique sont définis par l’ensemble des points de la surface de la pièce dont les rayons X associés sont localement tangents à la surface de la pièce. Pour déterminer ces points silhouettes du modèle informatique de la pièce associée à une position , on utilise une définition discrétisée de la silhouette (voir la figure 6) : à chaque côté C de triangle du modèle informatique, dit côté-silhouette, est associé un booléen spécifiant s’il appartient à la silhouette, par exemple 1, ou s’il n’appartient pas à la silhouette, par exemple 0). Le point est silhouette si et seulement si : . Puis chaque côté-silhouette est échantillonné en quelques points équidistants (sans perte de généralité), dits points-silhouette, de manière à pouvoir manipuler par la suite des points plutôt que des côtés. Chaque point-silhouette est associé au même booléen que le côté-silhouette auquel il appartient. Construction de la fonction à minimiser Dans un mode de réalisation 2D, on définit les éléments suivants. Soit l’ensemble des vecteurs paramètres régissant une transformation rigide Par exemple, le vecteur paramètre comprend les angles d’Euler et des coordonnées du vecteur translation : . Pour un point s, la transformation rigide donne une image du point s : où et où : , et . Soit l’ensemble des points-pixels appartenant à la bordure dans l’image de bordure, correspondant aux centres des pixels concernés. Soit un point de bordure de l’image radiographique, correspondant à un centre de pixel de l’image de bordures ; Soit l’ensemble des points appartenant à la silhouette, la silhouette dépendant de la position du modèle informatique. La position du modèle informatique est repérée par le vecteur paramètre . Soit un point-silhouette, c’est-à-dire, l’un des points-échantillons appartenant à un côté-silhouette. On définit une silhouette projetée comme une projection dans un plan image 2D de la silhouette 3D , la silhouette ayant subi une transformation rigide de paramètre . . Soit un point de la silhouette 3D tel que ; alors le projeté 2D de est donné par . Dans un mode de réalisation 4D, on définit de plus les éléments suivants. Soit le vecteur normal 3D en un point de la silhouette 3D, la silhouette projetée étant associée à un vecteur paramètre . On peut alors calculer un vecteur normal projeté 2D de de la manière suivante : Calcul de la pointe du vecteur normal 3D associé au point ; Calcul du vecteur unitaire 2D dans le plan image : Soit une nouvelle représentation 4D des points-bordures et des points-silhouettes projetée 2D. Dans cette représentation, deux dimensions localisent le point en 2D et deux dimensions définissent le vecteur normal 2D extérieur associé : A un point 2D de la bordure on associe un point 4D : où et sont les inverses des différences entre la valeur maximale et la valeur minimale des coordonnées en et en , respectivement. A un point 2D de la silhouette projetée 2D, on peut associer un point 4D : On définit un ensemble des points bordure 4D tel que ainsi qu’un ensemble des points silhouette projetée 4D tel que On obtient ainsi un ensemble des points bordure 4D ainsi qu’un ensemble des points silhouette projetés 4D. On définit une fonction comme une transformée en distance d’un point par rapport à un ensemble de point : où est une distance entre et . La distance peut, par exemple, être une distance euclidienne (norme de la différence), une distance de Tchebichev (norme de la différence) ou une distance cityblock (c’est-à-dire une norme de la différence). Par exemple, dans un mode de réalisation, peut être une transformée en distance 2D entre un point de la silhouette 2D projetée et l’ensemble des pixels bordures . Par exemple, dans le mode de réalisation 4D, peut être une transformée en distance 4D entre un point 4D de la silhouette 2D projetée et l’ensemble des pixels bordures 4D : : De manière complémentaire, la fonction f peut être robustifiée par un filtrage des adéquations anormales. Une fonction de robustification est définie et telle que pour tout point , un poids égal à 1 est affecté audit point si la valeur de sa transformée en distance est statistiquement normale ; dans le cas contraire, un poids égal à 0 lui est affecté. Dans ce cas, un pourcentage de points à prendre en compte et une probabilité maximale d’acceptabilité des points sont fixés. Un écart type des plus petites valeurs de transformée en distance des points tels que des poids sont positifs est aussi estimé (étape E413). Par exemple, dans le cas où la distance est une distance euclidienne, la pondération de robustification peut être définie de sorte que pour tout poids positif, on a où est la loi du chi-2 à 4 degrés de liberté et la fonction de répartition associée. Dans le mode de réalisation 2D, on utilise une loi du chi-2 à deux degrés de liberté. En variante, d’autres méthodes de robustification peuvent être utilisées telle que la fonction de Huber par exemple. Ensuite, on définit une fonction de pondération à valeur supérieure ou égale à . Dans un mode de réalisation, la fonction de pondération peut-être uniformément égale à . Dans un autre mode de réalisation, et notamment lors de l’utilisation d’un collimateur, la fonction de pondération peut être construite de sorte à améliorer la consistance de la fonction . En effet, lors de l’utilisation d’un collimateur, le collimateur protège une partie de l’écran-capteur hors cliché d’une exposition saturante. Dans ce cas, il n’est pas garanti que l’ensemble de la pièce soit projeté dans le champ de collimation (et donc certaines parties de la silhouette projetée peuvent être hors champ de collimation) ni que, durant le processus itératif de minimisation, les positions successives génèrent une silhouette projetée dans le champ de collimation. Sans précaution, la fonction peut donc être inconsistante dans le sens où elle peut atteindre un minimum local égal à en déplaçant la silhouette projetée en dehors du champ de collimation. Par exemple, afin d’améliorer la régularité de la fonction f et de contrôler le risque d’inconsistance de la fonction , il est possible de construire une fonction de pondération suivant les étapes suivantes. Dans un mode de réalisation, l’image 2D est décomposée en trois zones auxquelles sont affectés des poids non-normalisés : une première zone hors champ de collimation: si est localisé dans cette zone, alors . - une deuxième zone centrale du champ de collimation : si est localisé dans cette zone, alors . Une troisième zone qui est la couronne du champ de collimation : si est localisé dans cette zone, alors est déterminé par une fonction d’interpolation régulière (splines, polynomiale, etc.) entre les zones 1 et 2 et conduit à des valeurs entre et . A chaque tel que est affecté un poids en suivant la procédure de robustification. Cette étape permet de supprimer l’impact des points dont l’adéquation au point bordure correspondant est statistiquement anormale ; A chaque tel que est affecté un poids La fonction de pondération est ensuite normalisée : . L’évaluation de la fonction dépend du vecteur paramètre Dans le mode de réalisation 2D, la fonction f peut prendre la forme suivante : . Dans le mode de réalisation 4D, la fonction f peut prendre la forme suivante : où est le point 2D qui hérite de . Procédure d’évaluation rapide de la transformée en distance Dans un mode de réalisation, pour accélérer le temps d’évaluation de la transformée en distance associée à un vecteur paramètre de position du modèle informatique, une fonction d’interpolation sur est précalculée sur une grille 4D échantillonnée du domaine d’intérêt 4D avant le processus itératif d’optimisation : . La grille 4D d’échantillonnage du domaine peut, par exemple, avoir les dimensions dim suivantes : , où k est un paramètre d’échantillonnage. La valeur permet de paramétrer l’échantillonnage du vecteur normal en échantillonnant régulièrement et de -1 à 1 sur valeurs régulièrement espacées. Le vecteur est donc échantillonné sur valeurs vectorielles. Sans perte de généralité, le paramètre k peut par exemple être égal à 10. Puis, lors d’une étape d’optimisation de la fonction , l’évaluation de la transformée en distance en un point du domaine 4D peut être approximée par le calcul de la fonction d’interpolation au point . Toute famille d’interpolateurs peut être utilisée, par exemple linéaire, splines, ou polynomiale. Minimisation de la fonctionnelle (étape E42) Ensuite, on minimise la fonction (étape E42) afin d’estimer le vecteur paramètre optimal du modèle informatique permettant de maximiser l’adéquation entre l’image des bordures et l’image silhouette projetée. Autrement dit, tant que l’image de silhouette projetée et l’image des bordures ne coïncident pas, on déplace virtuellement la pièce dans le cadre d’un schéma d’optimisation jusqu’à ce que la silhouette projetée coïncide avec l’image des bordures. Afin de déterminer la position * , on minimise la fonction par un algorithme d’optimisation. Au cours du processus itératif de minimisation, la fonction est évaluée pour des vecteurs paramètres de position successifs . En particulier, à chaque itération , le vecteur position est modifier pour réduire la valeur de . Lorsque certains critères sont remplis (voir ces critères plus loin), l’itération est noté , on considère que la valeur de est optimale au sens de et on obtient donc la postition recherchée : . Ceci est résumé sur la . Types d’algorithmes d’optimisation Dans un premier mode de réalisation, la fonction peut être minimisée par un algorithme d’optimisation dit « schéma numérique de descente ». Le schéma numérique de descente peut, par exemple, être d’ordre 0 (algorithme du simplex). En variante, le schéma numérique de descente peut être d’ordre 1 ou 2, par exemple un algorithme du gradient ou un algorithme de Broyden-Fletcher-Goldfarb-Shannon (BFGS). Ceci n'est toutefois pas limitatif. Sans perte généralité, un critère d’arrêt de l’algorithme de schéma numérique de descente peut être que la décroissance de la fonction f entre deux itérations k et k+1 est non-significative. En variante ou de façon complémentaire, un critère d’arrêt de l’algorithme de schéma numérique de descente peut être que la variation de la variable entre deux itérations est non-significative : Ce mode de réalisation permet une optimisation fine de la fonction , mais nécessite une estimation initiale de la position de la pièce relativement proche de la position correspondant à la position 3D de la pièce pour laquelle l’image 2D a été acquise. Dans un deuxième mode de réalisation, la fonction peut être minimisée par un algorithme de minimisation globale tel qu’un algorithme stochastique, ou une approche combinatoire par exemple de type « branch and bound » en anglais. Dans un troisième mode de réalisation, la fonction peut être minimisée par l’utilisation séquentielle d’une optimisation globale (deuxième mode de réalisation), partant d’une valeur arbitraire et délivrant une valeur optimisée « grossièrement » , puis d’une optimisation locale (premier mode de réalisation) partant de la position est délivrant une estimation finale de la position. Le procédé de recalage comprend donc l’estimation du vecteur paramètre optimal du modèle informatique permettant de maximiser l’adéquation entre l’image des bordures et les points de la silhouette projetée 2D. La figure 7 illustre cette adéquation bordure image/silhouette projetée. Sur cette figure, le plus à gauche il s’agit d’une coupe de la pièce 1 (un cylindre). Les traits mettent en évidence des points de la coupe qui sont des points silhouettes et qui sont projetés sur un écran 2. Sur cette figure, l’image est représentée ainsi que l’image , dite image d’épaisseur, dont les valeurs des pixels correspondent à une fonction linéaire de l’épaisseur de matière traversée par les rayons X visant le pixel. Toujours sur cette figure, les points silhouettes projetés correspondant sont illustrés. Enfin, en bas l’image des bordures qui coïncide avec les points silhouettes projetés est illustrée. Procédé de recalage d’un modèle 3D d’une pièce, ledit modèle 3D étant un nuage de points 3D qui constituent un échantillonnage de la surface 3D de la pièce ; de préférence associé à une triangulation de Delaunay, par rapport à une image radiographique 2D de ladite pièce, ledit procédé comprenant les étapes suivantes, mises en œuvre dans une unité de traitement : a) acquisition (E0) d’une image radiographique d’une pièce (1) ; b) détermination (E11), à partir de l’image radiographique d’une image , dite image d’épaisseur, dont les valeurs des pixels correspondent à une fonction linéaire de l’épaisseur de matière traversée par les rayons X visant le pixel ; c) détermination (E12) d’une image 2D des bordures de l’image 2D radiographique de la pièce à partir de l’image d’épaisseur ; d) détermination (E4) d’une position 3D de la pièce pour laquelle l’image 2D a été acquise, ladite position étant celle pour laquelle l’image 2D des bordures et des points silhouette projetés 2D de la pièce pour cette position diverge le moins. Procédé de recalage selon la revendication 1, dans lequel, l’étape (E4) d) comprend une minimisation d’une fonction de mesure de la divergence entre l’image 2D des bordures et les points silhouettes 2D projetés pour une position du modèle 3D de la pièce de manière à déterminer une position * du modèle 3D minimisant ladite fonction , la position * correspondant à la position 3D de la pièce pour laquelle l’image 2D a été acquise. Procédé selon l’une des revendications 1 à 2, dans lequel l’image est obtenue (E11) au moyen de l’opération suivante : avec la densité moyenne de la pièce, une image transmittance définie par où est l’intensité reçue par le capteur et est l’intensité émise par la source peut être estimée, par exemple, selon l’une des deux techniques décrites plus haut. Procédé selon l’une des revendications 1 à 3, dans lequel l’image des bordures est obtenue en filtrant (E12) l’image I L au moyen d’un filtre de Canny appliqué sur l’image d’épaisseur . Procédé selon l’une des revendications 2 à 4, dans lequel la fonction est minimisée par un algorithme d’optimisation globale et/ou par un schéma numérique de descente. Procédé selon l’une des revendications 2 à 5, dans lequel, la détermination de la fonction comprend le calcul d’une transformée en distance DT d’un point silhouette àl’image des bordures. Procédé selon la revendication précédente, dans lequel la transformée en distance DT est définie à partir d’une distance, par exemple, une distance euclidienne, distance de Tchebichev ou une distance de type cityblock. Procédé selon l’une des revendications 6 à 7, dans lequel le calcul de la transformée en distance est effectué au moyen d’un interpolateur calculé à cet effet avant l’étape d’optimisation, conduisant à des évaluations rapides de la mesure de divergence f durant l’étape d’optimisation. Procédé selon l’une des revendications 2 à 8, dans lequel la fonction pour une position est définie par avec - un ensemble des points silhouettes 2D projetés; - s un point de la silhouette projetée, appartenant à l’ensemble ; - un point de l’image des bordures ; - un poids associé au point s ; - une transformée en distance du point s par rapport à l’image de bordure. Procédé de contrôle non destructif d’une pièce, de préférence aéronautique, comprenant les étapes suivantes : - acquisition (E0) par imagerie en 2D d’une image d’une pièce, l’image 2D comportant plusieurs pixels ; - traitement (REC3D2D) de recalage d’un modèle 3D sur l’image 2D acquise de manière à maitriser finement la géométrie d’acquisition de la pièce imagée, ledit traitement de recalage étant mis en œuvre au moyen d’un procédé selon l’une des revendications précédentes afin d’obtenir une estimation fine de la position de la pièce au moment de l’acquisition de l’image radiographique ; - analyse (EANA) de l’image 2D couplée à la géométrie d’acquisition affinée au moyen du traitement de recalage, de manière à détecter dans l’image et localiser finement dans le modèle des anomalies ou à dimensionner finement certains éléments de la pièce visibles dans les images.