FR2922316A1 - Wave field phase calculating method for measuring roughness or forms of e.g. microelectronics surface, involves calculating phase of wavefield on surface according to absence or existence of singularity - Google Patents

Abstract

The method involves determining values of intensity of wave field in points of a rectangular surface. Variance of the intensity of the wave field with respect to a coordinate is determined in the points of the surface, where the coordinate is defined along an axis perpendicular to the surface. A singularity function representative of the absence or the existence of a phase singularity is calculated. The singularity function is obtained by a second order Fredholm equation. A phase of the wavefield on the surface is calculated according to the absence or the existence of the singularity. Independent claims are also included for the following: (1) a device for calculating a phase of wave field (2) a computer programmable storage medium comprising a set of instructions for performing a wave field phase calculating method.

Description

 Method and device for calculating the total phase of a wave field from intensity images

TECHNICAL FIELD The present invention relates to a method and a device for calculating the phase of a wave field from intensity images. Such a device or method allows a user to build an image of the phase, and can thus make it possible to make a raised image of an object. Typically, this object consists of a micro object such as a MEMS (Micro-Electro-Mechanical Systems). Such a device or method can also make it possible, from the phase image, to produce an image representative of the topography of a surface from which the wave field originates, for example in order to measure the roughness of the surface. this surface. Such a device or method typically allows but not limited to a measurement of roughness and / or surface shapes in microelectronics, semiconductors, and in the paper or plastic industry. STATE OF THE PRIOR ART Document WO 2006/082327 describes a method and a device for reconstructing a phase image of a wave field from intensity images. According to WO 2006/082327, the reconstruction typically requires an illumination of an object, the translation of a camera in a vertical direction (z direction) and, for several vertical positions of the camera, an acquisition of an image of the camera. object on the sensitive elements (pixels) of the camera arranged in a plane (x, y) perpendicular to the vertical direction z. This technology is based on the characterization of the energy components of optical waves to reconstruct the vector components. According to this document, the phase is reconstructed by a simple integration of the longitudinal variation of the intensity of the wave field along the z direction. However, the result of such a phase reconstruction is not always satisfactory. The formula used for the reconstruction is not exact, and can lead to reconstruction errors, especially when the imaged object has important reliefs.

The aim of the invention is to propose a device and a method for calculating the phase of a wave field that is more efficient than those existing in the state of the art.

DISCLOSURE OF THE INVENTION This objective is achieved with a method of calculating the phase of a wave field, comprising: a determination of values of the intensity of the wave field at several points of a surface, a determination a variation, at the points of the surface, of the intensity of the wave field with respect to a defined coordinate along an axis perpendicular to the surface, preferably a first derivative of the intensity of the field of view; waveform relative to the coordinate perpendicular to the surface, characterized in that it further comprises: a calculation of a singularity function representative of an absence or existence of a phase singularity for each point of the surface, and - a calculation of the phase of the wave field on the surface which depends on the existence or the absence of singularities.

By definition, it is said that there is a singularity (or singular point) of the phase of the wave field at a point on the surface if the circulation of the gradient of the phase in the surface and around this point is not substantially zero but is substantially equal to a multiple of 2n. This circulation is called the value of singularity.

The surface is preferably finished, i.e. limited by edges. This surface is preferably flat. Thus, the surface can typically be rectangular. According to the invention, a wave field may contain at least one monochromatic wave.

The singularity function can be obtained by solving a second order Fredholm equation. The calculation of the singularity function may include a decomposition of the singularity function into a series of decomposition functions. The decomposition functions are preferably eigenfunctions of the one-dimensional Helmholtz equation, preferably for Neumann or Dirichlet boundary conditions in the surface. The decomposition of the function of singularities can be of the type: f (i -) = E fk F- (F) where f (r) is the function of singularities, k = {k, l} is a k =]

vector, meaning that the vector traverses all the values of k = k = 1 = (k = 1.1 = 1) at k = K = (k = K, l = L), K and T being two integers, Fk (r ") are the functions of decomposition, fi are coefficients of decomposition of the function of singularities.

The method according to the invention may comprise a calculation of the decomposition coefficients by solving a linear system of the type:

KE (Smk + smk) fk = b, ~ ,, m = 1, ..., M k = 1 where Fil = {m, n} is a vector, m = 1, ..., 111 meaning that the vector m runs all the values of m = 1 = (m = 1, n = 1) to m = M = (m = M, n = N), M and N being two integers, k = {k, 1} is a vector, K meaning that the vector runs all the values of k = 1 k = 1 = (k = l, l = 1) to ii = K = (k = K, l = L), K and L being two numbers integers, S is the surface, the surface being a rectangle whose sides have for lengths a and b, F and Fo are coordinates of points of the surface, are the coefficients of decomposition of the function of singularities, 1, szm =, b, -h = fJsb (T) F, ~, (r) dF, Fe, (r) and F (r) and Fk (r) 0, in other cases. are the decomposition functions, b (F) = $ fs2rz -, - 18I (io) [VlnIO) xVg (i-, io)] dio, X is the average wavelength of the wave field, I (F ) is the value of the intensity at the coordinate point F, z is a coordinate along the axis (Z) perpendicular to the surface, 8I (Fo) is the first derivative of the intensity in Fo with respect to the fk20 -4 coordinated along the axis (Z) perpendicular to the surface, g (F, io) is a function of Green used for the calculation of the function of singularities, smk = uk2 w, ff [VF-. (r) .VFk (Y)] Fm (r) a Uk, l ù k 4k = k7i / a,, 2 2 = lrr / b, the matrix w- being defined by Vlnl (i) = E w- • OFT (r).

The Green function used for calculating the singularity function is preferably a series of eigenfunctions of the two-dimensional Helmholtz differential equation, preferably for Neumann or Dirichlet boundary conditions in the surface. The method according to the invention may further comprise a calculation of the positions and signs of the singularities. The calculation of the positions and signs of the singularities can comprise: a search on the surface of regions where the absolute value of the singularities function is greater than a threshold, and a determination of the barycentres of these regions, each position of the barycentres as well as determined being the coordinates of a singularity. In the method according to the invention, the calculation of the phase can comprise: - a calculation, from the values and the variation of the intensity, of a first estimate of the phase, - a computation of a corrective term phase dependent on the existence of phase singularities, and - a combination of the first estimate of the phase and the phase correction term.

The calculation of the first estimate of the phase can comprise: a calculation of a gradient of the phase; an integration of the gradient of the phase. The calculation of the gradient of the phase may comprise a calculation of a first term which does not depend on the existence of singularities, this calculation comprising for Neumann or Dirichlet conditions respectively: a transform respectively of cosine or sinus applied on the -5 variation of the intensity, then, - on the result of this transform, an application of a filter in the spatial domain, then - on the result of this filtering, an application of a transform respectively of sinus or of cosine, then a division by the values of the intensity. The calculation of the gradient of the phase may furthermore comprise a calculation of a second term which depends on the existence of singularities and which is added to the first term, the calculation of the second term comprising for Neumann or Dirichlet conditions respectively: a respectively cosine or sine transform applied to the singularities function, then, on the result of this transform, an application of a filter in the spatial domain, then - on the result of this filtering, an application of a transform respectively of sine or cosine, then a division by the intensity values. The integration of the gradient of the phase can comprise, for Neumann or Dirichlet conditions respectively: a respectively sine or cosine transform applied on the gradient of the phase, then, on the result of this transform, an application of a filter in the spatial domain, then - on the result of this filtering, an application of a respectively cosine or sinus transform. The term phase correction may comprise a harmonic phase, calculated by means of a series, on all the singularities, of the value of the singularity multiplied by an integral of a gradient product of a Green function used for the calculation. of the harmonic phase. The Green function used to calculate the harmonic phase is preferably a series of eigenfunctions of the two-dimensional Helmholtz differential equation, preferably for Neumann or Dirichlet boundary conditions in the surface. The method according to the invention may furthermore comprise a grouping of the singularities into at least one pair of singularities of opposite signs, and for each pair a construction on the surface of a cutoff line starting from a singularity of the couple up to to the other singularity of the couple. The cutoff line of a couple preferably passes through a region of low intensity of the surface. For the at least one pair, the signs of the singularities of the pair are preferably opposite. For the at least one pair, the torque cut-off line preferably passes between the singularities of the pair by as many positive singularities as negative singularities. The term phase correction may comprise a hidden phase, and the method according to the invention may comprise a calculation of the hidden phase comprising: for each pair, the calculation of a torque function, and a sum, on all the pairs of singularities, of a torque operator applied to the functions of torque The calculation of the torque function can comprise a calculation of a series, on all the singularities of the pair, of the value of the singularity multiplied by: a difference between a first integral, according to a first coordinate, of a variation, according to a second coordinate, of a function of Green used for the calculation of the hidden phase, and a second integral, according to the second coordinate, of a variation, according to the first coordinate of the Green function used for the calculation of the hidden phase. The Green function used for the calculation of the hidden phase is preferably a series of eigenfunctions of the two-dimensional Helmholtz differential equation, preferably for Neumann or Dirichlet boundary conditions in the surface. The application of the torque operator to the torque function can comprise: - a calculation of a spatial gradient of the torque function, - an integration of the gradient of the torque function along an integration path by adding or removing 27r at the gradient of the torque function each time the integration path intersects the torque cutoff line The invention also relates to a device for calculating the phase of a -7 field of view. wave, comprising means for implementing a method according to the invention. For each of the steps of the method according to the invention, the device according to the invention may comprise means for implementing this step.

According to yet another aspect of the invention, there is provided a computer program for performing the method according to the invention. Finally, according to yet another aspect of the invention, there is provided a storage medium, such as a memory, a hard disk, a floppy disk, a CD or a USB key, storing a computer program for executing the method according to the invention. 'invention. DESCRIPTION OF THE FIGURES AND EMBODIMENTS Other advantages and particularities of the invention will appear on reading the detailed description of implementations and non-limiting embodiments, and the following appended drawings: FIG. 1 illustrates an outline method for calculating a value of a singularity at the point r ,, - Figure 2 is a schematic profile view of a preferred embodiment of the device according to the invention, - Figure 3 is a flowchart. of a preferred embodiment of the method according to the invention, - Figure 4 is an example of an object imaged by the device of Figure 2, and for which the method of Figure 3 is implemented, - the figures 5 to 10 are representations of data calculated according to the method of FIG. 3 by the device of FIG. 2 illustrating the object of FIG. 4; FIG. 5 illustrates the intensity I (i), FIG. the function f (r), FIG. 7 illustrates the first estimate of the phase 0, FIG. 8 illustrates the hidden phase x (F), FIG. 9 illustrates the harmonic phase x (1 °), and FIG. the total phase) = 01 (r) -x (x, Y) + x (r) • We will first make different notations and conventions used in this document.

The intensity of a wave field at a coordinate point F on a surface S at a z position along a Z axis will be designated I (r, z). The surface S is also called domain or support region. An intensity image includes all the intensity values known on the surface S. The phase of the wave field at point F will be designated by (t (r), after which we will consider the particular but not limiting case. where the surface S is planar, rectangular, the two-dimensional coordinates r = (x, y) of a point of S being expressed according to the x and y Cartesian coordinates contained in the surface S along axes X and X respectively. ordinates Y perpendicular to each other, and z being the position of the plane S along an axis of height Z perpendicular to said plane, the coordinates of some of these points may be provided with indices, for example Fo = (xo, yo) ), r, _ (xyy), r- _ (x, y,.) It is understood that the skilled person is able to adapt the teaching of the invention for other types of surface S and / or other types of coordinates.In particular, S can be plane and circular and r and all equations can be expressed with polar coordinates (r, 9) contained in the plane S, which is of interest especially if we consider intensity images in the form of a disk. To simplify the notations, we will omit the variable z in the list of coordinates of the intensity I. Nevertheless, when calculating the derivative of the intensity with respect to the variable z, this variable will be implicitly included in the list of contact information.

Moreover, the phase c1 is supposed to change very little during the propagation of the light, and consequently we will omit the variable z from the list of coordinates for the phase cl). The gradient operator V is computed with respect to the variable F, that is to say in Cartesian coordinates in the directions perpendicular to z i.e. V = [a / ax, a / ay] T. In particular, O, o is the gradient with respect to ro, and D; is the gradient with respect to F.

A vector or vector product between parentheses [] Z marked with an index z means the component of that vector or vector product along the Z axis.

The integral transforms kt (k) cosine (ICT) and sinus (IST) are defined for any function K (x) by: For ICT:

x (k) =? Jcos (kx) K (x) dx, 0 And, for the IST: k) = 12 Jsin (kx) K (x) dx, 7r o To numerically compute the cosine and sine integral transforms, we use the discrete cosine transforms of type 1 and sinus type 1.

The discrete transforms kc [k] cosine of type 1 (DCT1) and sine of type 1 (DST1), in one embodiment, can be defined for any discrete function K [j] by: For DCT1: nx [k] = EE [j] K [j] cos (j k7r / n), j = 0

for k = 0, ..., n, with:

0.5 if j = 0 or j = n, 1.0 in other cases. (1) (2) (3) (4) And, for the DST1: 292231 6 -10u i [k] = En-2 x [j] sin [(j + 1) (k + 1) z / n] , i = o (5) For k = 0, ..., nû2. Algorithms for making these transforms are known to those skilled in the art and will therefore not be precisely described. The vast majority of algorithms, written in the form of C / C ++ program libraries for the implementation of sine and cosine transforms, are based on the use of Fourier Transforms which, in turn, are implemented in the form of the Fast Fourier Transform (FFT) algorithm (see for example Ref [9] and Ref [10]). Indeed, it is sufficient to do a search on the Internet to find algorithms using many versions of the use of the FFT to calculate the 10 cosine and sinus transforms. The Fast Fourier Transform commonly used uses complex number arithmetic: to calculate the cosine and sine transforms of a real vector of 1V digits, the size of the Fast Fourier transform must be at least 2V. Lesser known are recursive algorithms, based on linear algebra calculation, which perform the computation of cosine and sinus transforms using real arithmetic: they do not use complex numbers. In the Ref. [11], a recursive algorithm (i.e., using a series of uniform and consecutive computations, whose input data is the output data of the previous equivalent computations performed with the number of elements in the divided processed data in two) is presented, which allows to calculate the cosine or sine transform, without using the Fourier Transform. To increase the speed of the calculation, the method according to the invention will preferably use, for calculating these cosine and sinus transforms, an algorithm that does not use fast Fourier transforms, such as for example a recursive algorithm; we will then speak of cosine transforms or fast sinuses. Finally, we use the convention: in =} m, n}, =} k, l}, (6) i =} 1, j}. 2922316 -11 In this document, an expression, formula, equation, ... given with some indices among in, Tc, T, m, n, k, l, i, j, ... remains valid re-expressed with d other clues.

We will first describe some equations and algorithms used

5 by the device and the method according to the invention of Figures 1 to 3.

The equations of light energy transport in paraxial approximation give a link between the speed of change of the flow of light energy in the direction of propagation of the wave field (direction z), and the rate of change of the flux in the plane S 10 perpendicular to the direction of propagation z. This link is described by the equation between the three components of the Poynting vector [I (Y, z) V (Y), 29LX II (r, z)] T

V. [I (i:, z) V a)] = - 2rci.-1 al (r'z) az where. is the average wavelength of the wave field for the propagation environment of the wave field. The environment of

Propagation notably comprises the reflection or transmission media of the wave field. According to Helmholtz's theorem, a vector field I (F) oc (r) can be decomposed into a potential part O (r) which satisfies the identity V x VO (F) = and in the part without divergence V x A ( @): 1 (i) v (F) = vo (F) + VxA (r). (8) Similarly, the gradient of the Ocl phase (r) can be decomposed into a potential part V (F) = o) and a rotational part I (F)

vxB (Y) = vxA (i ^) I (r) Vc (r) = Vyr (Y) + oXB (r) (9) The rotational part V x B (~ ^) of the field o (D (^) is particularly due to the presence of phase singularities. (7) -12 With reference to FIG. 1, the singularities (or singular points) of phase are located in the isolated points r, E R2, around which the circulation co (F ) (also called the singularity value) of o (D (F) along a contour C (@ 1) in S and with small radius is not zero but is equal to a multiple of 2i w ( FI) = ck C (@,) is traveled in the opposite direction of the needles of a clock, which we take as the positive direction, around the point rä and c is the unit vector tangent to the contour C (;) and oriented in The number F is called the order of the singularity In optics, the order of a singularity can only take integer values: F = 1,2,3,. ... In addition, the probability that the order of a singularity exceeds 1 is close to ze ro, because a singularity of the order F is decomposed into F singularities of the order 1 in the presence of an infinitely small perturbation during the propagation of

the light. For this reason, in practice, it is considered that there are only typically singularities of the order 1. The singularity is positive if the sign of w (F1) is positive in equation (10), otherwise she is negative. It follows from the principle of the conservation of energy that there are as many positive singularities as there are negative singularities, so that the sum

Signed orders of all singularities are zero. For a two-dimensional (I) (F) phase image (co-ordinates r) decomposed into pixels, it can be said that there is a singularity pixel singularity if it exists between this singularity pixel and a pixel adjacent to this singularity pixel, a difference in the value of the phase c (1 ^) substantially equal to a multiple of r radians, ie a difference in height along the Z axis substantially equal to a multiple of - in the case where the wave field has been transmitted by the imaged object, or 2

a difference in height along the Z axis substantially equal to a -13 multiple of 4 in the case where the wave field has been reflected by the imaged object. This is typically the case if the device according to the invention images an object having important reliefs with large slopes. In practice, it is noted that the intensity is generally substantially equal to zero in a singular point: I (r,) = 0, r, being the coordinates of a singular point. (11) According to the Stokes theorem, the vile rotor (r) is: vxvcI *.) = Vxvx130 = lco (Y,) 11, (Y,) Ô (YùY,), (12)! = 0 where 3 (FyY,) is the function of Dirac and n, (r,) is a unit vector perpendicular to the plane S and oriented in the positive direction of the Z axis, with as origin the coordinates r, of the lem singularity, co ( F ,,) being the value of the singularity at the point r,. In the practical situation where only samples of the phase on the S = (x, y) plane are known for pixels of size Ax along the X and Dy axis along the Y axis, a value of a singularity in r, is defined using the integration, on the contour C, of the components of the gradient of the phase around the point r, using the finite differences, as represented on the figure 1: vcD.di ° c (x ,,) +0 (D (x1 + Ax, y,) zW (x, ~ + 4y) 0I (x, y,) yy Ofi (x, yr) _ [fi (x, + Ax) , y,) ùcD (x, y,) Umm, (13) = ['v (x ,, y, + Ay) û (x1UAy Qy. We will show how to construct an image of the phase cb (r) d a wave field from the longitudinal component of the Poynting vector of the wave field (i.e., the intensity times the wavenumber 2i1 - ') and the variation of this vector component of Poynting along a z-direction of propagation of the wave field, taking into account the phase singularities -14 For the vector potential  (F) corresponding to the r Supporting vector S, there exists a function of Green g (Y ', ro) which reflects the geometry of S and the boundary conditions on the potential according to the equation: A (Y) ù ffsLvro xvro xA (YO)] g ( Y ~ YO) ~ o • (14) According to the Poisson equation, the potential part vyp (F) of the gradient of the phase v (D (r) can thus be rewritten in the form: vr 2) = I (r) 1 Ifs [((o Y, z) laz] vg (F,) dr0. (15) w (~ o Similarly, the rotational part V x B (r) of the gradient of the phase vcl) (F) can be rewritten as: vxB (r) = ùJ) ffs [vJ (YO) xvcp (F0)] xvg (Y, F0) dro (16) The vector product VI (F) xvc (F) is a vector oriented along the Z axis, and the rotor vx  (r) has its components in the plane S. Let the function G (F, ro) linked to the function of Green g (i ^, r''o) by the conditions of Cauchy-Riemann: aG ag ax ay '(17) DG ag ay ax The part The rotational V x B (r) of the gradient of the phase vc (F) can therefore be rewritten as: vxB (Y) = ùI (@) f f [VI (Y0, z) xvq) (YO)] ZvG ( Y YO) dFo.

We thus have: 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 where 8I (Fo) = al (F0, z) / az This equation represents the general solution at -15

calculation problem of the vector field \ 'D from measurements of the longitudinal component of the Poynting vector (i.e. intensity) and its variation along the axial direction z. The vector product of this equation by \ I (F) gives the following second-order Fredholm integral equation: f (r) = b (r) û JJS f (F) V (F, r0) ur0, (19) or .

f (j) = {VI (Y) xvcD (F) L 1 ~ -b (r) = JJ, S 27L / ~, 1 (SIfro) [V1n1 (;) xVgfr Fo)] d V (F, r 'o) = [V1n1 (F) xVG (F, ro) 17, where In is the natural logarithmic function. It is necessary to obtain the function f (r) for:

obtaining the α (F) gradient of the total phase according to formula (18); this gradient will include both the potential part and the rotational part; the calculation of the total phase. The total phase is calculated by integrating the gradient v (D (F) of the total phase, the integration of the potential part is straightforward, while the integration of the rotational part requires the placement of cleavage between singularities of opposite signs The positions and signs of the singularities are determined by means of the function f (r).

According to the invention, it can be shown from an integration of the formula (18) that the phase (D (F) is the sum of a potential part tg (F) and a rotational part also called hidden phase x (r)

(cr) = (r) + x (r), (21) the hidden phase x (F) being the function bound by the Cauchy-Riemann conditions with the function W (F) = w (r ',) g ( F, r,), that is to say in Cartesian 1 coordinates: (20) -16 ax`'J = 1 w (ri) gy (r, rl), ax ax (r) - E w ( rl) gx (F, rl), where a 1 signifies a sum over 1 with 1 varying from 0 to L-1, where L is the number of singularities, gx (F, Fo) being the derivative of g (F, ro) by relative to the component x of the vector r = (x, y), and gy (F, ro) being the derivative of g (F, io) with respect to the component y of the vector r = (x, y). The potential part tg (F) of the phase is equal to a first estimate of the phase 0,0 to which a harmonic phase x (x, y) is subtracted: W (r) = 01 (Y) û '(x, Y)

With: (F) = f J (Fo) • ~ TbW., R0) ui0 X (x, y) = JJS [\ x (F)) J.vrog (r, r0) dio = ù E rw (rl) J Referring to FIG. 1, the device 1 according to the invention embodying a method according to the invention comprises an optical system 2 (FIG. , a sensor 3 such as a CCD camera, a calculation unit 4 such as a central unit of a computer, a screen 5, and means (not shown) for illuminating an object 6.

The illumination means may be arranged to illuminate the object 6 in transmission or in reflection, so that the optical system 2 is arranged to project on the sensor 3 the light wave field 7 respectively transmitted or reflected by the object 6 and propagating parallel to the Z axis. The imaged object 6 is substantially in the plane

object focal 8 of the optical system 2.

The optical axis of the optical system is aligned with the Z axis. The device 1 further comprises means (not shown) for modifying the relative position of the sensor 3 along the axis Z with respect to the system (22) ( 23) (24) - 17 optical 2 secured to the object 6. The device may therefore comprise either means for moving the sensor 3, or means for moving the optical system 2 integrally with the object 6. The sensor 3 comprises a two-dimensional grid of rectangular pixels arranged in a plane perpendicular to the Z axis and respectively aligned along M lines parallel to the Y axis and N columns parallel to the X axis, the X and Y axes being perpendicular to each other and The sensor is thus arranged to acquire, for several possible positions of the sensor along the axis Z, an intensity image of the wave field, each intensity image being delimited by a flat surface. perpendicular to the Z axis, and comprising M times N poin ts or pixels. In practice, these M times N pixels may be only a subset of the total number of pixels of the sensor. The region S can be chosen by the user, and therefore it can contain a number M times N of pixels less than the total number of pixels that the sensor has. For example, the sensor 3 contains 1280 x 1024 pixels. However, the user can choose a calculation area S of arbitrary size, for example 268 x 873 pixels, without however exceeding the useful surface of the sensor, that is to say without going beyond the size of 1280 x 1024 pixels . The sensor 3 is connected to the computing unit. The calculation unit 4 stores software according to the invention and is arranged to calculate different data 1 (F), aI (F, z) / az, f (r), v11) (F), O, (Y) , x (F), x (x, y), (D (r) for each of these points, the screen 5 makes it possible to display these data in the form of graphs.

We will now describe, with reference to Figures 2 to 10, different steps of the method of calculating the phase according to the invention. These successive steps A to L comprise: A. An acquisition, by the sensor 3, of at least one pair of images of intensity 1, (r) = I (F, z) and I20 = I (F, z2) of the wave field coming from the object 6 (by reflection or by transmission), the images of the couple being acquired for positions 3a, 3b of the sensor 3 symmetrical with respect to the image focal plane 9 of the optical system 2, z , and z, being, along an axis Z, the respective positions of the planar surface respectively imaged S, -18 and S2 for each of the images, x and y being coordinates along respectively perpendicular axes X and Y between them and at the Z axis. S 1 and S 2 are plane surfaces perpendicular to the Z axis. All the acquired images have the same spatial scale and dimensions along the X and Y axes, which necessitates scaling of the raw data acquired by the sensor 3. The calculation unit 4 then calculates the value of the intensity I (FF) of the wave field for each point F of the plane surface S perpendicular to the axis Z and located at a position z along the Z axis, and calculating the first derivative SI (r) = al (F, z) / az of the intensity on the surface S with respect to the coordinate z for each point r of the S. S plane surface has the same spatial scale and dimensions along the X and Y axes as S, and S2. For the computation of the values of the intensity I (F), the intensity I (i) at the point 15 of coordinates F in S is typically equal to an average of the values at the point of coordinates r in S, and S2 of the different images. of intensity acquired respectively 1, (F) and 12 (F). Subsequently, we consider the case where: I (F) = [I, (r) + I2 (F) 1/2, (25) z = (z, + z2) / 2 But we could have chosen n ' any value for 1 (F) which is a more or less weighted average of acquired intensity images: I (F) = 1, (7) or I (r) = I2 (i-) or I (F) = [I, (i °) +212 (r) 1i3, etc. (26) S being between S, and S2, z being an average of z, and z2 according to this same weighting. Next we consider the non-limiting case for which the image area S is a rectangle delimited by four edges, each of the edges of this rectangle being aligned with one of the axes of the x and y coordinates used for the calculations. We have S = (0, a) x (O, b), that is to say that two edges of S aligned along the X axis have a length a, and two edges of S 292231 6 -19 aligned along the Y axis have a length b. Just like 11 (F) and 12 (F), the image I (F) comprises M times N points. The distance between two points of I (F) along the X axis is Ox, and the distance between two points of I (F) along the Y axis is Ay. A point of S imaged by a pixel located on the ith line (0 <û i <û M -1) and on the jth column (0 <û j <û N -1) of the sensor will have its coordinates designated by Fi _ ( x ;, y;), with x; = i & and y; = j Ay. Similarly, in the case of an acquisition of images, one could have chosen: I (r) = [I1 (r) + I2 (r) + ... + (r)] / n (27) The calculation of the first derivative SI (F) = al (ïF, z) / az is produced by an approximation. This approximation can be written as the difference: 31 (F) = [12 (F) -I1 (Y) 1 (Z2 -z1) (28) More generally, we could have chosen any value for 31 ( F) representative of a difference between intensity images acquired for several surfaces perpendicular to the Z axis, for example in the case of acquisition of four images 11 (F), 12 (F), 13 ( F) and 14 (F) whose respective positions of the planes along Z are z1, z2, z3 and z4: 31 (F) {[13 (r) +14 (r)] - [II (r) +12 (r)]] l [(Z3 + Z4) ù (Z1 + Z2)] (29) B. A calculation, by unit 4, and from intensity values 1 (F) and the derivative of the intensity 8I (i) = a (rF, z) / az, of the singularity function f (r) = [VI (F) xvi (F)], the value of which is significant for the presence or the absence in r of a singularity. This calculation is carried out by the direct, that is to say non-iterative, resolution of the following equation: f (Y) û b (,) û JjC .f (F0) V (F, F0) do, ( 30) which is a Fredholm equation of the second order, and where: -20

b (r) = 27r. IF (ro) [VlnI (F) xVrg (F, ro)] dFo V (F, Fo) = [V lnI (r) xV7.G (F, F0)] 1, being the average wavelength of the wave field, which is in practice almost monochromatic centered on a wavelength). with a bandwidth of a few nanometers. g (F, Fo) is the Green function that reflects the geometry of the support region S and the boundary conditions on the potential. The function of Green g (F, Fo) can be defined as the expression in the form of a series of eigenfunctions (F ,, (r) Fm (Fo)) of the two-dimensional Helmholtz differential equation of preferably for the boundary conditions of Neumann or Dirichlet in S: g (r, ro) = E v 2 F ~ n (r) Fm (ro), (32), where n = 0 u being the eigenvalues, with: = m7c / a 7 ~ n = n ir / b Um2 4m2 + fn2 where and / h = 0

. oo represents all the indices required to define a required eigenfunction, and the decomposition functions Fm (F) (or FF (F), F, - (F) ... according to the chosen index) form the orthogonal set. normal (see ref [1] Morse & Feshbach, page 820), and are eigenfunctions of the one-dimensional Helmholtz differential equation, preferably for the boundary conditions of Neumann or Dirichlet in S ff F i) Fk (F) =. (33) For Neumann boundary conditions, the value of a component of the spatial gradient of g (F, Fo) is fixed along the edges of S, this component being contained in the plane S and normal to the edge . Typically, this value is set equal to 0, in particular if the acquired images are images of an object 6 comprising a substantially planar surface imaged in particular by the edges of the images, and a (31) - 21 - object in relief centered on images. For the boundary conditions of Neumann, we take: g (F, T "o) = An cos (riny) cos (rinyo) n = 0 x COS (mX) COS (mx0) L ~ m + nnn 2 In practice , this calculation on the infinite of one or two sums is done using one of these two approaches: 1. By replacing the infinite (m = oo and n = oo) by integers m = M and n = N sufficiently large (typically 256, this value depending on the computational power of computation unit 4) to allow the required precision of the representation of the Green function Also, to allow the calculation of the Green function in using fast transforms (for example, Fourier, Sinus, Cosine transforms, etc.), the numbers M and N are preferably equal to powers of 2, M being a power of 2 greater than or equal to m and as close as possible to M, N being a power of 2 greater than or equal to N and as close as possible to N; 2. By calculating one of the two sums analytically (that is, expressing this sum in the form of an analytic function), and replacing the infinite (oo) of the remaining sum (with respect to m or n) by an integer sufficiently large to allow the required precision of the representation of the function of Green. In the implementation of the algorithm that follows, the second of these approaches is used, where the second sum on m in equation (34) is computed analytically, while the infinite (oo) of the first sum on n from equation (34) is replaced by N which is a power of 2 greater than or equal to N and as close as possible to N. It is possible to use comparable calculation approaches for calculating other sums in the process according to 'invention. For Dirichlet boundary conditions, the value of g (r, Fo) is set along (34) m = 0, m2 + n2SO - 22 - at each edge of S. Typically, this value is set equal to 0, in particular if the acquired images are images of an object comprising a substantially flat surface imaged in particular by the edges of the images, and an object in relief centered on the images ... DTD: For the boundary conditions of Dirichlet, take: llg (;, ro) = E Am ä sin (rin y) sin (rin y0) n = 0 xm = 0, m` + n2 m0 sin (mx) sin (mxo) ) m2 + n2 (35) Typically, the second sum over m in equation (35) is computed analytically, while the infinite (oo) of the first sum over n of equation (35) is replaced by N which is a power of 2 greater than or equal to N and as close as possible to N. The parameters 4, n, n ,,, Am, n will be detailed later. The computation of f (F) is based on a decomposition of f (F) into a series of decomposition functions FF (r): .f (F) = E fi Fi (F.). (36) k = 0 fi (or fm, J. ... according to the chosen index) being the decomposition coefficients of the function f (F). In formula (36), the infinite (oo) is replaced by K = (K, L), K and L being sufficiently high values (and not necessarily powers of two, and typically between 30 and 80 depending on the power of the unit 4) to represent the function f (F) with a required accuracy. In the process according to the invention, the procedure is typically the same for the calculation of other sums up to infinity. The function b (r) is decomposed into a series of decomposition functions F ,,, (i °):

CO b (F) = E b ,, Fm (r) • (37) m = 0 b ,,, (or bk, b; ... depending on the chosen index) being the coefficients of - 23 - decomposition of the function b (F). In formula (37) the infinite (x) is replaced by M = (M, N), M and N being the sufficiently high values

(Not necessarily powers of two, and typically between 30 and 80 depending on the power of the unit 4) to represent the function b (r) with a required accuracy. In the process according to the invention,

typically proceeds in the same way for the calculation of other sums up to infinity. According to the invention, the functions f (F) and b (F) are therefore decomposed

in series of decomposition functions F,;, (F) which are functions

of the one dimensional Helmholtz equation, preferably for the Neumann or Dirichlet boundary conditions in the S domain. For the Neumann boundary conditions, f (r) and

b (F) using the Discrete Cosine Transform (DCT), that is, taking: Fm (@) = Fm n (x, y) = A ,,,, n cos (mx) cos ( r ~ ny) (38) For the boundary conditions of Dirichlet, f (r) and b (F) are decomposed using Discrete Sinus Transform (DST), that is, taking: Fm (r) ) = Fm, n (x, y) = Am, n sin (mx) sin (77 n. @), (39) The computation of f (F) comprises the steps B.1. to B.6. following: B.1. A calculation of the parameters v, 1, 4m, Tb,, Am and Bn according to equation (40) and sm according to equation (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 ,,? if m> O has 1

a B ,, = ~ if n> 0 Bo 1 b if m = 0 m = if m> 0 (41) B.2.A calculation of the function b (r) = f fs27rl '81 (ro) [V1nI ( F) xVFg (r, Fo)] Z dFo For this, the calculation unit uses the calculation of: - the scalar function 810 [12 (F) -11 (F) J / Az, where Az = z2 - z, is the distance between the two shots (S, and S2) of two intensity images 11 (r) and I2 (r); and the vector function Oln1 (r), where I (rr) = [1, (r) +12 (x) 1/2 from intensity images I, (r) and 12 (r). To compute the function b (F), we calculate the first derivative gx (r, ro) of g (r, ro) with respect to the x component of the vector r = (x, y), and the first derivative gy ( r, ro) of g (r, ro) with respect to the y-component of the vector r "= (x, y): b (r) = f fs 22r) 'SI (ro) [VinI (r) x \ ## STR2 ## wherein R 1 (R 1), R 1, R 1, R 1, R 1, R 1, R 1, R 1, R 1, R 1, R 1, R 1, R 1, R 1, R 1, R 1, R 1, R 2, R 1, R 2 For the boundary conditions of Neumann, the derivatives of the function of Green with respect to the variables x and y are: ## EQU1 ##

-25-

t gx (ri0) = ù2a-1E sin (S ,,, x) cos (4, nxo) J (y, Yo) (43) m = 1 t gy (i ', r0) = - 2bu1 sin (r) ny) cos (nny0) / n (x, x0). (44) n = 1 These equations show us that fm and f 'are measurements of a representation of a function of Green g (F, i ^ o) in a space of eigenfunctions of the coordinate Helmholtz equation separable.

For Dirichlet boundary conditions, the derivatives of the Green function are those given by the two formulas above in which the sine function is replaced by the cosine function and vice versa. That is, for the boundary conditions of Neumann: b (x, y) = fis -13I (xo, y0) [01nI (x, Y) xV g (x, y; xo, YO)] _ dx0dy0 1 1 al (x 'y) / axle sin (riny) f f. (x, xo) f 31 (xo, Yo) cos (rinyo) dyo dxo (45) + 4zc ~ 1 a-1 aI (x'Y) / ay ° sin x 8I x, yo) cos (, nxo) dxo dyo I (x, y) "7, (~, n) fo _0 fm (y ~ Yo) f (o For the boundary conditions of Dirichlet, it suffices to rewrite the preceding equation by replacing the cosine function with the function sinus and vice versa For the boundary conditions of Neumann: fx, xo) û_ 1 cosh [rin (a û x)] cosh (rinxo), if x> xo sinh (yn7r) cosh [n ,, (a xo) j cosh (rinx), if x <xo 1 fcosh [n, (b û y)] coshyo), if y> 0 .i, n (Y'Yo) = sinh (ÿm7r) cosh [n, (b yy) )] cosh (n, y), if y <yo and, for the boundary conditions of Dirichlet, / nr (x, x0) _ 1 sinh [rhi (a x)) sinh (rinxo), if x> x0 Sinh (Y /! TC) Sinh [r) n (a xo)] S1n11 (nnx), if x <xo Y 1 sinh [n, (b y Y)] s (mYo) 'if y> Yo Î , ~ (Y ~ yo) = sinh () 7rnir) sinh [{,,, (b ù yo)] sinh (4n, y), if y <Yo The samples of an image I (x, y) being known on the discrete grid and -26 {x ,, y;}, or x, = iAx, i = O, ..., Mû1 and y, = jAy, j = O, ..., N 1, the integrals between 0 and a = MAx and between 0 and b = NAy in equation (45) are approximated, for the Neumann boundary conditions, by a sum: b (x;, Y,) _ _4- 1 al al N N ((((((((((((((((((((( 0 = 1 Jo = 1 N / (46) 47r ~, - 0y ôl x, y / yy 7c + '1 sin mùZ If:, (y, y, o) SI (x, 0, y, o) cos MI (x ,, Yi) m = 1 M ,, 0 = 1 10 = 1 M / For Dirichlet boundary conditions, it is sufficient to rewrite the previous equation by replacing the cosine function with the sine function and vice versa. B.3. A calculation of the coefficients bm according to the equation: bm û JJsb (r) F; n (@) dF. (47) To find the coefficients, the calculation unit 4 can use conventional methods of calculation such as the transform type cosine or sine, preferably cosine or fast sinus. For the Neumann boundary conditions, this calculation is done using the cosine transforms: ben = JJsb (r) F; n (r) dr ah = Am n Jcos (4mx) cos (rinY) b (x, y) dy 0 0 For Dirichlet boundary conditions, this calculation is done using sinus transforms: bm = J JSb (r) F, - (F) dr ah = Am n Jsin (4mx) Jsin (riny) b (x, y) dy dx 0 0 M- / ù Nù1 Amn Ax Ayl sin miù 1 sin nj = o Al, = o N The sums in equation (48) and (49) are computed by the dx (48) ## EQU1 ## ## EQU1 ## (b) (x, M Jo N / i)] (49)) b (x ;, 27 cosine or sine transforms, preferably cosine or fast sine. B.4 A calculation of coefficients smk, comprising the following steps B.4.1 and B.4.2: B.4.1 A calculation of the values wm and w ,, being part of the matrix wm _ wm = (50) The calculation of the wm and wm values is done using the Discrete Sinus Transform (DST) algorithm for Neumann boundary conditions or using the Neumann algorithm. the Discrete Cosine Transform (DCT) for the Dirichlet boundary conditions of the components Lx (r) and Ly (r) of the vector function i (F) = VlnI (r) respectively according to x and according to y. For the Neumann boundary conditions, let L (r) = \ 1nI (r) = lw,;, • VF ,, (r) Lx (r) = Ewm VF, (r) lax Ly (r) = wm VF, (F) / ay Lx (i -) = Ewm aFm (Y) / ax m wm aFm (Y) / ax MN = ù ~ m S1Il mx `4m nWm n cos (y) m = 1 n = 1 Ly (Y) ùwm), / Ôy ~ aa _ wm UFm (Y) / ay and NM y = nn Sln (1in.Y) 1 Am, nw m, n COS (mx) n = 1 m =] (51) ( 52) (53)

-28 For the Dirichlet boundary conditions, simply repeat these equations by replacing the sines with cosines and replacing the cosines with sines. The cosine and sinus transforms, preferably fast ones, make it possible to find the values wm and wy. For the boundary conditions of Neumann: ah wm, n = ù ml Jsin (mx) [L (xY) cos (ni, y) dy dx 0 0 Mi1 / A N1 (7t \ 1yl sin mii E Lx (x ;, y1) cos nj = o \. M / i = o \ N / m = 1, ..., n = 1, ..., N,

(54) Mu1 .7r - a w ,,; = û Jsin (riny) JLy, (x, Y) cos (mx) dx dy o _o

= AyE sin / nj N \ E Ly (x ,, yJ) cos / mi J j = o /, = o M m = 1, ..., M, n = 1, ..., N. For the Dirichlet boundary conditions, it is sufficient to take these

equations by replacing the sines with cosines and replacing the cosines with sinuses.

B.4.2. A smk = uk-2 wi JJS [vF (r) • vFk (F)] F,;, (r) di '. For the Neumann boundary conditions, the coefficients are: (55) coefficients of calculation -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 AB (k + m4kWk + m, l + n + fl1 + n1JiWk + m , l + n) EmEn mn Ak + m Ak B1-n B, K x Y l + (k + m kWk + m, ln + 11l-nnlWil + m, In) Em n Am Bk l Ak + m Ak Bn 1 Bl (x YJ + k + m kWk + m, n-1 -77n-1 ~ lWk + m, nI Em Am B ,, Ak_mAk BI + n BI x Y + ((k-mkWk-m, 1 + n + L En En En Ak Ak BI BI BI BI BI BI BI BI BI BI BI BI BI BI BI BI BI BI BI BI BI BI BI BI BI BI BI BI BI BI BI BI BI BI BI m4kWk-m, ln + nl-nnlWy km, ln EmEn ((km kWkx-m, n-1 ûrIn-111iWy km, n-1 Em xy ((mk kWm-k, l + n nl + nfllWm-k, l + n In (In-k kWmx -k, ln _11l-nnlWmy -k, ln In y (mk kWm-k, n-1 + nn-lnlWr-k, nI (56) Am Bn Am_kAk BI-n B, Am In sum (56), the components of the vector ï = (i, j) change as follows: i = 1, ... ,, j = 1, ..., N, the parameters uk2 ,, 4m, rin, Am and Bn and Em being previously described For the Dirichlet boundary conditions, the coefficients are: Smk 4Uk 2 Wi if [v) .VFk (Y)] F, (F) a 4 Uk 2 Wl, j {4i k ('8ù1-km + Si + k, m) (8j-1, n gj + l, n) (57) + 11j71 (Si-k, m ù (5i + k, m) (S j-1, n + (-j + I, n)} where the components of the vector i = (i, j) change as follows: i = 1, ..., M, j = 1, ..., N and where: ## EQU1 ## A calculation of the coefficients fk by solving the system of equations: E (swk + Sn ~ k) Jk = m = 1, ..., M (59) k = 1 (58) -30 where m = 1, .. ., M means that the vector m travels all the values of m = 1 = (m = 1, n-1) to m = M = (m = M, n = IV), where k = 1, ..., K means that the vector travels all values from = 1 = (k = 1,1 = 1) to k = K = (k = K, l = L), and where: 1sim = k,

Smk 0 in the other cases. Thus, the system of linear equations (59) contains MN linear equations to find the unknown K L which are the Ki coefficients fi of the decomposition of f (r) into a series of functions Fi (F) according to the formula (36). This system can be solved by means of a classical algorithm of

10 resolution of systems of linear equations, such as a decomposition LU (Low Upper Decomposition) or a decomposition SVD (Singular Value Decomposition). Bibliographic reference [9] gives a preferential list of system solving methods of linear equations. B.6. A calculation of the function f (F) from the coefficients fi obtained at the end of the preceding step, by the use of the equation: ## EQU1 ## Fk (F). (61) T To reconstruct the values of the samples of an image f (x, y) on the discrete grid {x ,, y1}, or x; = iAx, i = 0, ..., M1 and y. = jAy, j = 0, ..., N -1, for the Neumann boundary conditions, the sum (61) is effected by means of the cosine transform, preferably the fast cosine:

K \ (Z f (x ,, y1) = 1 coskiù Ak, / fk, lcos l ~ k = 1 mi N,

i = 0, ..., M1, j = 0, ..., N1. For the Dirichlet boundary conditions, the sum (61) is effected by means of the sinus transform, preferably the fast sinus: (60) (62) -31 K / 'L 7 7rf (x;, y;) = E sin kiu EAk, rfk ,, sin lj-k = 1 M. / l = 1 N ~ i = 0, ..., M -1, j = 0, ..., N -1. C. A determination, by the unit 4, of the presence of singularities in the phase comprising the steps C.1. at C.3. following: C.1. A formulation of a threshold fs; C.2. If there are points r in the support region S where f (F)> fs, singularities are present in the phase; C.3. If f (F) l <fs for every point r of S, the singularities are absent in the phase. D. If there are no singularities, a calculation, by the unit 4, of a first estimate of the phase comprising the steps D.I. and D.2. following: D.1. A calculation of the phase gradient according to equation (64), ie taking into account only the first term of the equation, the modified formula (18) assuming f (i) 0 for all 1 vee Y 21 (SI Y, z YY The components of the gradient are calculated according to the equation: -1 (i) 21 () If afro) g Y, Y 0 d: 0 cl) y (F) = 21 (r) fis 8I (r) gy (Y, Y0) dY 0 \ 0 I (With 1) X (7-) the component along the X axis of the gradient vcD (r) and component along the Y axis of the gradient V / c (@). The result of the insertion of equations (43) and (44) is, for the boundary conditions of Neumann: X (a sin x J) SI (x) cos (x) dx d ~ x, y) = ~ ( fm (y, y0 o ~ yo moa yo I (x, y) m = 1 oo _ (66) 47c2, a -b (D) (x, y) = I (x, y) E sin (nn y) (f ,, (x, xo) J SI (xo, yo) cos (i y0) dy0 dx0 For the Dirichlet boundary conditions, simply rewrite (63) (64) (65) (D), (r) the previous equation by replacing the cosine function with the sine function and vice versa, the samples of an image I (x, y) being known on the discrete grid {xi, yj} 1 or x; idx, i = 0, ..., Mû1 and yj = j 0y, j = 0, ..., Nû1, the integrals 5 between 0 and a = Mdx and between 0 and b = NAy in equation (66) are, for the boundary conditions of Neumann, approximated by a sum: (re Mii Mi1 / 7r X (xi, yj) = Y / M 1 sin m 7t / ij N1lf, n yj '@ ll jol & (xi, yjol / cos muiio -47LX 10 JI (xi, yj) m = 1 v M jo = 1-io = 1 M-47Ci -1 dx / N N-1 (7t \ M-1 () - N-1 l ~ Y (xi Y j) _ 1 sin nù jlf ,, \ xi 'xio l 51 (xio, yjo) cos Hi ° I (xi,) n = 1 N io = 1 _jo = 1 For the Dirichlet boundary conditions, simply rewrite the previous equation by replacing the cosine function by the sine function and vice versa. These equations show that, in order to calculate the gradient of the phase, the following operations are carried out for the boundary conditions respectively of Neumann or Dirichlet: a respectively cosine or sine transform applied on the derivative of the intensity 8I ( ro, z), then, on the result of this transform, an application of a filter fm (y, yo) or Mx, in the space domain, then on the result of this filtering, an application of a transformed respectively of sine or cosine, then - a division by the intensity I (r). 20 D.2. A calculation of the first estimate of the MF phase) by integration on the surface S of the gradient of the phase V (r) according to the equation: (F) = J J, 4 vcp (i0) vrog (r, F0) di0. (68) In terms of the phase V (xo, yo) components of the phase, the first estimate of the phase is calculated, for the Neumann boundary conditions, according to the equation: (67) (70) -33 - E COS (4mx) Jfä2 (y, yo) ~ EDxo (xo, y0) sin (4mxo) dx0 dy am = ~ oa - 41L / i -1 0o _ ba (69) - 'ab + V cos (rin y ) Jf, (x, xo) J x dx Yo (, y0) sin (rny0) dy0 0 oo _ With c> px, (xo, yo) the component along the X axis of the gradient of the phase VO (xo, y0) in r0 and (Dyo (x0, y0) the component along the Y axis of the gradient of the phase V'1 (x0, y0) in ro For the boundary conditions of Dirichlet, it suffices to rewrite the equation by replacing the cosine function with the sine function and vice versa, from the data cx (xi, y) and y (xi, y) obtained by equation (67), the first estimate of the phase is in calculated practice, for the boundary conditions of Neumann, according to the equation: For Dirichlet boundary conditions, simply rewrite 10 the previous equation replacing the cosine function with the sine function and vice versa. These equations show that, in order to calculate the first estimate of the phase, the following operations are carried out for the Neumann or Dirichlet boundary conditions respectively: a sinus or cosine transform respectively applied on the gradient of the phase, and then - on the result of this transform, an application of a filter f, n (y, yo) or ffy (x, x0) in the spatial domain, then - on the result of this filtering, an application of a transform 20 cosine or sine respectively. Those skilled in the art will recognize that step D corresponds to the method according to WO 2006/082327. U n = 1 4 ~ jl Mu1 (ù \ 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 n0 j0 jo = 1 N / - 34 - E. If there are no singularities, the total phase is defined by unit 4 as equal to the first estimate of the MF phase calculated in step DF. singularities, the first estimate of the phase is calculated by the unit 4 according to the following steps F.1 and F.2: F.1 A calculation of the gradient of the phase, using the function 8I (F) - [12 (F) û1, (r)] / Oz and the function f (F), obtained in step B, according to the equation: vo (r 27r 'SI rrr) d (71) I (r) $$ f (ro) VG (F, ro) d) To calculate the components (DX (x ,, yj) along the X axis and cy (x ,, yj) along the Y axis of the V (x, y) on the discrete 1 x ;, yj} e S grid, this calculation takes the following form for the Neumann boundary conditions: -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; ## EQU1 ## 0, yio) cos nNJoJ J a lo n = 1 \ N i = 1 = 1 ù 4e -'Ax / NI (x;, yi) û 2Ay / MM -1 I (Xi, yj) m = 1 Ni 7r sin ~ m M i ~ fm (Yi ~ Yio) E f (xio Yio) cos ~ mm io mi lo = 1 b = 1 For the boundary conditions of Dirichlet, it is sufficient to rewrite the previous equation by replacing the cosine function with the sine function and vice versa.

These equations show that, in order to calculate the gradient of the phase, the following operations are carried out: a calculation of a first term which does not depend on the existence of singularities, this calculation comprising for the boundary conditions respectively of Neumann or of Dirichlet: - 35 a cosine or sine transform respectively applied on the derivative of the intensity, then, on the result of this transform, an application of a filter fm (y, yo) or fy (x, xo) in the space domain, then on the result of this filtering, an application of a sine or cosine transform, respectively, and a division by the intensity values I (r). a calculation of a second term which depends on the existence of singularities, the calculation of the second term comprising, for the boundary conditions respectively of Neumann or Dirichlet, respectively a cosine or sine transform applied to the singularities function, and on the result of this transform, an application of a filter fm (y, yo) or fr (x, xo) in the spatial domain, then on the result of this filtering, an application of a transform respectively of sine or cosine, then a division by the I (F) values of the intensity. an addition of the first term with the second term. 20 F.2. A calculation of the first estimate of the phase çl (r) by integration on the surface S of the gradient of the phase ((D (F) according to the equation: 01 (F) - f JS VcD (r0). (F, F0) ~ 0. (73) This calculation of q1 (F) is carried out using the formulas given in step D.2., But with phase gradient values which are calculated by the equations of step F.1 and which therefore differ from the values of the gradient of the phase of the step DIG If there are singularities, a calculation by the unit 4 of the positions and signs of the singularities according to the steps GI to G.3: 30 GI A search on the surface S, that is to say the plane (x, y), of -36u regions where f (F) fs;

G.2. A determination of the centers of gravity of these regions. The centers of the regions thus determined are the coordinates of the singularities, r,,

G.3. A determination of the signs of f (F,). These signs are the signs of the singularities present at the points F, _ (x, y,); The singularities are localized in the points where f (F) is not zero and superior in absolute value to the threshold fS: the position of the local maximum or the local minimum, or the position of the center of gravity of f (F) can serve

as the indication of the position of a phase singularity. The singularity is positive [typically co (r,) = 27r], if f (r,)> 0, and negative [typically co (;) _ -2Tr], if f (F,) <0.

L is the number of singularities identified in S. H. If there are singularities, a calculation by unit 4 of the harmonic phase ex, y) using the positions and the signs of the singularities calculated at

step G, by means of the sum over all the singularities: x (x, Y) = ùE w (F,) JJs ° 'og (Fo, Fr) XO., g (F, Fo)] dFo The calculation of the harmonic phase x (F) at a point F therefore comprises

the computation of a series, on all the singularities, of the value of the singularity multiplied by an integral on the surface S of a product of gradients of the function of Green. Note that the integral transform and the calculation of x (x, y)

depends on the number and the signs of the singularities, but does not depend on the intensity values I (F) and the derivative of the intensity BI (F) = aI (F, z) / az The calculation of this sum is realized , in the case where the function of Green g (F, Fo) is a series of eigenfunctions (F, b (FF) Fm (FO)) of the equation

two-dimensional Helmholtz differential, by the following equation: L1 x (x, Y) = w (Fi) Tr [wY (Y1, x1) wY (x, Y) -wY (x1, y1) wY (Y, x)], (74) r = o -37 where Tr is the trace of a matrix, co (r,) = 27r is the value of the singularity in a point r ', = (x /, y /), L is the total number of singularities identified in the phase image: L / 2 positive singularities and L / 2 negative singularities.

The components W, nr n (v, u) of matrices W "(v, u) are, for the boundary conditions of 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 (iu) I}, (75) m = 1 , ..., M -1, n = 1, ..., Nû1, where m is the row index of the matrix WY (v, u), n being the column index of the matrix W '(v, u), v and u being able to take for value for example x ,, y ,, x, or y. By replacing in this equation y by) 7, one obtains the expression of the components of matrices W n (v, u). We have ay = a / b, ÿ = b / a, M = a / Ax, N = b / Ay, = m7r / a, Ax and Ay are the dimensions of a point in the image. boundary conditions of 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, In practice, under the Dirichlet boundary conditions, the phase Harmonic X (r) is posited as being substantially zero. There are singularities, a grouping of the singularity singularity singularities by the unit 4, the singularities of a pair having opposite sign values, according to the steps I.1. at I.7. following: 20 I.1. Determine the low intensity continuous regions (ie for which IO <IS, where I, is a positive intensity threshold) containing the same amount of positive singularities as negative singularities. This is possible because the low intensity lines are characteristic of the large jumps in the surface topography of the imaged object 6, resulting in almost-zero reflection in the direction of the camera. These jumps are also characterized by the presence for each singularity of its counterpart with the opposite sign; I.2. For one of these regions, construct an inking line (skeleton) contained in S and whose length is minimum, which does not leave the region, and which passes through all the singularities (both positive and negative) of the region; I.3. For the same region as in step I.2., From an initial singularity belonging to the region, construct on the surface S a cut line (also called cutoff path) connecting this initial singularity to a singularity final belonging to the region. By construction of the cut line is not necessarily meant for each pair a graphic display of this line, but rather preferably a determination of coordinates of a set of points forming on the surface S the cutoff line connecting the initial singularity to the final singularity. The singularities of a couple are of opposite signs, that is, they have opposite sign values. Break lines are breaks in continuity of the phase. The cut line is a part of the anchor line of this region. Along the cut line, between the initial and final singularities, intermediate singularities may be interspersed. A counter is defined as having an equal value, at the level of the initial singularity, at the order of the initial singularity, and incrementing along the cutoff line by a value +1 for each positive order singularity. 1 encountered and a value -1 for each negative singularity of order 1 encountered. The final singularity is the first sign singularity opposite to the initial singularity for which the counter value is zero. For example: if all the singularities are of order 1, and if the initial singularity is positive, and the two singularities neighboring the initial singularity along the inking line are two positive singularities, it is necessary that the line d Inking passes through the two positive singularities and then again through three negative singularities, the third of these negative singularities along the anchor line from the initial singularity being the final singularity; in other words, the breaking line of the couple passes between the singularities of the pair by as many positive singularities of order 1 as of negative singularities of order 1. 5 I.4. Mark the initial and final singularities as paired; I. S. Iterate steps I.3. and I.4. until all the singularities belonging to the region are connected in pairs; I.6. Repeat steps I.2. at I.S. for all regions of low intensity; I.7. Organize the links in pairs between the singularities in the form of a table including the signs, the positions, and the links between the connected singularities and store it in the memory of the unit 4. The singularities belonging to the same couple will be thereafter Numbered by indices i, where c = 0, ..., L / 2û1 is the number of the pair, i = {0,1} is the number of a singularity in the pair. The values of the singularities are represented by w (i;), where _ (x,;, y,) is the location of a singularity belonging to the pair c. By definition w (r, J = ùw (r,), because the singularities of opposite signs are related, the links between the singularities of opposite signs and the orientation of the cleavage path will be represented by the set R, _ { Rc i, Oc,}, where c = 0, ..., L / 2 -1 is the number of the pair, j = 0, ..., J. -1, where J is the number of coordinate points k , j in the support region S through which the cut-off path between the singularities (ie between the points r 1 and r 1) belonging to the pair with the number c passes in each point Rc belonging to the set Rc the orientation of the clipping path is represented as the vector Ô, which vector Ôc, is tangent to the clipping path and is oriented in the direction of the path along the path from the positive singularity to the negative singularity; If there are singularities, a calculation by the unit 4 of the hidden phase x (@) -40 according to the following equation: x (@) = LEl E 2 co (r, <; gv (u, y; rr <;) of û Jgx (x, v; @r <,;) dv c = o r = o 0 0 gx (x, v; ) being the first derivative of the function of Green g (F, r, with respect to the component x, with rv = (x, v), and gy (u, y; r,) being the first derivative of the function of Green g (ru, r,;) with respect to the component y, with ru = (u, y).

The calculation of the hidden phase at a point r of the surface therefore comprises: for each pair of index c, the computation of a torque function _x y Jgy (u, y; rr <;) of the û Jgx (x , v; rr <;) dv 0 0 - a sum, over all pairs of singularities, of a torque operator Q (, r,)) applies to the torque functions w (@@, ,, For the couple of index c, the computation of the torque function w (@, @, o, @, includes a calculation of a series, on all the singularities of the pair, of the value of the singularity w (@,) multiplied by : a difference between a first integral, according to the x-coordinate, of the first derivative, according to the second coordinate y, of the function of Green, and a second integral, according to the second coordinate y, of the first derivative, according to the first coordinate x, of the function of Green We note that the calculation of x (@) depends on the number and the signs

singularities, but does not depend on the intensity values I (F) and the intensity derivative 8I (@) = aI (@, z) / az. In the case where the Green function g (@, @ 0) is a series of eigenfunctions (F ,, (@) F ,, (Fo of the two-dimensional Helmholtz differential equation, x (FF) is computed according to the equation: (77) -41 L1 ~ 1 çà, co \ Yi ,,; ex, xi c = 0 i = 0 m> 0 n> 0 l + b-'yH (xùx,.,) ùa In this equation, therefore, i has only two possible values: 0 or 1. L is the total number of singularities identified in the phase image: typically L / 2 positive and L-singularities. / 2 negative singularities of order 1, H (x û x,) is a walk function (also called Heaviside function) 5 centered on = 0, and H (y û y,;) is a function on (also called function of Heaviside) centered on .yu y, = 0. In the case of Neumann boundary conditions: FyY (y, Yi ,,;) = û [sinh (mÿir) -1x [H (yyy,) sih (mÿ7c my my) cosh ((myrc ..,) û H (y, y) cosh (m) 77r m) 4 ,,) sinhmy)], i '(x, x,) = û [sinh (nyr) ] x [H (xxx,;) 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 where H (x - x,), H (x, xx), are functions of Heaviside centered on 10 = 0, H (y - y,), and H (y, yy) are Heaviside functions centered on y - y, = 0. For the boundary conditions of Dirichlet: Fm (Y, Y ,,,) = û [sinh (m @ 7r) x {H (yyY,) 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)],) Fy (Y, Y,) û Q n (Y, Y) Fr (x 'x,,;) (78 ) (79) Q> n (x, x1,) = 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,), where H (x ù x ,), H (x,; x), are Heaviside functions centered on x ù x,, = 0, Myu y,), and H (y, ù y) are Heaviside functions centered on y , = 0. In these equations, the operator S2 (w (F, r, o,)) represents an operation, applied to a function w (r, r, J, of placement of the cleavage paths between the two singularities belonging to the pair of index c which lie in the points r, o and i, This operation contains the following steps: 1. Calculate the spatial gradient modulo 27r of the torque function w (r, r, p, r,). modulo 27r means that each value of the gradient components of the function w (7 ,,) (a first derivative of w (F, r,, @,) with respect to x and a first derivative of w (@, @, 0 , F1) with respect to y) is reduced to the interval [û7r, 7r] by the addition or subtraction of 27r 2. Integrate along a line, starting from one and the same point of the surface S for all pairs indexed by c on which the operator S2, (w (r, r ,, @, is realized (preferably, from the origin of the coordinates (x = 0, y = o)) , going to point i-, and along paths of freely chosen integrations in the surface S (preferably comprising a segment parallel to the X axis and a segment parallel to the Y axis), the gradient field obtained in step 1, by adding or removing 27r each time the integration path intersects the clipping path of the index pair c according to the following rule: a) 27r is added if the orientation Ô ,, j of the cutoff path in a point Rc, of the path of crossing cutoff and the orientation of the integration path make an angle in the range (0.7r); (82) -43

b) 27r is subtracted if the orientation δ, of the cut-off path in a point k, j of the cut-off path traversed and the orientation of the integration path form an angle in the interval (0, û7r). K. If there are singularities, a computation by the unit 4 of the total phase c (F "): of the result ç, (@) of the step F, is subtracted the harmonic phase X (x, y ) obtained in step H and is added the hidden phase x (F) obtained in step J: c (F) = oi (F) + x (F) ûx (x, y) L. A storage of the total phase in the memory of the unit 4: the total phase is equal to the result of the step E if there are no singularities, otherwise the total phase is given at the end of step K. The image of the total phase is then displayed on the screen 5.

FIG. 4 is an example of an object 6 imaged by the device of FIG. 2, and for which the method of FIG. 3 is implemented. The object comprises a plane perpendicular to the axis Z and parallel to the axes X and Y, in the middle of which is disposed a rectangular plate 11 provided with four sides, including a first side 21 parallel to the axis Y and located in continuity with the surface 10, a second side 22 parallel to the Y axis and opposite to the first side 21, and third and fourth sides 23, 24 which progressively move away from a certain slope from the plane 10 and the side 21 to the side 22. The Z axis is expressed in radians, a radian corresponding to a height along the Z axis of D (27r). The length of the projections of the sides 21 to 24 on the plane 10 is equal to the width of 20 pixels of 5.2 micrometers of sides. The slope is such that the phase changes from 0 radians to 1.257r radians over these 20 pixels. FIGS. 5 to 10 are data images calculated according to the method of FIG. 3 by the device of FIG. 2 illustrating the object of FIG. 4. These images are represented in the plane S comprising the X and Y axes. pixel size of an image is 5.2 μm x 5.2 μm respectively along X and Y. For images 5, 7, 8, and 10, plus one point of the image is dark, and the value of the data represented is low at this point. FIG. 5 illustrates the intensity I (F) calculated in step A. It is noted that, at the location of the continuous edges 22, 23 for which there is a surface discontinuity, the intensity forms a line 12. Figure 6 illustrates the function f (r) calculated in step B. Note a region 25 for which f (F)> fs and f (,)> 0. There is therefore a positive singularity 26 at the position of the center of gravity of f (r) of the region 25. Similarly, note a region 27 for which f (r)> fs and f) <0. There is therefore a singularity negative 28 at the position of the center of gravity of f (r) of the region 27. The determination of the positions of the singularities 26 and 28 then makes it possible to connect these singularities on the intensity image of FIG. 5 by a cutoff line 29 passing through the dark line 12.

FIG. 7 illustrates the first estimate of the phase, (F) calculated in step F. FIG. 8 illustrates the hidden phase x (F) calculated in step J. FIG. 9 illustrates the harmonic phase x (r) calculated in step H. In this figure, there are three zones 30, 31, 32 separated by lines parallel to the Y axis, and corresponding to values of X (F) which are larger and larger in the direction croissants. Finally, FIG. 10 illustrates the total phase '(F) = ç, (r) - x (x, y) + x (F) calculated in step K. It is noted that the result of the total phase c1 (F ) in Figure 10 is closer to the reality of the imaged object of Figure 4 than the result of the first estimate of phase 01 (F) in Figure 7. Those skilled in the art will therefore understand that taking on account of the singularities is an obvious advantage of a method and a device for calculating the phase according to the invention compared to the state of the art, in particular with respect to what is described in the document WO 2006 / 082,327. Of course, the invention is not limited to the examples which have just been described and numerous adjustments can be made to these examples without departing from the scope of the invention. In particular, acquired intensity images could be circular. In this case, the method according to the invention is carried out as described above, but using all the formulas and equations just described expressed in polar coordinates, the polar coordinates being preferably centered on the center of the acquired circular images. In addition, when it is said that a quantity is calculated according to a given equation or formula, this quantity may more generally be calculated according to the invention substantially by this formula.

REFERENCES: [1] P.M.Morse, H.Feshbach, Methods of Theoretical Physics, McGraw-Hill, New York, 1953. [2] E.G. Abramochkin, V.G. Volostnikov, Relationship between twodimensional intensity and Fresnel phase diffraction zone, Optics Communications, vol. 74, no. 3,4, pages 144-148, December 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, pp. 921-925, November 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, pp. 4536-4540, July 1998. [5] Nugent et al, WO 00 26 622. [6] I. Lyuboshenko, Obtaining a phase image from an intensity image, patent application WO 2006/082327. [7] AP Prudnikov, YA Brychkov, IO Maichev, Integrals and Series, Gordon and Breach, New York, 1986. -46 [8] IS Gradshteyn and IM Ryzhik, Table of Integrals, Series and Products, Academic, New York, 1980 [9] W. Press, SA Teukolsky, WT Vetterling, BP 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).

Claims (25)

A method of calculating the phase (c (@)) of a wave field, comprising: -a determination (A) of intensity values (1 (F)) of the wave field at several points (F) a surface (S), - a determination (A) of a variation, at the points (F) of the surface, of the intensity of the wave field with respect to a coordinate (z) defined along an axis (Z) perpendicular to the surface, preferably a first derivative of the intensity of the wave field relative to the coordinate (z) perpendicular to the surface, characterized in that it further comprises a calculation (B) of a function of singularities (f (F)) representative of an absence or existence of a phase singularity for each point (F) of the surface (S), and - a calculation (D to K) of the phase (i) of the wave field on the surface (S) which depends on the existence or the absence of singularities. 2. Method according to claim 1, characterized in that the function of singularities is obtained by solving a second-order Fredholm equation. 3. Method according to claim 1 or 2, characterized in that the calculation (B) of the singularities function comprises a decomposition of the singularities function (f (F)) into a series of decomposition functions. 4. Method according to claim 3, characterized in that the decomposition functions are eigenfunctions of the one-dimensional Helmholtz differential equation, preferably for boundary conditions of Neumann or Dirichlet in the surface (S). 48 5. Method according to claim 3 or 4, characterized in that the decomposition of the function of singularities is of the type: f (F) = lf Fk (FF) where f (F) is the function of singularities, k = { k, 1} is a k = 1 K vector, E meaning that the vector runs all the values of k = 1 k = 1 = (k = 1,1 = 1) to k = K = (k = K, l = L), K and L being two integers, Fk (i-) are the functions of decomposition, fi are coefficients of decomposition of the function of singularities. 6. Method according to claim 5, characterized in that it comprises a calculation of decomposition coefficients by resolution of a linear system of the type: E m- + sm-) fk = b- & m = 1, ... , M k = 1 where m = {m, n} is a vector, m = 1, ..., M meaning that the vector m travels all the values of m = 1 = (m = 1, n = 1) to m = M = (m = M, n = N), M and N being two integers, k = {k, l} is a vector, E meaning that the vector runs all values of k = Î ii = 1 = (k = 1,1 = 1) at k = K = (k = K, l = L), K and L being two integers, S is the surface, the surface being a rectangle whose sides have for lengths a and h, F and Fo are coordinates of points of the surface, fi are the coefficients of decomposition of the function of singularities, 1sim = k,, b- = JJ s (F) F-r ,,, (F), Fm (F) ,,, (F) and F r, (F) and F 1 (1) br dr are 0 in other cases. the functions of decomposition, b (F) = IL 27r ~ -'SI (ro) [V lnl (F) x ~ g (r, ro) dio, is the average wavelength of the wave field, I ( F) is the value of the intensity at the point of coordinates F, z is a coordinate along the axis (Z) perpendicular to Smk-49 the surface, 31 (ro) is the first derivative of the intensity in Fo with respect to the coordinate along the axis (Z) perpendicular to the surface, g (F, FO) is a Green function used for calculating the singularity function, smk = uk z wf [vF- (F ) (F)], (F) di ', AND sVk,! - k n12 'k = kIr / a, n, = lzr / b, the matrix w- being defined by V1n1 (F) = E w- • VF (r). 7. Method according to claim 6, characterized in that the Green function used for calculating the singularities function is a series of eigenfunctions of the two-dimensional Helmholtz differential equation, preferably for Neumann boundary conditions. or Dirichlet in the surface (S). 8. Method according to any one of claims 1 to 7, characterized in that it further comprises a calculation (G) of the positions and signs of the singularities. 9. Method according to claim 8, characterized in that the calculation of the positions of the singularities comprises: a search on the surface (S) of regions where the absolute value of the singularities function is greater than a threshold, and a determination centers of these regions, each position of the centers of gravity thus determined being the coordinates (r) of a singularity. 10. Method according to any one of claims 1 to 9, characterized in that the calculation of the phase (q (Y)) comprises: - a calculation (D, F), from the values and the variation of the intensity, a first estimate of the phase (ç, (F)), - a calculation (H, I, J) of a phase correction term which depends on the existence of phase singularities, and a combination the first estimate of the phase and the phase correction term. 292231 6 - 50 11. Method according to claim 10, characterized in that the calculation (D, F) of the first estimate of the phase (O, ("")) comprises: a calculation of a gradient of the phase, an integration of the gradient of the phase. 12. Method according to claim 11, characterized in that the calculation (D, F) of the gradient of the phase comprises a calculation of a first term which does not depend on the existence of singularities, this calculation comprising for conditions respectively of Neumann or Dirichlet: a respectively cosine or sine transform applied to the variation of the intensity, then, on the result of this transform, an application of a filter in the spatial domain, then on the result of this filtering, an application of a sine or cosine transform respectively, then a division by the values of the intensity. 20 13. Method according to claim 12, characterized in that the calculation (F) of the gradient of the phase comprises a calculation of a second term which depends on the existence of singularities and which is added to the first term, the calculation of the second term. term comprising, for Neumann or Dirichlet conditions respectively: a respectively cosine or sine transform applied to the singularities function, then, on the result of this transform, an application of a filter in the spatial domain, then on the result of this filtering, an application of a sine or cosine transform respectively, then a division by the intensity values. 2922316 -51 14. Method according to any one of claims 11 to 13, characterized in that the integration of the gradient of the phase comprises, for Neumann or Dirichlet conditions respectively: a respectively sine or cosine transform applied on the gradient of the phase, then, on the result of this transform, an application of a filter in the spatial domain, then on the result of this filtering, an application of a respectively cosine or sinus transform. 10 15. Method according to any one of claims 10 to 14, characterized in that the phase correction term comprises a harmonic phase (x (x, y)), calculated by means of a series, over all the singularities, of the value of the singularity multiplied by an integral of a product of gradients of a Green function used for the calculation of the harmonic phase. 16. The method according to claim 15, characterized in that the Green function used for the computation of the harmonic phase is a series of eigenfunctions of the two-dimensional Helmholtz differential equation, preferably for Neumann boundary conditions. or Dirichlet in the surface (S). 17. Method according to any one of claims 10 to 16, characterized in that it comprises a grouping of the singularities in at least one pair of singularities of opposite signs, and for each pair a construction on the surface (S) d a cut line starting from a singularity of the couple until the other singularity of the couple. 30 18. The method of claim 17, characterized in that the cut line of a couple passes through a region of low intensity of the surface. 19. The method of claim 17 or 18, characterized in that for the at least one pair, the signs of the singularities of the pair are opposite, and in that the line of breaking of the torque passes, between the singularities of the couple, by as many positive singularities as negative singularities. 20. A method according to any of claims 17 to 19, characterized in that the phase correction term comprises a hidden phase (x (F)), and in that it comprises a calculation of the hidden phase comprising: each pair, the calculation of a torque function, and a sum, over all the pairs of singularities, of a torque operator applied to the torque functions 21. Method according to claim 20, characterized in that the calculation of the torque function comprises a calculation of a series, on all the singularities of the pair, of the value of the singularity multiplied by: a difference between a first integral, according to a first coordinate, a variation, according to a second coordinate, of a Green function used for the calculation of the hidden phase, and a second integral, according to the second coordinate, of a variation, according to the first coordinate, of the Green function used for the calculation of the hidden phase. 22. Method according to claim 21, characterized in that the Green function used for the calculation of the hidden phase is a series of eigenfunctions of the two-dimensional Helmholtz differential equation, preferably for Neumann or Neumann boundary conditions. Dirichlet in the surface (S). 23. Method according to any one of claims 20 to 22, characterized in that the application of the torque operator to the torque function comprises: a calculation of a spatial gradient of the torque function, integrating the gradient of the torque function along an integration path, adding or removing 27r at the gradient of the torque function each time the integration path intersects the torque cutoff line 24. Device for calculating the phase (cD (F)) of a wave field, characterized in that it comprises means for implementing the steps of the method according to any one of claims 1 to 23. Storage medium storing a computer program for performing the method according to any one of claims 1 to 23.15.