Brevet FR2881520A1 — publié le 20060804

OBTENTION D'UNE IMAGE DE PHASE À PARTIR D'UNE IMAGE D'INTENSITÉ La présente invention se rapporte au domaine du 5 traitement d'image. La présente invention se rapporte plus particulièrement à un procédé d'obtention d'une image de phase à partir d'une image d'intensité. On s'intéresse depuis une vingtaine d'années à l'obtention d'une image de phase à partir d'une image d'intensité. Dans la publication Relationship between twodimensional intensity and phase in Fresnel Diffraction zone , Abramoshkin et al. (1989) fournissent la valeur du gradient de la phase en fonction de l'intensité, du gradient de l'intensité et des fonctions de Green dans l'espace infini. Plus précisément, dans la publication de Fornaro et al. Interferometric SAR phase unwrapping Using Green's Formulation , on divulgue l'utilisation de transformées de Fourier rapides pour calculer la phase en fonction du gradient de phase par les fonctions de Green. Cette formulation par transformée de Fourier est utilisée dans la demande PCT WO 00/26622 (Nugent et Al.) qui décrit un procédé d'obtention d'une telle image de phase. Ce procédé concerne la récupération de phase de champ d'ondes de rayonnement par résolution de l'équation de transfert d'intensité. On détermine d'abord le taux de variation de l'intensité de l'image, orthogonalement par rapport à une surface qui s'étend dans le champ d'ondes (c'est-à-dire en mesurant l'intensité dans deux plans séparés). Ce taux est ensuite soumis au processus de calcul suivant. établissement d'une transformée intégrale, multiplication par un filtre correspondant à l'inversion d'un opérateur différentiel, et établissement d'une transformée intégrale inverse. On multiplie le résultat par une fonction de l'intensité par rapport à la surface. Les filtres ont une forme fondée sur les caractéristiques du système optique utilisé pour l'acquisition des données d'intensité, du type ouvertures numériques et fréquences spatiales. Plus spécifiquement, la transformée de Fourier utilisée dans ce document de l'art antérieur est une transformée de Fourier à deux dimensions. Or, un des problèmes majeurs dans le domaine du traitement d'image, est le temps de calcul des applications représentant les résolutions des équations caractéristiques de l'image. Dans le présent cas, les temps de calculs proposés par le document de l'art antérieur ne sont pas satisfaisants pour l'homme du métier désirant effectuer ces calculs sur un grand nombre d'images. La présente invention résout donc ce problème technique en fournissant une nouvelle méthode d'obtention d'une image de phase à partir d'une image d'intensité. Par ailleurs, on notera que la demande PCT précitée ne 30 mentionne pas l'intervention de contraintes aux limites, en supposant seulement que la phase est nulle à l'infini. Un objet de la présente invention est donc de prendre en compte ces contraintes, tout en gardant un temps de calcul efficace. Pour ce faire, la présente invention est du type décrit ci-dessus et elle est remarquable, dans son acception la plus large, en ce qu'elle concerne un procédé de récupération de phase de champ d'ondes de rayonnement comprenant les étapes suivantes: (a) calculer une mesure représentative de la variation d'intensité dudit champ d'ondes de rayonnement sur une surface sélectionnée s'étendant globalement d'un bout à l'autre du champ d'ondes; (b) calculer une mesure représentative de l'intensité 15 dudit champ d'ondes de rayonnement sur ladite surface sélectionnée; (c) transformer ladite mesure de la variation d'intensité pour produire une première représentation de transformée intégrale et appliquer à ladite première représentation de transformée intégrale un premier filtre, pour produire une première représentation de transformée intégrale modifiée; ( d)appliquer un inverse de ladite première transformée intégrale à ladite première représentation de transformée intégrale modifiée pour produire une représentation non transformée; (e) appliquer une correction basée sur ladite mesure d'intensité sur ladite surface sélectionnée à ladite représentation non transformée; (f) transformer ladite représentation non transformée corrigée pour produire une seconde représentation de transformée intégrale et appliquer à ladite seconde représentation de transformée intégrale un second filtre; (g) appliquer un inverse de ladite seconde transformée intégrale à ladite seconde représentation de transformée intégrale modifiée pour produire une mesure de phase dudit champ d'ondes de rayonnement d'un bout à l'autre dudit plan sélectionné, caractérisé en ce que lesdits premier et second filtres sont des mesures d'une représentation d'une fonction de Green dans un espace de fonctions propres de l'équation de Helmholtz à coordonnées séparables, de sorte de que lesdites première et seconde transformées intégrales sont des transformées intégrales en une dimension de l'espace dudit champ d'ondes de rayonnement. Avantageusement, lesdites première et seconde 15 transformées intégrales sont produites en utilisant une transformée de Fourier, par exemple une transformée rapide. Selon une variante, lesdits premier et second filtre sont sensiblement les mêmes. Selon un mode de réalisation, lesdites surfaces sont sensiblement planes. L'invention concerne également un programme 25 informatique, éventuellement stocké sur un support, pour exécuter les différentes étapes du procédé selon l'invention. L'invention concerne également un appareil pour la 30 récupération de la phase de champ d'ondes de rayonnement comprenant les étapes suivantes: un moyen d'acquisition d'une mesure représentative de la variation d'intensité dudit champ d'ondes de rayonnement sur une surface sélectionnée s'étendant globalement d'un bout à l'autre du champ d'ondes; - un moyen d'acquisition d'une mesure représentative de l'intensité dudit champ d'ondes de rayonnement sur ladite 5 surface sélectionnée; - un moyen de traitement pour séquentiellement réaliser les étapes suivantes: (i) transformer ladite mesure de la variation d'intensité pour produire une première représentation de transformée intégrale et appliquer à ladite première représentation de transformée intégrale un premier filtre, pour produire une première représentation de transformée intégrale modifiée; (ii) appliquer un inverse de ladite première transformée intégrale à ladite première représentation de transformée intégrale modifiée pour produire une représentation non transformée; (iii) appliquer une correction basée sur ladite mesure d'intensité sur ladite surface sélectionnée à 20 ladite représentation non transformée; (iv) transformer ladite représentation non transformée corrigée pour produire une seconde représentation de transformée intégrale et appliquer à ladite seconde représentation de transformée intégrale un second filtre; (v) appliquer un inverse de ladite seconde transformée intégrale à ladite seconde représentation de transformée intégrale modifiée pour produire une mesure de phase dudit champ d'ondes de rayonnement d'un bout à l'autre dudit plan sélectionné, caractérisé en ce que lesdits premier et second filtres sont des mesures d'une représentation d'une fonction de Green dans un espace de fonctions propres de l'équation de Helmholtz à coordonnées séparables, de sorte de que lesdites première et seconde transformées intégrales sont des transformées intégrales en une dimension de l'espace dudit champ d'ondes de rayonnement. Selon un premier mode de réalisation, ledit moyen d'acquisition d'une mesure représentative de la variation d'intensité dudit champ d'ondes de rayonnement comprend au moins une plateforme mobile. Selon un second mode de réalisation, ledit moyen d'acquisition d'une mesure représentative de la variation d'intensité dudit champ d'ondes de rayonnement comprend au moins une lentille à focale variable. On comprendra mieux l'invention à l'aide de la description, faite ciaprès à titre purement explicatif, d'un mode de réalisation de l'invention, en référence aux figures annexées où : la figure 1 représente les données d'image selon 20 l'invention; la figure 2 illustre la comparaison des multiplications complexes entre la présente invention et les solutions de l'art antérieur, en particulier celle du document WO 00/26622 pour une image traitée allant de 64*64 pixels à 128*128 pixels; la figure 3 illustre la comparaison des multiplications complexes entre la présente invention et les solutions de l'art antérieur, en particulier celle du document WO 00/26622 pour une image traitée allant de 1024*1024 pixels à 2047*2047 pixels; - la figure 4 illustre les données d'images initiales à l'origine de l'algorithme selon la présente invention; - les figures 5A, 5B et 5C, qui doivent se lire de façon successive, représentent un mode d'implémentation exemplaire associé à la présente invention; - les figures 6 à 9 représentent des dispositifs 5 pour la mise en oeuvre de l'invention. La présente invention est décrite en détail ci-dessous dans un système de coordonnées cartésiennes et une image de type rectangulaire. Il est entendu que l'homme du métier est apte à adapter l'enseignement de la présente demande pour d'autres types d'images et/ou d'autres types de coordonnées, à la lecture de la présente description. NOTATIONS ET PRINCIPE GÉNÉRAL DE L'INVENTION La présente invention concerne de façon générale une méthode de reconstruction de la fonction e) de phase en un point r=ky) du plan 2D à partir d'échantillons de la différence du premier ordre aux points (x, y)= (nix,mAy). Nous calculons la phase seulement dans les points pour lesquels x = nAx et y = moy où les lignes d'une grille bi-dimensionnelle (2D) se croisent, c'est-à-dire dans les centres des points où se situent les éléments (pixels) d'un capteur de l'énergie lumineuse, par exemple une caméra CCD. Les différences de phase du premier ordre aux points (x, y)= (nAx,mAy) sont obtenues par le calcul du gradient de la phase à partir de la dérivée de l'intensité aI/az selon l'axe optique z. Ce calcul sera détaillé ultérieurement. Afin de faciliter la transition de la notation discrète vers la notation continue, le descriptif des termes utilisés dans la description du procédé selon l'invention est présenté dans le Tableau I annexé à la présente demande. Dans le document mentionné ci-dessus, Fornaro et al. ont utilisé la première identité de Green pour l'estimation de la phase à partir des différences de premier ordre présentées sous forme continue comme le vecteur du gradient F(r). La phase est alors calculée comme la somme d'une intégrale sur la région du support S de l'image de phase et sur la frontière C de cette région: Équation (1) . 0(r')= f)'(J).vg(r,r')dr _ f (r)[ns Évg(r,r')]dc où ns est un vecteur unitaire, perpendiculaire à C, orienté vers la surface extérieure de la région S comme illustré figure 1, est la fonction de Green. Dans l'équation (1), la phase 0 entre tant dans la partie gauche que dans l'intégrale sur le contour C. Cette double introduction de la même valeur dans deux parties de l'équation (1) ne permet pas d'obtenir la phase si elle n'est pas égale à zéro ou connue a priori sur le contour C. Ainsi, l'intégration de la phase sur le contour C signifie que la phase doit soit être connue, soit être égale à zéro, soit être pré-calculée par un autre dispositif optique, par exemple des réseaux de lentilles de Shack- Hartmann avant la reconstruction de la phase à l'intérieur de la surface S. L'intégration sur le contour dans l'équation (1) peut également être supprimée si la fonction de Green choisie satisfait les conditions aux limites de Neumann sur la frontière C de la surface S: Équation 2: ns - Vg(r, r')= 0 Or, l'équation (2) est valable pour la fonction de Green, donc la phase reconstruite selon l'équation (1) satisfait les conditions aux limites de Neumann. L'intégration sur le contour dans l'équation (1) peut aussi être supprimée si la fonction de Green choisie satisfait sur le contour C les conditions aux limites de Dirichlet et s'il est connu que la phase 00 est nulle sur C, ce qui s'écrit: Équation 3: g(Y,Y')= O, r EC, r'ES, 00= 0, r EC. Or, comme précédemment, l'équation (3) est valable pour la fonction de Green et pour la phase sur le contour C, donc la phase calculée selon (1) en utilisant la fonction de Green de l'équation (3) satisfait les conditions aux limites de Dirichlet. Par ailleurs, il est connu de l'homme du métier que la fonction de Green satisfait l'équation: Équation 4: vzg(@,@')= -b(@,@') où S(r,r') est la fonction de Dirac. Selon un aspect de l'invention, pour satisfaire les équations précédentes, on choisit une fonction de Green exprimée sous forme d'une série de fonctions propres de l'équation de Helmholtz comme suit. L'homme du métier remarquera qu'elle satisfera automatiquement les conditions aux limites (2) ou (3) selon la forme choisie ainsi que l'identité (4). Pour une région de support rectangulaire S, ces fonctions spéciales sont: Équation 5 ri Um n (Y)= Am n COS(nX)COS nmy m,n = 0,..., o0 pour les conditions aux limites de Neumann, et 15 Équation 6 Umn(É)= Amn Slll(nX) Sffi(11mi m,n = 0,...,o0 pour les conditions aux limites de Dirichlet. Pour un support S par exemple de dimension axb1 la constante Am n dans (5) et (6) est égale à a-yb_y,_ si m ou n sont nuls et 2a-b-Y' dans les autres cas, où a est la longueur, b est la hauteur de la région rectangulaire de support S comme sur la figure 1. On utilise aussi selon l'invention, les valeurs propres associées aux fonctions (5) et (6) qui sont: Équation 7 2 nit mn Qm n = Sn + nm 9 n = a 9 nm = b -, m, n = 0,...,00 2881520 11 La fonction de Green pour un domaine S s'écrit donc comme une série infinie aux fonctions propres de l'équation de Helmholtz: Équation 8: (r r,l_ _2 rl r,l gr, r / Qm,n Um, n (r l Um,n (r / m,n=O, a' ,, s0 La dérivation de l'expression générale (8) et la vérification des conditions aux limites (2) et (3) sont facilement réalisables par l'homme du métier pour une telle fonction. Ainsi, si la fonction de Green (8) est utilisée à la place de g(r, r') dans l'équation (1) l'intégration sur le contour dans l'équation (1) peut être supprimée et la phase peut être exprimé sous forme compacte comme: Équation 9: )=f F'(J).og(r,r)dS dans laquelle on notera simplement que le gradient est calculé par rapport à la variable r c'est-à-dire que CALCULS NUMÉRIQUES PAR TRANSFORMÉE DE FOURIER SELON L'INVENTION Un des objets de la présente invention est de calculer numériquement l'équation (9), à l'aide des fonctions de Green telles que définies pour obtenir la phase cherchée, tout en obtenant un temps de calcul acceptable. La solution est obtenue par la possibilité d'utilisation d'une transformée de Fourier rapide uni-dimensionnelle (1DFFT ou One Dimensional Fast Fourier Transform en anglais). L'utilisation de cette transformée en une dimension ne s'impose pas du tout à l'homme du métier au regard des équations générales ci-dessus, et nous fournissons ci-dessous les différentes étapes nécessaires à la compréhension du fait qu'il est possible d'utiliser une telle transformée en une dimension. Plusieurs transformations préliminaires sont nécessaires. Afin de faciliter les calculs de l'homme du métier désirant réaliser l'invention, nous mentionnons différentes étapes intermédiaires du calcul qui ne doivent pas être considérées comme des composantes essentielles de l'invention, mais simplement des moyens d'obtenir une forme d'équation sur laquelle il est possible de mettre en oeuvre les transformations en une dimension adéquate. Pour plus de clarté, on considère d'abord les conditions aux limites de Neumann, avec y=a/b et la fonction de Green dans l'expression de l'équation (8) avec les fonctions spéciales (5) : Équation 10: COS(e x) COS n.x (jr/a)2 +(2r/b)2 g(r, r,)= oo Am2 COS(tJmY)COS(YlmY m=0, m +n2eO oo 15 On sépare alors les cas n = 0, m e 0 et n 0, m =0 pour obtenir: Équation 11 g(r,r')= 2n-2y n-2 cos(e x)cos( xi) n= Or, on sait que: Équation 12 cos(nu) _ r cosh[a (ir u)] 1 1n2 + a2 2a sinh(aur Donc, Équation 13 1 -1 COS((nx)COS((nx') .?L y 1 a-b fm(xx) 2, m#0 n=1 (n/a + (m/bY 2m 2ym avec Équation 14 cosh[r1m (a x)]cosh(r)mx' if x > x lcosh[rim (a x')1 cosh(11mx if x (n/af + (m/bY m=1 2a2, Osus2ar, fm (x, x . sinh Or, Équation 15: n-2 cos(nu)cos(nv)= n 1 1 3u2+3(v JrY-?L2, if0susv 12 3v2 + 3(u -nf -7L2, if 0 s u s v Donc la fonction de Green est donnée par selon cette 5 variante de l'invention par: Équation 16: g(r,, r,)= 2z-' L m-' COS(flmY) COS(rimY/)fm (x, x') mal + y 3x2+3(x' af a2, if xsx 6a2 3(x'y + 3(x - ay - a2, if x' y' fn ' Y) sinh(Y nz) coshkn (b - Y) ' cosh(nY) if y En particulier, nous choisissons ces deux nombres comme les puissances de deux, plus grands que M +l et N +l, respectivement, M et N représentant les dimensions de 25 l'image à traiter. gx (r, r l)= -2a -1 n= Or, la présence des facteurs sin(r),ny) et sin( ,x) dans les équations (19) et (20), l'intégration (par rapport à y dans la première somme de la partie droite de l'équation (21) pour obtenir DX(th), et par rapport à x dans la seconde somme de la partie droite pour obtenir Dy( n)) permet de calculer les séquences DX(rO et Dyn) comme les parties imaginaires (prises avec le signe négatif) de la transformée de Fourier (FT) des dérivés (par rapport à y et par rapport à x, respectivement) de la phase. Ces deux séquences peuvent être calculées numériquement en utilisant l'algorithme de la transformée de Fourier Discrète, par l'algorithme de la transformée de Fourier rapide (1DFFT ou FFT 1D indifféremment) appliqué aux différences de phase de premier ordre Fy et FX, respectivement: Équations 22 et 23: Dx (r)m) = s 4fFT[F(n,m n = const j m Dy( )= {FFT [F (n, m) m = const ]n] J On notera que n=const signifie que le premier indice (i.e. n) est fixé dans Fy(n,m) tandis que la FFT 1D de la séquence résultante, dans laquelle le second indice (i.e. m) varie, est calculé. Inversement, m=const signifie que le second indice (i.e. m) est fixé dans F (n,m) tandis que la FFT 1D de la séquence résultante, dans laquelle le premier indice (i.e. n) varie, est calculé. Les indices 7m et n indiquent les fréquences discrètes pour lesquelles les FFT 1D doivent être évaluées. L'homme du métier comprendra par ailleurs facilement au regard des calculs ci-dessus que dans le cas des conditions aux limites de Dirichlet, la phase est donnée de 5 la même façon par Équation 24 M'-1 a 0(x', .y') = 2b -1 Sln('ImY') f { J m (x, xl)Dx \T1m)dx m=1 0 N'-1 b + 2a -1 Sln( nx') f fm (y,Y')Dy( )dy avec les parties réelles des transformées de Fourier 1D. Équation 25 Dx (re) = FT [Fy (n, m n = const Dy n)=NIFFT[Fx(n,m m= const CALCUL À PARTIR DE L'INTENSITÉ En référence encore une fois à la publication de Abramoshkin et Al. précédemment mentionnée, on connaît le 20 gradient de la phase en fonction de la fonction de Green Équation 26 v p(' = aIr)vgr, r')dr ) f az où la fonction de Green utilisée est la fonction de Green pour l'espace infini: Équation 27 g(r, r,)= - 1 lnlr -r, 2n Donc, en réutilisant l'équation (1), Équation 28 r, ai(") r r r (r") = kf I(r;z) vgr,r)dr +f (r) [rîs É vg(r, r') ]dr. En utilisant les fonctions de Green selon l'invention 15 (qui annulent la deuxième intégrale sur le contour), il vient donc: Équation 29: 0(') _ k Kr. zl-' aI(r'")og(r, r)dr É vg(r, r')dr S fs aZ avec, pour les conditions de Neumann: Équation 30 M'-1 a Mx, l l y,(x', y)= 2b -1 COS(ilmy') f xl)Px 1Îm ldx, m=1 o /, N'-1 b (px, (7C I, y . 2a -1 COS( n x' ) J m (y, y' ) Py) dy, n= 0 avec Équation 31: Px) _ Zs FFT al (n, m), n = const az Py (fin) = FFT al (n, m m = const az en On obtient également les fonctions équivalentes pour les conditions aux limites de Dirichlet. L'homme du métier comprendra que l'équation 28 est substantiellement celle utilisée par Nugent et al. dans la demande PCT WO 0026622. Dans cette demande, le deuxième terme est annulé en supposant que la phase est nulle sur le contour, alors que dans la présente demande, la phase peut être non nulle, le calcul est quand même réalisable par les conditions aux limites de Neumann et de Dirichlet annulant ce terme. Par ailleurs, dans cette demande, l'équation 28 est exprimée dans le domaine de Fourier et calculée par des transformées de Fourier en deux dimensions. Au contraire, la présente invention entend effectuer ce calcul de phase en utilisant des transformées de Fourier à une dimension, grâce à l'usage de fonctions de Green spécifiques, ce qui diminue 1m fortement le nombre d'opérations complexes, et donc le temps de calcul des algorithmes associés au procédé selon l'invention. Illustré figure 2 et 3, on a calculé le nombre de multiplications complexes entre la solution de la présente invention par transformée de Fourier 1D, et la solution de l'art antérieur du document WO 00/26622 en fonction de la taille de l'image traitée. Il est clair pour l'homme du métier que la solution de la présente invention économise grandement le temps de traitement de l'image pour la reconstruction de la phase. ALGORITHME DE TRAITEMENT POUR LA RÉCUPÉRATION DE LA 15 PHASE SELON L'INVENTION. Nous détaillons maintenant un exemple d'implémentation de la présente invention afin de compléter la description pour l'homme du métier. Pour cela, on utilise les notations suivantes: Abréviations: - CL Conditions aux limites; - 1D FFT Transformée de Fourier Directe unidimensionnelle; - 1D IFFT Transformée de Fourier Inverse unidimensionnelle; - f[...] Partie réelle d'un nombre complexe entre parenthèses; - s[...] Partie imaginaire d'un nombre complexe entre parenthèses. Les paramètres d'images en entrée suivants: - M x N, où M est le nombre de lignes de l'image, indexées par i = 1,..., M; et N est le nombre de colonnes, 5 indexées par j = 1, ..., N. Constantes: - a = N' = 21 a N + 1; - b = M' = 2' z M + 1. - Hmax = 225 (pour un ordinateur 32 bit). - y= a/b, Avy = 2r/b, Mmax = min{M', Hmax/y}; (3 = b/a, Avx = Jr/a, Nmax = min{N', Hmax/f3} La figure 4 illustre une première étape de l'algorithme selon l'invention dans lequel, à partir d'une image, on calcule ai/az par différence entre les images de l'intensité dans les plans z1 et 22. Cette différence est bien sûr calculée pixel par pixel sur les images I2(m,n) et I1(m, n) . L'homme du métier comprendra que, dans le cas de coordonnées quelconques, cette mesure de aI/az correspond en fait à toute mesure de la variation d'intensité sur une surface. De façon générale, l'homme du métier remarquera, en lien avec la description de l'invention ci-dessus, que l'algorithme donne la possibilité à chaque fois que cela est nécessaire, de choisir entre les conditions aux limites de Neumann et de Dirichlet. Illustré figure 5A, 5B et 5C, l'homme du métier reconnaîtra une implémentation des équations données plus haut selon la présente invention. On notera aussi que l'aspect discret des données réelles sur les images (les valeurs en x et en y étant obtenues par des données de pixels), permettent un calcul discret des intégrales données dans les équations et une interversion des calculs dans les sommes discrètes finies. Les modes de calculs de ces sommes discrètes sont données dans le présent algorithme à titre d'exemple et peuvent bien sûr varier selon les choix de l'homme du métier sans limiter la portée de la présente invention. Illustré figure 5A, les hachures obliques de l'image initiale représente la fenêtre de FFT utilisée pour les calculs. Cette fenêtre dépasse l'image utile et peut être remplie par des 0, selon la méthode connue du zero padding , ou remplissage par des zéros , ou par tout autre moyen bien connu dans le domaine du traitement des images. Les lignes verticales de la figure 1 indiquent alors le sens de la FFT 1D utilisée selon l'invention dans cette partie illustrée de l'algorithme. Dans ce mode de réalisation de l'invention, on calcule la FFT en fixant une colonne et en faisant varier l'indice de ligne. On notera par ailleurs figure 5A l'introduction du paramètre Mmax utilisé de préférence selon ce mode d'implémentation pour éviter la divergence des termes en sinus hyperbolique et un calcul spécifique du coefficient cy(m,n) au cas où l'indice représentant les coordonnées dépasse cette valeur seuil. On reconnaîtra à l'opération A l'obtention de la dérivée de la phase selon y tel qu'il est décrit dans la première partie de l'équation 30. On notera que l'opération A nécessite la division par une mesure de l'intensité qui peut faire diverger la solution si cette intensité est nulle. Dans ce cas, il est possible d'effectuer préalablement une détermination des zéros de l'intensité pour corriger cette divergence, et/ou d'affecter une valeur minimale à l'intensité en cas de nullité en un point de l'espace considéré. Selon l'équation 22, on applique ensuite le procédé de calcul selon l'invention utilisant une transformée de Fourier 1D au gradient de la phase selon y pour obtenir le premier membre de la phase comme dans l'équation 21. Ce premier membre cp2(m,n) est obtenu à l'opération B de la figure 5C. De la même façon (non illustré), on obtient la dérivée de la phase selon x tel qu'elle est décrite dans la première partie de l'équation 30, et selon l'équation 23, on applique ensuite le procédé de calcul selon l'invention utilisant une transformée de Fourier 1D au gradient de la phase selon y pour obtenir le second membre de la phase comme dans l'équation 21, cp1(m,n). La phase reconstruite selon l'invention est alors calculée par cp(m,n) = cp1(m,n)+ cp2(m,n). De façon générale, illustré figure 5A, 5B et 5C, l'algorithme selon l'invention comprend donc, suite au calcul d'une mesure de la variation d'intensité sur une surface (aI/az dans le mode de réalisation illustré), l'application d'une première transformation intégrale, de préférence sous la forme d'une transformée de Fourier rapide, pour obtenir la variable nommée (rr(m,n) dépendant des conditions aux limites utilisées. Cette variable est ensuite multipliée (dans l'espace de Fourier ou plus généralement dans l'espace de projection) par un filtre correspondant aux fonctions fY n et Ln, modifiées éventuellement par les contraintes de calcul de type MmaxÉ Conformément aux équations 19 et 20, ces filtres et fY sont des représentations d'une fonction de Green dans un espace de fonctions propres de l'équation de Helmholtz, 10 par exemple les fonctions associées aux équations 5 et 6. On applique ensuite une transformée intégrale inverse, par exemple de type transformée de Fourier inverse, en prenant les parties réelles ou imaginaires selon les types de conditions aux limites. Comme l'indique par exemple l'équation 26, la division par une mesure de l'intensité, par exemple celle prise sur un plan de calcul de la variation de l'intensité, permet d'obtenir le gradient de la phase par rapport à une des coordonnées. En conformité avec les équations 21, 22, et 23, on réapplique ensuite une transformée intégrale à une dimension, par exemple de type Fourier, on réapplique un second filtre à la variable obtenue (DX(m,n) dans le mode de réalisation proposé), et on reprend une transformée intégrale inverse pour obtenir sensiblement un des termes de l'expression définissant la phase recherchée. Il est entendu qu'au regard de l'exemple de mise en oeuvre de l'algorithme donné dans les figures et décrit ci- f Y dessus, l'homme du métier est apte à apporter des modifications d'implémentation, par exemple selon les capacités de calcul des machines utilisées, ou selon une caractéristique spécifique de l'image, sans toutefois sortir du cadre de l'invention. AUTRES MODES DE RÉALISATION ASSOCIÉS À L'INVENTION. Au regard de la description faite ci-dessus, d'une part des équations nécessaires à la mise en oeuvre de l'algorithme de l'invention, et d'autre part de l'algorithme de traitement lui-même, l'homme du métier comprendra que le mode de réalisation décrit n'est qu'un exemple non limitatif de la mise en oeuvre de l'invention. En particulier, on comprendra que dans le cas d'un support non rectangulaire, l'homme du métier est apte à déterminer les fonctions propres de l'équation de Helmholtz pour une forme de support spéciale. À titre indicatif seulement, on donne la forme de ces fonctions dans le cas d'un support circulaire Équation 5' : Uni,n(p,cp)=[cos(mcp)+sin(m(p)] m(ramnp/a) m,n=0,K zoo où Jm est la fonction de Bessel d'ordre m avec J, \rama)= 0, m, n = O,K, oo, i. e. na,nn est le n-ième zéro de la dérivée de la fonction de Bessel d'ordre m. L'équation 8 permet alors d'obtenir la fonction de Green associée. L'utilisation de la fonction de Green associée à cette fonction propre est particulièrement utile dans le cas où la reconstruction de la phase dans une région circulaire donne les renseignements suffisants sur la topographie de l'objet de géométrie circulaire, comme par exemple dans le cas de la visualisation et de la mesure de la face avant d'une fibre optique. L'homme du métier est enfin apte à déterminer ces fonctions propres dans le cas de supports triangulaires ou semi-infinis. Par ailleurs, il est entendu que nous avons décrit un mode de réalisation faisant intervenir des transformées de Fourier discrètes en une dimension, mais que l'homme du métier est apte à mettre en oeuvre des transformées de Bessel discrètes en une dimension, en tant qu'analogue curviligne des transformées de Fourier. L'homme du métier peut également utiliser tout type de transformation intégrale, en tant que projection d'une fonction sur un espace de fonction, par exemple déterminé par le type de fonction utilisé dans les équations 5 ou 5' sans pour autant sortir du cadre de la présente invention. Dans tous les cas, l'homme du métier aura compris que si les fonctions propres de l'équation de Helmholtz s'expriment en des variables séparées, en particulier si le domaine de l'image n'est pas trop distordu, on pourra utiliser la présente invention et n'effectuer des transformations intégrales qu'en une seule dimension. Pour des domaines très distordus, on pourra utiliser des transformations conformes pour ramener ces domaines à des domaines plus simples, ou bien définir ce domaine comme une composition de domaine plus simples. De plus, il est entendu que le mode de réalisation décrit pour une image est également applicable pour tout type d'onde de rayonnement selon l'ensemble du spectre lumineux. Afin d'appliquer l'enseignement de la présente invention aux différents spectres, l'homme du métier est apte à modifier les accessoires de captation tels que décrits dans la partie suivante consacrée aux accessoires associées à l'invention. DISPOSITIFS DE MESURE POUR LA MISE EN UVRE DE L'INVENTION Afin de mettre en oeuvre l'invention, et conformément par exemple à la figure 4, il convient d'utiliser des dispositifs de mesure de l'intensité dans des plans distincts. Selon un mode de réalisation, ces mesures sont réalisées dans des plans perpendiculaires à l'axe optique. 20 Pour ce faire, illustré figure 6, on réalise une prise d'image sur un échantillon transparent. L'homme du métier comprendra que pour effectuer des mesures sur deux plans perpendiculaires, on peut soit modifier la position de l'échantillon par rapport au dispositif optique fixe, soit déplacer le dispositif optique. Sur la figure 6, pour un échantillon transparent, le dispositif selon un mode de réalisation de l'invention comprend donc un capteur d'images 61, par exemple sous la forme d'un capteur CCD, un système optique de projection 62 pour projeter la lumière issue de l'échantillon 63 vers le capteur. La lumière est émise par un système optique d'éclairage 64. L'échantillon est mobile selon l'axe des z grâce à un dispositif de déplacement (non représenté). Illustré figure 7, selon un autre mode de réalisation de l'invention, l'échantillon est fixé à une plateforme immobile 75. Le dispositif optique comprend toujours un capteur 71, un système de projection 72 et un système d'éclairage 73, et il comprend en outre une plateforme 76 déplaçable à laquelle sont fixés les éléments 71, 72 et 73. Cette plateforme est actionnée par un système de déplacement 77 sous la forme par exemple d'un moteur ou d'un élément piézo-électrique. Illustré figure 8, l'invention comprend également un dispositif de mesure utilisant une ou plusieurs lentilles à focale variable. Dans ce cas, on dispose par exemple une lentille à focale variable 88 au niveau du système de projection et selon une variante une lentille à focale variable 89 au niveau du système d'éclairage pour des mesures sur un échantillon transparent fixe. L'homme du métier comprendra que dans ce cas le déplacement des mesures d'intensité correspond à la différence de focale, c'est-à-dire Oz= F2-F1. Dans le cas de deux lentilles, les modifications de la focale se feront simultanément pour les deux lentilles. On comprendra également que ce système est également réalisable avec une seule lentille à focale variable. Illustré figure 9, pour un échantillon diffusant et non plus transparent, on peut, selon un mode de réalisation, utiliser une séparatrice du faisceau 90 destiné à envoyer un rayon lumineux issu du système d'éclairage 94. On peut alors utiliser l'ensemble des éléments précédemment décrits tels que le capteur, le système de projection, l'ensemble étant soit fixe soit mobile, ainsi qu'une ou plusieurs lentilles à focale variable tel que précédemment décrit. Il est entendu que dans tous ces modes de réalisation, le capteur d'image est relié à un dispositif de traitement de l'image comprenant un logiciel pour la mise en oeuvre de l'algorithme selon l'invention. Les images de phase obtenues par ce logiciel peuvent 10 alors être visualisées sur un écran d'ordinateur. L'invention est décrite dans ce qui précède à titre d'exemple. Il est entendu que l'homme du métier est à même de réaliser différentes variantes de l'invention sans pour autant sortir du cadre du brevet. Ternie dans la Description du ternie dans Relation à la notation discrète Relation à la géométrie notation continue la notation continue d'un capteur de type CCD t' = , t et Vecteurs depuis r origine Réseau de points avec les Coordonnées du centre vers les points dans coordonnées (in. n) pour chaque d-un élément sur le -= (x l'espace 2D nombre entier ira et n_ capteur CCD S Région du support de la Réseau de points avec les Tous les éléments actifs phase coordonnées (In. n) où du capteur CCD nl =1 M et n = L.... N _ X = n Ligne de la frontière de la a cc 1 ie du réseau des points La colonne des éléments région de support S de la avec l'indice ta = + passifs du capteur CCD phase. en direction de r. situés la plus à droite de tous les éléments actifs du capteur CCD. y = b Ligne de la frontière de la Ligne du réseau des points avec La ligne des éléments région de support S de la l'indice In =M+ 1. passifs du capteur CCD phase, en direction de x. situés au-dessus de tous les éléments actifs du capteur CCD. C Frontière de la région de Deux lignes du réseau des points Tous les éléments support S de la phase. avec les coordonnées (o. n) , et passifs du capteur (ou sa (Ill + 1.n) où r1= 1 v. frontière) qui entourent ainsi que deux colonnes avec les le domaine rectangulaire coordonnées (m.0) et des pixels actifs. tin.N+1) où in= 1,...,M. Tableau I REVENDICATIONS 1. Procédé de récupération de phase de champ d'ondes de rayonnement comprenant les étapes suivantes: - calculer une mesure représentative de la variation d'intensité dudit champ d'ondes de rayonnement sur une surface sélectionnée s'étendant globalement d'un bout à l'autre du champ d'ondes; - calculer une mesure représentative de l'intensité 10 dudit champ d'ondes de rayonnement sur ladite surface sélectionnée; transformer ladite mesure de la variation d'intensité pour produire une première représentation de transformée intégrale et appliquer à ladite première représentation de transformée intégrale un premier filtre, pour produire une première représentation de transformée intégrale modifiée; appliquer un inverse de ladite première transformée intégrale à ladite première représentation de transformée intégrale modifiée pour produire une représentation non transformée; - appliquer une correction basée sur ladite mesure d'intensité sur ladite surface sélectionnée à ladite représentation non transformée; - transformer ladite représentation non transformée corrigée pour produire une seconde représentation de transformée intégrale et appliquer à ladite seconde représentation de transformée intégrale un second filtre; - appliquer un inverse de ladite seconde transformée intégrale à ladite seconde représentation de transformée intégrale modifiée pour produire une mesure de phase dudit champ d'ondes de rayonnement d'un bout à l'autre dudit plan sélectionné, caractérisé en ce que lesdits premier et second filtres sont des mesures d'une représentation d'une fonction de Green dans un espace de fonctions propres de l'équation de Helmholtz à coordonnées séparables, de sorte de que lesdites première et seconde transformées intégrales sont des transformées intégrales en une dimension de l'espace dudit champ d'ondes de rayonnement. 2. Procédé selon la revendication 1, caractérisé en 10 ce que lesdites première et seconde transformées intégrales sont produites en utilisant une transformée de Fourier. 3. Procédé selon la revendication 2, caractérisé en ce que ladite transformée de Fourier est une transformée de 15 Fourier rapide. 4. Procédé selon la revendication 1, caractérisé en ce que lesdits premier et second filtres sont sensiblement les mêmes. 5. Procédé selon la revendication 1, caractérisé en ce que lesdites surfaces sont sensiblement planes. 6. Programme informatique pour exécuter les étapes 25 selon l'une quelconques des revendications 1 à 5. 7. Programme informatique stocké sur au moins un support de stockage lisible par ordinateur, caractérisé en ce qu'il comprend un moyen pour exécuter les étapes selon l'une des revendications 1 à 5. 8. Appareil pour la récupération de la phase de champ d'ondes de rayonnement comprenant les étapes suivantes: - un moyen d'acquisition d'une mesure représentative de la variation d'intensité dudit champ d'ondes de rayonnement sur une surface sélectionnée s'étendant globalement d'un bout à l'autre du champ d'ondes; - un moyen d'acquisition d'une mesure représentative de l'intensité dudit champ d'ondes de rayonnement sur ladite surface sélectionnée; - un moyen de traitement pour séquentiellement réaliser les étapes suivantes: (i) transformer ladite mesure de la variation d'intensité pour produire une première représentation de transformée intégrale et appliquer à ladite première représentation de transformée intégrale un premier filtre, pour produire une première représentation de transformée intégrale modifiée; (ii) appliquer un inverse de ladite première transformée intégrale à ladite première représentation de transformée intégrale modifiée pour produire une représentation non transformée; (iii) appliquer une correction basée sur ladite mesure d'intensité sur ladite surface sélectionnée à ladite représentation non transformée; (iv) transformer ladite représentation non transformée corrigée pour produire une seconde représentation de transformée intégrale et appliquer à ladite seconde représentation de transformée intégrale un second filtre; (v) appliquer un inverse de ladite seconde transformée intégrale à ladite seconde représentation de transformée intégrale modifiée pour produire une mesure de phase dudit champ d'ondes de rayonnement d'un bout à l'autre dudit plan sélectionné, caractérisé en ce que lesdits premier et second filtres sont des mesures d'une représentation d'une fonction de Green dans un espace de fonctions propres de l'équation de Helmholtz à coordonnées séparables, de sorte de que lesdites premières et seconde transformées intégrales sont des transformées intégrales en une dimension de l'espace dudit champ d'ondes de rayonnement. 9. Appareil selon la revendication 8, caractérisé en ce que ledit moyen d'acquisition d'une mesure représentative de la variation d'intensité dudit champ d'ondes de rayonnement comprend au moins une plateforme mobile. 10. Appareil selon la revendication 8, caractérisé en ce que ledit moyen d'acquisition d'une mesure représentative de la variation d'intensité dudit champ d'ondes de rayonnement comprend au moins une lentille à focale variable.