Holographic imaging with a Shack-Hartmann wavefront sensor

A high-resolution Shack-Hartmann wavefront sensor has been used for coherent holographic imaging, by computer reconstruction and propagation of the complex field in a lensless imaging setup. The resolution of the images obtained with the experimental data is in a good agreement with the diffraction theory. Although a proper calibration with a reference beam improves the image quality, the method has a potential for reference-less holographic imaging with spatially coherent monochromatic and narrowband polychromatic sources in microscopy and imaging through turbulence. © 2016 Optical Society of America OCIS codes: (090.0090) Holography; (090.1995) Digital holography; (100.3190) Inverse problems; (110.0110) Imaging systems; (010.7350) Wave-front sensing References and links 1. D. Gabor, “A new microscopic principle,” Nature 161, 777–778 (1948). 2. M. K. Kim, Digital Holographic Microscopy (Springer, 2011). 3. U. Schnars, C. Falldorf, J. Watson, and W. Jüptner, Digital Holography and Wavefront Sensing (Springer, 2015). 4. J. W. Goodman and R. Lawrence, “Digital image formation from electronically detected holograms,” Appl. Phys. Lett. 11, 77–79 (1967). 5. J. W. Goodman, Introduction to Fourier Optics (Roberts and Company Publishers, 2005). 6. P. Ferraro, S. De Nicola, A. Finizio, G. Coppola, S. Grilli, C. Magro, and G. Pierattini, “Compensation of the inherent wave front curvature in digital holographic coherent microscopy for quantitative phase-contrast imaging,” Appl. Opt. 42, 1938–1946 (2003). 7. G. Vdovin, H. Gong, O. Soloviev, P. Pozzi, and M. Verhaegen, “Lensless coherent imaging by sampling of the optical field with digital micromirror device,” J. Opt. 17, 122001 (2015). 8. J. R. Fienup, “Phase retrieval algorithms: a comparison,” Appl. Opt. 21, 2758–2769 (1982). 9. R. A. Gonsalves, “Phase retrieval and diversity in adaptive optics,” Opt. Eng. 21, 829 (1982). 10. Y. Shechtman, Y. C. Eldar, O. Cohen, H. N. Chapman, J. Miao, and M. Segev, “Phase retrieval with application to optical imaging,” IEEE Signal Process Mag. 32(3), 1–25 (2014). 11. R. W. Gerchberg and W. O. Saxton, “A practical algorithm for the determination of phase from image and diffraction plane pictures,” Optik 35, 237 (1972). 12. F. Meng, D. Zhang, X. Wu, and H. Liu, “A comparison of iterative algorithms and a mixed approach for in-line x-ray phase retrieval,” Opt. Commun. 282, 3392–3396 (2009). 13. C. Guo, C. Wei, J. Tan, K. Chen, S. Liu, Q. Wu, and Z. Liu, “A review of iterative phase retrieval for measurement and encryption,” Opt. Lasers Eng., http://www.sciencedirect.com/science/article/pii/S0143816616300197 (2016). 14. B. C. Platt and R. Shack, “History and principles of Shack-Hartmann wavefront sensing,” J. Refract. Surg. 17, 573–577 (2001). 15. S. A. Alexandrov, T. R. Hillman, T. Gutzler, and D. D. Sampson, “Synthetic aperture fourier holographic optical microscopy,” Phys. Rev. Lett. 97, 168102 (2006). #264688 Received 6 May 2016; revised 1 Jun 2016; accepted 2 Jun 2016; published 13 Jun 2016 © 2016 OSA 27 Jun 2016 | Vol. 24, No. 13 | DOI:10.1364/OE.24.013729 | OPTICS EXPRESS 13729 16. J. Holloway, M. S. Asif, M. K. Sharma, N. Matsuda, R. Horstmeyer, O. Cossairt, and A. Veeraraghavan, “Toward long distance, sub-diffraction imaging using coherent camera arrays,” http://arxiv.org/abs/1510. 08470 (2015). 17. W. Southwell, “Wave-front estimation from wave-front slope measurements,” J. Opt. Soc. Am. 70, 998–1006 (1980). 18. Y. Carmon and E. N. Ribak, “Phase retrieval by demodulation of a Hartmann-Shack sensor,” Opt. Commun. 215, 285–288 (2003). 19. A. Talmi and E. N. Ribak, “Direct demodulation of Hartmann-Shack patterns,” J. Opt. Soc. Am. A 21, 632– 639 (2004). 20. C. Canovas and E. N. Ribak, “Comparison of Hartmann analysis methods,” Appl. Opt. 46, 1830–1835 (2007). 21. K. R. Freischlad and C. L. Koliopoulos, “Modal estimation of a wave front from difference measurements using the discrete fourier transform,” J. Opt. Soc. Am. A 3, 1852–1861 (1986). 22. F. Roddier and C. Roddier, “Wavefront reconstruction using iterative Fourier transforms,” Appl. Opt. 30, 1325– 1327 (1991). 23. G. Dai, Wavefront optics for vision correction, Vol. 179 (SPIE Press, 2008). 24. A. Talmi and E. N. Ribak, “Wavefront reconstruction from its gradients,” J. Opt. Soc. Am. A 23, 288–297 (2006). 25. D. C. Ghiglia and M. D. Pritt, Two-dimensional phase unwrapping: theory, algorithms, and software, vol. 4 (Wiley New York, 1998). 26. R. M. Goldstein, H. A. Zebker, and C. L. Werner, “Satellite radar interferometry: two-dimensional phase unwrapping,” Radio Sci. 23, 713–720 (1988). 27. R. Horstmeyer, R. Heintzmann, G. Popescu, L. Waller, and C. Yang, “Standardizing the resolution claims for coherent microscopy,” Nat. Photon. 10, 71 (2016). 28. M. Loktev, O. Soloviev, S. Savenko, and G. Vdovin, “Speckle imaging through turbulent atmosphere based on adaptable pupil segmentation,” Opt. Lett. 36, 2656–2658 (2011). 29. D. de Lima Monteiro, G. Vdovin, and P. Sarro, “High-speed wavefront sensor compatible with standard CMOS technology,” Sens. Actuat. A 109, 220 – 230 (2004).


Introduction
Recording of the complex electromagnetic field with the method of holography allows to reconstruct both the phase and the amplitude at any propagation point [1,2].In general, optical holograms are created by recording the interference patterns created by the object and the reference waves.The object wave field is restored then by scattering the reference beam on the hologram [3].
Methods of digital holography register the wave field in the computer [4], to facilitate the object field reconstruction by computational wave propagation using physical optics models of diffraction and interference [5].Digital holography has potential advantages of higher speed of the hologram acquisition, simple experimental setups, and most notably, the availability of the complete amplitude and phase information about the optical field.Once the complex amplitude is known, it is possible to manipulate the wavefront to correct aberrations retrospectively [6,7], extending the methods of adaptive wavefront correction to virtual domain.
Since the phase information is not directly available, different methods of phase reconstruction from intensity measurements have been developed.Phase retrieval and diversity methods have been used to derive the computer description of the optical field [8][9][10].These methods reconstruct the complex amplitude by solving an inverse source problem, using numerical methods, such as Gerchberg-Saxton algorithm [7,11,12].In these cases reconstruction of the complex field from recorded intensities alone represents an ill-posed problem [13].
The Shack-Hartmann (SH) sensor is a simple tool commonly used in adaptive optics to register the arrays of local wavefront tilts [14].Each subaperture builds a spot-like point-spread function (PSF) image, where the deviation of the centroid of the light spot is proportional to the local wavefront derivative.The wavefront can be efficiently reconstructed from the array of the centroid deviations, in a single shot.
In this work we describe experimental realization of holographic image reconstruction from the phase and intensity obtained directly with a high-resolution SH sensor.

Method
The scheme of a holographic imaging setup based on a SH sensor is shown in Fig. 1.A transparent object is illuminated by a spatially coherent light beam and the scattered field is registered by a high resolution SH sensor.Registered intensity contains sampled information about the phase gradients and the intensity field, which can be used to reconstruct the complex field in the sensor plane, and to obtain a coherent holographic image of the object by back-propagating the complex field to the object plane.Though the method is applied to a transparent object, a similar setup could be used for imaging of reflective and scattering objects.For a SH sensor having full aperture of A and pitch p, the resolution of the sensor can be estimated from the analysis of interference produced by two coherent point sources S 1 , S 2 belonging to the object plane.If the distance between these two sources is equal to B, then the interference fringe period δ in the sensor plane, at a distance L from the object, is given by: To avoid any information loss, the fringe period should be smaller than the sensor aperture: δ < A. From this condition, we can derive the expression for minimum distance between two point sources that is still resolved by the sensor: where object features smaller than B min will not be resolved by the sensor.On the other hand, from the Nyquist criterion, the minimum fringe period in the sensor plane should cover at least 2 micro-lenses: δ = 2p, from which we can derive the expression for the maximum field of view in the object plane: Condition (3) has a simple physical meaning: the whole object should be small enough to remain unresolved by a single lenslet.This condition defines the difference between our method and the approaches described in Refs.[15,16], where each subaperture resolves the object.If the object is larger than B max , the reconstructed field will have parasitic low-frequency modulation, caused by aliasing.In some cases, when the low frequency modulation is not important, the field of view can be chosen larger than that defined by the condition (3).
The total number of resolved elements, in the ideal aliasing-free case, is given by: The number of resolved elements is one half of the number of micro-lenses along the chosen coordinate, and the total number of pixels resolved in the object plane is one quarter of the number of micro-lenses in the sensor matrix.
Reconstruction of the complex optical field from a high-resolution SH pattern requires several steps.The phase is obtained by integrating the wavefront gradients [17].Since the number of resolved elements in the hologram is proportional to the number of lenslets, large lenslet arrays should be used for high-resolution imaging.Special methods of phase reconstruction, suitable for processing of large arrays of spot images produced by a SH sensor have been developed in Refs.[18][19][20].In these methods, the SH intensity pattern is represented as a composition of a series of interference fringes.The pattern is Fourier transformed, then the first side lobes along the frequency axes f x and f y in the Fourier domain are moved to the origin.Then, the phase derivatives φ x and φ y are obtained as the argument of the inverse Fourier transform of the shifted distribution of intensity.Then, the gradients W x and W y can be derived as: where F is the focal length of the micro-lens array.
Once the gradient information has been obtained, the wavefront can be reconstructed by a Fourier-based modal reconstruction [21][22][23].
Here the Fourier series is seen as a set of basis function of the wavefront.Ŵ is used to denote the angular spectrum of W .In Ref. [24] it was shown that replacement of f by 2 sin( f /2) in Eq. ( 6) gives a better noise performance.We have tried both reconstructors, but found no noticeable difference.Before the Fourier based wavefront reconstruction can take place, it is essential to preprocess the gradient fields.The φ x , φ y derived from the arguments of exponential functions range in [0, 2π], possessing numerous discontinuities.Depending on the method, phase unwrapping [25] might be required to obtain a smooth phase function that can be used for building the physical optics model of propagation.
The intensity field can be obtained by smoothing the spot pattern with a suitable lowfrequency filter.The final complex field in the sensor plane is obtained by combining the phase and intensity fields.Finally, the object amplitude and phase distributions are reconstructed by back propagation of the complex field to the object plane.

Simulation
For the numerical simulation we have adopted the parameters of the experimental setup, described in section 4, shown in Fig. 1: a collimated beam with the wavelength of λ = 633 nm is scattered by positive transparent resolution test chart, shown in Fig. 2. To represent the complex field, we used a square array with sampling dimensions of n × n = 2048 × 2048 pixels with pixel pitch µ = 5.5 µm.After passing through the sample, the wave is propagated to the SH sensor at a distance of L = 0.5 m.The SH array is formed by N ×N = 140×140 micro-lenses in orthogonal geometry with a pitch of p = 63 µm.The focal length of each micro-lens is 2 mm, and the total size of the array is A = 8.82 mm.According to Eq. ( 2) the resolution of our setup in the object plane is equal to: and the total field-of-view B max = B min N/2 = 2.45 mm.The simulation was initiated with filtering the coherent light wavefield through the intensity mask.Then the filtered wave was propagated to the SH sensor, filtered through the phase mask corresponding to the SH array, and propagated to the image plane, where the intensity distribution was registered by a simulated camera with linear intensity response with pixel pitch of µ = 5.5µm.The simulated intensity pattern I 0 (x, y), registered by the SH sensor, is shown in the middle of Fig. 2.
The right part of Fig. 2 shows the result of discrete Fourier transform (DFT) of the registered SH sensor intensity pattern I 0 (x, y).The four symmetrical sidelobes at distance d = µN/p ≈ 179 pixels from the center contain the phase information.The scale of the DFT image is defined by the zero padding of the input array.The sidelobes have been extracted in the window of 140×140 pixels, corresponding to the number of lenslets N, and translated to the origin.The wavefront gradients were obtained by applying inverse DFT to the centered sidelobe.Figure 3 shows the wrapped distributions φ x and φ y .Two-dimensional Goldstein branch cut unwrapping algorithm [25,26] was used to unwrap the gradient fields.The wavefront W (x, y), reconstructed using Eqs.(5,6) is shown in the right part of Fig. 3.As expected, the wavefront has low-amplitude high-frequency modulation, with some large phase jumps localized in the areas of low intensity, where the phase reconstruction is ill-defined.
The SH pattern I 0 (x, y, L) was smoothed by a Gaussian filter to obtain the intensity field Ĩ0 (x, y, L), shown in Fig. 4. The complex field registered in the SH sensor plane was composed as U(x, y, L) = Ĩ0 (x, y, L) exp(iθ (x, y, L)), where θ (x, y, L) = kW (x, y, L) and k = 2π/λ .The complex field in the resolution test chart plane U(x, y, 0) can be reconstructed by back propagation: Figure 4 shows the filtered intensity in the sensor plane and the numerical reconstructed image of the resolution test chart in the object plane, at a distance of L = −0.5 m.The resolution of the image corresponds to the theoretical limit.The results of numerical experiment clearly demonstrate the validity of the method.As with any other SH sensor, the reconstruction is expected to be not very sensitive to the degree of temporal coherence in the illumination.We have simulated the hologram reconstruction obtained with broadband sources composed of three monochromatic lines with the bandwidth ∆λ = 20 nm, formed by three lines at λ 1 = 623, λ 2 = 633, λ 3 = 643 nm; ∆λ = 50 nm, formed by three lines at λ 1 = 598, λ 2 = 633, λ 3 = 658 nm; and ∆λ = 200 nm, formed by three lines at λ 1 = 533, λ 2 = 633, λ 3 = 733 nm.For each polychromatic source, we simulated the SH spot pattern by incoherent summation of three intensity patterns obtained for three wavelengths The resulting polychromatic spot pattern I ∑ has been used to coherently reconstruct the wavefront and the object, using the central wavelength of λ 2 = 633 nm.The reconstruction results, shown in Fig. 5, demonstrate robust reconstruction for source bandwidth of up to 50 nm (up to 200 nm with significant resolution loss), proving the method conditional usability with polychromatic spatially coherent sources.

Experiment
For experimental validation we have built a setup according to the scheme shown in Fig. 1, with all parameters matching the simulation described in the previous section.Positive USAF 1951 test target (R1DS1P, Thorlabs, U.S), shown in Fig. 6  In order to eliminate the systematic aberrations, the reference wavefront has been registered without the object.In principle, this step can be excluded, but then some virtual adaptive optics needs to be used for the compensation of systematic system aberrations.Calibration with the reference beam is a simple step, which needs to be done only once, then the calibration data can be used as long as the illumination beam is stable.After filtering and correction for the reference phase, we retrieved the intensity and phase distributions in the plane of SH sensor.Figure 6 illustrates the intensity corresponding to the object field, registered by the SH sensor.The image of the resolution test target, obtained by back propagation to the distance of L = −0.5 m, is shown in Fig. 7.
Robust reconstruction of the wavefront is possible only when spots from all subapertures are present in the image registered by the SH sensor.Theoretically, the high harmonics of the scattered light should enter all subapertures, even those in the geometrical shadow.However, registration of these harmonics requires a large dynamic range and high signal to noise ratio of the image sensor.In our experiment, with standard USAF 1951 resolution chart, we have resolved the second element in the third group.The resolution is defined by a standard formula R = 2 G+(E/6) ≈ 8.95 lp/mm, where G = 3 and E = 2 are the group and element numbers.This corresponds to the resolved line pitch of 56 µm, which is close to the theoretical limit of 35 µm.We attribute the resolution loss to the measurement noise, limited dynamic range of the camera, and aliasing caused by the too large object size.Also, the resolution of the method can not be ; reported exactly, due to some ambiguity in the definition of coherent resolution, caused by the intrinsic nonlinearity of coherent imaging.See [27] for further discussion on the definition of coherent resolution.

Discussion
We have demonstrated holographic imaging with a SH sensor, experimentally realizing a transmissive lensless imaging setup with close to the diffraction limited resolution on a high-contrast object.
This technique holds promise for microscopy and general coherent imaging with extended depth of field in applications that require fast registration and fast processing of a large number of relatively low-resolution holograms, such as dynamic holographic interferometry, flow cytometry, dynamic microscopic imaging.In addition, direct reconstruction of the complex field allows to create a virtual adaptive optical system, computationally correcting for the aberrations in the imaging path.This approach can be used for coherent multi-aperture imaging of remote objects through atmospheric turbulence, in systems similar to suggested in [16,28].
The experimental setup can be further optimized for speed and simplicity.In particular, the number of camera pixels used for hologram registration, can be optimized.In our setup we used more than 100 pixels to register the light spot under each single lenslet.Essentially, a quad cell with 4 pixels is sufficient for the registration of both the intensity and the center of gravity of the light spot.This brings us to the minimum requirement to the number of pixels in the registration camera: 4 imaging sensor pixels are needed per micro-lens, and 4 micro-lenses are needed per reconstructed hologram pixel, resulting in the minimum requirement of 16 camera pixels per reconstructed pixel in the final image.This conclusion sets theoretical limit to the number of pixels, in the assumption of linear response of the quad cell to the spot coordinate.In practice, a system with large number of quad cells is expected to be non-linear and difficult to align, but nonetheless a quad-cell based Hartmann sensor proved possible in the previous work from the authors [29].Coherent back propagation extends the information contents beyond the single intensity distribution, by using the phase to reconstruct the intensity in a number of different planes along the propagation path.
There is no principal requirement for a reference beam, though the quality of reconstruction is better when the reference beam is used to account for systematic aberrations.To achieve the theoretical resolution, the imaging setup requires the incoming light to have a high degree of spatial coherence over the whole lenslet array.However, simulation proves that the instrument is not very sensitive to the degree of temporal coherence in the illumination and can be used with spatially coherent polychromatic sources.
In this work we have investigated the bare-bone proof of concept.It can be further developed for a higher resolution and higher numerical aperture by using denser and larger micro-lens arrays, and by introducing additional optics for proper optical coupling of the illumination beam to the object and the sensor, towards usable coherent imaging instrument.

Fig. 2 .
Fig. 2. Numerical model: resolution test chart (left).Intensity pattern retrieved from the SH sensor at a distance of 0.5 m from the chart (middle).The central part of the Fourier transform of the intensity pattern with sidelobes (inset) used to reconstruct the x and y components of local WF tilts (right).

Fig. 3 .
Fig. 3. Gradients φ x , φ y , corresponding to the diffraction on USAF test chart, reconstructed from the inverse Fourier transform (left, middle), and the wavefront reconstruction (right).

Fig. 4 .
Fig. 4. Sensor intensity obtained in simulation, by filtering the SH pattern (left), and the reconstruction of the resolution test chart by back propagation of the reconstructed complex field to the object plane (right).

Fig. 5 .
Fig. 5. Simulated image reconstruction obtained with spatially coherent 633 nm monochromatic source (left) and polychromatic sources with bandwidth of 20, 50 and 200 nm, centered at 633 nm (images 2 to 4, counted from left to right).
was used as the test object.The object was placed at a L = −0.5 m distance in front of the SH mask.The SH array (OKO Tech) is formed by N × N = 140 × 140 micro-lenses in orthogonal geometry with a pitch of p = 63 µm, focal length of F = 2 mm and the total size of the square array of A = 8.82 mm.High resolution camera UI-3370 (IDS corporation) has been used for image registration.

Fig. 7 .
Fig. 7. Filtered experimentally registered intensity (left), and object reconstruction, obtained by back propagation of the reconstructed wave to -0.5 m (right).