RECONSTRUCTION OF A WAVEFRONT OF A LIGHT BEAM CONTAINING OPTICAL VORTICES

20230296444 · 2023-09-21

    Inventors

    Cpc classification

    International classification

    Abstract

    A method for reconstructing the wavefront of a light beam by analyzing wavefront-gradient data of said light beam, the light beam containing at least one optical vortex, considering the contribution of the optical vortices to the wavefront. The method including providing a phase-gradient map g of the wavefront of the light beam, generating a Laplacian of a vector potential based on the phase gradient map g, the resulting Laplacian of the vector potential map, called “Laplacian map”, exhibiting peaks, the location of each peak corresponding to the location of an optical vortex and the integral of the peak being proportional to a topological charge n of said optical vortex, computing a singular phase map φ.sub.s based on the topological charge n and location of each optical vortex, the singular phase φ.sub.s map being representative of the contribution of the optical vortex.

    Claims

    1. A method for reconstructing the wavefront of a light beam by analyzing wavefront-gradient data of said light beam, the light beam containing at least one optical vortex (V), considering the contribution of the optical vortices (V) to the wavefront, the method comprising: a) providing a phase-gradient map g of the wavefront of the light beam, b) generating a Laplacian of a vector potential based on the phase gradient map g, the resulting Laplacian of the vector potential map (−ΔA.sub.z), called “Laplacian map”, exhibiting peaks (P), the Laplacian map being obtained by applying a curl operator to the phase-gradient map g, the location of each peak (P) corresponding to the location (u) of an optical vortex (V) and the integral of the peak being proportional to a topological charge n of said optical vortex (V), c) applying an image processing method to the Laplacian map to extract the location (u) and the charge n of optical vortices, the image processing method comprising a segmentation of the Laplacian map into a plurality of segments (Ω) to obtain a segmented Laplacian map, segments (Ω) being determined so that different peaks of the Laplacian map belong to different segments; d) for these different segments (Ω) of the segmented Laplacian map (−ΔA.sub.z): i) determining the location (u) of the corresponding peak (P) in the segment (Ω) of the Laplacian map (−ΔA.sub.z) to locate the corresponding optical vortex (V), ii) computing an approximate value of the integral of the corresponding peak (P) over the segment (Ω), and determining the topological charge n of the corresponding vortex (V) based on said approximate value of the integral of said peak (P); e) computing a singular phase map φ.sub.s based on the topological charge n and location (u) of each optical vortex (V) computed in steps i) and ii), the singular phase φ.sub.s map being representative of the contribution of the optical vortex and being computed by: iii) generating a vortex-location map, where map values are equal to the corresponding topological charges n of optical vortices at optical vortices locations (u), and iv) convoluting said vortex-location map with a single +1 spiral phase profile; f) reconstructing the wavefront of the light beam with: computing a regular phase φ.sub.R map of the light beam, the regular phase corresponding to a vortex-free phase computed from the irrotational component, adding the regular phase φ.sub.R to the singular phase φ.sub.s to obtain the said reconstructed wavefront of the light beam).

    2. The method according to claim 1, wherein in step ii) the topological charge n of an optical vortex is obtained by: dividing the approximate value of the integral of the corresponding peak by 2kπ, with k an integer rounding the result of the former operation to the closest integer to get the topological charge n.

    3. The method according to claim 2, wherein the integer k is equal to −1 or +1.

    4. The method according to claim 1, wherein the segmentation at step c) is realized with a numerical tool among the list: the watershed algorithm, the normalized cut algorithm, the SLIC algorithm, or convolutional neuronal networks, and/or a deep learning algorithm.

    5. The method according to claim 1, step i) comprising computing the weighted centroid of each peak (P) in each segment (Ω) corresponding to said peak, determining the location (u) of the said peak based on the weighted centroid.

    6. The method according to claim 1, wherein in step c) the image processing method comprises fitting the peaks (P) with bell-shaped functions, for example with Gaussian functions.

    7. The method according to claim 6, wherein in step ii), the integral of the peak is estimated by computing the area under the corresponding fitting bell-shaped curve.

    8. The method according to claim 1, wherein in step c) the image processing method comprising solving an inverse problem to estimate the distribution of the peaks (P), by using tools such as a compressed sensing algorithm, a Markov Chain Monte Carlo Algorithm (MCMC) methods or neuronal networks.

    9. The method according to claim 1, the image processing method comprising identifying the extrema values of the Laplacian, notably using a threshold to ignore noise, wherein in step i) the location of the peaks is determined based on the location of the extrema.

    10. The method according to claim 1, wherein in step ii) the integral of the peaks is determined from the extrema magnitude only.

    11. The method according to claim 1, wherein the convolution in step e)iv) to obtain the singular phase φ.sub.s is: achieved numerically by adding as many spiral phase profiles as the number of detected optical vortices, each spiral phase profile being centered at the same location and having the same topological charge as the detected optical vortices (V).

    12. The method according to claim 1, wherein in the step f), the regular phase φ.sub.R map is computed by applying a divergence operator to the phase gradient map g, followed by a double spatial integration.

    13. The method according to claim 1, wherein the phase-gradient g is decomposed according to Helmholtz decomposition as follows: g=∇φr+∇×A where φr is the regular phase, the gradient ∇φr corresponding to the irrotational component of g and where A is the potential vector corresponding to the solenoidal component of g.

    14. A computer-readable storage medium comprising instructions which, when executed by a computer, causes the computer to carry out the steps of the method according to claim 1.

    15. A device for reconstructing the wavefront of a light beam by analyzing wavefront-gradient data of the said light beam, the light beam, the light beam containing at least one optical vortex (V), in order to estimate the contribution of the optical vortices (V) to the wavefront, comprising: a wavefront sensor to measure the phase-gradient of the wavefront of the light beam; coupled to the wavefront sensor, means for carrying out the steps of the method of claim 1, some calculating means comprising a memory storage unit, a set of instructions and a processor for executing said set of instructions, these instructions being: a) providing a phase-gradient map g of the wavefront of the light beam; b) generating a Laplacian of a vector potential based on the phase gradient map g, the resulting Laplacian of the vector potential map (−ΔA.sub.z), called “Laplacian map”, exhibiting peaks (P), the Laplacian map being obtained by applying a curl operator to the phase-gradient map g, the location of each peak (P) corresponding to the location (u) of an optical vortex (V) and the integral of the peak being proportional to a topological charge n of said optical vortex (V); c) applying an image processing method to the Laplacian map to extract the location (u) and the charge n of optical vortices, the image processing method comprising a segmentation of the Laplacian map into a plurality of segments (Ω) to obtain a segmented Laplacian map, segments (Ω) being determined so that different segments contain different peaks of the Laplacian map; d) for these different segments (Ω) of the segmented Laplacian map (−ΔA.sub.z): i) determining the location (u) of the corresponding peak (P) in the segment (Ω) of the Laplacian map (−ΔA.sub.z) to locate the corresponding optical vortex (V), ii) computing an approximate value of the integral of the corresponding peak (P) over the segment (Ω), and determining the topological charge n of the corresponding vortex (V) based on the said approximate value of the integral of the said peak (P); e) computing a singular phase map φ.sub.s based on the topological charge n and location (u) of each optical vortex (V) computed in steps i) and ii), the singular phase φ.sub.s map being representative of the contribution of the optical vortex and being computed by: iii) generating a vortex-location map, where map values are equal to the corresponding topological charges n of optical vortices at optical vortices locations (u), and iv) convoluting said vortex-location map with a single +1 spiral phase profile; f) reconstructing the wavefront of the light beam with: computing a regular phase φ.sub.R map of the light beam, the regular phase corresponding to a vortex-free phase computed from the irrotational component, adding the regular phase φ.sub.R to the singular phase φ.sub.s to obtain the said reconstructed wavefront of the light beam.

    16. The device according to claim 15, configured to reconstruct the wavefront of the light beam, by calculating the regular phase φ.sub.R which is added to the singular phase φs.

    17. The device according to claim 15, wherein the calculating means are configured to: pre-compensate the aberrations or compensate a posteriori the aberrations, of the wavefront of the light beam which result of the emission sources and of the properties of the materials the light beam travels across.

    18. The device according to claim 15, wherein the wavefront sensor is chosen among the list: Shack-Hartmann wavefront sensor Multi-wave interferometer such as Quadri-wave lateral shearing interferometer; Pyramid Wavefront sensor.

    19. (canceled)

    Description

    BRIEF DESCRIPTION OF THE DRAWINGS

    [0087] The accompanying drawings, which are incorporated in and constitute a part of the specification, illustrate presently preferred embodiments of the invention, and together with the description, serve to explain the principles of the invention. In the drawings:

    [0088] FIGS. 1 and 2 previously described, illustrate the consequences of neglecting the divergence-free component in adaptive optics and in imaging, respectively,

    [0089] FIG. 3 to 5a are schematic of one embodiment of a method according to the invention,

    [0090] FIG. 5b shows detail IV of FIG. 5b,

    [0091] FIGS. 6a to 7d show examples of reconstruction of the shapes of the wavefront using a method according to the invention, and

    [0092] FIGS. 8a to 8c is a schematic of one embodiment of a system according to the invention.

    [0093] Whenever possible, the same reference numbers will be used throughout the drawings to refer to the same or like parts.

    [0094] In accordance with the invention, and as broadly embodied, a method according to the present invention is provided, aiming at analyzing a wavefront of a light beam 15 containing at least one optical vortex V to estimate the contribution of the optical vortices to the wavefront.

    [0095] First, at step 101, a phase-gradient map g of the wavefront of the light beam 15 is provided. This phase-gradient map g may be generated using a WFS 50.

    [0096] For example, the WFS 50 may be a Shack-Hartmann WFS, a Multi-wave interferometer such as Quadri-wave lateral shearing interferometer, or a Pyramid WFS.

    [0097] Alternatively, the WFS 50 comprising a diffuser may be used. Such a WFS is described in EP 19305382 and in (Berto, Pascal, Hervé Rigneault, and Marc Guillon. “Wavefront sensing with a thin diffuser.” Optics Letters 42.24 (2017): 5117-5120.).

    [0098] Then at step 107, a Laplacian of a vector potential based on the phase gradient map g is generated.

    [0099] To generate the Laplacian map, the method may comprise prior to step 107, a step 103 consisting in performing a Helmholtz decomposition of the phase-gradient map g to split it into an irrotational component and a solenoidal component, as illustrated in FIG. 4.

    [0100] According to Helmholtz decomposition, the phase-gradient g may be expressed as follows: g=∇φ.sub.r+∇×A where φr is the regular phase, the gradient ∇φr corresponding to the irrotational component and A is the potential vector resulting in the solenoidal contribution.

    [0101] In order to characterize the solenoidal contribution, the method may comprise at step 105 applying a curl operator to the phase-gradient map g. The curl of the phase-gradient map g is given by the following equation: ∇× g=∇(∇.Math.A)−ΔA.

    [0102] In the following, a Coulomb gauge is used requiring ∇.Math.A=0. However, the method can be generalized to any other choice of gauge. Using the Coulomb gauge yields ∇×g=−ΔA.

    [0103] Since g is a two-dimensional vector field, in say a x-y plane, the potential A may be written as A=A.sub.ze.sub.z, where e.sub.z is an unit vector perpendicular to the x-y plane.

    [0104] Hence, applying the curl operator to the phase map g yields obtaining the Laplacian map—ΔA.sub.z.

    [0105] Examples of Laplacian maps—ΔA.sub.z are shown in FIGS. 4, 5a and 5b. In the illustrated embodiments, the Laplacian map—ΔA.sub.z exhibits peaks P related to optical vortices. In the illustrated embodiments, the peaks—which are considered—correspond to the useful signal. Peaks due to noise are not taken into account.

    [0106] The location of each peak P corresponds to the location u of an optical vortex.

    [0107] The integral of the peak being proportional to a topological charge n of said optical vortex.

    [0108] After determining the Laplacian map—ΔA.sub.z, the next steps 111 and 113 consist in calculating the location of the peaks P on the Laplacian map—ΔA.sub.z and in computing an approximate value of the integral of each peak in the Laplacian map, respectively.

    [0109] In step 113, one may use the location of the peaks P on the Laplacian map—ΔA.sub.z to locate the corresponding optical vortex.

    [0110] In step 115, based on Stokes theorem, the approximate value of the integral of each peak in the Laplacian map makes it possible to evaluate each topological charge n of each corresponding optical vortex.

    [0111] On the one hand, the integral of g over a close contour 1 surrounding the vortex is custom-character.Math.dl=2πn.

    [0112] On the other hand, according to Stokes theorem, the integral of gcustom-character.Math.dl is equal to the surface integral of its curl ∫∫∇×g.Math.dS, S being the surface enclosed by the close contour 1.

    [0113] As mentioned above, the curl of the phase gradient map ∇×g is equal to the Laplacian map. Therefore, the value of the integral of each peak in the Laplacian map ∫∫−ΔA.sub.z.Math.dS advantageously evaluates the topological charge n of the corresponding optical vortex. More precisely, and as shown in FIG. 6b, when calculating the approximate value of the integral, the topological charge n of an optical vortex may be obtained by dividing the approximate value of the integral of the corresponding peak by 2kπ, with k an integer being, in general, equal to −1 or +1. The sign of the integer k provides information relative to the direction of phase circulation, clockwise or counterclockwise, of the vortex.

    [0114] It is worth mentioning that in the method according to the invention, the estimation of the vortex charge through the integration of −ΔA.sub.z/2kπ is not or very slightly affected by the limited resolution of the wavefront sensor or by the diverging phase gradient magnitude in the vicinity of the vortex. This is mainly related to the use of the Stokes theorem for computing the topological charges. Since at large-enough distances from the vortex location, the phase-gradient estimate is accurate, the contour integration at this distance yields the correct charge value. Therefore, in virtue of the equivalence ensured by the Stokes theorem, surface integration of −ΔA.sub.z/(2kπ) over the surface enclosed by this close contour yields the same accurate charge value.

    [0115] In the embodiment of FIG. 5a, the surface integral of −ΔA.sub.z/(2kπ) is rounded to the closest integer to get the topological charge n of the corresponding optical vortex.

    [0116] In some applications, the Laplacian can be affected by noise and/or by filtering caused by the WFS, notably by a transfer function of the WFS. In this case, it is preferable to apply an image processing method to the Laplacian map to extract the location and the charge of optical vortices from the computed Laplacian map −ΔA.sub.Z.

    [0117] In this case, the method may comprise an additional step 109 preceding step 111, consisting in applying the image processing method to the Laplacian map −ΔA.sub.z.

    [0118] In the embodiments of FIGS. 4, 5a and 5b, the image processing method comprises a segmentation of the Laplacian map −ΔA.sub.z into a plurality of segments 2 to obtain a segmented Laplacian map.

    [0119] Advantageously, segments 2 are determined so that different peaks P of the Laplacian map belong to different segments.

    [0120] In the illustrated embodiments of FIGS. 4, 5a and 5b, the segmentation is performed using a watershed algorithm.

    [0121] Other segmentation algorithms may be used, for example, normalized cut algorithm, SLIC algorithm, convolutional neural networks and deep learning algorithm.

    [0122] Once the Laplacian map −ΔA.sub.z is segmented, peaks locations u are then calculated in step 111.

    [0123] For this purpose, step 111 comprises first computing the weighted centroid of each peak P in each segment Ω corresponding to said peak P, and second determining the location u of the said peak P based on the weighted centroid, as illustrated in FIGS. 4 and 5a.

    [0124] As shown in FIGS. 4, 5a and 5b, in step 113, the approximate value of the integral of each peak P is determined by computing the surface integral of the Laplacian −ΔA.sub.z over the corresponding segment.

    [0125] The topological charge n of the corresponding vortex V may then be determined in a same way as described above.

    [0126] At step 115, the singular phase map φ.sub.s is calculated based on the topological charge n and location u of each optical vortex V computed in steps 111 and 113.

    [0127] As embodied in FIGS. 4 and 5a, the singular phase map φ.sub.s may be computed by: [0128] iii) generating a vortex-location map, where map values are equal to the corresponding topological charges of optical vortices at optical vortices locations, and [0129] iv) convoluting said vortex-location map with a single +1 spiral phase profile.

    [0130] FIG. 5a relates to the case where a plurality of peaks P is detected.

    [0131] In this case, the singular phase φ.sub.s is advantageously obtained by adding as many spiral phase profiles as the number of detected optical vortices, each spiral phase profile being centered at the same location and having the same topological charge as the corresponding detected optical vortex.

    [0132] The method may comprise as shown in FIG. 4, a step 121 consisting in computing the regular phase φ.sub.r map of the light beam. As mentioned above, the regular phase corresponds to the vortex-free phase computed from the irrotational component.

    [0133] To that end, a divergence operator may be applied to the phase gradient map g at step 117, followed by a double spatial integration at step 119. The regular phase is then obtained at step 121.

    [0134] A suitable method to compute the regular phase may be found in (L. Huang, M. Idir, C. Zuo, K. Kaznatcheev, L. Zhou, and A. Asundi, Opt. Lasers Eng. 64, 1 (2015).)

    [0135] At step 123, the regular phase may be added to the singular phase φ.sub.s to reconstruct the wavefront of the light beam with adding the regular phase φ.sub.R, as embodied in FIG. 4.

    [0136] FIGS. 6A and 6B illustrate the possibility to rebuild optical vortices of charges larger than 1. The first line represents the Laplacian map and the second line the wavefront. In the example illustrated herein, the value of the topological charge n varies from − 10 to +10. This is of typical interest for information multiplexing in telecommunications through air or optical fibers. One specific difficulty in this case is the increasing gradient at the vortex location, leading to an increased filtering by the optical transfer function of the WFS. As a consequence, the −ΔA.sub.z maps exhibit peaks P whose width increases with the charge of the vortex n, as shown in line 1.

    [0137] As can be seen in FIGS. 6B, despite those difficulties, the integration of −ΔA.sub.z/(2π) over segmented regions yields the topological charge n with a 3% accuracy. This slight discrepancy can be attributed to the calibration uncertainties of the WFS 50 and of the spatial light modulator 40 used for generating the vortices, and can be dealt with by rounding the integrated value to the closest integer.

    [0138] FIGS. 7A to 7D illustrate the possibility to retrieve the phase of complex random wavefields which is of special interest for imaging in depth through or inside scattering media like biological tissues. FIG. 7A displays a random wavefield exhibiting a high density of optical vortices of charge +1 and −1. These vortices exhibit elliptical phase-increase and intensity profiles along the azimuthal coordinate which may alter the ability to detect them.

    [0139] Furthermore, as can be seen in FIG. 7A, the separation distances between the vortices may be much smaller than a speckle grain size, especially when close to creation or annihilation events of pairs of vortices. Despite these specific additional difficulties, FIG. 7B shows that such a complex wavefield can be efficiently rebuilt using the Laplacian map displayed in FIG. 7D. The difference between the rebuilt and the actual phase profiles is shown in FIG. 7C demonstrating the almost perfect identification of optical vortices.

    [0140] FIGS. 8a to 8c displays a further example of reconstructing a phase relying on the method of the present invention. In the embodied example, a light beam 10 generated by a laser 5 is spatially filtered through a single mode fiber (SMF) 20. The light beam 10 is then collimated by a lens 30 and polarized by a polarizing element 33. Here, a Polarizing beam splitter (PBS) 33 is used. The resulting spatially filtered and polarized laser beam 13 is then sent onto a phase-only spatial light modulator (SLM) 40.

    [0141] A phase pattern 60, wrapped between 0 and 2π is then addressed to the SLM 40. As illustrated, the phase pattern 60 exhibits both smooth local lenses, for the eyes and the face contour for instance, as well as left- and right-handed optical vortices at the tips of the wiggling mustache. The resulting beam 15 contains thus a plurality of vortices V. The beam 15 is then imaged onto a WFS 50 with a 4-f Galileo telescope 53, 55.

    [0142] FIG. 8b shows the phase rebuilt when neglecting the vortices. As illustrated, the rebuilt phase fails to reproduce a large part of the features of the phase pattern 60. This demonstrates the significant loss of information when vortices contribution is neglected, which, depending on the considered application, can lead to undesired effects and computation errors.

    [0143] By contrast, in FIG. 8c the phase profile 60 addressed onto the SLM is rebuilt thanks to a Helmholtz decomposition of the phase gradient map recorded by the WFS, in a same way as described above. As shown, taking the divergence-free component into account allows the proper reconstruction of spiral phase profiles drawing the mustaches of the profile 60.

    [0144] The invention is not limited to the described embodiments, and various variations and modifications may be made without departing from its scope.

    [0145] For example, at step 109 the image processing method may comprise fitting the peaks P with bell-shaped functions, for example with Gaussian functions.

    [0146] In this case, at step 111 one may compute the center of the bell-shaped functions and determine the location of the corresponding peaks P based on the center of the bell-shaped functions.

    [0147] In step ii), the integral of the peak P may be then estimated by computing the area under the corresponding fitting bell-shaped curve.

    [0148] In a variant, the image processing method may comprise solving an inverse problem to estimate the distribution of the peaks, by using tools such as a compressed sensing algorithm, a Markov Chain Monte Carlo Algorithm (MCMC) methods or neural networks.

    [0149] According to a further embodiment, the image processing method may comprise identifying the extrema values of the Laplacian, notably using a threshold to ignore noise, wherein in step i) the location of the peaks is determined based on the location of the extrema.

    [0150] In this case, in step ii), the integral of the peaks may be determined from the extrema magnitude only.

    Additional Example of Invention

    [0151] In an embodiment, a method according to the present invention relies on the use of the Helmholtz's theorem, according to which the vector field g can be split into an irrotational (curl-free) component and a solenoidal (divergence-free) component.

    [0152] All current integration techniques for WFS basically consists in computing ∇.Math.g, so implicitly cancelling the solenoidal component of the vector field. When processing expectedly-smooth wavefront, the solenoidal term typically results from noise in low-intensity regions. However, ignoring optical vortices in wavefronts detected in the far field for instance is usually not possible. Neglecting a single optical vortex is then indeed equivalent to adding a complementary spiral phase mask which has been described as performing a two-dimensional Hilbert transform of the intensity pattern. For complex wavefields, such a single-spiral transform induces a major change in patterns since resulting in an inversion of intensity contrasts. The Helmholtz decomposition has been applied to many fields in physics, mostly in fluid mechanics. Optical phase vortices appearing in freely propagating beam have the specificity to be quantized and are multiples of 27. Here, using a high resolution WFS described in EP 19305382, it has been demonstrated the possibility to rebuild complex wavefronts containing optical vortices.

    [0153] The experimental system used is shown in FIG. 8a. A phaseonly spatial light modulator (SLM) (Hamamatsu, LCOS-X104 468-01), illuminated by a spatially filtered and polarized laser beam is imaged onto a high-resolution wavefront sensor (WFS).

    [0154] A phase pattern, wrapped between 0 and 21 is then addressed to the SLM 50, (FIG. 8a) exhibiting both smooth local lenses (for the eyes and the face contour for instance) as well as left- and right-handed optical vortices at the tips of the wiggling mustache.

    [0155] The Helmholtz decomposition of the phase-gradient vector-field measured by the WFS (FIG. 8c) according to:


    g=∇φ.sub.r+∇×A

    [0156] then allows splitting apart the contribution of a regular phase φ.sub.r and of a vector potential A. The sought-for phase profile φ, corresponding to the integral of g, is then φ=φ.sub.r+φ.sub.s, where the singular phase contribution φ.sub.s, defined over a multiply-connected domain, satisfies ∇φ.sub.p=∇×A. Noteworthy, the unicity of the Helmholtz decomposition both requires setting boundary conditions and taking into account the contribution of a so-called additional “harmonic” (or translation) term h, which is both curl-free and divergence-free. In practice, the translation term h, accounting for tilted wavefronts, may be included in the curl-free component ∇φ.sub.r, which it is achieved here thanks to symetrizations operations on the phase gradient map. Unicity of the solution is further ensured by implicit periodic boundary conditions applied by achieving derivation and integration computations through discrete Fourier transforms. The Helmholtz decomposition can typically be achieved by computing ∇.Math.g=∇φ.sub.r and ∇× g=∇(∇.Math.A)−ΔA. A double spatial integration of ∇.Math.g allows retrieving T.sub.r. Determining the potential vector A requires fixing gauge. The Coulomb gauge ∇.Math.A=0 is chosen here for obvious convenience.

    [0157] Since g is a two-dimensional vector field (in say x, y plane), A can be written as A=A.sub.ze.sub.z. Then, Helmholtz decomposition can be achieved according to the scheme shown in FIG. 4.

    [0158] Noticing that gx+igy=(∂x+i∂y)(φ.sub.r, +iA.sub.z), the Helmholtz decomposition may be conveniently achieved thanks to a single computation step by projecting vectors onto the circular unit vector σ+=e.sub.x+ie.sub.y:

    [00003] φ r - iA z = F - 1 { F [ g . σ + ] ib σ + }

    [0159] where b stands for the two-dimensional coordinate vector in the reciprocal Fourier space. The regular phase component φ.sub.r can then be directly recovered in a similar was as usually performed (see FIG. 4). The divergence-free component may require further processing steps to get the phase with spirals φ.sub.s from the potential-vector component A.sub.z.

    [0160] Considering an optical vortex of topological charge n, by definition of the topological charge, the integral over a close contour surrounding the singularity is: custom-character.Math.dl=2πn.

    [0161] Applying the Stokes' theorem then yields ∫.sub.s−ΔA.sub.ZA.Math.dS=2πn

    [0162] Reducing the contour length to zero, it then appears that −ΔA.sub.z/2πn=δ(r.sub.⊥) is the Green function of the Laplace equation (see FIG. 4) making it easy to identify. Nevertheless, rebuilding the corresponding sought-for singular phase component φs is more difficult and no systematic procedure has ever been proposed. In principle, φs could be obtained by simply convoluting −ΔA.sub.z/2πn by a single +1 optical vortex.

    [0163] However, in practice, the −ΔA.sub.z map may not be a perfect Dirac distribution, or a single-pixeled non-zero data map, but may be affected by noise and more importantly, by filtering by the optical transfer function.

    [0164] As detailed in (Berto, Pascal, Hervé Rigneault, and Marc Guillon. “Wavefront sensing with a thin diffuser.” Optics Letters 42.24 (2017): 5117-5120), the optical transfer function may be dominated by the non-overlapping condition of WFS, which limits the maximum local value of the eigenvalues of the Jacobi matrix of g (i.e. the Hessian matrix of φ). As a first consequence, the curvatures of the phase component φ.sub.r (i.e. its second derivatives) are filtered. For an optical vortex, the magnitude of g diverges as 1/r at the vortex center, and so the Hessian coefficients ∂.sub.xg.sub.y and ∂.sub.yg.sub.x. The vector potential A.sub.z is thus similarly filtered. An image processing step may thus be required to extract the location and the charge of optical vortices from the computed ΔA.sub.z map.

    [0165] First, noise is removed by applying a Gaussian filter to ΔA.sub.z (FIG. 4), whose width is set according to a rough minimal estimate of the WFS resolution (10 camera-pixels here). The resulting image is then segmented using a watershed operation. Integration and weighted-centroid computation over each segmented regions yields the charge of each vortex. Finally, a Dirac-like vortex map is rebuilt according to the location and the charge of vortices, and convoluted by a +1 spiral phase mask to yield φs (see FIG. 4). The complete phase reconstruction φ=1φr+φs is then computed and wrapped between 0 and 2π (FIG. 4).

    [0166] To demonstrate the possibility to rebuild optical vortices of charge larger than 1, a phase spirals from −10 to +10 (FIGS. 6a and 6b) is addressed. One specific difficulty in this case is the increasing gradient at the vortex location, leading to an increased filtering by the optical transfer function of the WFS. As a first consequence, the ΔA z maps exhibit peaks whose width increases with the charge of the vortex n (FIG. 6a). Nevertheless, integration of −ΔA.sub.z/(2π) over segmented regions yields n with a 3% accuracy which it may attributed to the calibration uncertainties of the SLM and WFS (FIG. 6b). Rounding the integrated value to the closest integer allowed proper reconstruction of the optical vortices (FIG. 3a). Differences of rebuilt phase profiles with perfect ones addressed to the SLM is due to the noisy contribution of the regular phase φr. In principle, φr should be perfectly flat but imperfections of the SLM and in the imaging system leads to wavefront distortions.

    [0167] The possibility to retrieve the phase of complex random wavefields is demonstrated. Random wavefields exhibit a high density of optical vortices of charge +1 and −1. These vortices exhibit elliptical phase and intensity profiles along the azimuthal coordinate which may alter the ability to detect them.

    [0168] Furthermore, their separation distances may be much smaller than the speckle grain size, especially when close to creation or anihilation events of pairs of vortices. Despite these specific additional difficulties, such a complex wavefields addressed to the SLM (FIG. 7a) could be efficiently rebuilt (FIG. 7b). Some of the 0-2π equiphase lines of the input phase are materialized. These equiphase lines are easy to identify because of the abrupt black to-white drop. The difference between the rebuilt and the input phase profiles is shown in FIG. 7c demonstrating the almost perfect identification of optical vortices. Differences mostly appear on the φr contribution on the edges of the SLM, where the SLM reliability degrades. The dense experimental map of vortices distribution and charges is shown in FIG. 7d.