FR 2922316 A1 20090417 FR 0707246 A 20071016 �Procédé et dispositif de calcul de la phase totale d'un champ d'ondes à partir d'images d'intensité � Domaine technique La présente invention concerne un procédé et un dispositif de calcul de la phase d'un champ d'ondes à partir d'images d'intensité. Un tel dispositif ou procédé permet à un utilisateur de construire une image de la phase, et peut donc permettre de réaliser une image en relief d'un objet. Typiquement, cet objet consiste en un micro objet tel un MEMS (� Micro-Electro-Mechanical Systems �). Un tel dispositif ou procédé peut aussi permettre, à partir de l'image de phase, de réaliser une image représentative de la topographie d'une surface d'où provient le champ d'ondes, par exemple dans le but de mesurer la rugosité de cette surface. Un tel dispositif ou procédé permet typiquement mais de façon non limitative une mesure de rugosité et/ou de formes de surfaces en micro-électronique, en semi-conducteurs, et dans l'industrie du papier ou du plastique. Etat de la technique antérieure On connaît le document WO 2006/082327, qui décrit un procédé et un dispositif de reconstruction d'une image de phase d'un champ d'ondes à partir d'images d'intensité. Selon WO 2006/082327, la reconstruction requiert typiquement une illumination d'un objet, la translation d'une caméra dans une direction verticale (direction z) et, pour plusieurs positions verticales de la caméra, une acquisition d'une image de l'objet sur les éléments sensibles (pixels) de la caméra arrangés dans un plan (x,y) perpendiculaire à la direction verticale z. Cette technologie repose sur la caractérisation des composants énergétiques des ondes optiques pour en reconstruire les composants vectoriels. Selon ce document, la phase est reconstruite par une simple intégration de la variation longitudinale de l'intensité du champ d'ondes le long de la direction z. Cependant, le résultat d'une telle reconstruction de phase n'est pas toujours satisfaisant. La formule utilisée pour la reconstruction n'est pas exacte, et peut conduire à des erreurs de reconstruction en particulier lorsque l'objet imagé présente des reliefs importants. Le but de l'invention est de proposer un dispositif et un procédé de -2 calcul de la phase d'un champ d'ondes plus performant que ceux existant dans l'état de la technique. Exposé de l'invention Cet objectif est atteint avec un procédé de calcul de la phase d'un champ d'ondes, comprenant : une détermination de valeurs de l'intensité du champ d'ondes en plusieurs points d'une surface, une détermination d'une variation, en les points de la surface, de l'intensité du champ d'ondes par rapport à une coordonnée définie le long d'un axe perpendiculaire à la surface, de préférence une dérivée première de l'intensité du champ d'ondes par rapport à la coordonnée perpendiculaire à la surface, caractérisé en ce qu'il comprend en outre : un calcul d'une fonction de singularités représentative d'une absence ou d'une existence d'une singularité de phase pour chaque point de la surface, et - un calcul de la phase du champ d'ondes sur la surface qui dépend de l'existence ou de l'absence de singularités. Par définition, on dit qu'il existe une singularité (ou point singulier) de la phase du champ d'ondes en un point de la surface si la circulation du gradient de la phase dans la surface et autour de ce point n'est pas sensiblement zéro mais est sensiblement égale à un multiple de �2n. Cette circulation est appelée � valeur de la singularité �. La surface est de préférence finie, c'est-à-dire limitée par des bords. Cette surface est de préférence plane. Ainsi, la surface peut typiquement être rectangulaire. Selon l'invention, un champ d'ondes peut contenir au minimum une onde monochromatique. La fonction de singularités peut être obtenue par résolution d'une équation de Fredholm du second ordre. Le calcul de la fonction de singularités peut comprendre une décomposition de la fonction de singularités en une série de fonctions de décomposition. Les fonctions de décomposition sont de préférence des -3 fonctions propres de l'équation de Helmholtz à une dimension, de préférence pour des conditions limites de Neumann ou de Dirichlet dans la surface. La décomposition de la fonction de singularités peut être du type: f(i-)= E fk F-(F) où f(r) est la fonction de singularités, k = { k, l } est un k=] vecteur, signifiant que le vecteur parcours toutes les valeurs de k=Î k = 1 = (k =1,1=1) à k = K = (k = K,l = L ), K et T étant deux nombres entiers, Fk(r") sont les fonctions de décomposition, fi sont des coefficients de décomposition de la fonction de singularités. Le procédé selon l'invention peut comprendre un calcul des coefficients de décomposition par résolution d'un système linéaire du type : K E (Smk +smk)fk =b,~,, m = 1,...,M k=1 où Fil = { m, n } est un vecteur, m = 1,...,111 signifiant que le vecteur m parcours toutes les valeurs de m = 1 = (m =1,n =1) à m = M = (m = M,n = N), M et N étant deux nombres entiers, k = { k, 1 } est un vecteur, K signifiant que le vecteur parcours toutes les valeurs de k=1 k = 1 = (k = l,l =1) à ii= K = (k = K, l = L ), K et L étant deux nombres entiers, S est la surface, la surface étant un rectangle dont les côtés ont pour longueurs a et b , F et Fo sont des cordonnées de points de la surface, sont les coefficients de décomposition de la fonction de singularités, 1,szm= , b,-h= fJsb(T)F,~,(r)dF, Fe,(r) et F(r) et Fk(r) 0, dans les autres cas. sont les fonctions de décomposition, b(F)=$fs2rz-,-18I(io)[VlnIO)xVg(î-,io)]dio, X est la longueur d'onde moyenne du champ d'ondes, I(F) est la valeur de l'intensité au point de coordonnées F, z est une coordonnée le long de l'axe (Z) perpendiculaire à 25 la surface, 8I(Fo) est la dérivée première de l'intensité en Fo par rapport à la fk20 -4 coordonnée le long de l'axe (Z) perpendiculaire à la surface, g(F,i-o) est une fonction de Green utilisée pour le calcul de la fonction de singularités, smk =uk2 w, f f [ VF-.(r).VFk(Y)]Fm(r) a Uk,l ù k 4k =k7i/a, , 2 2 =lrr/b, la matrice w- étant définie par Vlnl(i)=E w-•OFT(r). La fonction de Green utilisée pour le calcul de la fonction de singularités est de préférence une série de fonctions propres de l'équation différentielle de Helmholtz à deux dimensions, de préférence pour des conditions limites de Neumann ou de Dirichlet dans la surface. Le procédé selon l'invention peut comprendre en outre un calcul des positions et des signes des singularités. Le calcul des positions et des signes des singularités peut comprendre : - une recherche sur la surface de régions où la valeur absolue de la fonction de singularités est supérieure à un seuil, et - une détermination des barycentres de ces régions, chaque position des barycentres ainsi déterminée étant les coordonnées d'une singularité. Dans le procédé selon l'invention, le calcul de la phase peut comprendre : - un calcul, à partir des valeurs et de la variation de l'intensité, d'une première estimation de la phase, - un calcul d'un terme correctif de phase qui dépend de l'existence de singularités de phase, et - une combinaison de la première estimation de la phase et du terme correctif de phase. Le calcul de la première estimation de la phase peut comprendre : - un calcul d'un gradient de la phase, - une intégration du gradient de la phase. Le calcul du gradient de la phase peut comprendre un calcul d'un premier terme qui ne dépend pas de l'existence de singularités, ce calcul 30 comprenant pour des conditions respectivement de Neumann ou de Dirichlet: - une transformée respectivement de cosinus ou sinus appliquée sur la -5 variation de l'intensité, puis, - sur le résultat de cette transformée, une application d'un filtre dans le domaine spatial, puis -sur le résultat de ce filtrage, une application d'une transformée respectivement de sinus ou de cosinus, puis une division par les valeurs de l'intensité. Le calcul du gradient de la phase peut comprendre en outre un calcul d'un deuxième terme qui dépend de l'existence de singularités et qui est additionné au premier terme, le calcul du deuxième terme comprenant pour des conditions respectivement de Neumann ou de Dirichlet: une transformée respectivement de cosinus ou sinus appliquée sur la fonction de singularités, puis, sur le résultat de cette transformée, une application d'un filtre dans le domaine spatial, puis - sur le résultat de ce filtrage, une application d'une transformée respectivement de sinus ou de cosinus, puis une division par les valeurs de l'intensité. L'intégration du gradient de la phase peut comprendre, pour des conditions respectivement de Neumann ou de Dirichlet : - une transformée respectivement de sinus ou de cosinus appliquée sur le gradient de la phase, puis, - sur le résultat de cette transformée, une application d'un filtre dans le domaine spatial, puis - sur le résultat de ce filtrage, une application d'une transformée respectivement de cosinus ou de sinus. Le terme correctif de phase peut comprendre une phase harmonique, calculée au moyen d'une série, sur toutes les singularités, de la valeur de la singularité multipliée par une intégrale d'un produit de gradients d'une fonction de Green utilisée pour le calcul de la phase harmonique. La fonction de Green utilisée pour le calcul de la phase harmonique est de préférence une série de fonctions propres de l'équation différentielle de Helmholtz à deux dimensions, de préférence pour des conditions limites de Neumann ou de Dirichlet dans la surface. -6 Le procédé selon l'invention peut comprendre en outre un regroupage des singularités en au moins un couple de singularités de signes opposés, et pour chaque couple une construction sur la surface d'une ligne de coupure partant d'une singularité du couple jusqu'à l'autre singularité du couple. La ligne de coupure d'un couple passe de préférence par une région de faible intensité de la surface. Pour l'au moins un couple, les signes des singularités du couple sont de préférence opposés. Pour l'au moins un couple, la ligne de coupure du couple passe de préférence, entre les singularités du couple, par autant de singularités positives que de singularités négatives. Le terme correctif de phase peut comprendre une phase cachée, et le procédé selon l'invention peut comprendre un calcul de la phase cachée comprenant : - pour chaque couple, le calcul d'une fonction de couple, et - une somme, sur tous les couples de singularités, d'un opérateur de couple appliqué aux fonctions de couple Le calcul de la fonction de couple peut comprendre un calcul d'une série, sur toutes les singularités du couple, de la valeur de la singularité multipliée par: une différence entre une première intégrale, selon une première coordonnée, d'une variation, selon une deuxième coordonnée, d'une fonction de Green utilisée pour le calcul de la phase cachée, et une deuxième intégrale, selon la deuxième coordonnée, d'une variation, selon la première coordonnée, de la fonction de Green utilisée pour le calcul de la phase cachée. La fonction de Green utilisée pour le calcul de la phase cachée est de préférence une série de fonctions propres de l'équation différentielle de Helmholtz à deux dimensions, de préférence pour des conditions limites de Neumann ou de Dirichlet dans la surface. L'application de l'opérateur de couple sur la fonction de couple peut comprendre : - un calcul d'un gradient spatial de la fonction de couple, - une intégration du gradient de la fonction de couple le long d'un chemin d'intégration, en ajoutant ou en supprimant 27r au gradient de la fonction de couple à chaque fois que le chemin d'intégration coupe la ligne de coupure du couple L'invention concerne aussi un dispositif de calcul de la phase d'un -7 champ d'ondes, comprenant des moyens pour mettre en oeuvre un procédé selon l'invention. Pour chacune des étapes du procédé selon l'invention, le dispositif selon l'invention peut comprendre des moyens pour mettre en oeuvre cette étape. Suivant encore un autre aspect de l'invention, il est proposé un programme informatique pour exécuter le procédé selon l'invention. Enfin, suivant encore un autre aspect de l'invention, il est proposé un support de stockage, comme par exemple une mémoire, un disque dur, une disquette, un CD ou une clef USB, stockant un programme informatique pour exécuter le procédé selon l'invention. Description des figures et modes de réalisation D'autres avantages et particularités de l'invention apparaîtront à la lecture de la description détaillée de mises en oeuvre et de modes de réalisation nullement limitatifs, et des dessins annexés suivants : - la figure 1 illustre un contour d'intégration pour un calcul d'une valeur d'une singularité au point r,, - la figure 2 est une vue schématique de profil d'un mode de réalisation préférentiel de dispositif selon l'invention, - la figure 3 est un organigramme d'un mode de réalisation préférentiel de procédé selon l'invention, - la figure 4 est un exemple d'objet imagé par le dispositif de la figure 2, et pour lequel le procédé de la figure 3 est mis en oeuvre, - les figures 5 à 10 sont des représentations de données calculées selon le procédé de la figure 3 par le dispositif de la figure 2 imageant l'objet de la figure 4 : ^ la figure 5 illustre l'intensité I(i), • la figure 6 illustre la fonction f(r), ^ la figure 7 illustre la première estimation de la phase �,0, • la figure 8 illustre la phase cachée x(F), ^ la figure 9 illustre la phase harmonique x(î°), et -8- ^ la figure 10 illustre la phase totale �)=01(r)-x(x,Y)+x(r)• Nous allons tout d'abord poser différentes notations et conventions utilisées dans ce document. L'intensité d'un champ d'ondes à un point de coordonnées F sur une surface S située à une position z le long d'un axe Z sera désignée par I(r,z). La surface S est aussi appelée domaine ou région de support. Une image d'intensité comprend l'ensemble des valeurs d'intensité connues sur la surface S. La phase du champ d'ondes au point F sera désignée par (t(r). On considèrera par la suite le cas particulier mais nullement limitatif où la surface S est plane, rectangulaire, les coordonnées bidimensionnelles r = (x, y) d'un point de S étant exprimées selon les coordonnées cartésiennes x et y contenues dans la surface S le long d'axes respectivement des abscisses X et des ordonnées Y perpendiculaires entre eux, et z étant la position du plan S le long d'un axe de hauteur Z perpendiculaire audit plan. Les coordonnées de certains de ces points peuvent être munies d'indices, comme par exemple Fo = (xo,yo), r, _ (xäy, ), r- _ (x,,y, ). Il est entendu que l'homme du métier est apte à adapter l'enseignement de l'invention pour d'autres types de surface S et/ou d'autres types de coordonnées. En particulier, S peut être plane et circulaire et r et toutes les équations peuvent être exprimées avec des coordonnées polaires (r,9) contenues dans le plan S, ce qui présente un intérêt notamment si on considère des images d'intensité ayant la forme d'un disque. Pour simplifier les notations, nous allons omettre la variable z dans la liste des coordonnées de l'intensité I. Néanmoins, lors du calcul de la dérivée de l'intensité par rapport à la variable z , cette variable sera implicitement incluse dans la liste des coordonnées. De plus, la phase c1 est supposée changer très peu lors de la propagation de la lumière, et par conséquent nous omettrons la variable z de la liste des coordonnées pour la phase cl). -9 L'opérateur gradient V est calculé par rapport à la variable F, c'est à dire en coordonnées cartésiennes dans les directions perpendiculaires à z i.e. V = [ a/ax, a/ay ] T . En particulier, O,o est le gradient par rapport à ro, et D; est le gradient par rapport à F . Un vecteur ou un produit vectoriel entre des parenthèses [ ]Z marquées d'un indice z signifie la composante de ce vecteur ou produit vectoriel le long de l'axe Z. Les transformées intégrales kt(k) cosinus (ICT) et sinus (IST) sont définies pour n'importe quelle fonction K(x) par: Pour la ICT : x (k) = ? Jcos(kx) K(x) dx, 0 Et, pour la IST : �k)= 12 Jsin(kx) K(x) dx, 7r o Pour calculer numériquement les transformées intégrales cosinus et sinus, nous utilisons les transformées discrètes cosinus de type 1 et sinus de type 1. Les transformées discrètes kc[k] cosinus de type 1 (DCT1) et sinus de type 1 (DST1), dans un mode de réalisation, peuvent être définies pour n'importe quelle fonction discrète K[j] par: Pour la DCT1 : n x[k]=E E[j]K[j]cos( j k7r/n ), j=0 pour k = 0,...,n, avec : 0,5 si j=0ou j = n, 1,0 dans les autres cas. (1) (2) (3) (4) Et, pour la DST1 : 292231 6 -10û i[k]=En-2 x[j]sin[ (j+1)(k+1)z/n ], i=o (5) Pour k = 0,...,nû2. Des algorithmes pour réaliser ces transformées sont connues par l'homme du métier et ne seront donc pas précisément décrites. La grande majorité des algorithmes, écrits sous forme de bibliothèques de programmes en C/C++ pour l'implémentation des transformées de sinus et 5 de cosinus est basée sur l'utilisation des Transformées de Fourier qui, à leurs tours, sont implémentées sous forme de l'algorithme de Transformée de Fourier Rapide (FFT) (cf. par exemple Réf. [9] et Réf. [10]). En effet, il suffit de faire une recherche sur Internet pour trouver des algorithmes utilisant de nombreuses versions de l'utilisation de la FFT pour calculer les 10 transformées de cosinus et de sinus. La Transformée de Fourier Rapide couramment employée utilise l'arithmétique des nombres complexes : pour calculer les transformées de cosinus et de sinus d'un vecteur réel de 1V chiffres, la taille de la transformée de Fourier Rapide doit être au moins 2ÎV . Moins connus sont des algorithmes récursifs, basés sur le calcul d'algèbre 15 linéaire, qui effectuent le calcul des transformées de cosinus et de sinus en utilisant l'arithmétique réel : ils ne font pas appel à des nombres complexes. Dans la Réf. [11], un algorithme récursif (c'est à dire utilisant une série de calculs uniformes et consécutifs, dont les données en entrées sont les données en sortie des calculs équivalents précédents, effectués 20 avec le nombre d'éléments dans les données traitées divisées en deux) est présenté, qui permet de calculer la transformée de cosinus ou sinus, sans faire appel à la Transformée de Fourier. Pour augmenter la vitesse du calcul, le procédé selon l'invention utilisera de préférence, pour calculer ces transformées de cosinus et de sinus, un algorithme qui ne fait pas appel aux transformées de Fourier rapides comme par exemple un algorithme récursif ; on parlera alors de transformées cosinus ou sinus � rapides �. Enfin, nous utilisons la convention : in=}m,n }, =}k,l }, (6) i=}1,j}. 2922316 -11 Dans ce document, une expression, formule, équation,... donnée avec certains indices parmi in, Tc, T, m, n, k , l , i , j, ... reste valable ré-exprimée avec d'autres indices. On va tout d'abord décrire quelques équations et algorithmes utilisés 5 par le dispositif et le procédé selon l'invention des figures 1 à 3. Les équations de transport de l'énergie lumineuse en approximation paraxiale donnent un lien entre la vitesse de changement du flux de l'énergie lumineuse dans la direction de propagation du champ d'ondes (direction z ), et la vitesse de changement du flux dans le plan S 10 perpendiculaire à la direction de propagation z . Ce lien est décrit par l'équation entre les trois composantes du vecteur de Poynting [I(Y,z)V (Y),29LX II(r,z)]T V. [ I(i:,z)V a)]=-2rci.-1 al(r'z) az où . est la longueur d'onde moyenne du champ d'ondes pour l'environnement de propagation du champ d'ondes. L'environnement de 15 propagation comprend notamment les milieux de réflexion ou de transmission du champ d'ondes. Selon le théorème de Helmholtz, un champ vectoriel I(F)oc(r) peut être décomposé en une partie potentielle O�(r) qui satisfait l'identité V x VO(F) = et en la partie sans divergence V x A(@) : 1(i)v (F)=vo(F)+VxA(r). (8) 20 De même, le gradient de la phase Ocl)(r) peut être décomposé en une partie potentielle V (F)= o� ) et en une partie rotationnelle I(F) vxB(Y)=vxA(i^) I(r) Vc(r)=Vyr(Y)+oXB(r) (9) La partie rotationnelle V x B(~^) du champ o(D(î^) est notamment due à la présence de singularités de phase. (7) -12 En référence à la figure 1, les singularités (ou points singuliers) de phase sont localisées dans les points isolés r, E R2, autour desquels la circulation co(F) (aussi appelée � valeur de la singularité �) de o(D(F) le long d'un contour C(@1) compris dans S et'avec petit rayon n'est pas zéro mais est égale à un multiple de �2i w(FI)=ck C(@,) est parcouru dans la direction opposée des aiguilles d'une montre, que nous prenons comme la direction positive, autour du point rä et c est le vecteur unitaire tangent au contour C(;) et orienté dans le sens inverse des aiguilles d'une montre. Le nombre F s'appelle l'ordre de la singularité. En optique, l'ordre d'une singularité ne peut prendre que des valeurs entières : F =1,2,3,.... En plus, la probabilité que l'ordre d'une singularité dépasse 1 est proche de zéro, car une singularité de l'ordre F se décompose en F singularités de l'ordre 1 en présence d'une perturbation infiniment petite au cours de la propagation de la lumière. Pour cette raison, en pratique, on considère qu'il y a seulement en général des singularités de l'ordre 1. La lèm, singularité est positive si le signe de w(F1) est positif dans l'équation (10), sinon elle est négative. Il suit du principe de la conservation de l'énergie qu'il y a autant de singularités positives que de singularités négatives, de sorte que la somme des ordres signés de toutes singularités est zéro. Pour une image de phase (I)(F) bidimensionnelle (coordonnées r ) décomposée en pixels, on peut dire qu'il existe une singularité au niveau d'un pixel de singularité s'il existe, entre ce pixel de singularité et un pixel adjacent à ce pixel de singularité, une différence de la valeur de la phase c(1^) sensiblement égale à un multiple de r radians, c'est à dire une différence de hauteur le long de l'axe Z sensiblement égale à un multiple de - dans le cas où le champ d'ondes a été transmis par l'objet imagé, ou 2 une différence de hauteur le long de l'axe Z sensiblement égale à un -13 multiple de 4 dans le cas où le champ d'ondes a été réfléchi par l'objet imagé. Cela est typiquement le cas si le dispositif selon l'invention image un objet présentant des reliefs importants avec de grandes pentes. En pratique, on remarque que l'intensité est en général sensiblement égale à zéro dans un point singulier : I(r,)=0, r, étant les coordonnées d'un point singulier. (11) Selon le théorème de Stokes, le rotor de vil(r) est : vxvcI*.)=vxvx130=lco(Y,)11,(Y,)Ô(YùY, ), (12) !=0 où 3(FùY,) est la fonction de Dirac et n,(r,) est un vecteur unitaire perpendiculaire au plan S et orienté dans la direction positive de l'axe Z, avec pour origine les coordonnées r, de la lème singularité, co(F,,) étant la valeur de la singularité au point r, . Dans la situation pratique où seulement des échantillons de la phase sur le plan S= (x, y) sont connus pour des pixels de taille Ax le long de l'axe X et Dy le long de l'axe Y, une valeur d'une singularité en r, est définie à l'aide de l'intégration, sur le contour C, des composantes du gradient de la phase autour du point r, en utilisant les différences finies, comme représenté sur la figure 1 : vcD.di °c(x,, )+0(D(x1+Ax,y,)zW(x,~+4y)0I(x,,y,) y y Ofi(x,,yr) _ [ fi(x, +Ax,y,)ùcD(x,,y,) Umm, (13) = [ 'v(x,,y, +Ay) û(�x1UAy Qy . Nous allons montrer comment construire une image de la phase cb(r) d'un champ d'ondes à partir de la composante longitudinale du vecteur de Poynting du champ d'onde (c'est à dire de l'intensité fois le nombre d'onde 2iri -') et de la variation de cette composante du vecteur de Poynting le long d'une direction de propagation z du champ d'ondes, en tenant compte des singularités de phase. -14 Pour le potentiel vectoriel Â(F) correspondant à la région de support S, il existe une fonction de Green g(Y',ro ) qui reflète la géométrie de S et les conditions limites sur le potentiel selon l'équation : A(Y)ù ffsLvro xvro xA(YO)]g(Y~YO)~o • (14) D'après l'équation de Poisson, la partie potentielle vyp(F) du gradient de la phase v(D(r) peut donc se réécrire sous la forme : v r 2 ) = I(r) 1 Ifs [ al(o Y ,z)laz ] vg(F, ) dr0 . (15) w (~ o De même, la partie rotationnelle V x B(r) du gradient de la phase vcl)(F) peut se réécrire sous la forme : vxB(r)=ùJ ) ffs [vJ(YO)xvcp(F0)]xvg(Y,F0)dro (16) Le produit vectoriel VI(F)xvc(F) est un vecteur orienté le long de l'axe Z, et le rotor v x Â(r) a ses composantes dans le plan S. Soit la fonction G(F,ro) liée à la fonction de Green g(i^,r''o) par les conditions de Cauchy- Riemann : aG ag ax ay ' (17) DG ag ay ax La partie rotationnelle V x B(r) du gradient de la phase vc(F) peut donc se réécrire sous la forme : vxB(Y)=ùI(@) f fÇ [VI(YO,z)xvq)(YO)]ZvG(Y,YO)dFo. On a donc : vcp(F) = ) ' Ya)vg (Y, Fo ) d@ I(;. ffs CSI( o 1 (18) fJ, [VI(Fo)xv(D(Y0) ] VG(Y,F0)dFo, 1(~) s où 8I(Fo)=al(F0,z)/az. Cette équation représente la solution générale au -15 problème de calcul du champ vectoriel \'D à partir des mesures de la composante longitudinale du vecteur de Poynting (i.e. intensité) et de sa variation le long de la direction axiale z. Le produit vectoriel de cette équation par \I(F) donne l'équation intégrale de Fredholm de second ordre suivante: f(r)=b(r)û JJS f(F )V(F,r0)ur0 , (19) où . f(j.)={VI(Y)xvcD(F)L 1 ~-b(r)= JJ,S 27L/~, 1(SIfro)[V1n1(;)xVgfr Fo)] d V(F,r'o)= [ V 1n1(F)xVG(F,ro) 17, In étant la fonction logarithme népérien. Il est nécessaire d'obtenir la fonction f(r) pour : l'obtention du gradient oc(F) de la phase totale selon la formule (18) ; ce gradient comportera tant la partie potentielle que la partie rotationnelle ; le calcul de la phase totale. Le calcul de la phase totale s'effectue en intégrant le gradient v(D(F) de la phase totale. L'intégration de la partie potentielle est directe. Par contre, l'intégration de la partie rotationnelle nécessite le placement de lignes de coupure entre singularités de signes opposés. Les positions et les signes des singularités sont déterminés à l'aide de la fonction f(r). Selon l'invention, on peut montrer à partir d'une intégration de la formule (18) que la phase (D (F ) est la somme d'une partie potentielle tg(F) et d'une partie rotationnelle aussi appelée phase cachée x(r) 'c�r)= (r)+x(r), (21) la phase cachée x(F) étant la fonction liée par les conditions de Cauchy-Riemann avec la fonction W(F)= w(r',)g(F,r,), c'est-à-dire en coordonnées 1 cartésiennes: (20) -16 ax`'J=1 w(ri)gy(r,rl), ax ax(r)--E w(rl)gx(F,rl), a 1 signifiant une somme sur 1 avec 1 variant de 0 à L-1, L étant le nombre de singularités, gx(F,Fo) étant la dérivée de g(F,ro) par rapport à la composante x du vecteur r = (x,y), et gy(F,ro) étant la dérivée de g(F,i-o) par rapport à la composante y du vecteur r = (x,y). La partie potentielle tg(F)de la phase est égale à une première estimation de la phase 0,0 à laquelle est soustraite une phase harmonique x(x,y): W(r)=01(Y)û '(x,Y) Avec : (F )= f J (Fo)•~TbW.,r0)ui0 X(x,y)= JJS[\ x (F)) J.vrog(r,r0)dio =ùE rw(rl) J Js [ V ro g(ro , rr) X V ro g(r , F0) J di-0 En référence à la figure 1, le dispositif 1 selon l'invention mettant en oeuvre un procédé selon l'invention comprend un système optique 2, un capteur 3 tel une caméra CCD, une unité de calcul 4 telle une unité centrale d'un ordinateur, un écran 5, et des moyens (non représentés) pour illuminer un objet 6. Les moyens d'illumination peuvent être agencés pour illuminer l'objet 6 en transmission ou en réflexion, de sorte que le système optique 2 soit agencé pour projeter sur le capteur 3 le champ d'ondes lumineuses 7 respectivement transmis ou réfléchi par l'objet 6 et se propageant parallèlement à l'axe Z. L'objet imagé 6 se trouve sensiblement dans le plan focal objet 8 du système optique 2. L'axe optique du système optique est aligné avec l'axe Z. Le dispositif 1 comprend en outre des moyens (non représentés) pour modifier la position relative du capteur 3 le long de l'axe Z par rapport au système (22) (23) (24) - 17 optique 2 solidaire de l'objet 6. Le dispositif peut donc comprendre soit des moyens pour déplacer le capteur 3, soit des moyens pour déplacer le système optique 2 solidairement avec l'objet 6. Le capteur 3 comprend une grille bi-dimensionnelle de pixels rectangulaires arrangés dans un plan perpendiculaire à l'axe Z et alignés respectivement selon M lignes parallèles à l'axe Y et N colonnes parallèles à l'axe X, les axes X et Y étant perpendiculaires entre eux et à l'axe Z. Le capteur est donc agencé pour acquérir, pour plusieurs positions possibles du capteur le long de l'axe Z, une image d'intensité du champ d'ondes, chaque image d'intensité étant délimitée par une surface plane perpendiculaire à l'axe Z, et comprenant M fois N points ou pixels. En pratique, ces M fois N pixels peuvent n'être qu'un sous ensemble du nombre total de pixels du capteur. La région S peut être choisie par l'utilisateur, et donc elle peut contenir un nombre M fois N de pixels moindre que le nombre total de pixels que compte le capteur. Par exemple, le capteur 3 contient 1280 x 1024 pixels. Or, l'utilisateur peut choisir une zone de calcul S de taille arbitraire, par exemple 268 x 873 pixels, sans toutefois dépasser la surface utile du capteur, c'est à dire sans aller au-delà de la taille de 1280 x 1024 pixels. Le capteur 3 est relié à l'unité de calcul. L'unité de calcul 4 stocke un logiciel selon l'invention et est agencée pour calculer différentes données 1(F), aI(F ,z)/az, f(r), v11)(F), O,(Y), x(F), x(x,y), (D(r) pour chacun de ces points. L'écran 5 permet d'afficher ces données sous la forme de graphiques. On va maintenant décrire, en référence aux figures 2 à 10, différentes étapes du procédé de calcul de la phase selon l'invention. Ces étapes successives A à L comprennent : A. Une acquisition, par le capteur 3, d'au moins un couple d'images d'intensité 1,(r)=I(F,z,) et I20=I(F,z2) du champ d'ondes provenant de l'objet 6 (par réflexion ou par transmission), les images du couple étant acquises pour des positions 3a, 3b du capteur 3 symétriques par rapport au plan focal image 9 du système optique 2, z, et z, étant, le long d'un axe Z, les positions respectives de la surface plane imagée respectivement S, -18 et S2 pour chacune des images, x et y étant des coordonnées le long d'axes respectivement X et Y perpendiculaires entre eux et à l'axe Z . S, et S2 sont des surfaces planes perpendiculaires à l'axe Z. Toutes les images acquises ont une même échelle spatiale et des mêmes dimensions le long 5 des axes X et Y, ce qui nécessite une mise à cette échelle des données brutes acquises par le capteur 3. L'unité de calcul 4 procède alors au calcul de la valeur de l'intensité I(FF) du champ d'ondes pour chaque point F de la surface plane S perpendiculaire à l'axe Z et située à une position z le long de l'axe Z, et au 10 calcul de la dérivée première SI(r )= al(F ,z)/az de l'intensité sur la surface S par rapport à la coordonnée z pour chaque point r de la surface plane S. S a la même échelle spatiale et les mêmes dimensions le long des axes X et Y que S, et S2. Pour le calcul des valeurs de l'intensité I(F), l'intensité I(î) au point 15 de coordonnées F dans S est typiquement égale à une moyenne des valeurs au point de coordonnées r dans S, et S2 des différentes images d'intensité acquises respectivement 1,(F) et 12 (F) . Par la suite, on considèrera le cas où : I(F)= [ I,(r)+ I2(F) 1/2, (25) z=(z, +z2)/2 Mais on aurait pu choisir n'importe quelle valeur pour 1(F) qui est une 20 moyenne plus ou moins pondérée des images d'intensité acquises : I(F)=1,(7) ou I(r)=I2(i-) ou I(F)=[I,(i°)+212(r) 1i3 etc... (26) S étant comprise entre S, et S2, z étant une moyenne de z, et z2 selon cette même pondération. On considère par la suite le cas nullement limitatif pour lequel la surface imagée S est un rectangle délimité par quatre bords, chacun des 25 bords de ce rectangle étant aligné avec un des axes des coordonnées x et y utilisés pour les calculs. On a S = (0, a)x (O,b), c'est-à-dire que deux bords de S alignés le long de l'axe X ont une longueur a, et deux bords de S 292231 6 -19 alignés le long de l'axe Y ont une longueur b. Tout comme 11(F) et 12(F), l'image I(F) comprend M fois N points. La distance entre deux points de I(F) le long de l'axe X est Ox, et la distance entre deux points de I(F) le long de l'axe Y est Ay . Un point de S imagé par un pixel situé à la ième ligne (0 b(r)= fi 27rÀ. 'SI(ro)[VlnI(F)xVrg(F,ro)] dFo V(F,Fo)= [ V lnI(r)xV7.G(F,F0) ] 1, étant la longueur d'onde moyenne du champ d'ondes, qui est en pratique quasi-monochromatique centré sur une longueur d'onde ). avec une bande passante de quelques nanomètres. g(F,Fo) est la fonction de Green qui reflète la géométrie de la région de support S et les conditions limites sur le potentiel. La fonction de Green g(F,Fo) peut être définie comme l'expression sous forme d'une série de fonctions propres (F,,(r)Fm(Fo)) de l'équation différentielle de Helmholtz à deux dimensions, de préférence pour les conditions limites de Neumann ou de Dirichlet dans S : g(r,ro)=E v 2 F~n(r)Fm(ro), (32) ,n=0 u étant les valeurs propres, avec : = m7c/a 7~n = n ir/b Um2 4m2 +fn2 ù et où /h= 0 . oo représente tous les indices requis pour définir une fonction propre requise, et les fonctions de décomposition Fm(F) (ou FF(F), F,-(F)... selon l'indice choisi) forment l'ensemble ortho-normal (voir réf.[1] Morse & Feshbach, page 820), et sont des fonctions propres de l'équation différentielle de Helmholtz à une dimension, de préférence pour les conditions limites de Neumann ou de Dirichlet dans S ff F�i )Fk(F)=. (33) Pour des conditions aux limites de Neumann, on fixe, le long des bords de S, la valeur d'une composante du gradient spatial de g(F,Fo), cette composante étant contenue dans le plan S et normale au bord. Typiquement, on fixe cette valeur égale à 0, en particulier si les images acquises sont des images d'un objet 6 comprenant une surface sensiblement plane imagée notamment par les bords des images, et un (31) - 21 û objet en relief centré sur les images. Pour les conditions aux limites de Neumann, on prend : g(F,T"o) = Aän cos(riny)cos(rinyo) n=0 x COS(mX) COS( m x0 ) L ~m +nnn 2 En pratique, ce calcul sur l'infini d'une ou de deux sommes s'effectue en utilisant une de ces deux approches : 1. En remplaçant l'infini (m = oo et n = oo) par des nombres entiers m= M et n = N suffisamment grands (typiquement 256, cette valeur dépendant de la puissance de calcul de l'unité de calcul 4) pour permettre la précision requise de la représentation de la fonction de Green. Aussi, pour permettre le calcul de la fonction de Green en utilisant des transformées rapides (par exemple, les transformées de Fourier, de Sinus, de Cosinus, etc.) ; les nombres M et N sont de préférence égaux à des puissances de 2, M étant une puissance de 2 supérieure ou égale à m et la plus proche possible de M, N étant une puissance de 2 supérieure ou égale à N et la plus proche possible de N; 2. En calculant une des deux sommes analytiquement (c'est-à-dire en exprimant cette somme sous la forme d'une fonction analytique), et en remplaçant l'infini (oo) de la somme restante (par rapport à m ou n) par un nombre entier suffisamment grand pour permettre la précision requise de la représentation de la fonction de Green. Dans l'implémentation de l'algorithme qui suit, la deuxième de ces approches est utilisée, où la deuxième somme sur m dans l'équation (34) est calculée analytiquement, tandis que l'infini (oo) de la première somme sur n de l'équation (34) est remplacé par N qui est une puissance de 2 supérieure ou égale à N et la plus proche possible de N. On peut utiliser des approches de calcul comparables pour un calcul d'autres sommes dans le procédé selon l'invention. Pour des conditions aux limites de Dirichlet, on fixe, le long de (34) m=0, m2+n2SO - 22 û chaque bord de S, la valeur de g(r,Fo). Typiquement, on fixe cette valeur égale à 0, en particulier si les images acquises sont des images d'un objet comprenant une surface sensiblement plane imagée notamment par les bords des images, et un objet en relief centré sur les images...DTD: Pour les conditions aux limites de Dirichlet, on prend : l l g(;, ro) = E Am ä sin(rin y) sin(rin y0 ) n=0 x m=0, m`+n2 m0 sin(mx) sin(mxo ) m2 + n2 (35) Typiquement, la deuxième somme sur m dans l'équation (35) est calculée analytiquement, tandis que l'infini (oo) de la première somme sur n de l'équation (35) est remplacé par N qui est une puissance de 2 supérieure ou égale à N et la plus proche possible de N. 10 Les paramètres 4,n, n,,, Am,n seront détaillés par la suite. Le calcul de f(F) est basé sur une décomposition de f(F) en une série de fonctions de décomposition FF(r) : .f(F)=E fi Fi (F.). (36) k=0 fi (ou fm , J. ... selon l'indice choisi) étant les coefficients de décomposition de la fonction f(F). Dans la formule (36), l'infini (oo) est 15 remplacé par K = (K,L), K et L étant les valeurs suffisamment élevées (et pas nécessairement des puissances de deux, et typiquement comprises entre 30 et 80 selon la puissance de l'unité 4) pour représenter la fonction f(F) avec une précision requise. Dans le procédé selon l'invention, on procède typiquement de la même manière pour le calcul d'autres sommes 20 allant jusqu'à l'infini. La fonction b(r) est décomposée en une série des fonctions de décomposition F,,,(i°) : CO b(F)=E b,, Fm(r)• (37) m=0 b,,, (ou bk , b;... selon l'indice choisi) étant les coefficients de - 23 û décomposition de la fonction b(F). Dans la formule (37) l'infini (x) est remplacé par M = (M,N), M et N étant les valeurs suffisamment élevées (pas nécessairement des puissances de deux, et typiquement comprises entre 30 et 80 selon la puissance de l'unité 4) pour représenter la fonction b(r) avec une précision requise. Dans le procédé selon l'invention, on procède typiquement de la même manière pour le calcul d'autres sommes allant jusqu'à l'infini. Selon l'invention, les fonctions f(F) et b(F) sont donc décomposées en séries des fonctions de décomposition F,;,(F) qui sont des fonctions propres de l'équation de Helmholtz à une dimension, de préférence pour les conditions limites de Neumann ou de Dirichlet dans le domaine S. Pour les conditions aux limites de Neumann, on décompose f(r) et b(F) en utilisant la Transformée de Cosinus Discrète (DCT), c'est-à-dire en prenant : Fm (@) = Fm n (x, y) = A,,,,n cos( mx ) cos(r~ny ) (38) Pour les conditions aux limites de Dirichlet, on décompose f(r) et b(F) en utilisant la Transformée de Sinus Discrète (DST), c'est-à-dire en prenant : Fm (r) = Fm,n (x, y) = Am,n sin( mx ) sin( 77 n.@ ) , (39) Le calcul de f(F) comprend les étapes B.1. à B.6. suivantes : B.1. Un calcul des paramètres v, 1 , 4m, Tb, , Am et Bn selon l'équation (40) et sm selon l'équation (41) : k = k z/a =lir/b `gym = m7r/a (40) i7, =nn/b 2 2 2 uk,! = k + ni -24- Am n = A,,, B,, ? si m>O a 1 a B,, =~ si n>0 Bo 1 b si m=0 �m = si m>0 (41) B.2.Un calcul de la fonction b(r)= f fs27rl '81(ro)[V1nI(F)xVFg(r,Fo)]Z dFo Pour cela, l'unité de calcul utilise le calcul : - de la fonction scalaire 810 [ 12(F)-11(F) J/Az, où Az = z2 - z, est la distance entre les deux plans (S, et S2) de prise de vue de deux images d'intensité 11(r) et I2(r) ; et de la fonction vectorielle Oln1(r), où I(rr)= [ 1,(r)+12(x)1/2 à partir des images d'intensité I,(r) et 12(r). Pour calculer la fonction b(F), sont calculées la dérivée première gx(r,ro) de g(r,ro) par rapport à la composante selon x du vecteur r =(x,y), et la dérivée première gy(r,ro) de g(r,ro) par rapport à la composante selon y du vecteur r" = (x, y) : b(r) = f fs 22r) 'SI(ro)[VinI(r)x\rg(r,ro)]Z dio = ôl(r)/ I(rff 2zl '8I(iô)g(r,ro)(lii) (42) ,~-_ a1(r)/ ff 27r%l -'8I(ro)gX(r,ro) Fo I(r) s Pour les conditions aux limites de Neumann, les dérivés de la fonction de Green par rapport aux variables x et y sont : -25- t gx(ri0)=ù2a-1E sin(S,,,x)cos(4,nxo)J (y,Yo) (43) m=1 t gy(i',r0)=-2bù1 sin(r)ny)cos(nny0)/n (x,x0). (44) n=1 Ces équations nous montrent que fm et f' sont des mesures d'une représentation d'une fonction de Green g(F,î^o) dans un espace de fonctions propres de l'équation de Helmholtz à coordonnées séparables. Pour les conditions aux limites de Dirichlet, les dérivés de la fonction de Green sont celles données par les deux formules ci-dessus dans laquelle la fonction sinus est remplacée par la fonction cosinus et inversement. C'est-à-dire, pour les conditions aux limites de Neumann : b(x,y)= fis -13I(xo,y0)[01nI(x,Y)xV g(x,y;xo,YO)]_ dx0dy0 1 1 al (x' y) /ax `° e sin(riny) f f. (x,xo) f 31(xo,Yo)cos(rinyo)dyo dxo (45) +4zc~ 1 a-1 aI(x'Y)/ây ° sin x 8I x , yo)cos(,nxo) dxo dyo I (x, y) "7, (~,n ) f o _0 fm (y~ Yo) f (o Pour les conditions aux limites de Dirichlet, il suffit de réécrire l'équation précédente en remplaçant la fonction cosinus par la fonction sinus et inversement. Pour les conditions aux limites de Neumann : f �x, xo) û_ 1 cosh[rin (a û x)] cosh(rinxo ), si x > xo sinh(yn7r) cosh[n,, (a û xo )j cosh(rinx), si x 0 .Î,n (Y'Yo) = sinh(ÿm7r) cosh[ n, (b û yo )] cosh( n,y), si y x0 ' Sinh(Y/!TC) Sinh[r)n (a ù xo )] S1n11(nnx), si x Yo Î,~ (Y~ yo) = sinh()7rnir) sinh[{,,, (b ù yo )] sinh(4n, y), si y -28 Pour les conditions aux limites de Dirichlet, il suffit de reprendre ces équations en remplaçant les sinus par des cosinus et en replaçant les cosinus par des sinus. Les transformées cosinus et sinus, de préférence � rapides �, permettent de retrouver les valeurs wm et wy . Pour les conditions aux limites de Neumann : a h wm,n = ù ml Jsin( mx) [L(xY) cos(ni,y) dy dx 0 0 Mù1 / A Nù1 ( 7t \ 1yl sin miù E Lx(x;,y1)cos nj =o \. M/ i=o \ N/ m =1,..., n =1,..., N , (54) Mù1 .7r - a w,,;, = û Jsin(riny) JLy,(x,Y) cos(mx) dx dy o _o = AyE sin/nj N\ E Ly (x,, yJ)co- s/mi J j=o / ,=o M m =1,..., M, n =1,..., N . Pour les conditions aux limites de Dirichlet, il suffit de reprendre ces équations en remplaçant les sinus par des cosinus et en replaçant les cosinus par des sinus. B.4.2. Un smk = uk-2 wi JJS [vF (r)• vFk (F) ] F,;, (r) di' . Pour les conditions aux limites de Neumann, les coefficients sont : (55) des coefficients calcul -29 s;nk = v - 2 W- J Js [ OF: (Y) . OFk (Y) ]F,, i (r) di -2 Ak+m Ak B,+n B, x Y l = 4 vk,l A B ( k+m4kWk+m,l+n +fl1+n1JiWk+m,l+n )EmEn m n Ak+m Ak B1-n B, K x Y l + (k+m kWk+m,l-n +11l-nnlWil+m,I-n)Em�n Am Bk l Ak+m Ak Bn 1 Bl ( x Y J + k+m kWk+m,n-1 -77n-l~lWk+m,n-I Em Am B,, Ak_mAk BI+n BI x Y + ((k-mkWk-m,1+n +17l+n'llwk-m,l+n)EmEn Am Bn + Ak_m Ak BI_, B1 Am Bn Ak-m Ak Bn l Bl + Am B,, Am-k Ak B,+n B, ((x k-m4kWk-m,l-n +nl-nnlWy k-m,l-n EmEn ((k-m kWkx-m,n-1 ûrIn-111iWy k-m,n-1 Em x y ((m-k kWm-k,l+n nl+nfllWm-k,l+n En (In-k kWmx -k,l-n _11l-nnlWmy -k,l-n En y ( m-k kWm-k,n-1 +nn-lnlWr-k,n-I (56) Am Bn Am_kAk BI-n B, Am Bn Am_kAk B,_,B1 Am Bn Dans la somme (56), les composants du vecteur ï = (i, j) changent comme suit : i =1,...,, j =1,...,N, les paramètres uk2,, 4m, rin, Am et Bn et Em étant décrits précédemment. Pour les conditions aux limites de Dirichlet, les coefficients sont : Smk 4Uk 2 Wi if [v).VFk(Y) ]F, (F) a 4 Uk 2 Wl,j { 4i k ('8ù1-k m + Si+k,m )(8j-1,n gj+l,n) (57) + 11j71 (Si-k,m ù (5i+k,m ) (Sj-1,n + (-j+I,n) } où les composants du vecteur î = (i, j) changent comme suit : i =1,...,M, j =1,...,N , et où: 27r Cl 1ù(ù1)m-(i-k) +1ù(ù1)m+(i-k) mù(iùk) m+(iùk ) = Am A Ak B.5. Un calcul des coefficients fk par résolution du système d'équations : E (swk +Sn~k)Jk =m=1,...,M (59) k=1 (58) -30 où m =1,...,M signifie que le vecteur m parcours toutes les valeurs de m =1 = (m =1,n-1) à m = M =(m = M,n =IV), où k = 1,...,K signifie que le vecteur parcours toutes les valeurs de = 1 = (k =1,1=1) à k=K=(k=K,l=L),etoù: 1sim=k, Smk 0 dans les autres cas. Ainsi, le système d'équations linéaires (59) contient MN équations linéaires pour trouver les K L inconnus qui sont les Ki coefficients fi de la décomposition de f(r) en une série de fonctions Fi(F) selon la formule (36). Ce système peut être résolu au moyen d'un algorithme classique de 10 résolution des systèmes d'équations linéaires, telle une décomposition LU (� Low Upper Decomposition�) ou une décomposition SVD (� Singular Value Decomposition �). La référence bibliographique [9] donne une liste préférentielle de méthodes de résolution de système d'équations linéaires. B.6. Un calcul de la fonction f(F) à partir des coefficients fi obtenus à 15 l'issue de l'étape précédente, par l'utilisation de l'équation : ÎC f(r)=l fk Fk(F). (61) T Pour reconstruire les valeurs des échantillons d'une image f(x, y) sur la grille discrète {x,,y1 }, ou x; =iAx, i=0,...,Mû1 et y. = jAy, j = 0,...,N -1, pour les conditions limites de Neumann, la somme (61) s'effectue au moyen de la transformée cosinus, de préférence cosinus 20 � rapide � : K \ ( îZ\ f(x,,y1)=1 coskiù Ak,/fk,lcos l~ k=1 mi N, i = 0,...,Mû1, j = 0,...,Nû1. Pour les conditions limites de Dirichlet, la somme (61) s'effectue au moyen de la transformée sinus, de préférence sinus � rapide � : (60) (62) -31 K / ' L 7 7rf(x;,y;)=E sin kiù EAk,rfk,,sin lj- k=1 M./ l=1 N ~ i = 0,...,M -1, j = 0,...,N -1. C. Une détermination, par l'unité 4, de la présence de singularités dans la phase comprenant les étapes C.1. à C.3. suivantes : C.1. Une formulation d'un seuil fs ; C.2. S'il existent des points r dans la région de support S où f(F) > fs, des singularités sont présentes dans la phase ; C.3. Si f(F)l px,(xo,yo) la composante selon l'axe X du gradient de la phase VO(xo,y0) en r0 et (Dyo(x0,y0) la composante selon l'axe Y du gradient de la phase V'1(x0,y0) en ro. Pour les conditions aux limites de Dirichlet, il suffit de réécrire l'équation précédente en remplaçant la fonction cosinus par la 5 fonction sinus et inversement. A partir des données cx (xi, y;) et y (xi, y,) obtenues par l'équation (67), la première estimation de la phase est en pratique calculée, pour les conditions aux limites de Neumann, selon l'équation : Pour les conditions aux limites de Dirichlet, il suffit de réécrire 10 l'équation précédente en remplaçant la fonction cosinus par la fonction sinus et inversement. Ces équations montrent que, pour calculer la première estimation de la phase, on effectue les opérations suivantes pour les conditions aux limites respectivement de Neumann ou de Dirichlet : 15 - une transformée respectivement de sinus ou de cosinus appliquée sur le gradient de la phase, puis, - sur le résultat de cette transformée, une application d'un filtre f,n (y,yo) ou ffy(x,x0)dans le domaine spatial, puis - sur le résultat de ce filtrage, une application d'une transformée 20 respectivement de cosinus ou de sinus. L'homme de l'art reconnaîtra que l'étape D correspond au procédé selon WO 2006/082327. U n=1 4~j l Mù1 ( ù \ Nù1 y 01 \xi,yj ) sin mùi fm yj yjo M m=1 M / Jo=1 47i%, -1 Ax Nù1 ( \ Mù1 + E sin nùj Efy(xi,xio) N n=1 \ N / io=1 V CI) xo (xi° , Yi COS m ù i0 0=1 M N-1 ( 7 E ~vo (xi0,yjo )cos nù j0 jo=1 N / - 34 û E. S'il n'y a pas de singularités, la phase totale est définie par l'unité 4 comme étant égale à la première estimation de la phase MF) calculée à l'étape D. F. S'il y a des singularités, la première estimation de la phase est calculée par l'unité 4 selon les étapes F.1. et F.2. suivantes : F.1. Un calcul du gradient de la phase, en utilisant la fonction 8I(F)- [ 12(F)û1,(r) ]/Oz et la fonction f(F), obtenue à l'étape B, selon l'équation : vo(r 27r ' SI r r r )d (71) I(r) $$ f(ro) VG(F,ro)d . Pour calculer les composantes (DX(x,, yj) selon l'axe X et c y(x,, yj) selon l'axe Y du gradient V (x, y), sur la grille discrète 1 x;, yj }e S, ce calcul prend la forme suivante pour les conditions aux limites de Neumann: -47m. ay MM1 / \ N1 (DX(x,,yj)= E sin mùi EfY(yi yj0) Ix, , Yi m=1 \ M / ia=1E SI(xio,yjo) cos 10=1( ir m ù 10 M _ Nù1 ( Mù1 ?~/N 1 sin n _ J EfY(x1,xj I x;,Y; n=1 \ N io=1 Nù1 E f(xb,yjo)cos n Njo _la (72) (Dy(Xi,Y Nù1 Nù1 E sin n- j1 f fY(xl,x,o) Y 6I(x,0,yio)cos nNJoJ J a lo n=1 \ N i =1 =1 ù 4e -'Ax/N I(x;,yi) û 2Ay/M M -1 I(Xi , yj) m=1 N-i 7r sin ~m M i~ fm (Yi ~Yio) E f(xio Yio ) cos ~m m io m-i lo =1 b =1 Pour les conditions aux limites de Dirichlet, il suffit de réécrire l'équation précédente en remplaçant la fonction cosinus par la fonction sinus et inversement. Ces équations montrent que, pour calculer le gradient de la phase, on effectue les opérations suivantes: - un calcul d'un premier terme qui ne dépend pas de l'existence de singularités, ce calcul comprenant pour les conditions aux limites respectivement de Neumann ou de Dirichlet : - 35 une transformée respectivement de cosinus ou sinus appliquée sur la dérivée de l'intensité, puis, ^ sur le résultat de cette transformée, une application d'un filtre fm (y,yo) ou fy(x,xo) dans le domaine spatial, puis ^ sur le résultat de ce filtrage, une application d'une transformée respectivement de sinus ou de cosinus, puis une division par les valeurs I(r) de l'intensité. - un calcul d'un deuxième terme qui dépend de l'existence de singularités, le calcul du deuxième terme comprenant pour les conditions aux limites respectivement de Neumann ou de Dirichlet ^ une transformée respectivement de cosinus ou sinus appliquée sur la fonction de singularités, puis, ^ sur le résultat de cette transformée, une application d'un filtre fm (y, yo) ou f r (x,xo) dans le domaine spatial, puis ^ sur le résultat de ce filtrage, une application d'une transformée respectivement de sinus ou de cosinus, puis une division par les valeurs I(F) de l'intensité. - une addition du premier terme avec le deuxième terme. 20 F.2. Un calcul de la première estimation de la phase çl(r) par intégration sur la surface S du gradient de le phase \(D(F) selon l'équation : 01 (F.)- f JS VcD(r0).\7g(F,F0)~0 . (73) Ce calcul de ql(F) est réalisé au moyen des formules données à 25 l'étape D.2., avec toutefois des valeurs du gradient de la phase qui sont calculées par les équations de l'étape F.1. et qui diffèrent donc des valeurs du gradient de la phase de l'étape D.I. G. S'il y a des singularités, un calcul par l'unité 4 des positions et des signes des singularités selon les étapes G.I. à G.3. suivantes : 30 G.I. Une recherche sur la surface S, c'est à dire le plan (x, y), des -36ùrégions où f(F) fs ; G.2. Une détermination des barycentres de ces régions. Les centres des régions ainsi déterminés sont les coordonnées des singularités, r, , G.3. Une détermination des signes des f(F,). Ces signes sont les signes des singularités présentes aux points F, _ (x,, y,) ; Les singularités sont localisées dans les points où f(F) n'est pas zéro et supérieur en valeur absolue au seuil fS : la position du maximum local ou du minimum local, ou la position du centre de gravité de f(F) peut servir comme l'indication de la position d'une singularité de phase. La singularité est positive [typiquement co(r,)=27r], si f(r,)>0, et négative [typiquement co(;) _ -2Tr ], si f (F,) L est le nombre de singularités identifiées dans S. H. S'il y a des singularités, un calcul par l'unité 4 de la phase harmonique ex, y) en utilisant les positions et les signes des singularités calculés à l'étape G, au moyen de la somme sur toutes les singularités : x(x,Y)=ùE w(F,) JJs °'og(Fo,Fr)XO.,g(F,Fo)] dFo Le calcul de la phase harmonique x(F) en un point F comprend donc le calcul d'une série, sur toutes les singularités, de la valeur de la singularité multipliée par une intégrale sur la surface S d'un produit de gradients de la fonction de Green. On remarque que la transformée intégrale et le calcul de x(x, y) dépend du nombre et des signes des singularités, mais ne dépend pas des valeurs d'intensité I (F) et de la dérivée de l'intensité BI(F )= aI(F ,z)/az Le calcul de cette somme est réalisé, dans le cas où la fonction de Green g(F, Fo) est une série de fonctions propres (F,b(FF) Fm (FO )) de l'équation différentielle de Helmholtz à deux dimensions, par l'équation suivante : Lù1 x(x,Y)= w(Fi)Tr [wY(Y1,x1)wY(x,Y)-wY(x1,y1)wY(Y,x)], (74) r=o -37 où Tr est la trace d'une matrice, co(r,)=�27r est la valeur de la singularité dans un point r', =(x/,y/), L est le nombre total des singularités identifiés dans l'image de phase : L/2 singularités positives et L/2 singularités négatives. Les composantes W,nr n (v,u) des matrices W" (v,u) sont, pour les conditions aux limites de Neumann : 1 W', n (v,u) = Am n cos(rinv) [ ((m +in )sinh(nyir) 1-I x{~n si.nh(ny7r)sin( mu)+4m [COSh(U_fly)_(_1)m cosh(îu) I }, (75) m =1,...,M -1, n =1,...,Nû1, m étant l'indice de ligne de la matrice WY (v,u), n étant l'indice de colonne de la matriceW'(v,u), v et u pouvant prendre pour valeur par exemple x,, y,, x, ou y. En remplaçant dans cette équation y par )7, on obtient l'expression des composantes des matrices W n (v, u) . On a y = a/b , ÿ = b/a , M = a/Ax , N = b/Ay, = m7r/a , Ax et Ay sont les dimensions d'un point dans l'image. De même, pour les conditions aux limites de Dirichlet : Wm n (v,u)= Am n sin(rinv) [ (4m +11,2)sinh(ny7r) x {Tin sinh(ny7r) cos(mu)+ m [ sinh(rinu û nyir)û (û1)m sinh(rinu) J }, (76) m =1,...,Mû1, n =1,...,Nû1, En pratique, dans les conditions aux limites de Dirichlet, la phase 15 harmonique X(r) est posée comme étant sensiblement nulle. I. S'il y a des singularités, un regroupage des singularités en couples de singularités par l'unité 4, les singularités d'un couple ayant des valeurs de signes opposés, selon les étapes I.1. à I.7. suivantes : 20 I.1. Déterminer les régions continues de faible intensité (c'est-à-dire pour lesquelles IO Le calcul de la phase cachée en un point r de la surface comprend donc : pour chaque couple d'indice c, le calcul d'une fonction de couple _x y Jgy (u, y ; rr des singularités, mais ne dépend pas des valeurs d'intensité I (F) et de la dérivée de l'intensité 8I(@ )= aI(@ ,z)/az . Dans le cas où la fonction de Green g(@,@0) est une série de fonctions propres (F,, (@) F,, (Fo� de l'équation différentielle de Helmholtz à deux dimensions, x(FF) est calculée selon l'équation : (77) -41 L1~1 çà, co\Yi,,; / ex, xi c=0 i=0 m>0 n>0 l +b-'yH(xùx,.,)ùa lx H(Yùy1c,,)}) Dans cette équation, i n'a donc que deux valeurs possibles : 0 ou 1. L est le nombre total des singularités identifiés dans l'image de phase : typiquement L/2 singularités positives et L/2 singularités négatives d'ordre 1, H(x û x, ) est une fonction marche (aussi appelée fonction de Heaviside) 5 centrée sur = 0, et H(y û y, ;) est une fonction marche (aussi appelée fonction de Heaviside) centrée sur .yû y, = 0 . Dans le cas des conditions aux limites de Neumann : FäY (y, Yi,,; ) = û [ sinh(mÿir) -1 x [H(yûy, )sih(mÿ7c û my) cosh((myrc.., ) û H(y,; û y)cosh(m)77r m)4,, )sinhmy) ], î' (x,x, )=û[sinh(nyr)] x [ H(xûx, ; )sinh(nyir ûrinx)cosh(rinx, ) û H(x, ûx)cosh(ny7r ûr)nxl ; )sinh(rhix) ]. Q' (Y, Y/, )= 2 (n7r) -1 sin(r1ny) cos(rinyr ), (80) n où H(x û x, ), H(x, ûx), sont des fonctions de Heaviside centrées sur 10 = 0, H(y - y, ;), et H(y, ; ûy) sont des fonctions de Heaviside centrées sur y û y, = 0 . Pour les conditions aux limites de Dirichlet : Fm (Y, Y,, , ) = û [ sinh(m@7r) x { H(yûY, )cosh(mÿrc ûmY)sinh(mY,., ) û H(y, û y)sinh(mgir )cosh(my) ], Fy (x' ) _ û [ sinh(nyir) ] -' x [ H(x )cosh(nyir ûr)nx) sinh(ri,,x,, ) û H(x, ; ûx)sinh(nyir ûnnx, )cosh(nnx) ], ) F y (Y,Y, , ) û Q n (Y,Y, ) Fr (x' x, ,; ) (78) (79) Q>n(x,xl,; )= 2 (mir)-'sin( mx)cos( mx, )' (81) 2922316 -42 Qm x, )= 2 (m7r) 1 COS(mx) SjI1(mxl r ), Qn (y, y, )= 2 (na)-1 cos(riny) sin(riny1,), où H(x ù x, ), H(x, ; ù x), sont des fonctions de Heaviside centrées sur x ù x, , =0, Myù y,), et H(y,, ù y) sont des fonctions de Heaviside centrées sur yù y, = 0 . Dans ces équations, l'opérateur S2�( w(F,r, o,)) représente une opération, appliquée à une fonction w(r,r, J, de placement des chemins de coupure entre les deux singularités appartenant au couple d'indice c qui se situent dans les points r, o et i, . Cette opération contient les étapes suivantes : 1. Calculer le gradient spatial modulo 27r de la fonction de couple w(r,r, p,r, ). L'opération modulo 27r signifie que chaque valeur des composantes du gradient de la fonction w(7,,) (une dérivée première de w(F,r, ,@, ) par rapport à x et une dérivée première de w(@,@, 0,F1 ) par rapport à y) est ramenée à l'intervalle [û7r,7r] par l'addition ou la soustraction de 27r ; 2. Intégrer le long d'une ligne, en partant d'un seul et même point de la surface S pour tous les couples indicés par c sur lesquels l'opérateur S2,( w(r,r,,@, � est réalisé (de préférence, à partir de l'origine des coordonnées (x = 0,y = o)), en allant jusqu'au point i- , et selon des chemins d'intégrations librement choisis dans la surface S (comprenant de préférence un segment parallèle à l'axe X et un segment parallèle à l'axe Y), le champ de gradient obtenu à l'étape 1, en lui ajoutant ou en lui supprimant 27r à chaque fois que le chemin d'intégration coupe le chemin de coupure du couple d'indice c selon la règle suivante : a) 27r est ajouté si l'orientation Ô,,j du chemin de coupure dans un point Rc, du chemin de coupure traversé et l'orientation du chemin d'intégration font un angle compris dans l'intervalle (0,7r) ; (82) -43 b) 27r est soustrait si l'orientation Ô, du chemin de coupure dans un point k ,j du chemin de coupure traversé et l'orientation du chemin d'intégration font un angle compris dans l'intervalle (0,û7r). K. S'il y a des singularités, un calcul par l'unité 4 de la phase totale c(F") : du résultat ç,(@) de l'étape F, est soustraite la phase harmonique X(x,y) obtenue à l'étape H et est ajoutée la phase cachée x(F) obtenue à l'étape J : c(F)=oi(F)+x(F)ûx(x,y) L. Un stockage de la phase totale dans le mémoire de l'unité 4 : la phase totale est égale au résultat de l'étape E s'il n'y a pas des singularités, sinon la phase totale est donnée à l'issue de l'étape K. L'image de la phase totale est alors visualisée sur l'écran 5. La figure 4 est un exemple d'objet 6 imagé par le dispositif de la figure 2, et pour lequel le procédé de la figure 3 est mis en oeuvre. L'objet comprend un plan 10 perpendiculaire à l'axe Z et parallèle aux axes X et Y, au milieu duquel est disposé un plateau rectangulaire 11 muni de quatre côtés, dont un premier côté 21 parallèle avec l'axe Y et situé en continuité avec la surface 10, un deuxième côté 22 parallèle avec l'axe Y et opposé au premier côté 21, et un troisième et un quatrième côtés 23, 24 qui s'éloignent progressivement selon une certaine pente à partir du plan 10 et du côté 21 vers le côté 22. L'axe Z est exprimé en radians, un radian correspondant à une hauteur le long de l'axe Z de D (27r). La longueur des projections des côtés 21 à 24 sur le plan 10 est égale à la largeur de 20 pixels de 5,2 micromètres de côtés. La pente est telle que la phase change de 0 radians à 1,257r radians au cours de ces 20 pixels . Les figures 5 à 10 sont des images de données calculées selon le procédé de la figure 3 par le dispositif de la figure 2 imageant l'objet de la figure 4. Ces images sont représentées dans le plan S comprenant les axes X et Y. La taille d'un pixel d'une image est de 5,2 pm x 5,2 pm respectivement le long de X et Y. Pour les images 5, 7, 8, et 10, plus un point de l'image est foncé, et plus la valeur de la donnée représentée est -44 faible en ce point. La figure 5 illustre l'intensité I(F) calculée à l'étape A. On remarque que, à l'emplacement des bords continus entre eux 22, 23 24 pour lesquels il existe une discontinuité de surface, l'intensité forme une ligne d'encrage sombre 12. La figure 6 illustre la fonction f(r) calculée à l'étape B. On remarque une région 25 pour laquelle f(F) > fs et f(,)>0. Il existe donc une singularité positive 26 à la position du centre de gravité de f(r) de la région 25. De même, remarque une région 27 pour laquelle f(r) > fs et f ) La figure 7 illustre la première estimation de la phase �,(F) calculée à l'étape F. La figure 8 illustre la phase cachée x(F) calculée à l'étape J. La figure 9 illustre la phase harmonique x(r) calculée à l'étape H. Sur cette figure, on distingue trois zones 30, 31, 32 séparées par des lignes 20 parallèles à l'axe Y, et correspondant à des valeurs de X(F) de plus en plus grande dans la direction des x croissants. Enfin, la figure 10 illustre la phase totale '(F) = ç, (r)- x(x, y)+ x(F) calculée à l'étape K. On remarque que le résultat de la phase totale c1(F) sur la figure 10 est plus proche de la réalité de l'objet imagé de la figure 4 25 que le résultat de la première estimation de la phase 01(F) sur la figure 7. L'homme du métier comprendra donc que la prise en compte des singularités est un avantage évident d'un procédé et d'un dispositif de calcul de phase selon l'invention par rapport à l'état de la technique, en particulier par rapport à ce qui est décrit dans le document WO 2006/082327. 30 Bien sûr, l'invention n'est pas limitée aux exemples qui viennent -45 d'être décrits et de nombreux aménagements peuvent être apportés à ces exemples sans sortir du cadre de l'invention. En particulier les images d'intensité acquises pourraient être circulaires. Dans ce cas, le procédé selon l'invention est réalisé comme décrit précédemment, mais en utilisant toutes les formules et équations venant d'être décrites exprimées en coordonnées polaires, les coordonnées polaires étant de préférence centrées sur le centre des images circulaires acquises. De plus, quand il est dit qu'une grandeur est calculée selon une équation ou formule donnée, cette grandeur peut plus généralement être 10 calculée selon l'invention sensiblement par cette formule. REFERENCES : [1] P.M.Morse, H.Feshbach, Methods of Theoretical Physics, McGraw-Hill, New-York, 1953. 15 [2] E.G. Abramochkin, V.G. Volostnikov, � Relationship between twodimensional intensity and phase in a Fresnel diffraction zone, � Optics Communications, vol. 74, no. 3,4, pages 144-148, Décembre 1989. [3] V. P. Aksenov, V.A. Banakh, O.V. Tikhomirova, � Potential and vortex parts of optical speckle fields, � Atmospheric and Oceanic Optics, vol. 9, 20 no. 11, pages 921-925, Novembre 1996. [4] V. Aksenov, V. Banakh, O. Tikhomirova, � Potential and vortex features of optical speckle fields and visualization of wave-front singularities, � Applied Optics, vol. 37, no. 21, pages 4536-4540, July 1998. [5] Nugent et al, WO 00 26 622. 25 [6] I. Lyuboshenko, � Obtention d'une image de phase à partir d'une image d'intensité �, demande de brevet WO 2006/082327. [7] A. P. Prudnikov, Y. A. Brychkov, I. O. Maichev, Integrals and Series, Gordon and Breach, New York, 1986. -46 [8] I. S. Gradshteyn and I. M. Ryzhik, Table of Integrals, Series and Products, Academic, New York, 1980. [9] W. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery, Numerical Recipes in C, Cambridge University Press, 1992. [10] www.fftw.orq [11] G. Baszenski, M. Tasche, � Fast polynomial multiplication and convolutions related to the discrete cosine transform �, Linear Algebra and its Applications, vol. 252, pp 1-25 (1997). REVENDICATIONS 1. Procédé de calcul de la phase (c(@)) d'un champ d'ondes, comprenant : -une détermination (A) de valeurs de l'intensité (1(F)) du champ d'ondes en plusieurs points (F) d'une surface (S), - une détermination (A) d'une variation, en les points (F) de la surface, de l'intensité du champ d'ondes par rapport à une coordonnée (z) définie le long d'un axe (Z) perpendiculaire à la surface, de préférence une dérivée première de l'intensité du champ d'ondes par rapport à la coordonnée (z) perpendiculaire à la surface, caractérisé en ce qu'il comprend en outre : - un calcul (B) d'une fonction de singularités (f(F)) représentative d'une absence ou d'une existence d'une singularité de phase pour chaque point (F) de la surface (S), et - un calcul (D à K) de la phase (�i)) du champ d'ondes sur la surface (S) qui dépend de l'existence ou de l'absence de singularités. 2. Procédé selon la revendication 1, caractérisé en ce que la fonction de singularités est obtenue par résolution d'une équation de Fredholm du second ordre. 3. Procédé selon la revendication 1 ou 2, caractérisé en ce que le calcul (B) de la fonction de singularités comprend une décomposition de la fonction de singularités (f(F)) en une série de fonctions de décomposition. 4. Procédé selon la revendication 3, caractérisé en ce que les fonctions de décomposition sont des fonctions propres de l'équation différentielle de Helmholtz à une dimension, de préférence pour des conditions limites de Neumann ou de Dirichlet dans la surface (S).-48 5. Procédé selon la revendication 3 ou 4, caractérisé en ce que la décomposition de la fonction de singularités est du type: f(F)=l fi Fk(FF) où f(F) est la fonction de singularités, k = { k, 1 } est un k=1 K vecteur, E signifiant que le vecteur parcours toutes les valeurs de k=1 k =1 = (k =1,1=1) à k = K = (k = K,l = L), K et L étant deux nombres entiers, Fk (i-) sont les fonctions de décomposition, fi sont des coefficients de décomposition de la fonction de singularités. 6. Procédé selon la revendication 5, caractérisé en ce qu'il comprend un calcul des coefficients de décomposition par résolution d'un système linéaire du type : E m- +sm-)fk =b- & m = 1,...,M k=1 où m = { m, n } est un vecteur, m = 1,...,M signifiant que le vecteur m parcours toutes les valeurs de m = 1 = (m =1,n =1) à m = M = (m = M,n = N), M et N étant deux nombres entiers, k = { k, l } est un vecteur, E signifiant que le vecteur parcours toutes les valeurs de k=Î ii= 1 =(k=1,1=1) à k = K = (k = K,l = L), K et L étant deux nombres entiers, S est la surface, la surface étant un rectangle dont les côtés ont pour longueurs a et h, F et Fo sont des cordonnées de points de la surface, fi sont les coefficients de décomposition de la fonction de singularités, 1sim=k, , b- = JJ s (F) F- r ,,,(F) , Fm (F) ,,,(F) et F r ,(F) et Fi-(î) b r dr sont 0 dans les autres cas. les fonctions de décomposition, b(F) = IL 27r~ -'SI(ro) [ V lnl(F)x ~g(r, ro) dio , est la longueur d'onde moyenne du champ d'ondes, I(F) est la valeur de l'intensité au point de coordonnées F, z est une coordonnée le long de l'axe (Z) perpendiculaire à Smk- 49 la surface, 31(ro) est la dérivée première de l'intensité en Fo par rapport à la coordonnée le long de l'axe (Z) perpendiculaire à la surface, g(F,FO) est une fonction de Green utilisée pour le calcul de la fonction de singularités, smk =uk z w-f [vF-(F)•\7] (F) ]F;.,(F)di' , ET sVk,! - k n12 ' k =kIr/a, n, =lzr/b, la matrice w- étant définie par V1n1(F)=E w- •VF (r). 7. Procédé selon la revendication 6, caractérisé en ce que la fonction de Green utilisée pour le calcul de la fonction de singularités est une série de fonctions propres de l'équation différentielle de Helmholtz à deux dimensions, de préférence pour des conditions limites de Neumann ou de Dirichlet dans la surface (S). 8. Procédé selon l'une quelconque des revendications 1 à 7, caractérisé en ce qu'il comprend en outre un calcul (G) des positions et des signes des singularités. 9. Procédé selon la revendication 8, caractérisé en ce que le calcul des positions des singularités comprend : - une recherche sur la surface (S) de régions où la valeur absolue de la fonction de singularités est supérieure à un seuil, et - une détermination des barycentres de ces régions, chaque position des barycentres ainsi déterminée étant les coordonnées (r,) d'une singularité. 10. Procédé selon l'une quelconque des revendications 1 à 9, caractérisé en ce que le calcul de la phase (q(Y)) comprend : - un calcul (D, F), à partir des valeurs et de la variation de l'intensité, d'une première estimation de la phase (ç,(F)), - un calcul (H, I, J) d'un terme correctif de phase qui dépend de l'existence de singularités de phase, et une combinaison de la première estimation de la phase et du terme correctif de phase. 292231 6 - 50 11. Procédé selon la revendication 10, caractérisé en ce que le calcul (D, F) de la première estimation de la phase (O,(Î")) comprend : un calcul d'un gradient de la phase, 5 une intégration du gradient de la phase. 12. Procédé selon la revendication 11, caractérisé en ce que le calcul (D, F) du gradient de la phase comprend un calcul d'un premier terme qui ne dépend pas de l'existence de singularités, ce calcul 10 comprenant pour des conditions respectivement de Neumann ou de Dirichlet: une transformée respectivement de cosinus ou sinus appliquée sur la variation de l'intensité, puis, sur le résultat de cette transformée, une application d'un filtre dans le 15 domaine spatial, puis sur le résultat de ce filtrage, une application d'une transformée respectivement de sinus ou de cosinus, puis une division par les valeurs de l'intensité. 20 13. Procédé selon la revendication 12, caractérisé en ce que le calcul (F) du gradient de la phase comprend un calcul d'un deuxième terme qui dépend de l'existence de singularités et qui est additionné au premier terme, le calcul du deuxième terme comprenant pour des conditions respectivement de Neumann ou de Dirichlet: 25 une transformée respectivement de cosinus ou sinus appliquée sur la fonction de singularités, puis, sur le résultat de cette transformée, une application d'un filtre dans le domaine spatial, puis - sur le résultat de ce filtrage, une application d'une transformée 30 respectivement de sinus ou de cosinus, puis une division par les valeurs de l'intensité. 2922316 -51 14. Procédé selon l'une quelconque des revendications 11 à 13, caractérisé en ce que l'intégration du gradient de la phase comprend, pour des conditions respectivement de Neumann ou de Dirichlet : une transformée respectivement de sinus ou de cosinus appliquée sur 5 le gradient de la phase, puis, sur le résultat de cette transformée, une application d'un filtre dans le domaine spatial, puis sur le résultat de ce filtrage, une application d'une transformée respectivement de cosinus ou de sinus. 10 15. Procédé selon l'une quelconque des revendications 10 à 14, caractérisé en ce que le terme correctif de phase comprend une phase harmonique (x(x, y) ), calculée au moyen d'une série, sur toutes les singularités, de la valeur de la singularité multipliée par 15 une intégrale d'un produit de gradients d'une fonction de Green utilisée pour le calcul de la phase harmonique. 16. Procédé selon la revendication 15, caractérisé en ce que la fonction de Green utilisée pour le calcul de la phase harmonique est une série 20 de fonctions propres de l'équation différentielle de Helmholtz à deux dimensions, de préférence pour des conditions limites de Neumann ou de Dirichlet dans la surface (S). 17. Procédé selon l'une quelconque des revendications 10 à 16, 25 caractérisé en ce qu'il comprend un regroupage des singularités en au moins un couple de singularités de signes opposés, et pour chaque couple une construction sur la surface (S) d'une ligne de coupure partant d'une singularité du couple jusqu'à l'autre singularité du couple. 30 18. Procédé selon la revendication 17, caractérisé en ce que la ligne de coupure d'un couple passe par une région de faible intensité de la surface.- 52 û 19. Procédé selon la revendication 17 ou 18, caractérisé en ce que pour l'au moins un couple, les signes des singularités du couple sont opposés, et en ce que la ligne de coupure du couple passe, entre les singularités du couple, par autant de singularités positives que de singularités négatives. 20. Procédé selon l'une quelconque des revendications 17 à 19, caractérisé en ce que le terme correctif de phase comprend une phase cachée (x(F)), et en ce qu'il comprend un calcul de la phase cachée comprenant : pour chaque couple, le calcul d'une fonction de couple, et une somme, sur tous les couples de singularités, d'un opérateur de couple appliqué aux fonctions de couple 21. Procédé selon la revendication 20, caractérisé en ce que le calcul de la fonction de couple comprend un calcul d'une série, sur toutes les singularités du couple, de la valeur de la singularité multipliée par: une différence entre une première intégrale, selon une première coordonnée, d'une variation, selon une deuxième coordonnée, d'une fonction de Green utilisée pour le calcul de la phase cachée, et une deuxième intégrale, selon la deuxième coordonnée, d'une variation, selon la première coordonnée, de la fonction de Green utilisée pour le calcul de la phase cachée. 22. Procédé selon la revendication 21, caractérisé en ce que la fonction de Green utilisée pour le calcul de la phase cachée est une série de fonctions propres de l'équation différentielle de Helmholtz à deux dimensions, de préférence pour des conditions limites de Neumann ou de Dirichlet dans la surface (S). 23. Procédé selon l'une quelconque des revendications 20 à 22, caractérisé en ce que l'application de l'opérateur de couple sur la fonction de couple comprend :-53 - un calcul d'un gradient spatial de la fonction de couple, une intégration du gradient de la fonction de couple le long d'un chemin d'intégration, en ajoutant ou en supprimant 27r au gradient de la fonction de couple à chaque fois que le chemin d'intégration coupe la ligne de coupure du couple 24. Dispositif de calcul de la phase (cD(F)) d'un champ d'ondes, caractérisé en ce qu'il comprend des moyens pour mettre en oeuvre les étapes du procédé selon l'une quelconque des revendications 1 à 23. 25. Support de stockage, stockant un programme informatique pour exécuter le procédé selon l'une quelconque des revendications 1 à 23.15