Phase imaging by spatial wavefront sampling

Phase imaging techniques extract the optical path-length information of a scene, whereas wavefront sensors provide the shape of an optical wavefront. Since these two applications have different technical requirements, they have developed their own specific technology. Here we show how to perform phase imaging combining wavefront sampling using a reconfigurable spatial light modulator with a beam position detector. The result is a time-multiplexed detection scheme, capable of being shortened considerably by compressive sensing. This robust reference-less method does not require the phase unwrapping algorithms demanded by conventional interferometry, and its lenslet-free nature removes tradeoffs usually found in Shack-Hartmann sensors.


Introduction
Even though the physical nature of light has been fully understood for more than a century, there are still no available detectors capable of directly imaging both the amplitude and phase information of a wavefront (Fig. 1a). Information about those two quantities is of capital interest when trying to perform biomedical imaging [1,2], aberration measurement and correction in visual optics [3,4], and threedimensional imaging [5], among other applications. The limitation in performing optical measurements ultimately arises from the extremely fast oscillations of the optical fields, which current electronics are unable to resolve. As detectors only capture light irradiance, several approaches have been proposed over time to tackle the phase problem. Historically, Gabor suggested in 1949 the first quantitative technique, which used interferometric information to recover the complex optical field [6], stablishing the basis of modern holography and paving the way for applications such as quantitative phase microscopy [7]. In parallel, phase information plays a fundamental role in adaptive optics, a technique that intends to measure and correct optical aberrations in real time [8]. Although adaptive optics was initially conceived for circumventing atmospheric turbulences in Astronomy, its simple operation principle has found applications in other areas, such as visual optics [3,4], microscopy [9,10], and biomedical imaging [11,12].
Two main groups of techniques have emerged to tackle the problem of recovering the complex amplitude of the light field. The first one includes interferometric approaches, which measure the interference between light coming from the object and a reference beam (Fig. 1b). Although they are extremely powerful for conducting precise phase measurements, their high sensitivity to environmental perturbations (such as mechanical vibrations and changes in temperature) and the need of a reference beam (not always attainable) hinder the implementation of portable and compact interferometric imaging systems in many applications. As an alternative, a second group of techniques has emerged, whose objective is to recover the same information without the need of a reference beam. These noninterferometric approaches rely on several assumptions about the object beam, and use mathematical algorithms to infer the wavefront information [13]. Examples of this kind of techniques are Fourier ptychography [14], coherent diffractive imaging [15,16], phase imaging based on the transport-of-intensity-equation [17], and phase imaging with randomized illumination [18]. By eliminating the need of a reference beam, simpler and more robust devices can be designed. Furthermore, computational approaches also enable to extend the physical capabilities of imaging systems, providing increased field of view [19], optical sectioning [20] or superresolution [21,22]. However, the recovery algorithms used in those techniques usually entail high post-processing times or multiple data acquisitions to recover one image.
Among the non-interferometric techniques, Shack-Hartmann (SH) wavefront sensors [23] are currently the most employed for measuring optical aberrations. SH wavefront sensors combine a lenslet array with a pixelated detector, like a CCD or a CMOS camera. By placing the detector at the focal plane of the lenslet array, the positions of the foci generated by each lenslet can be linked to the phase of the wavefront in each region (Fig. 1c). In other words, each lenslet samples the phase information of a region of the wavefront. This method provides a very compact system, easy to use and align, and with no complex postprocessing stages. For these reasons, SH sensors are used in manifold scientific fields [4,[24][25][26][27]. However, the number of lenslets, their focal length, and their diameter limit the spatial resolution, dynamic range, and sensibility of SH wavefront sensing. Although multiple solutions have been proposed to manage the several tradeoffs between the above magnitudes [28][29][30][31][32][33][34][35][36], manufacturing processes still place a boundary on the size and curvature of available lenslets, constraining the attainable spatial resolution of commercial SH sensors. Quadriwave lateral shearing interferometry [37], also known as a modified Hartmann mask technique, makes it possible to increase the SH spatial resolution by a factor four [38]. Alternatively, phase imaging can be performed using pyramid wavefront sensors, which employs a four-sided refractive pyramid to split the Fourier plane in four parts, each one generating a phase gradient image on a pixelated detector [39]. Further developments replace the pyramid by a quadri-foil lens [40] or a liquid crystal display [41]. Despite the benefits of this approach in terms of spatial resolution, the illumination numerical aperture determines the detectable dynamic range, so in practice only relatively smooth phase gradients can be precisely recorded [42].
Here, we perform phase imaging using an alternative noninterferometric approach to measure the complex amplitude of a wavefront. We overcome the inherent limitations in spatial resolution, optical efficiency and dynamic range that are found in SH wavefront sensing. We sample the wavefront by using a set of binary amplitude masks generated by a high-speed spatial light modulator. A single focusing lens forms a time-dependent light distribution on its focal plane, where an irradiance detector is placed. Measuring the changes of the total irradiance and the centroid position of that distribution, both the amplitude and phase of the wavefront can be recovered (Fig. 1d). One advantage of this time-multiplexed detection scheme is that it can be applied using compressive sensing (CS), allowing us to go beyond the Shannon-Nyquist condition, with the subsequent reduction in acquisition time [43,44]. Unlike previous phase retrieval methods based on SH wavefront sensors or wavefront modulation [45][46][47], our approach is lenslet-free and does not rely on any kind of iterative algorithms. Moreover, the use of phase unwrapping algorithms is not required, such as those employed in interferometric approaches due to their limited phase unambiguity range [48]. As a proof of concept, we perform aberration sensing, as well as phase imaging of a low-contrast sample with variable optical thickness. The corresponding results are compared, respectively, with Shack-Hartmann wavefront sensing and phase-shifting interferometry.

Operation principle
The basic idea of our technique is to use the relationship between the amplitude and phase of a wavefront and the statistical moments of its Fourier irradiance distribution. This idea is also at the heart of current SH sensors, as they sample a wavefront by means of an array of lenslets and measure the centroid position (the first-order statistical moment) of the Fourier distribution created by each lenslet, allowing one to derive the local slope of the wavefront. By placing a pixelated detector (typically a CCD) at the focal plane of the lenslet array, SH sensors fully map the local slopes of the wavefront. After that measurement, numerical integration provides the phase of the wavefront at each lenslet position [49] (see Fig. 2a).
Several tradeoffs arise from the SH configuration. First, the spatial resolution of the measurement is determined by the total number of lenslets and their size. In order to increase the spatial resolution, arrays with higher density of lenslets can be built. However, reducing the size of each lenslet also decreases the amount of light at each region of the detector. If the signal-to-noise ratio of the measurement is low, the accuracy of the centroid position detection is greatly reduced, which hinders a reliable wavefront reconstruction. Even if the detector can work in low-light-level scenarios, building large arrays of lenslets with a diameter below a hundred of microns and an accurate curvature is technologically challenging [50]. Moreover, centroid displacement and wavefront slope are related via the lenslet focal length, so, for a given slope, the higher the focal length is, the larger the centroid displacements that are measured. This brings about a second tradeoff, since a strongly aberrated wavefront may produce centroid displacements larger than the size of the detection area allocated to every lenslet on the sensor. The result is a crosstalk between nearby spots that produces errors in the reconstructed phase, limiting the attainable dynamic range of the sensor. This problem can be circumvented (without sacrificing sensitivity) by a number of techniques that track the real spot location or infer it by using computational techniques [29]. Instead of trying to expand the dynamic range of a sensor by increasing its hardware and/or software complexity, an alternative approach is to implement a reconfigurable (adaptive) SH sensor that includes an array of diffractive lenslets programmed onto a spatial light modulator [28,33]. The lenslet characteristics can be then chosen to match the requirements of a specific application optimally, albeit the tradeoff between sensitivity and dynamic range still remains. Another aspect of SH sensors that is usually not considered is the fact that the total amount of light arriving at each detector region provides a measurement of the light power at the corresponding lenslet position. Using that information, one can map the irradiance of the light coming from an object. However, due to the relatively poor spatial resolution of SH sensors (usually around a thousand of lenslets), they barely can compete with interferometric systems to measure the complex amplitude of an object with high spatial frequency content.
Another approach for measuring the local slope of a wavefront is the use of a small scanning aperture and a single lens rather than an array of lenslets (Fig. 2b). For each position of the aperture, only light from a region of the wavefront will be focused by the lens. This light will generate a focal spot on the Fourier plane of the lens, where the detector is placed. For a small enough aperture, the wavefront can be locally approximated by a plane wave, and it can be demonstrated (see Methods) that the position of the centroid of the distribution is related to the gradient of the phase inside the aperture by the equation where ∇ ⃗ ⃗ = ( ⁄, ⁄) represents the gradient operator in two dimensions, λ is the light wavelength, and is the focal length of the focusing lens. One can see the combination of this small aperture and the single lens as one of the lenslets of a SH array. By scanning the wavefront with the aperture, the phase can be easily recovered by numerical integration. Furthermore, calculating the zero-order moment of the light distribution (i.e., its total irradiance) provides a measurement of the wavefront amplitude in the position of the aperture, leading to the recovery of an amplitude image. Although this approach would solve the spatial resolution problem of SH sensors, given that it is easier to generate small apertures than small lenslets, it still presents a tradeoff between the spatial resolution and the amount of light at each area of the detector plane, unless the scanning is performed with a laser combined with a galvanometric mirror. Our technique resolves the lack of efficiency of the scanning-based approach and makes it possible a compressed sensing process. As can be seen in Fig. 2c, instead of just selecting a small region of the wavefront, a time-variable mask (a sequence of programmed patterns) selects simultaneously several regions of the wavefront. This procedure resembles those used in structured illumination microscopy or in single-pixel imaging [44,51,52]. It must be noted that, as in single-pixel imaging, the technique presented here can work in two different configurations: either the object (or wavefront) under study can be imaged onto the spatial light modulator or the codified patterns can be projected onto the object plane [53][54][55][56]. The key factor that has to be accomplished is that both the mask and the object wavefront overlap in one plane. Even though the irradiance distribution generated by the focusing lens results from the contribution of light coming from different areas of the sampling plane, the amplitude and phase information can be still retrieved, provided that the sequence of sampling patterns is known (see Methods). Now, the amount of light arriving at the detector for every pattern has been greatly increased in comparison to the scanning-aperture scheme. This fact is commonly known as the multiplex or Fellgett's advantage, and has been extensively used to increase the signal-to-noise ratio of optical measurements over the last decades [52,[57][58][59]. In our proposal, the set of masks used are a shifted and rescaled version of the Walsh-Hadamard functions [60], since this election provides several benefits. First, Walsh-Hadamard functions are binary, so they are a good choice when working with binary amplitude modulators, as is the case of digital micromirror devices (DMDs), which reach very high refresh rates (around 22 kHz). Second, the sampling patterns have the same number of absorbing and transmitting (or reflecting) areas, no matter the size of those areas. This fact is important, because once the window size on the modulator has been fixed, increasing or decreasing the spatial resolution of the masks does not change the total amount of light transmitted or reflected by the device. Then, irrespective of the chosen spatial resolution, each measurement always works with the same number of incident photons, unlike the wavefront sensing approaches described above. Furthermore, since in every measurement only a single spatial light distribution is detected, the cross-talk errors that limit the dynamic range in a SH sensor are not present here.
In order to recover the amplitude and phase information of a wavefront, our method requires the measurement of the zero and first order statistical moments of the spatial distribution of light, i.e., its total A wavefront coming from an object passes through an array of lenslets, which produce a distribution of focal spots on the detector. The position of each spot is linked to the local slope of the wavefront on every lenslet. Numerical integration of the slope data provides the phase of the wavefront. (b) Raster scanning wavefront sensing. A small amplitude aperture is moved over the wavefront plane and, for each consecutive scanning position, the Fourier distribution of the emerging light is recorded. The total irradiance of that distribution at each aperture position can be used to recover the amplitude image of the object. Additionally, the centroid relative location of each Fourier distribution provides the local slope of the wavefront at each scanning position, which allows one to obtain the phase image of the object. (c) Spatial wavefront sampling. Instead of using a small scanning aperture, the wavefront coming from the object is sampled by a set of amplitude masks (reconstruction basis). Now, from the irradiance of each Fourier distribution generated on the detector plane, one can measure the mathematical projection of the object amplitude into the reconstruction basis. By demultiplexing that information, the object amplitude is spatially resolved. In the same way, demultiplexing the data concerning the centroid locations provides the slope map of the wavefront and then, after numerical integration, the object phase image.
irradiance and the position of its centroid. The classical approach to measure those quantities is to place a pixelated sensor in the Fourier plane of the lens to acquire a digital image, and then calculate them computationally. However, other experimental approaches can be explored. Here we have used a lateral position detector instead of a camera. The benefit of using this detector lies in its capability to provide the information about the power of the beam and the position of its centroid at speeds in the order of several kHz, without need of computational procedures. Despite digital cameras offer better sensitivity and accuracy in the detection of the centroid, we have opted for a simpler design to demonstrate the feasibility of our approach.

A. Experimental verification
To demonstrate our method, we present the experimental device shown in Fig. 3a. Light coming from a laser (Oxxius slim-532) emitting at 532 nm is collimated with a lens (L1) and impinges onto a DMD (DLP Discovery 4100 V-7000 from ViALUX). By using a 4-f system, formed by lenses L2 and L3, light is projected onto an object, which modifies the wavefront phase. After going through it, light passes through a condensing lens (CL), with a focal length of 150 mm. In its Fourier plane, the lateral position detector (Thorlabs PDP90A) measures the irradiance of the beam and the position of its centroid. To retrieve the amplitude and phase information of the object, a full set of Hadamard patterns must be projected. The size of this set depends on the pixel size of the image one wants to obtain. Then, for an × image, 2 patterns must be sent (see Methods). This sequential acquisition represents a limitation of our technique. Whereas a SH sensor captures all the information in a single shot, our technique relies on the sequential sampling of the wavefront to work.
As a first example, we show the results for a photoresist plate representing a typical coma aberration ( Fig. 3c-d). In this case, the spatial resolution was set to 64 × 64 pixels, with a pixel pitch of 82.08 µm. For this size, the total acquisition time was roughly one second. This is not caused by the refresh rate of our DMD (that would need roughly 200 ms to project 64 2 patterns), but by the limited bandwidth of our detector (Fig. 3b). Faster lateral position sensing detectors, which are commercially available, would speed up the acquisition process. A comparison of this result with those obtained by other wavefront sensing techniques (using a SH sensor and through an interferometric measurement) can be found in Supplement 1.

B. Comparison with Shack-Hartmann wavefront sensing
To test the capabilities of the proposed technique, we have compared it with a commercial SH wavefront sensor (Thorlabs WFS150-5C) in several scenarios. As a first object, we employed a spherical lens (Fig. 4a) with a focal length of 500 mm aligned with the optical axis. This simple choice facilitates to compare both methods quantitatively, since we expect to recover the phase of a spherical wavefront. The object plane was imaged onto the lenslet array plane of the SH wavefront sensor by using a 4-f system. The focal spots produced by the lenslet array can be seen in Fig. 4b. Using the displacements of each focal spot respect to a reference position (previously determined by a calibration measurement), the 3D view of the wavefront can be recovered by direct integration. To ease the visualization, we also present a wrapped view of the wavefront phase (Fig. 4c).
Using the optical setup shown in Fig. 3a, we reconstructed the wavefront phase by using spatial wavefront sampling. A small region of the full measured 64×64 map of centroid displacements is shown in Fig.  4d. The modulo-2π phase reconstructed from those displacements can be observed in Fig. 4e. The total phase change over the aperture is 9.67λ, a value that differs in 0.06λ from that obtained using the SH sensor.
As a second experiment, we moved the lens away from the optical axis, thus producing a larger phase gradient. Given the technical specifications of our SH sensor, the maximum measurable total phase change over the aperture of the sensor is around 100λ. When the phase gradient produced by the off-axis lens produces a higher phase variation, the focal spots move so much in one direction that cross-talk appears and the sensor is not able to provide reliable results.
In the top row of Fig. 5, we show the result obtained by using the SH sensor at the limit of its dynamic range (here we show a wrapped view of the phase to ease visualization). This limit can be considerably overcome using the setup shown in Fig. 3a. As can be seen in the bottom row of Fig. 5, spatial wavefront sampling allows us to measure a total gradient phase of 217. The above improvement in the dynamic range is due to the relative increase in the size of the detector area attained by our technique. In a SH sensor, the detector is divided into a number of regions equal to the number of lenslets in the array. Each region, containing a small number of pixels, is used to map the position of each focal spot and sets a limit to the maximum spot displacement. As was mentioned before, walking off beyond that limit makes several focal spots to share the same region of the detector, leading to significant errors in the reconstructed phase. Since in our technique we sequentially detect the light distribution as a whole, the full size of the sensor can be used in every consecutive measurement. As a result, the cross-talk penalty typically found in a SH sensor is not present anymore. In the same way as in the design stage of SH sensors, fine tuning of the focal length of the condenser lens used and the size of the detector will provide bigger or smaller dynamic ranges and sensitivities in the phase detection.

C. Use of Compressive Sensing
Given the similarities between our technique and single-pixel imaging, some advanced undersampling approaches can be used to improve both the acquisition speed and the light efficiency of our method.
CS provides a recovery framework for signals sampled below the limit imposed by the Shannon-Nyquist theorem [43]. For a sequential signal sampling, this implies a reduction in measurement time. Once the signal has been sampled, an optimization algorithm provides an estimation of the original signal. This estimation is usually based on sparsity constraints and provides a high-fidelity reconstruction (see Methods).
In our experiments, we have used standard CS algorithms ( 1 -magic package [61]) to recover the coma aberration presented in Fig. 3. The results are shown in Fig. 6.
As can be seen, for relatively smooth aberrations, estimations with fidelity higher than 90% can be recovered by using 10% of the total number of measurements established by the Shannon-Nyquist criterion. Fig. 6 includes two insets that show individual reconstructions for compression ratios corresponding to correlation coefficients over 0.9 and 0.99, respectively (points A and B in the fidelity plot). As in our commercial Shack-Hartmann sensor, the recovered wavefronts are expressed in the Zernike basis and high-resolution images are generated from that expansion (see Supplement 1 for details). By taking this into account, high-resolution acquisitions in 20 ms could be achieved operating at the full frame rate of the DMD (~20 kHz), allowing the capability to perform real time aberrometry. This could be very useful in ophtalmic scenarios, where the patient needs to stand still while the eye aberrations are measured. Moreover, even though half of the source light is not used when operating with Hadamard functions, the reduced number of projections makes the system more efficient in global terms of photon usage that a raster scan approach (even when a galvanometric mirror is employed). This could be convenient in biological scenarios, where photodamage thresholds limit the amount of light that can be used at each region of the sample. By using this CSbased procedure, one would be able to increase the measured signal level without causing photodamage when compared with the raster scan approach.

D. Complex amplitude retrieval and comparison with phase-shifting interferometry
Our technique offers the possibility of reconstructing not only the phase distribution of a wavefront coming from an object but also its amplitude. For the set of illumination patterns, the light power measured by the detector provides the projections of the wavefront amplitude into the basis of Hadamard functions. Measuring all the successive projections, the object amplitude image can be recovered offline, following the operation principle of single-pixel imaging [44] (see Methods). A simple example with an object composed of an amplitude mask attached to the lens used in the previous section can be seen in Supplement 1.
In principle, the complex amplitude of a wavefront can be retrieved by means of a SH sensor, but with a relatively poor spatial resolution (given by the lenslet diameter and the lens fill factor). As a consequence, spatial features smaller than a few hundreds of microns are lost in the   Fig. 4a is placed in an off-axis position. The maximum phase gradient measurable by the commercial sensor is 100. A zoomed region, marked in red, is shown in the right part of the figure. In the bottom row, we show the results for the same object obtained with the system shown in Fig. 3a. In this case, the lens has been displaced a bigger distance from the optical axis. The phase gradient measured is 217. A zoomed region, marked in blue, is shown in the right part of the figure. recovery process. Even though high resolution SH sensors are being developed and can be used to perform phase imaging, the tradeoffs between dynamic range and sensitivity are still present [62]. In our technique, the limit in the spatial resolution is ultimately given by the modulator pixel size (typically 10 microns) and the magnification of the projecting system. In our setup, this supposes an increment in the resolution of one order of magnitude respect to a typical SH sensor.
To demonstrate this improvement, we use as a sample a thin layer of a photoresist material. We placed the photoresist over a transparent plate, creating an object with different regions with and without material, as can be seen in Fig. 7a. Due to the difference of refractive index between the photoresist and the air, an impinging wavefront acquires a spatial phase distribution when it passes through the object. By using the setup shown in Fig. 3a, the amplitude and phase information of the sample is recovered (Fig. 7b). A 3D visualization of the phase information is shown in Fig. 7c. The images have a spatial resolution of 128 × 128 pixels, with a pixel pitch of 41.04 µm. Since both the photoresist and the transparent plate present a similar absorption in the green region of the spectra, the amplitude image presents low quality. However, fine details of the sample are recovered in the phase image, such as the small holes (with a diameter of around 80 µm) present in the photoresist stripes, produced by air bubbles in the manufacturing process. These small spatial features cannot be retrieved by our SH sensor, as the diameter of each lenslet is 150 µm. Additionally, due to the low number of lenslets (39 × 31 = 1209) the SH sensor would only provide an image clearly "pixelated".
In order to check the validity of our results, we resort to phaseshifting interferometry. This technique uses a pixelated detector, which makes it possible to exploit the high spatial resolution of a commercial CCDs (in our experiment, the pixel size is 6.5 µm). The recovered phase is shown in Fig. 7d. Unlike our wavefront sensing technique, phaseshifting interferometry only provides modulo-2 mappings (see the vertical scale in Fig. 6d), requiring the use of unwrapping algorithms to reconstruct a continuous phase like the one shown in Fig. 7b. We focus our attention on the small region between the two air bubbles in one of the photoresist stripes. After zooming in, it can be seen that the phase difference between the top of the strip and the substrate is 22.7 radians. In our measurements with the structured illumination setup, this difference is 23 radians, showing a very good agreement with the previous value. We also performed a thickness measurement of the region of interest with a mechanical profilometer (Dektak 6, Veeco). By using the strip height provided by the profilometer (3.8 µm), the refractive index of the photoresist layer can be estimated from the optical path length, ( = λΔ 2πΔ ⁄ ). In both phase shifting holography and our technique, the index estimation is 1.51, which is in good agreement with the nominal value of the photoresist material.

D. Discussion
We have introduced a wavefront sensor based on patterned illumination produced by an amplitude spatial light modulator. Instead of resorting to an array of lenslets, as in SH wavefront sensing, the spatial information is captured by illuminating the object with binary amplitude masks generated with a DMD. The sensor is a simple position sensitive photodetector. This provides several benefits. First, we eliminate the spatial resolution constraints of traditional SH wavefront sensors. Second, by using the full size of the detector in each measurement, the dynamic range is greatly extended and the cross-talk problem is eliminated. Last, once the window size of the DMD has been fixed, changing the spatial resolution of the masks does not change the total amount of light, thus increasing optical efficiency when compared to SH wavefront sensing. This, in combination with the multiplex advantage, makes us think that the technique is well suited to work at low-light level scenarios, where similar approaches have been already proposed to obtain amplitude information [63]. Additionally, we have exploited the high spatial resolution offered by commercial DMDs to perform phase imaging with a pixel size of tens of microns. To this end, we have imaged a photoresist sample and calculated its local refraction index. These results have been demonstrated to be comparable with those provided by phase shifting interferometry, a well-established interferometric technique.
A distinctive feature of our technique is the fact that the spatial resolution is fixed by the spatial light modulator and the projecting optics, while the other optical elements included in the sensor, the focusing lens and the light detector, set the attainable sensitivity and dynamic range. By calculating several statistical moments of the light distribution at the detector plane (its irradiance and the position of its centroid), the method provides both amplitude and phase information. This enables the technique to produce results that are closer to the conventional phase imaging approaches, but still retaining some of the benefits of the SH approach. These benefits are the reference-less nature of SH wavefront sensing and the fact that the recovery process does not require phase unwrapping or iterative algorithms. In addition, the easy implementation of our system also offers the possibility of being used as an add-on module of a conventional microscope [20,64], thus allowing one to carry out quantitative phase imaging with sub-micrometric spatial resolutions.
Due to the temporal multiplexing nature of the technique, there is a tradeoff between spatial resolution and image acquisition time. Here, we have shown results with resolutions between 64 × 64 and 128 × 128 pixels, which take seconds to be acquired. Using faster, already available detectors, will allow image capturing at video frame rates. However, the greater the resolution, the slower the acquisition will be. This drawback can be considerably alleviated using CS [65] or adaptive approaches [66,67]. The feasibility of applying CS has been shown in Section 3, where a typical ocular aberration has been reconstructed with high compression ratios. In particular, acquisitions times of around 20 ms could be achieved operating at the full frame rate of the DMD, which would open up the door to perform real-time aberrometry with our system. Furthermore, this approach increases the global light efficiency of the method when compared to the raster scan procedure. Last, given the technological challenge of manufacturing both arrays of lenslets and pixelated sensors that work outside the visible spectrum, the technique proposed here is a good candidate to operate in regions such as IR and THz, where similar approaches have been already used to obtain amplitude images [52,68].

A. Phase measurement from centroid position
Let us assume a light distribution, at position = 0, with the form: . If a thin lens, with a focal length , is placed at the position = , the complex amplitude distribution at the position = 2 , using Fresnel propagation, will be 2 ( , ) = 1̃1 ( , ), where ̃1 denotes the Fourier transform of 1 . The irradiance of this light distribution will be given by Our technique is based on the calculus of several statistical moments of 2 . In particular, we will use the energy, , and the centroid position, Δ ⃗ ⃗ = (Δ ; Δ ). Those quantities will be given by: .

(3)
By using Parseval's Theorem, it is easy to prove that = | 2 ( , )| 2 = | 1 ( , )| 2 , and the first measurement provides the energy of the wavefront. Using the Moment's Theorem of Fourier transformations and some algebra (see Supplement 1), it can be proved that the position of the centroid at the measurement plane and the phase of the object at the origin position are related by: For a small square aperture with lateral size , placed at position ( , ), the amplitude of the field can be described as: where K is a constant. Introducing this into Eq. (4), we get:

(6)
If the square aperture is small enough, the gradients of the phase can be considered constant over the integration domain. Then, the centroid positions results: where ∇ ⃗ ⃗ , = ( ⁄; ⁄)| = , = represents the bidimensional gradient evaluated at point = , = . By combining Eq. (5) and Eq.
(2), it is easy to prove that = 2 2 , and thus the centroid of the distribution and the phase of the object are related by:

(8)
It is clear that, for a small region of the wavefront, the centroid of the resulting distribution provides information about the gradient of the phase in that region. If the small square aperture is displaced, one can measure the correspondent centroid positions and then estimate the gradient map of the wavefront. After that, numerical integration provides the wavefront at the original position. If multiple square masks are used at the same time, it is also possible to relate the measurements with the phase at each one of the mask positions. Then, it is possible to understand Hadamard illumination as a particular case of this procedure, and the phase recovery can be performed while illuminating the full scene with a set of Hadamard patterns (see Supplement 1).

B. Object recovery by multiplexing with Hadamard patterns
Every measurement can be mathematically described by the equation = , where the object to be measured is described by the object vector , the measurements made are represented by the vector , and the measurement process is carried away by the sensing operator, represented by the matrix . Once the measurements have been done, one just needs to solve the algebraic problem presented by Eq. (9) to recover the object. The structure of the sensing operator will be determined by both the nature of the object under study and the system used to recover an image. For a traditional camera setup, matrix is just the identity matrix, and each element of is proportional to the energy of a given area of the object. By placing a detector at the image plane, the energies of each zone are measured and the object is recovered. Usually this is done with a sensor array, like a CCD camera, and all the measurements are done at the same time. However, there are scenarios, such as confocal microscopy, where this process is carried away by a bucket detector, usually a photomultiplier tube. In this case, the measurement process is sequential. It is possible to exploit Eq. (9) to increase the capabilities of an optical setup. In our approach, instead of just using the identity matrix as the sensing operator, each row of the matrix contains a Hadamard function. Now, the measurement process consists of making the superposition between the object and each one of the Hadamard functions. After that, the energy of that superposition is measured by the detector. This process can be made by projecting the functions onto the object with a spatial light modulator, as is done in our phase imaging system. When using this approach, each measurement contains light coming from several regions of the object. Doing this increases the amount of light at each measurement, and thus the signal-to-noise ratio. As Hadamard functions conform an orthonormal basis, Eq. (9) can be easily solved.
Once this way of measuring has been introduced, it is easy to see that it is possible to recover an image of different physical parameters, provided that the measurement process can be described by Eq. (9). It can be noted that our centroid measurements can be arranged in matrix form as: Here, Δ ⃗ ⃗ is the vector containing the measured positions of the centroids, and now the object is the spatial gradient distribution, expressed as a vector, ∇ ⃗ ⃗ ⃗⃗⃗⃗⃗⃗ (each element of this vector is a gradient for a given illumination pattern). Again, in our experiments we used the Hadamard functions as rows of the sensing matrix. For each Hadamard function projected onto the object, the energy of the distribution and the position of its centroid are measured with the position sensing detector. Then, Eq. (10) is solved and the gradient of the phase is obtained. After that, numerical integration provides the phase of the wavefront, ( , ).

C. Compressive Sensing
Solving the problem given by Eq. (9) is easy when sampling the object at the Shannon-Nyquist ratio. However, for < measurements the equation system presents infinite solutions. The fundamental idea of CS is to use the fact that most of the signals found in nature are sparse in some basis of functions, i.e. they have a representation in some basis where most of the expansion coefficients are zero, so only a small number of them contain most of the relevant information. This fact can be used to provide a solution for the above underdetermined equation system. We can express the measurements as = = Ψ = Φ , where we have represented the object in some basis of functions ( = Ψ ). Then, CS algorithms provide a solution to the 1 -norm minimization problem: = ‖ ′ ⃗⃗⃗ ‖ 1 such that ′ ⃗⃗⃗ = .

(11)
In order to find a good solution, some premises need to be fulfilled [69]. First, the number of measurements cannot be arbitrarily low. Depending on the sparsity of the object (which depends on the chosen recovery basis), higher or lower number of measurements are needed to obtain good results. Second, the basis pair (measurement and recovery) need to be incoherent. In practice, it is usually easy to find pairs of basis that fulfill this principle. Usually, the measurement basis is chosen for its ease of generation in a SLM, and once one has been picked, the user selects a recovery basis that fulfills the incoherence constraint and where the object is sparse. Last, a sensing strategy needs to be defined. Initially, it was proposed that randomly choosing a subset of elements in the measurement basis was enough to obtain good results [43]. Soon after that, other strategies were defined to reduce the number of measurements maintaining the image quality. For example, it is known that natural scenes tend to have a power spectrum centered around low frequencies, so mixed sampling strategies have been used with very good results [58].
In our experiments, represents the shifted and rescaled Hadamard basis (which elements are either 0 or 1), easily generated on a DMD, and Ψ is chosen to be the Haar wavelet basis, where most images tend to be sparse. The measurements have been performed following the mixed approach: fixed the number of samples, , we chose the lowest frequency /2 Hadamard patterns, and then we select /2 random Hadamard patterns from the remaining elements of the complete basis.

D. Phase shifting holography measurements
For the measurements using phase shifting digital holography, we used an interferometer in a Mach-Zehnder configuration. A collimated laser (Oxxius slim-532) was divided by a beam splitter into the object and the reference beam. The first one illuminated the object while the second traveled directly towards the camera (Allied Stingray F-145). The object, a layer of photoresist spin coated onto a transparent plate, was imaged onto the camera by an optical system in a 4-f configuration. The phase shifts were generated by shifting a grating codified onto a DMD (DLP Discovery 4100) and filtering the first diffraction order in the Fourier plane of a 4-f optical configuration located before the object plane [70]. The camera recorded four interferograms with a phase shift interval equal to π/2 and a standard phase-shifting algebraic operation [71] was used to measure the amplitude and phase at the object plane.

Measurement of an optical aberration and comparison with Shack-Hartman wavefront sensing and phase shifting interferometry
Fig. S1 shows the results for the plate that generates a coma aberration. Using our approach (experimental setup of Fig. 3a in the manuscript), the patterns projected onto the plate had a resolution of 64×64 pixels. As a way to smooth the retrieved phase distribution, the Zernike decomposition of the wavefront was calculated. Then, an image with high spatial resolution (512×512) was generated with that information. This process made it possible to compare that phase image with the one obtained with the Shack-Hartmann sensor used in our experiments, since the commercial software of this sensor directly calculates the Zernike decomposition and also provides a smoothed high-resolution image. Phase shifting interferometry results were obtained by applying four different phase shifts between the two arms of a Mach-Zehnder interferometer (see methods). The phase distributions provided by the three methods present high resemblance, as can be observed in Fig. S1. For our approach, the peak-tovalley phase value is 6.96 λ. The Shack-Hartmann sensor provides a value of 7.24 λ, while the corresponding value for phase measured by phase-shifting interferometry is 7.20 λ. The discrepancy between the above quantities is less than 5% of the total peak-to-valley phase value obtained with phase shifting interferometry, confirming the validity of the approach.

Complex amplitude retrieval of an object
Here we show a simple example of the retrieval of an amplitude and phase object with our technique. We used an object composed of an amplitude mask attached to the lens used in section 3 of the manuscript (Fig. S2a). The reconstruction basis is formed by Hadamard functions of 128 × 128 pixels, with a pixel pitch of 41.04 µm. The retrieved complex amplitude is shown in Fig. S2. As explained in the Methods section of the manuscript, after the position detector provides the total irradiance and the centroid position of the resulting light distribution for the full set of Hadamard patterns projected onto the sample. Using the irradiances provides the amplitude image (left inset of Fig. S2b), and the position of the centroids can be used to obtain the phase information (right inset of Fig. S2b).  study, a positive lens whose central part is covered by an amplitude mask. (b) Spatial distribution of the reconstructed complex amplitude. Given that the mask has opaque regions where the phase is not defined, we arbitrarily assign a null phase value to those object zones.

Relation between the phase of an object and the highorder statistical moments of its intensity Fourier distribution
Here we will discuss the relation between the phase of a wavefront in a given plane and the high-order statistical moments of the intensity distribution at the Fourier plane of a lens. This configuration is represented in Fig. S3. In order to work out the relation, Parseval's Theorem and the Moment's Theorem of Fourier Transforms will be invoked.
We assume the initial field to have the following form: 1 ( , ) = ( , ) ( , ) . By using Fresnel's propagation, we can link the field in the first plane, 1 ( , ), with the field in the second plane, 2 ( , ), in the following way: where ̃1 denotes the Fourier transform of 1 . As the planes we are working with are both at a distance from the lens, there is no global phase involved.
Given that we measure the light irradiance, the relevant quantity is the square modulus of this distribution: (2) Our detector will provide the total energy of this distribution, which we will call , and the centroid of the distribution, which we will call Δ ⃗ ⃗ = (Δ ; Δ ). Those physical measurements can be mathematically described with the following operations: .

(4)
Using Parseval's Theorem: that is, S is just the total energy at the first plane. Now, let us make the calculation of the centroid of the distribution. By using its definition, we have, for Δ : After a variable change ( ′ = ⁄ , ′ = ⁄ ) we obtain: For the sake of simplicity, from now on the marks in the spatial variables inside the integral will be removed. In order to solve the integral, we will and and are related via Fourier transform. Now, we see that Eq. (7) involves the first moment, 1,0 ( = 1, = 0). In order to solve the integral, we will also need the Fourier transform of |̃1( , )| 2 . It is easy to prove that [1]: In other words, the Fourier transform of the square modulus of a function is the complex cross correlation of its Fourier transform. Combining Eq. (7) and Eq.
Introducing the derivative into the integral: The derivative will be, after the variable change = + , = : So, what we get is the partial derivative of the complex field in one direction. In this way, Δ can be expressed as: In a similar way, Δ can be expressed as: Introducing the analytic form of 1 :  5) and (22). The field distribution at the first plane is related to the field distribution at the second plane through a thin lens. Both planes are separated a distance twice the focal length, f, of the lens.

(18)
The first term is zero for two main reasons. First, for the centroid position to be real. Second, the integrand is a perfect differential that integrates to zero, given that the amplitude goes to zero at infinity. It can easily be seen that: Then, for the first integral, one has: And Δ will be: The same procedure can be used to calculate Δ , and introducing the gradient operator, ⃗ = ( ; ) = (̂⁄ ;̂⁄ ), the centroid position can be described as: So, the centroid position of the intensity distribution in the detector plane is related to the gradient of the phase of the wavefront in the initial plane. By using Eq. (22), it is possible to relate different centroid measurements to the phase of the wavefront, as we will see in the next section.

Using a square mask to measure the local phase of the wavefront
Let us describe what happens if we sample the complex field with an amplitude square mask that only allows the light to pass through a square window of size at a position ( , ). The configuration is shown in Fig. S4.
Then, the field can be described by using the rectangular functions in the following way: Introducing this amplitude into the centroid expression, we get, for : Now, we can consider that if is small enough, the gradient of the phase of the object will be constant over the whole integration domain (i.e., the phase will be locally flat on that region). By doing this, the partial differentiation term goes out of the integral and we have: Doing the same for , we obtain: In this case, it is easy to prove that , the energy of the wavefront in the first plane, will be equal to 2 2 : then, by introducing the operator ⃗ ( , ) = (̂;̂)| ( , ) , we arrive to: = ( ; ) = 2 ⃗ ( , ) ( , ).

(28)
So, for a square aperture placed in the position ( , ), the value of the measured centroid will be proportional to the gradient of the object phase within the area of the aperture.

Spatial multiplexing of the phase information of a wavefront
Now we will see what happens if we have more than one aperture at the same time. In this case, the problem will be as depicted in Fig. S5. Now, the field will be given by: The energy at the object plane will be: The crossed term will be always zero if both masks do not overlap. This is the case we will always consider: several square masks acting onto the field at different positions. Then the energy will be: which is the sum of the energy of the two individual masks, as one can expect. Now, the centroid will be: We can see that this this result is similar to that obtained for the singleaperture one, but now with two integrals. By solving both integrals, Δ can be expressed as: If we identify 2 = as the energy of each square aperture, and 2 ( 1 2 + 1 2 ) as the energy of all the apertures combined, we have: and for Δ : which is the averaged sum of the two centroids, corresponding to the results of each individual aperture. We can generalize for an arbitrary number, , of non-overlapping square apertures: If all the apertures transmit the same energy, the first factor inside the summatory is simply 1/N, which leads to: So, for a combination of square apertures, the centroid of the resulting distribution is related to the gradient of the phase at each one of the positions of the apertures. If the wavefront is divided into square regions (or pixels), Eq. (39) can be used to recover the gradient of the phase at each region by generating N different distributions of squares. This is exactly what we do when we illuminate the object with Hadamard patterns, as each pattern can be viewed as a different combination of squares that either transmit or block the light. Fig. S5. Masking the wavefront with multiple square masks at the same time. Now, two regions of the wavefront will be transmitted through the system. We will see how the centroid of the distribution in the second plane is related to the phase of the wavefront in the first plane.