Phase-Error Correction for Coherent Array Imaging Systems

It is often possible to reduce the requirements on an imaging system by placing greater demands either on an illumination system or on post-detection processing of the data collected by the system. An extreme example of this is a system with no receiver optics whatsoever. By illuminating an object or scene with coherent light having a shaped illumination pattern, the receiver can be a simple detector array with no imaging optics, detecting the speckle intensity pattern reflected from the object; an image of the object can be reconstructed by a phase retrieval algorithm.


Introduction
This paper describes a concept for an active imaging approach that forms images without imaging optics, allowing the receiver to be wide, allowing fine spatial resolution, but is thin in the axial direction, making it compact and light weight.It reduces the requirements on the optical hardware by increasing the requirements on the post-detection processing.An illustration of the system concept is shown in Fig. 1.If an object is flood-illuminated by a coherent laser beam, the reflected optical field at the object propagates to produce a speckle pattern in a plane a distance from the object.If the coherent field were measured in that plane, by heterodyne detection, then it would be possible to reconstruct an image in the computer by digitally propagating the field back to the object plane using a Fresnel or Fourier transform.However, heterodyne detection over a large array of detectors, needed to form an image of a complicated object, is expensive and beyond the current state of the art at visible and near-IR wavelengths.If just intensity measurements are made, then it is also possible to reconstruct an image under some circumstances.One such circumstance is holography [1].If near the plane of the object there is a suitably placed corner reflector satisfying the holographic separation condition, then the measured intensity pattern is a Fourier transform hologram of the object, and an image of the object can be reconstructed by simply inspecting the autocorrelation function obtained by Fourier transforming the detected intensity pattern [2].An appropriate holographic reference (often referred to as a beacon) is, unfortunately, not available for most scenarios.A second approach is imaging correlography [3,4].If one can collect a large number of independent realizations of the speckle intensity patterns, and compute the autocovariance (the ensemble averaged autocorrelation, after subtracting the mean) of the speckle patterns, then one arrives at the squared magnitude of the Fourier transform of the incoherent object (the object intensity reflectivity had it been incoherently illuminated), in a manner equivalent to Hanbury Brown-Twiss intensity interferometry.Then a phase retrieval algorithm using a nonnegativity constraint and a support constraint, based on the illumination pattern, can reconstruct an image.Phase retrieval algorithms for real-valued, nonnegative images are usually successful even if the support constraint is fairly loose.A problem with imaging correlography is that the signal-to-noise ratio (SNR) of this process is very low for higher spatial frequencies, yielding images of poor resolution unless a very large number of speckle realizations (which can be had by translating the detector or rotating the object) are used.A third approach is PROCLAIM, which involves illuminating a 3-D object with a series of laser wavelengths or optical frequencies [5] to collect a "cube" of data.The opacity of 3-D objects makes reconstruction of 3-D images from the cube of data fairly robust, but the object is not allowed to move appreciably during the entire data-collection time, which is problematic for many scenarios.A fourth approach is to have an object that exists on a dark background, such as an airborne target or a satellite orbiting the earth.If the object has abrupt (sharp) edges, and the shape of the object is favorable, then it is possible to reconstruct an image of the object with a phase retrieval algorithm [6,7].One must, however, make an accurate estimate of the support of the object (the set of points outside of which the object is zero) [8,9], and not all objects have favorable supports for phase retrieval.

Conformal Array of Detectors Coherent Illuminator
In this paper we discuss an approach to lensless imaging that can get around the problems mentioned above.If we illuminate the object or scene with a laser beam having a known, favorable shape, then we can use that shape as a support constraint for a phase retrieval algorithm to reconstruct an image.It is important to realize that the illumination pattern should come from a projection system that is much smaller in diameter than the array of speckle intensity measurements at the receiver.If the transmitter (illumination) optics had a large, diffraction-limited aperture, then one could just use that as the receiver optics as well, and lensless imaging would not be needed.Hence we restrict our attention to transmitter diameters small compared with the diameter of the receiving detector array.For this reason, the illumination pattern would have soft, tapered edges.While phase retrieval of complexvalued objects can be robust for sharp-edged patterns, it has been shown to work poorly for soft-edged support constraints [10].In this paper we show that it is possible to reconstruct images with such soft-edged illumination patterns if one employs special illumination patterns and advanced forms of the phase retrieval algorithm.Such an imaging system has several advantages over conventional imaging systems.First, the system is thin in the axial direction; whereas a conventional f/1 optical system has a depth comparable to its width, this system can have a depth many times smaller than its width.Similarly, the weight of the receiver can be many times less than a conventional telescope.The array of detectors can be on a flat or curved surface, allowing it to be conformal to, say, the fuselage of an airplane.The detector array does not have to be steered, omitting the large gimbal system used conventionally, and that allows the aperture to be wider, for a given platform, improving the resolution of the system.With the laser illumination, it is also possible for the detector array to be one-dimensional, sweeping out the second dimension in time as a synthetic aperture, further reducing the size and weight of the detector array.Disadvantages of this approach include the need for a moderately high-power, coherent laser illuminator, substantial computational requirements, and a coherent image, which suffers from speckle.
Section 2 of this paper describes the support constraint formed by the illumination pattern.Section 3 describes the advanced phase retrieval algorithm: the hybrid input-output algorithm using an expanding Fourier modulus [9], and image refinement using an improved form of a "patching" algorithm [11].Section 4 shows some digital simulation results, including the effects of noise on the image reconstruction, and Section 5 draws conclusions.

Support constraint
Since a main advantage of the imaging concept is to have a compact, thin, and possibly conformal sensor, it would be disadvantageous to have a large-aperture optic projecting the illumination pattern.If, for example, the diameter of the projection optics were 1/4 that of the receiver array, then the volume of the projection optics would be approximately the cube of this, or 1/64 the volume that the optics would be if it were of the full diameter.For this reason, we consider projection optics having a diameter 1/4 that of the receiver optics to be small enough.Since atmospheric turbulence would be greatest nearer the ground, the projection of the illumination pattern should be minimally affected by turbulence (although if one were projecting in the opposite direction, it would be a problem).Furthermore, since the illumination pattern has soft edges, minor blurring of it should not have a large effect on the phase retrieval algorithm, which does not assume a tight support constraint.
It has long been known that certain classes of support constraints are more favorable than others.For example, support constraints that are highly asymmetric work better for phase retrieval than do symmetric supports [6], partially because phase retrieval is "more unique" in that case [12,13] and partly because it avoids the "twin image" mode of stagnation [11].Similarly, support constraints with separated parts are favorable.While phase retrieval is ordinarily unique for 2-D objects, it is usually highly ambiguous for 1-D objects [14].Nevertheless, phase retrieval for 1-D objects with supports having separated parts is usually unique, despite the general non-uniqueness of 1-D phase retrieval [15].Furthermore, phase retrieval has been shown to work better in practice for 2-D objects with supports having separated parts [6].For these reasons we choose an illumination pattern that is both highly asymmetric and has separated parts, as will be shown in Section 4.

Phase retrieval algorithm
Many types of phase retrieval algorithms have been developed, the two most prominent being the hybrid-input-output (HIO) version of the iterative transform algorithm and gradient-search nonlinear optimization algorithms [16,17].For this study we employed HIO, but it does not work well for this application in its simplest form, since it usually stagnates when one has a complex-valued object and only a soft-edged support constraint [10].
We here employ a more robust version of the algorithm using an expanding, weighted Fourier magnitude [9].Since that version of the algorithms is not in wide use, it will be briefly described here.It begins by multiplying the measured Fourier magnitude (the square root of the measured intensity) by a narrow weighting function (effectively using an aperture of narrower diameter), which reduces the resolution of the reconstructed image.Additionally, by virtue of the tapered nature of the weighting function, it reduces the sidelobes of the impulse response, thereby minimizing the energy due to diffraction that would fall far outside the support constraint.During the early iterations, the narrow width of the weighting function effectively allows only the phases over a narrower aperture to be retrieved.Since at this point only a low-resolution image of the object is being reconstructed, it finds a solution many times faster than it would for an image with the full space-bandwidth product.Also, for such low resolution images there are fewer local minima in which the algorithm can become trapped.After the algorithm has converged for this low-resolution image, the weighting function is replaced by a modestly wider one, over a wider effective aperture.At this stage we keep the previously retrieved phases over the narrower aperture, and the new phase values being retrieved by the iterative algorithm are those in the wider aperture just beyond that narrower aperture.After a few iterations, we also allow the previously retrieved phases to change as well.This process is repeated with weighting functions of successively larger diameters, as illustrated in Fig. 2, until the entire aperture is included.Effectively one bootstraps from phases over smaller apertures to the phase over the entire aperture.
A second algorithm enhancement employed here is an improved version of the patching algorithm.A given reconstructed image can have some residual artifacts due to the failure of the phase retrieval algorithm in certain areas of the Fourier plane [11].By reconstructing the image multiple times, each time from a different starting guess, one can get multiple reconstructed images, each with different residual artifacts.By combining information from these multiple images one can reconstruct an improved composite image.The patching algorithm described previously [11] can be generalized to employ more than two images, as follows [18].Run the HIO algorithm with multiple starting guesses to arrive at K different initial reconstructed images, each with different residual artifacts.If one of them is artifact free (which will be known if the reconstruction is in agreement with the Fourier data and support constraints to within the measurement noise), then no further work is required.If there is no artifact-free image, continue as follows.By finding the location of the maximum of the cross-correlation of the images, register (to within a fraction of a pixel) the images to one another.Zero out each of the images in their region of support, leaving the energy outside the support constraint.Fourier transform each of these.The location of the energy in the Fourier transforms are where the Fourier phase is in error for each reconstruction.Produce a smooth version of the magnitude of these Fourier transforms by convolving each with a kernel of several pixels in diameter, producing a smoothed map of the Fourier phase error.At each pixel in the Fourier domain, select the phase of the Fourier transform (of the entire initial reconstructed image) having the smallest Fourier phase error as computed above.This results in a composite Fourier transform.Inverse Fourier transforming that result gives a composite reconstructed image.That image is input to the HIO algorithm for further iterations.This extended patching algorithm is meant only as a refinement when the initial reconstructed images are of fairly good quality so that each has regions of its Fourier transform with a good approximation to the true Fourier phase.For a given illumination pattern, only a limited area of the scene is illuminated and imaged.To arrive at an image of the entire scene, one would use multiple lasers pulses with the illumination pattern translated to cover different, but partially overlapping, portions of the scene, and collect the far-field speckle patterns and perform image reconstruction for each.The partial overlap is important in order to more accurately register the multiple images and further to provide speckle reduction when the overlapping images are averaged in intensity.With a sufficient number of appropriately placed illumination patterns, one can arrive at a speckle-reduced, mosaicked image of the entire scene.

Digital simulation and reconstruction results
Figure 3 shows computer-simulated data used for digital reconstruction experiments.Figure 3(a) shows the magnitude of a portion of a complex-valued synthetic-aperture radar (SAR) image of Michigan Stadium taken with ERIM's (now, General Dynamics Advanced Information Systems') DCS SAR.Shown is the upper left quadrant of the 384x384 image.
Square roots of the SAR image magnitudes are shown here because of their large dynamic range.Figure 3(b) shows the upper left quadrant of a simulated laser illumination pattern.Unless otherwise indicated, subsequent SAR images show just the upper left quadrant, since that is all that was illuminated.Figure 3(c) shows the illuminated scene, which is the product of the object and the illumination pattern.The width of the illumination pattern, being less than 1/2 the width of the 384-width array, ensures that the intensity of the Fourier transform of the illuminated scene is Nyquist sampled.
The illumination pattern consists of three patches formed from three square areas of uniform intensity, each of width 38 pixels, convolved with the impulse response due to the projection-optics aperture.That illumination aperture was a circle of diameter 1/4 the width of the detector array, and it was weighted with a radial version of the Hanning (raised cosine) weighting function to reduce sidelobes in the illumination pattern.Three separated parts to the illumination pattern can be accomplished by a number of approaches, e.g., using a diffractive optical element.Figure 4 shows a horizontal cut through the center of the amplitude of the two horizontally displaced patches of the illumination pattern.The upper curve is the same cut multiplied by 50, to better show the sidelobes.This illustrates that the illumination pattern does not go to zero between the patches or outside of the patches.Yet the phase retrieval algorithm, employing a support constraint, assumes zero values in those areas.Despite the fact that this is not strictly true, the phase retrieval algorithm reconstructs an image by finding one agreeing with the measured Fourier intensity pattern that has the minimum energy outside the support constraint.x (pixels) Field Amplitude Fig. 4. Cut through illumination pattern amplitude (lower curve).Upper (dashed) curve is 50x.
Figure 5 shows the 384x384-pixel magnitude of the Fourier transform of the illuminated scene, which is the data collected at the detector array.Photon noise was simulated in the Fourier intensity data with, for this case, an average of 100,000 photons per pixel, for an SNR of 316.(Lower SNR cases are shown later.)Shown here, it includes a separable Hanning weighting on the magnitude, apodizing the values at the edges of the detected data in order to reduce sidelobes in the reconstructed image, thereby making the support constraint work better.Fine fringes in the horizontal, vertical, and 45-degree directions can be seen throughout the pattern, a result of the interference between the fields from the separated parts of the illuminated object.This pattern is not a hologram, since none of the separated parts of the object is a delta function.Nevertheless it is conceivable that these fringe patterns help to encode the phase of the field into the intensity (perhaps the local deviation of a fringe is approximately related to the phase, as it would be for an interferogram), thereby making phase retrieval more robust for objects like this with separated supports than for objects with contiguous supports.The initial estimate of the image was a field of complex-valued random numbers multiplied by a version of the illumination pattern magnitude thresholded at 0.2 times its maximum magnitude.Sequences of iterations were performed, each with 40 iterations of HIO followed by 5 iterations of the error-reduction algorithm (i.e., just satisfying constraints in each domain).For the first six sequences of iterations, the support constraint was the illumination pattern magnitude thresholded at 0.2 times its maximum.For the next seven sequences of iterations, the threshold was decreased to 0.05, and for the final subsequent sequences of iterations, the threshold was decreased to 0.01.Hence, the support constraint started narrower and was relaxed to a wider support constraint for later iterations.The Fourier taper function used was the autocorrelation of a circular aperture having a radius that was some fraction of the width of the sensor aperture.That fraction was 0.1 for the first sequence of iterations, 0.15 for the second, 0.2, for the third, and was increased by 0.1 to 1.0, which is the sequence shown in Fig. 2.After that, no additional Fourier tapering was performed aside from the Hanning weighting.Figure 6 shows a partially reconstructed image as the iterations progressed, starting with a low-resolution image with a narrow Fourier weighting function for early iterations and finishing with a fine resolution image for late iterations using the entire Fourier magnitude.This reconstruction sequence was repeated two more times, with different random-number starting guesses.The three resulting partially reconstructed images were combined using the extended patching method, and then two more sequences of iterations were performed on that result to arrive at the final reconstructed image.Figure 7 illustrates the steps in the extended patching method.In this case, the entire 384x384 array is shown.As discussed earlier, a composite Fourier transform is constructed by patching together the phase from the Fourier transforms of multiple partially reconstructed images.The selection of the best Fourier transform to use for a given spatial frequency (for a given detector pixel) is determined as follows.First, set to zero the field inside the region of support, leaving just the field outside the support, as shown in Fig. 7 (g)-(i).This is the component of the field we wish to suppress.Next, compute the Fourier transform each of those fields, as shown in Fig. 7 (j)-(l).Next, smooth the intensity of that pattern, and at each pixel note which one had the minimum value.The composite Fourier transform is then formed from the Fourier transforms of the entire partially reconstructed images according to which one had the minimum value at each pixel as computed above.For this particular case, each initial reconstruction was fairly good, with only minor artifacts, and the patching method served to clean up the image.It may be more important for other cases, in which the initial reconstructed images are of lower quality.If the initial reconstructed images are very poor, then patching together their Fourier transforms does little good.
Figure 8 shows a second imaging example for which the illumination pattern has the shape of a triangle, but with soft edges owing to the effect of the diameter of the illumination optics.Ordinarily a triangle would be considered to be a shape very favorable to phase retrieval, since it is highly asymmetric and (if it has nonzero values at the three apexes) always has a unique solution in the absence of noise [6,12].Hard-edged triangular-shaped complex-valued images reconstruct well [10].The reconstructed image is significantly blurred, and the quality of the image is substantially inferior to the case of an illumination pattern having separated parts as shown earlier.In this case the image quality is so poor that the patching method was  ineffective because the bad areas of the Fourier domain were so widespread.This result emphasizes the benefits of the separated parts within the illumination pattern.
Figure 9 shows the effect of SNR of the measured data on the quality of the reconstructed image.Figure 9(a) shows, for reference, the illuminated object. .This shows that the quality of the reconstructed image degrades slowly and gracefully as the SNR decreases, and the phase retrieval algorithm can reconstruct good-quality images under relatively low light levels (100 photons/pixel).

Conclusions
We have shown that it is possible to reconstruct a coherent image from a far-field speckle intensity pattern using a known illumination pattern, formed from a small-aperture illumination system.This allows one to image with a detector array with no imaging optics.Illumination patterns particularly favorable to phase retrieval are required, such as asymmetric patterns with separated parts.While straightforward applications of the iterative transform algorithm were not capable of reconstructing an image under these circumstances, more advanced forms were.Especially important was the use of an expanding weighted Fourier magnitude.An extended version of the patching algorithm for image refinement from multiple reconstructed images was also demonstrated.The phase retrieval algorithm was shown to be tolerant of noise in the data.
Portions of this paper were presented in [19,20]

Fig. 1 .
Fig. 1.Lensless imaging with a detector array, sensing the intensity of a laser speckle pattern backscattered from the object.[Figure modified from a figure from Brad Tousley (DARPA/TTO)]

Fig. 2 .
Fig. 2. Radial cut through expanding weighting function on the Fourier magnitude.Lower curves are used for earlier iterations and upper curves for later iterations.The diffraction cutoff frequency is at pixel 192.

Fig. 6 .
Fig. 6.Partially reconstructed images with increasing iterations and increasing resolution owing to the use of the expanding Fourier magnitude.(a) through (f) are for the first six weighting functions shown in Fig. 2.

Fig. 7 .
Fig. 7. Illustration of extended patching method.(a)-(c) Three partially reconstructed images from different starting guesses; (d)-(f) same, but overexposed by a factor of 20; (g)-(i) same, but with region of object support masked out; (j)-(l), the magnitudes of the Fourier transforms of (g)-(i) respectively, showing the regions in which the Fourier transforms of the images in (a)-(c) are in error.