## Abstract

Wavefront sensors and more general phase retrieval methods have recently attracted a lot of attention in a host of application domains, ranging from astronomy to scientific imaging and microscopy. In this paper, we introduce a new class of sensor, the *Coded Wavefront Sensor*, which provides high spatio-temporal resolution using a simple masked sensor under white light illumination. Specifically, we demonstrate megapixel spatial resolution and phase accuracy better than 0.1 wavelengths at reconstruction rates of 50 Hz or more, thus opening up many new applications from high-resolution adaptive optics to real-time phase retrieval in microscopy.

© 2017 Optical Society of America

## 1. Introduction

Wavefront sensing and phase retrieval are two closely related tasks that find applications in a wide range of fields. The classical Shack-Hartmann Wavefront Sensor [1] is commonly used in astronomy and ophthalmology. As a first-order method that tracks the 2D motion of focus spots generated by a lenslet array, the Shack-Hartmann sensor offers high frame rates, but suffers from low spatial resolution, corresponding to the number of sub-apertures rather than the number of pixels in the image sensor. While the spatial resolution can be improved by reducing the number of pixels per sub-aperture, e.g. to 2 × 2 pixels per lenslet as in the Altair adaptive optics system for Gemini [2], this comes at the cost of reduced range, such that large distortions can not be measured. Similar tradeoffs exist for other designs, including the pyramid wavefront sensor [3] and its variant [4], and the Hartmann wavefront sensor [5].

Better accuracy and higher spatial resolution can be achieved with Curvature Wavefront Sensors [6], and closely related phase retrieval methods developed recently in microscopy [7,8]. These methods are based on inverting the Transport of Intensity Equation (TIE) [9], which requires phase diversity, i.e. multiple image measurements along the optical axis. TIE has been extensively investigated [10], and proves to be effective and accurate approach for wavefront sensing. However, the requirement for phase diversity makes it difficult to design snapshot-capable systems that can image fast-moving dynamic phenomena.

In this paper, we propose a novel first-order wavefront sensor named Coded Wavefront Sensor. By placing a binary mask in close proximity of a camera sensor, one can numerically decode the wavefront from the apparent motion of the diffraction patterns. Our work is closely related to the sampling field sensor [11] and the quadri-wave lateral shearing interferometer wavefront sensor [12], but offers full sensor resolution under both coherent and incoherent illumination. Compared to the classical Shack-Hartmann wavefront sensor, the Coded Wavefront Sensor offers much higher pixel utilization and full sensor resolution reconstruction, overcoming the trade-off between spatial resolution and range. Specifically, the Coded Wavefront Sensor offers the ability to measure large distortions with high spatial resolution. Compared to TIE-based methods, the Coded Wavefront Sensor does not require phase diversity, works with both monochromatic and polychromatic illumination, is simple to build, and offers much faster reconstruction than video-rates. It is therefore highly suitable for imaging dynamic deformations of the wavefront.

## 2. Principle

A schematic of Coded Wavefront Sensor is shown in Fig. 1. A uniform random binary mask is placed closely (for instance, distance *z* ≈ 1.5 mm) in front of a bare image sensor. For calibration, the Coded Wavefront Sensor is illuminated by a planar wavefront, and the corresponding diffraction pattern of the mask on the image sensor is recorded. During measurement, distortions of the wavefront result in localized diffraction pattern displacements (small arrows in Fig. 1(b)), which can be tracked and used to measure the wavefront at the full resolution of the image sensor.

To model the principle of the Coded Wavefront Sensor, we denote the scalar field at the mask plane and the sensor plane as *u*_{0}(**r**) and *u _{z}* (

**r**) respectively, where

**r**= (

*x, y*) denotes a coordinate point on the 2D plane. The original scalar field

*u*

_{0}(

**r**) is the product of the mask transfer function

*p*

_{0}(

**r**) and the scalar optical field exp[j

*ϕ*(

**r**)], where

*ϕ*(

**r**) is the distorted wavefront under investigation.

For mask functions *p*_{0}(**r**) with small feature sizes, the effect of diffraction must be considered, whereas for low frequency wavefront *ϕ*(**r**), one may simply employ ray optics. Under proper approximation conditions, in the Appendix we show:

*k*= 2

*π/λ*is the wave number for wavelength

*λ*, and

*p*(

_{z}**r**) is the diffraction field of

*p*

_{0}(

**r**) under collimated illumination.

Denoting the calibration (reference) image and the measurement image as *I*_{0}(**r**) and *I*(**r**), given *I*_{0}(**r**) = |*p _{z}* (

**r**)|

^{2}, the measurement image

*I*(

**r**) is shifted relative to

*I*

_{0}(

**r**) by a point-wise apparent motion proportional to the wavefront slope ∇

*ϕ*(

**r**):

*k*if we consider the distorted wavefront

*ϕ*(

**r**) =

*ko*(

**r**) with optical path

*o*(

**r**), which means the Coded Wavefront Sensor allows for broadband illumination.

However, Eq. (2) is nonlinear. To retrieve *ϕ*(**r**) from *I*_{0}(**r**) and *I*(**r**), one may iteratively solve a linearized version of Eq. (2). At each step, the linearized version of Eq. (2) leads to the following formula, which is the basis for the so-called *optical flow* methods in computer vision [13]:

Compared to TIE that requires axial phase diversity, our model needs no *z*-direction displacements, because the distorted wavefront results in local transverse motion of the diffraction pattern. Our model Eq. (2) is nonlinear, involving first-order terms of the wavefront, whereas TIE is a linear model that employs second-order terms. Finally, in terms of optical path, Eq. (2) permits white light illumination, whereas TIE only allows for coherent or partially coherent illumination.

#### 2.1. Sensor performance

The wavefront slope is determined by the following relationship:

where*d*is the sensor pixel size, and

_{s}*a*is the amount of apparent motion measured in numbers of pixels. For a given tracking algorithm to solve Eq. (2)

*a*limits itself to a certain range, e.g.

*a*∈ [

*a*

_{min}

*a*

_{max}]. The sensitivity (or, the accuracy) of the sensor depends on

*a*

_{min}, whereas the dynamic range depends on

*a*

_{max}.

Denoting the local wavefront radius as *R*, one of the approximation made in deriving Eq. (1) requires (see Appendix):

*R*≫

*z*. For a given distance

*z*, Eq. (5) can be used to estimate the maximum curvature 1/

*R*that can be measured by our sensor.

#### 2.2. Numerical solver

While the purpose of the Coded Wavefront Sensor is to estimate *ϕ*(**r**) from *I*_{0}(**r**) and *I*(**r**), standard optical flow methods seek to reconstruct the per-pixel apparent motion vectors between the reference image and measurement image, i.e. the gradient of the phase function. Instead of using standard optical flow algorithms like Horn-Schunck [13], we therefore devise our own reconstruction model that directly solves for the phase function itself and regularizes the apparent motion to be curl free.

Define the unknown wavefront as ** ϕ**, the image gradient field as (

**g**

*,*

_{x}**g**

*) = ∇*

_{y}*I*

_{0}(

**r**), and a “time” derivative

**g**

*=*

_{t}*I*(

**r**) −

*I*

_{0}(

**r**), at each linearisation, the reconstruction problem is formulated as:

**G**= [diag(

**g**

*) diag(*

_{x}**g**

*)] is a concatenated diagonal matrix with the image derivatives on the diagonal, and*

_{y}**M**is a binary diagonal matrix that selects only the

*M*visible pixels from the

*N*wavefront samples that affect the measurement, for direction

*x*and

*y*respectively.

Like existing optical flow methods, Eq. (6) contains both a “photoconsistency” term (first term) that describes the apparent motion of the pattern in the image plane, and a smoothness regularization term (second term), which is controlled with an additional parameter *α* > 0. With a smoothness regularizer added directly to the phase function, Alternating Direction Method of Multipliers (ADMM) [15] is used to solve this optimization problem, where each updating step involves either element-wise operations or fast Fourier transforms (FFT), and hence is naturally parallelizable.

## 3. Implementation

#### 3.1. Hardware

Our prototype Coded Wavefront Sensor (see Fig. 4) exploits a uniformly random binary mask pattern, and a GS3-U3-14S5M-C PointGrey monochromatic 2/3″ CCD sensor with a pixel pitch of 6.45 µm. The mask is placed on top of the bare sensor, with a spacing of approximately 1.5 mm. The binary mask is fabricated using photolithography in a chrome layer deposited on a 4″ Fused Silica wafer, with a pixel pitch of 12.9 µm.

#### 3.2. Software

Our GPU version is implemented in
`CUDA`, and the
`CuFFT` library is used for all the involved FFT operations, which are the most time-consuming parts. To make the best utilization of
`CuFFT`, the unknown width and height are set to be the powers of small primes, e.g. two or three.

## 4. Results

In this section the results are presented, including the synthetic ones in simulation, the experimental validation using a Spatial Light Modulator (SLM), and the realistic wavefront visualization of heat flow and defocusing.

All numerical experiments are run with a fixed regularization parameter. Both the CPU and GPU version of our algorithm are run on a workstation that has 125.8 GB RAM, exploits Intel Xeon E5-2697V3 @2.60 GHz×16 as CPU, and GeForce GTX TITAN X (Pascal) as GPU, with a Ubuntu 14.04 Linux as operating system.

#### 4.1. Simulation

We have conducted two simulations to investigate the sensitivity and accuracy of the Coded Wavefront Sensor. In the simulations, the illumination is monochromatic (as of *λ* = 550 nm). The overall aperture size equals 6.6 mm × 6.6 mm with sensor and mask pixel pitch be 6.45 µm and 12.9 µm respectively. The scalar field of interest is sampled with 1.29 µm. Gaussian noise is added and the image signal-to-noise-ratio (SNR) equals 40 dB. The wave propagation is simulated using the angular spectrum method [16] with filtering [17] to suppress high frequency artifacts.

The first numerical experiment evaluates the dynamic range of our sensor. A planar wave (i.e. the reference), and sixteen different scales of spherical waves are simulated at the mask plane, respectively, for five different distances *z*. The reference image, and sixteen measurement images are consequently recorded at the sensor plane. Figure 2 shows the wavefront reconstruction error (in terms of root mean square, RMS) our sensor can attain, providing the fixed curvature wavefronts that progressively violate Eq. (2). With the increase of wavefront range, the decrease of accuracy can be partially explained by the approximations made to derive Eq. (1).

Figure 3 shows the second numerical experiment, where we evaluate the performance of our sensor by sensing typical atmospheric turbulence. The same turbulence is evaluated at different scales. The synthetic atmospheric turbulence respects the Kolmogorov’s theory, and is implemented using the sub-harmonic method [18]. The outer scale and inner scale of the base turbulence are set to be 4 m and 1 mm respectively. The mask-to-sensor distance *z* = 1.5 mm. The result indicates the possibility to apply our sensor for atmospheric turbulence measurement.

#### 4.2. Quantitative experimental demonstration

We also evaluated the accuracy of Coded Wavefront Sensor by using it to measure known wavefronts. The experimental setup is shown in Fig. 4. The different target wavefronts are generated using a reflective LCoS-based SLM, the PLUTO SLM by HOLOEYE, whose pixel pitch is 8.0 µm, and the maximum phase retardation of 2*π*. We use phase wrapping to simulate larger phase retardations. This introduces higher diffraction orders in the synthesized waveforms, which are detected by the Coded Wavefront Sensor and show up as noisy structures in Fig. 5. This limits the experiments to small and medium distortions. We used different scales of cubic wavefronts, spherical wavefronts, customized wavefronts, and Zernike single-mode wavefronts as ground-truth targets.

A distant, white point light source is used as collimated illumination. A polarizer is placed in the optical path to produce linearly polarized illumination, which is required for the SLM to operate as a pure phase modulator. Behind the beamsplitter, a Kepler telescope structure is imposed to ensure the SLM and our wavefront sensor are in conjugate planes, so the distorted wavefront measured by Coded Wavefront Sensor is the one produced by the SLM. In our setup, the Kepler telescope is composed of two lenses, with focus length *f*_{1} = 100 mm and *f*_{2} = 75 mm respectively. The sensor exposure time is set to be 5 ms at both the calibration step and the measurement step.

Selected experimental results for several different megapixel (1024 × 1024) wavefront shapes are visualized in Fig. 5 as synthetic interferograms (one interference ring equals *λ* = 632.8 nm). The typical mean absolute phase error is below 50 nm, and is slightly larger than the ones in the simulation. We believe that the quantization and phase wrapping in the SLM that we use to generate the distortion are introducing extra errors. The ground truth wavefront is not really the actual wavefront distortion for that reason, and the isocontours in the recovered wavefront can be clearly seen. The RMS for all the experimental wavefronts are listed in Table 1, where for most wavefronts the reconstruction error is below 0.1 wavelengths.

We achieve a reconstruction speed of up to 50 frames per second, i.e. 50 million individual phase measurements per second. It is possible to increase the frame rate by sacrificing spatial resolution. Further improvements in speed should be possible through numerical strategies such as “warm starting”, i.e. initializing the solution for one frame with the solution from the previous frame.

For reference, the GPU timing performance for different wavefront resolutions is in Table 2. All the timing experiments were run repeatedly for 1000 times, and the average is chosen as the final timing statistics. Our
`CUDA` GPU version greatly improves the timing performance, compared to the CPU version that is implemented in
`MATLAB`. For example, the average computing time for *N* = 1024 × 1024 is 2.00 s in the CPU version, over which our GPU implementation achieves ∼ 102 times speedup. For one megapixel wavefront resolution, we have achieved wavefront reconstruction speed as 51 Hz, which is much beyond real-time performance.

#### 4.3. Realistic wavefront imaging

We also visualized two realistic wavefronts, the ones created by heat flow and defocusing, using the Coded Wavefront Sensor. The heat flow is generated using a lighter, and defocus is achieved by manually moving a convex lens back and forth. To increase the field-of-view for better visualization, we employed a telescope system for wavefront magnification, with a ratio of two. All the measurement images were captured and the wavefronts were reconstructed on GPU in real-time. Please see Visualization 1 for the experimental setup, the whole captured data and the reconstructed wavefronts. Here, two frames are chosen to be shown in Fig. 6, with their reconstructed wavefronts in interference fringes respectively.

## 5. Discussion

The Coded Wavefront Sensor is related to a number of other imaging systems and designs. Just like the Shack-Hartmann sensor can be interpreted as a combination of a lenslet-based light field camera [19] and 2D spot tracking software, our sensor can be seen as a combination of a mask-based light field camera [20, 21] with a more sophisticated dense 2D motion tracking method. The Coded Wavefront Sensor also bears similarity to Background Oriented Schlieren (BOS) imaging [22], but with the patterned “background” moved into the camera for a compact form factor.

The design of the Coded Wavefront Senors allows it to be used as a drop-in replacement for any optical system currently using a Shack-Hartmann sensor, with an immediate gain in spatial resolution. In addition, we believe it can also be incorporated into optical systems configured for phase retrieval problems, including in microscopy. Since there is no need for phase diversity and coherent illumination, such an adaptation should in fact be easier than many existing phase retrieval setups. We will explore these applications in future works.

Several other extensions in both applications and usage are also conceivable in the future. The calibration or reference image does not always need to correspond to plane wave illumination, but could be a pre-distorted wavefront. For example when characterizing freeform lenses, a known ground-truth lens can be used to form the reference image, and the Coded Wavefront Sensor can then be used to characterize the difference between the reference lens and another lens.

For the accuracy of the Coded Wavefront Sensors it is necessary that the mask produces a locally distinctive diffraction pattern on the image sensor. To facilitate this process, the mask could be custom-designed (instead of random) to produce a specific diffraction pattern such as wavelet noise [23]. It should also be possible to use grayscale masks or even random phase gratings as an alternative to the binary masks employed in this work.

## 6. Conclusion

In conclusion, we introduce the Coded Wavefront Sensor, a novel sensor design that is physically implemented by a single binary masked sensor to encode the incoming wavefront, and is numerically implemented by an efficient optimization decoding algorithm, such that wavefront reconstruction with high spatio-temporal resolution is achieved within sub-wavelength accuracy. The theoretical principle behind Coded Wavefront Sensors offers a new approach to the wavefront sensing problem, namely the direct 2D tracking of diffraction patterns.

## Appendix

Our mathematical model is based on modifications and simplification of the Green’s function model by Teague [9] (Appendix A.3). Consider a general scalar field *u*_{0}(**r**) (where **r** = (*x, y*) denotes a coordinate point on a 2D plane of interest) with wave number *k*, on the mask plane:

*f*

_{0}(

**r**) = exp[j

*ϕ*(

**r**)] is the scalar field corresponding to the wavefront

*ϕ*(

**r**) that we want to measure, and

*p*

_{0}(

**r**) is the transfer function of the mask. We assume:

- The mask
*p*_{0}(**r**) is of high frequency (in our case, it is uniformly random binary), whose Fourier transform*P*_{0}() (where*ρ*is the Fourier dual of*ρ***r**) is broadband. - The wavefront
*ϕ*(**r**) is of low-frequency, and is second-order differentiable, and hence the distorted scalar field*f*_{0}(**r**) is smooth enough such that its spectrum is bandlimited, and decays sufficiently close to zero in high frequency regions.

We now consider the wave propagation problem for small distance *z*. For a general mask *p*_{0}(**r**), the Fresnel approximation may not be valid, and the more general Huygens-Fresnel based principle must be applied. Using the compact form of the Rayleigh-Sommerfeld diffraction formula, and expand it in the Fourier domain, the diffractive scalar field *u _{z}* (

**r**) is:

**″ =**

*ρ***–**

*ρ***′. And the approximation step comes from:**

*ρ**F*

_{0}(

**) with ${\lambda}^{2}{\Vert \mathit{\rho}\u2033\Vert}_{2}^{2}\ll 1$. The last term in Eq. (8) is the diffraction of**

*ρ*″*f*

_{0}(

**r**–

*λz*

**′), for which the**

*ρ**n*th diffraction order [(

*z*/

*k*)∇

^{2}

*ϕ*(

**r**)]

*= (*

^{n}*z*/

*R*)

*≈ 0 when local radius*

^{n}*R*≫

*z.*Also, a Taylor expansion of

*ϕ*(

**r**–

*λz*

**′) up to second-order suggests:**

*ρ**d*, by the Nyquist Theorem that ${\Vert \mathit{\rho}\prime \Vert}_{2}^{2}\le 1/\left(4{d}_{m}^{2}\right)$: In our prototype,

_{m}*d*= 12.9 µm and ${d}_{m}^{*}\approx 1$ µm for

_{m}*R*= 1 m. Given

*R*≫

*z*and the valid linear approximation of

*ϕ*(

**r**), we approximate the last term in Eq. (8) as:

*p*(

_{z}**r**− (

*z/k*)∇

*ϕ*(

**r**)) is the diffractive scalar field of

*p*

_{0}(

**r**) at the sensor plane.

## Funding

King Abdullah University of Science and Technology (KAUST) baseline funding; KAUST Advanced Nanofabrication Imaging and Characterization Core Lab.

## Acknowledgments

The authors would like to thank Yifan Peng for his help in the initial optical setup and valuable discussions.

## References and links

**1. **B. C. Platt and R. Shack, “History and principles of Shack-Hartmann wavefront sensing,” Journal of Refractive Surgery **17**, S573–S577 (2001). [PubMed]

**2. **L. K. Saddlemyer, G. Herriot, J.-P. Véran, and J. M. Fletcher, “Design aspects of the reconstructor for the Gemini adaptive optics system (Altair),” Proc. SPIE **3353**, 150–159 (1998). [CrossRef]

**3. **R. Ragazzoni, “Pupil plane wavefront sensing with an oscillating prism,” J. Mod. Opt. **43**, 289–293 (1996). [CrossRef]

**4. **O. Fauvarque, B. Neichel, T. Fusco, S. Thierry, and J.-F. Sauvage, “Variation around a pyramid theme: optical recombination and optimal use of photons,” Opt. Lett. **40**, 3528–3531 (2015). [CrossRef] [PubMed]

**5. **A. F. Brooks, T.-L. Kelly, P. J. Veitch, and J. Munch, “Ultra-sensitive wavefront measurement using a Hartmann sensor,” Opt. Express **15**, 10370–10375 (2007). [CrossRef] [PubMed]

**6. **F. Roddier, “Curvature sensing and compensation: a new concept in adaptive optics,” Appl. Opt. **27**, 1223–1225 (1988). [CrossRef] [PubMed]

**7. **L. Waller, L. Tian, and G. Barbastathis, “Transport of intensity phase-amplitude imaging with higher order intensity derivatives,” Opt. Express **18**, 12552–12561 (2010). [CrossRef] [PubMed]

**8. **Z. Jingshan, R. A. Claus, J. Dauwels, L. Tian, and L. Waller, “Transport of intensity phase imaging by intensity spectrum fitting of exponentially spaced defocus planes,” Opt. Express **22**, 10661–10674 (2014). [CrossRef] [PubMed]

**9. **M. R. Teague, “Deterministic phase retrieval: a green’s function solution,” J. Opt. Soc. Am. **73**, 1434–1441 (1983). [CrossRef]

**10. **T. E. Gureyev, A. Pogany, D. M. Paganin, and S.W. Wilkins, “Linear algorithms for phase retrieval in the Fresnel region,” Opt. Commun. **231**, 53–70 (2004). [CrossRef]

**11. **R. Tumbar, R. A. Stack, and D. J. Brady, “Wave-front sensing with a sampling field sensor,” Appl. Opt. **39**, 72–84 (2000). [CrossRef]

**12. **J. C. Chanteloup, “Multiple-wave lateral shearing interferometry for wave-front sensing,” Appl. Opt. **44**, 1559–1571 (2005). [CrossRef] [PubMed]

**13. **B. K. Horn and B. G. Schunck, “Determining optical flow,” Artificial intelligence **17**, 185–203 (1981). [CrossRef]

**14. **T. Brox, A. Bruhn, N. Papenberg, and J. Weickert, “High accuracy optical flow estimation based on a theory for warping,” in ECCV’04, 8th European Conference on Computer Vision (2004) pp. 25–36.

**15. **S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed optimization and statistical learning via the alternating direction method of multipliers,” Foundations and Trends in Machine Learning **3**, 1–122 (2011). [CrossRef]

**16. **J. W. Goodman, *Introduction to Fourier Optics* (Roberts and Company Publishers, 2005).

**17. **K. Matsushima and T. Shimobaba, “Band-limited angular spectrum method for numerical simulation of free-space propagation in far and near fields,” Opt. Express **17**, 19662–19673 (2009). [CrossRef] [PubMed]

**18. **R.G. Lane, A. Glindemann, and J.C. Dainty, “Simulation of a Kolmogorov phase screen,” Waves in random media **2**, 3, 209–224 (1992). [CrossRef]

**19. **R. Ng, M. Levoy, M. Brédif, G. Duval, M. Horowitz, and P. Hanrahan, “Light field photography with a hand-held plenoptic camera,” Computer Science Technical Report CSTR **2**, 1–11 (2005).

**20. **A. Veeraraghavan, R. Raskar, A. Agrawal, A. Mohan, and J. Tumblin, “Dappled photography: Mask enhanced cameras for heterodyned light fields and coded aperture refocusing,” ACM Trans. Graph. **26**, 69 (2007). [CrossRef]

**21. **K. Marwah, G. Wetzstein, Y. Bando, and R. Raskar, “Compressive light field photography using overcomplete dictionaries and optimized projections,” ACM Trans. Graph. **32**, 46 (2013). [CrossRef]

**22. **H. Richard and M. Raffel, “Principle and applications of the background oriented schlieren (bos) method,” Measurement Science and Technology **12**, 1576 (2001). [CrossRef]

**23. **B. Atcheson, W. Heidrich, and I. Ihrke, “An evaluation of optical flow algorithms for background oriented schlieren imaging,” Experiments in Fluids **46**, 467–476 (2009). [CrossRef]