High-resolution synthetic-aperture digital holography with digital phase and pupil correction

A 218 mega-pixel synthetic aperture was collected by raster scanning a CCD detector in a digital holography imaging experiment. Frames were mosaicked together using a two-step cross-correlation registration. Phase correction using sharpness metrics were utilized to achieve diffraction-limited resolution. ©2011 Optical Society of America OCIS codes: (090.1995) Digital holography; (100.3010) Image reconstruction techniques. References and links 1. B. Hecht, J. Colineau, and J. C. Lehureau, “Short-range synthetic aperture imaging at 633 nm by digital holography,” Appl. Opt. 41(23), 4775–4782 (2002). 2. J. H. Massig, “Digital off-axis holography with a synthetic aperture,” Opt. Lett. 27(24), 2179–2181 (2002). 3. D. Claus, “High resolution digital holographic synthetic aperture applied to deformation measurement and extended depth of field method,” Appl. Opt. 49(16), 3187–3198 (2010). 4. L. D. Weaver, J. S. Fender, and C. R. Dehainaut, “Design considerations for multiple telescope imaging arrays,” Opt. Eng. 27, 730–735 (1988). 5. P. Nisenson, “Speckle Imaging with the PAPA Detector and the Knox-Thompson Algorithm,” in DiffractionLimited Imaging, D. M. Alloin and J.-M. Mariotti, eds. (Kluwer Academic Publishers, 1989), pp. 157–169. 6. M. Guizar-Sicairos, S. T. Thurman, and J. R. Fienup, “Efficient subpixel image registration algorithms,” Opt. Lett. 33(2), 156–158 (2008). 7. R. A. Muller and A. Buffington, “Real-time correction of atmospherically degraded telescope images through image sharpening,” J. Opt. Soc. Am. 64(9), 1200–1210 (1974). 8. S. T. Thurman and J. R. Fienup, “Correction of anisoplanatic phase errors in digital holography,” J. Opt. Soc. Am. A 25(4), 995–999 (2008). 9. S. T. Thurman and J. R. Fienup, “Phase-error correction in digital holography,” J. Opt. Soc. Am. A 25(4), 983– 994 (2008). 10. A. E. Tippie and J. R. Fienup, “Phase-error correction for multiple planes using a sharpness metric,” Opt. Lett. 34(5), 701–703 (2009). 11. U. Grenander, Abstract Inference (Wiley, New York, 1981). 12. A. E. Tippie and J. R. Fienup, “Multiple-plane anisoplanatic phase correction in a laboratory digital holography experiment,” Opt. Lett. 35(19), 3291–3293 (2010). 13. A. Kozma and C. R. Christensen, “Effects of speckle on resolution,” J. Opt. Soc. Am. 66(11), 1257–1260 (1976). 14. A. J. Page, L. Ahrenberg, and T. J. Naughton, “Low memory distributed reconstruction of large digital holograms,” Opt. Express 16(3), 1990–1995 (2008).


Introduction
For many imaging applications, one desires to achieve high resolution while maintaining a wide field-of-view.Typically, with conventional optics, resolution or field-of-view must be sacrificed to achieve the other.This tradeoff is due primarily to the increase in aberrations as the lens diameter and field-of-view increase.For optical systems that have demanding specifications on both the resolution and field-of-view, such as for lithographic lenses, the number of optical elements in the system increases to compensate and correct for these effects.Furthermore, for focal-plane imaging applications, the number of pixels in the detector array can limit the space-bandwidth product.In this paper, we discuss a laboratory experiment in which we were able to achieve both high resolution and a wide field-of-view simultaneously, all without imaging lenses in the system, with digital holographic aperture synthesis by using a computational imaging approach and extensive post-detection correction algorithms.
Large synthetic apertures can be formed by combining together multiple frames of digital holography data, thereby increasing the resolution and the space-bandwidth product of the system [1][2][3].As the synthetic aperture becomes larger, the system tolerance requirements become more demanding.Without proper knowledge and alignment of frames within the synthetic aperture, residual wavefront errors will occur due to lateral and longitudinal pupil geometry mismatches [4].Claus [3] used an iterative approach to estimate frame-to-frame phase offsets.To avoid aberrations due to a drift in the reference beam (local oscillator) phase relative to the object beam, phase-error correction algorithms are needed.To correct for the axial position of the plane of the object, which may vary in depth across its surface, dynamic focusing must be performed when the object depth exceeds the depth of field.Massig corrected for defocus [2].Binet et al. [1] showed that optical aberrations such as defocus and astigmatism were present after forming an image from a synthetic aperture.Compensation for aberrations in [1] were found manually with user input.We describe a laboratory experiment in which we assembled a large, 218 mega-pixel synthetic-aperture digital hologram.Our synthetic aperture reconstruction represents significant increase in pixel count compared to recent synthetic apertures such as Ref [3].This demonstrates a non-trival advancement in synthetic aperture digital holography, given the demands on the system tolerance and stability of the experiment, as well as the greater effects of wde-angle anisoplanatic phase errors.We developed phase-error correction algorithms (which require no user intervention) to correct frame-by-frame piston errors due to reference beam drift, reference beam location errors, local focusing errors, and other higher-order errors.We show quantitative results, demonstrating the diffraction-limited images capabilities upon implementing these techniques.These algorithms correct for pupil mapping errors which would otherwise cause space-variant aberrations over the aperture for a wide field-of-view and alleviate the tight data collection requirements, making the synthesis of large apertures more practical and affordable.

Experimental setup and data acquisition
The experimental setup was configured in an off-axis, lensless Fourier transform geometry, as illustrated in Fig. 1.Light from an argon-ion laser with an output power of 450 mW, operating at a wavelength of 514 nm was spilt into two arms, a reference beam and an object illumination beam.In the reference arm, neutral density filters reduced the intensity of the beam and a 5X microscope objective focused the light into the end of a 3.5 μm single-mode fiber.The output of the fiber was placed in the plane of the object and served as a diffraction- limited reference point source.Through a collection of mirrors, a second beam was positioned to flood illuminate the object at an off-axis angle of 18 degrees.The object used in this experiment was a USAF chrome-on-glass bar target.Flat white paint was applied on the chrome side of the target as a means of providing diffuse reflection, and that side of the target faced away from the illumination beam and the camera, so the object was seen through the back of the glass plate.A portion of the USAF target was masked, such that the effective object area was 4.7 cm by 6.3 cm.The USAF target was positioned normal to the optical axis of the system, making it parallel to the detector; mirror reflections from the glass substrate and the mirror-like chrome surface reflected away from the detector due to the off-axis illumination angle.The distance from the object to the detector was z = 3.32 m.The detector was a Retiga 2000R Q-Imaging CCD camera with a pixel pitch of 7.2 μm in both dimensions and an area of 1200 x 1600 pixels (8.64 mm x 11.52 mm).The CCD was mounted on two translation stages.The horizontal translation stage was mounted to the table, while the vertical stage held the camera.
In order to acquire all the frames necessary for the full-aperture image, the camera was raster scanned, left to right, to gather a row of frames starting at the bottom left-hand corner of the aperture array, and then subsequent rows of frames were gathered bottom to top.We always acquired frames in the same direction to minimize hysteresis and mechanical backlash of the translation stages.The translation between frames was about 833 pixels horizontally and about 555 pixels vertically giving an approximately 50% overlap between frames; the redundant data from the overlapping areas served to help align and register frames since the accuracy and precision of the translation stages used in this experiment proved to unacceptable for alignment requirements.
The exposure time for each frame was 80 ms.The total acquisition time for 21 x 21 frames (over 800 Mpixels, of which ¾ were redundant), included a two second wait and settle time for the translation stages after moving to each location, a dark frame exposure (for calibration), an object-only frame and a reference-only frame (for diagnostics), and a holography frame for each position, was 5 hours and 58 minutes under automatic computer control.

Registration and mosaicking of a large synthetic aperture
The hologram intensity in the detector plane is where R is the reference field and H is the object field.A holographic image can be reconstructed by recovering the field in the detector plane by standard techniques (Fourier transform, zero out all but the desired term and inverse Fourier transform) and then digitally propagating the resulting field to the image plane.show zoomed-in regions of the image to emphasize the finer resolution elements of the USAF target.For a single frame, the finest resolved bars were (group, element) (1,1) in the horizontal direction and (1,3) in the vertical direction, corresponding to 2.00 line-pairs (lp) per mm and 2.52 lp/mm, respectively.Aperture synthesis was performed by first extracting the detector-plane complex field of each individual frame and then mosaicking together those fields.We found it important to register and mosaic the extracted fields rather than the hologram intensities because the latter's quasi-periodic fringe pattern, along with piston phase changes between frames, made frame-to-frame registration of holograms very challenging.Frame-to-frame registration must be done with great accuracy, and Appendix A shows that to keep space-variant phase errors due to frame translation to 1/10 wave, one must register the frames accurately to within 1/7 pixel (about one micron for our CCD).We found that the speckles in the object beam at the detector changed somewhat over the long data collection time, contributing to errors when registering each frame within the larger aperture.To minimize the accumulation of registration errors across the synthetic aperture in our initial mosaic, we registered frames in a spiral-like pattern, starting in the center and spiraling out counterclockwise (even though the frames were gathered in a raster scan order), analogous to what is done with phase reconstructors [5].Our early preliminary results showed that a spiral registration was more effective than starting in one corner of the array and registering in a raster scan order.Further improvements in the registration approach should be possible.
The mosaicking of the large synthetic aperture was a two-step process, consisting of an initial frame registration based on the cross-correlation of the new frame with a single adjacent frame, followed by the cross-correlation of the new frame with the overlap regions from several frames.Using the a priori knowledge that there is ~50% overlap between adjacent frames, the registrations were performed using an efficient subpixel cross-correlation [6] estimating the x and y translations between the frames to within a tenth of a pixel and the piston phase.The set of images formed using these initial positions and phases for the mosaic are shown in Fig. 2 The initial mosaic formed by averaging all overlapping frames, but only taking into account a single adjacent frame in the registration process, is shown in Fig. 3.This full aperture takes into account the translation and phase information from the initial crosscorrelations between adjacent frames.The overall mosaic had a blocky, patchwork-like appearance instead of the expected uniform speckle pattern.The darker regions indicate destructive interference between fields from different frames and reveal errors in the piston phase estimates and/or registration errors.
To improve the uniformity and accuracy of the mosaic when averaging areas of overlap between frames, we performed a second cross-correlation.Now, instead of instead of crosscorrelating each subsequent individual frame with one adjacent neighbor, we performed a cross-correlation in the area of overlap between the single frame and the averaged sum of previous frames.This second step takes into account multiple frames of data and allows for adjustments in the translations and piston phase of frames for better agreement within the total synthetic aperture.
Figure 4 shows the resulting full aperture after both the initial mosaicking and the refinement, clearly showing improved uniformity and demonstrating the need for this refinement.The image after this refined mosaicking is shown in Fig. 2(i) -(l).We noted a further increase in resolution after this secondary alignment of frames: (3,2) is resolved for horizontal bars and (3,1) is resolved for vertical bars.

Phase correction and image reconstruction
The full-aperture complex field from the object in the detector plane was used to reconstruct the fine-resolution, wide-field-of-view image.We applied phase correction techniques to the field to compensate for residual phase errors that remain after the piston phase correction performed during mosaicking.Three phase correction algorithms we used were correcting (i) the location of the reference point source, (ii) residual higher-order phase errors, and (iii) local region phase errors due to variations in axial depth.

Reference correction
The complex field retrieved in the detector plane after windowing out other hologram terms is * .
, exp , where z is the distance from the detector plane to the image plane and   F denotes a Fourier transform.We found that our measurement of the reference source location was insufficiently accurate, leading to residual aberrations.To determine and correct that error, we maximized the sharpness of the image [7,8]  as a function of the reference source location, where   , I  is the intensity of the image and  is a scalar quantity.For this work, we used a value of  = 0.6 which is favorable for the case of a single speckle realization [9].(If  > 1, then we would maximize .S  ) We implemented this using a conjugate gradient nonlinear optimization routine using analytic gradients of S with respect to the three-dimensional coordinate position of the reference point source.Derivations of these analytic gradients can be found in Appendix B. Figure 2(m) -(p) shows the reconstructed image after optimizing and correcting for the reference point source location.The resolvable bars after the appropriate reference wavefront were (3,2) for horizontal bars and (3,3) for vertical bars.

Higher-order phase-error correction
In addition to applying a phase correction that specifically corrects for wavefront errors related to the reference source, we found it necessary to include an additional step to correct for other phase errors in the system.This may correct any residual errors from the mosaicking process, focus correction to locate the plane of best focus for the image, etc.We model the complex field in the detector plane after reference correction, as where 1 H is the complex field with phase errors, H is the ideal field and   , e xy  is the residual phase error.We found the e  that maximizes the sharpness of the image using a conjugate gradient search over e  using the analytic gradient of e S   [10].For this residual phase correction, we used the method of sieves [11,12], in which we convolve the point-by-point version of the analytic gradient with a Gaussian kernel.Important for dealing with very large arrays, the method of sieves requires less computational time and memory storage than optimizing over coefficients of polynomial expansions of the phase estimates, which need storage of polynomial basis functions for efficiency.Using a bootstrapping method, we started with a kernel having a standard deviation of 2500 pixels and reduced the size of the Gaussian kernel every five iterations (i.e.2500, 2000, 1500, 1000, 500, 250), ending with a final kernel with a standard deviation of 250 pixels.This allows the algorithm to first optimize over larger, lower-frequency global phase errors and then eventually correct for phase errors that are more local on a frame-by-frame basis, allowing for only smooth phase maps and avoiding stagnation during the optimization process.Visible in Fig. 5, which shows the phase error estimate after maximizing the sharpness of the image, are low-order focus correction and individual-frame higher-order correction.An important thing to note is that the sieves method of sharpness optimization had no explicit knowledge of the acquisition or mosaicking procedures; by optimizing the sharpness in the image plane, the algorithm was able to estimate a phase error that shows characteristics that are consistent with the frame-by-frame nature of the data in the detector plane.Figure 2(q) -(t) shows the further improved image after this phase correction is applied; horizontal bars (3,5) were resolved, as were vertical bars (3,5).

Region-of-interest phase correction
A phase-error correction applied over the entire field of view sharpens the image to the plane of overall best focus correction.For the experimental case presented here, varied defocus over a wide field-of-view object is likely since the depth of field is μm, so even a planar object having a small tilt can have appreciable depth.The depth of field would not be an issue for long-range imaging, but it was a problem in our laboratory experiment on account of the relatively large angle (2.3 degrees) subtended by the synthetic aperture.It is possible, however, to apply a phase correction that is optimized for a specific object depth by selecting a region of interest (ROI).
For ROI phase correction, the field in the detector plane was first numerically propagated to the nominal image plane.A mask was applied in the image plane to window out the specific ROI and the windowed field was then propagated back to the detector plane.This field, , ROI H can be written in terms of the original field as where

 
; zH A is the angular spectrum propagation of the field H by a distance z and n M is a mask defining the n th ROI.ROI H was then used in the sharpness minimization routine described above.The ROI selected for this example was a 600 x 600 pixel region and the mask was a square with a weighted cosine tapering function of 24 pixels applied to the edge of the region.M 1 included two horizontal bars from Group -1, Element 1, and M 2 encompassed Group 2 and higher as shown in Fig. 6.After finding the phase that minimizes the sharpness of a ROI, we corrected the entire synthetic aperture for that phase error, and the corrected field was then propagated using angular spectrum to the image plane.
Figure 7 shows the phase-error maps for the two ROI sharpness optimizations.Compared with Fig. 5, where the phase correction exhibited both global phase errors as well as frameby-frame correction, the primary phase correction seen by the ROI M 1 phase estimate are smooth higher-order errors.For the image reconstructions, this ROI phase estimate replaces the phase estimate found by utilizing the entire FOV.A comparison between the images resulting from these two different phase corrections is shown in Fig. 8. From Fig. 8(b) and (g), we see that the ROI phase estimates do improve and sharpen the image within the masked region.From Fig. 8(d) and (f), we also see that areas far from the ROI mask were blurred due to space-variant errors, presumably due to a tilted plane or residual errors in the reference location.
The image reconstruction using M 2 as a selective ROI is shown in Fig. 2(u) -(x).The resolution increase is substantial in comparison to our previous, intermediate images.For the   reference-only phase correction, the finest resolved bars in Fig. 2(m) -(p) were (Group, Element) (3, 2) for horizontal bars, a modest improvement over the (3,5) horizontally resolved bars that were seen after entire-image sharpness correction.However, for the ROI phase correction, the resolution increased an entire group (by a factor of two), with the finest resolved horizontal bars now (4,2) and vertical bars (4,5).This corresponds to 17.9 line-pairs (lp) per mm and 25.4 lp/mm, respectively, equivalent to a resolution of 56 μm and 39 μm, respectively.
The final size of the synthetic aperture was 12,100 by 18,000 pixels (87 mm by 129 mm), about 200 megapixels.The theoretical image resolution, given by zD  , where D is the synthetic aperture width, is 19.5 μm horizontally and 13.1 μm vertically.Since no speckle averaging has been performed, the expected resolution for a speckled image of a diffuse object will be reduced [13].For a USAF bar target, we found that the resolution of a speckled image of three-bar target degrades by about a factor of three.Therefore, the expected resolution of the synthetic aperture used in this experiment is 59 μm for horizontal bars and 39 μm for vertical bars, consistent with our experimental.Further experiments are underway to produce a gigapixel image using these techniques.

Computational requirements
The computational requirements are quite demanding for mosaicking, phase-error correction and image formation for such a large synthetic aperture.An ordinary desktop computer having a few GB of RAM cannot handle everything in fast memory necessary for such large arrays of data.For the computations performed in this experiment, an IBM x3755 server with four AMD Opteron 8224 SE processors and 128 GB of memory was used.The computation time for calculating initial cross-correlations between adjacent frames was 17.9 minutes.The mosaicking of the synthetic aperture took 9 hours and 19 minutes and the phase-error estimation took roughly 3 days of computational processing.Parallel implementation using a distributed system for the propagation kernel [14] would be one way to reduce the computation time for the reconstruction process.

Conclusion
We have demonstrated, to these authors' knowledge, the largest diffraction-limited synthetic aperture of digital holography data published to date.The synthetic aperture was formed by mosaicking together 21 x 21 frames of digital holography data.Steps needed to obtain accurate mosaicking included cross-correlation and intensity optimization of regions of overlap between adjacent frames.Several phase correction techniques using sharpness metrics were developed and implemented including the automatic correction of higher-order phase errors.The resulting image agrees with the theoretical diffraction-limited resolution of a speckle object.

Appendix A: Requirements for pixel registration
The discrete Fourier transform of an object mn f .
where N is the total number of pixels in the array  rad.This is consistent with the analogous situation of pupil mapping errors in multiple-telescope arrays [4].Based on the Marechal criterion, the phase accuracy needed for diffraction-limited resolution can be approximated as 1/14 = 0.072 waves rms for all phase errors.For the standard deviation of the phase to be within the Marechal criterion, the maximum allowable standard deviation of the displacement is given by max 1 14 , op mN   or max 0 (14 ).pixel.Therefore, one must register the frames accurately within a 1/7 pixel to meet the phase accuracy required to meet this diffraction-limited criterion for the worst case.These theoretical calculations emphasize the requirements for very accurate sub-pixel registration on a frame-to-frame basis across the entire synthetic aperture in order to maintain diffraction-limited resolution over the entire field of view.Note, however, that for points in the image near the optical axis, o m is small and image quality is relatively insensitive to translation errors.

Appendix B: Derivation of analytic gradient of S
The field obtained by angular spectrum propagation of  

Figure 2 (
a) -(d) shows the image from a single CCD frame; Fig. 2(a) shows the entire field of view of the image, while Fig. 2(b) -(d)

Fig. 3 .
Fig. 3. Amplitude of the mosaic of 21 x 21 (12100 by 18000 pixels) frames after initial crosscorrelation of single adjacent frame and averaging overlap regions.
(e) -(h).The improvement over a single frame in resolved bars of the USAF target was from (1,1) to (2,4) horizontally and from (1,3) to (2,5) vertically.Note that an improvement by (n bars, m elements) is equivalent to a resolution improvement by a linear factor of /6 2. nm  We wish to divide out the complex conjugate of the reference field, due to the reference source location applied at the detector plane, is the position of the reference point source, with the origin   , , 0 x y z  at the center of the detector plane and k = 2.  Because a Fresnel transform is insufficiently accurate for this wide-aperture, wide field-of-view image, we used angular spectrum propagation of to the image plane is given by

Fig. 5 .
Fig. 5. Higher order phase correction over the entire FOV.Colorbar units in radians.

Fig. 8 .
Fig. 8. Image reconstruction with ROI higher-order phase correction (a) -(d) using mask M1 and (e) -(h) using mask M2.Subimages (b) and (f) are 160 pixel areas of (a) and (e), respectively and (d) and (h) are 160 pixel subsets of the central areas of (c) and (g).
to the image plane is given by Eq. (3).The partial derivative of the sharpness metric S given in Eq. (4), with respect to the x-component of the reference point source location r x , is given by