Abstract
This paper explores the use of single-shot digital holography data and a novel algorithm, referred to as multiplane iterative reconstruction (MIR), for imaging through distributed-volume aberrations. Such aberrations result in a linear, shift-varying or “anisoplanatic” physical process, where multiple-look angles give rise to different point spread functions within the field of view of the imaging system. The MIR algorithm jointly computes the maximum a posteriori estimates of the anisoplanatic phase errors and the speckle-free object reflectance from the single-shot digital holography data. Using both simulations and experiments, we show that the MIR algorithm outperforms the leading multiplane image-sharpening algorithm over a wide range of anisoplanatic conditions.
© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. INTRODUCTION
Digital holography (DH) uses coherent illumination and spatial-heterodyne detection to sense the amplitude and phase of light scattered off an object’s surface [1]. The resulting data are sensitive to phase errors caused by index-of-refraction perturbations in, for example, the atmosphere or optical systems. With this in mind, we can estimate these phase errors directly from the DH data for wavefront sensing purposes or to digitally correct aberrated images [2–6].
Image-sharpening (IS) algorithms can estimate the phase errors from DH data by maximizing an image-sharpness metric. In practice, the sharpened images result from estimates of the complex-valued reflection coefficient, , which is the complex-valued ratio of the reflected field to the incident field. For surfaces that are rough relative to the illumination wavelength, this estimation process leads to images with high-spatial-frequency variations known as speckle. Therefore, IS algorithms require incoherent averaging of multiple speckle realizations to accurately estimate the phase errors—a process referred to here as multishot DH [2].
Recently, we developed a model-based iterative reconstruction (MBIR) algorithm for jointly computing the maximum a posteriori (MAP) estimates of the phase errors, , and the real-valued reflectance, , from a single data realization—a process referred to here as single-shot DH [5]. We define the reflectance as , where indicates the expected value [7]. In general, is smoother and has higher spatial correlation when compared to the reflection coefficient, . Reconstructing , rather than , allows us to leverage its higher spatial correlation to better constrain the estimation process and produce more accurate estimates of the phase errors with fewer data and less signal, compared to IS algorithms [5,6].
Both the IS algorithm in Ref. [2] and the MBIR algorithm in Ref. [5] were designed for the estimation and correction of phase errors that are isoplanatic. In general, isoplanatic phase errors (IPEs) result in a shift-invariant point spread function (PSF) in the image domain. We typically model IPEs using a single phase screen located near the entrance pupil of the imaging system [8]. This model accurately represents scenarios where the phase errors are concentrated near the observer, such as the imaging of space-based objects through atmospheric turbulence along mostly vertical propagation paths.
For other scenarios, the phase errors are often distributed along the propagation path, such as the imaging of Earth-based objects through atmospheric turbulence along mostly horizontal propagation paths. In these cases, the phase errors are anisoplanatic, meaning they result in a shift-variant PSF in the image domain. We typically model anisoplanatic phase errors (APEs) using a discrete representation consisting of multiple phase screens distributed along the propagation path [8]. To help quantify APE effects, we use the isoplanatic angle, [9]. This atmospheric turbulence parameter, in practice, gives us a gauge for the image-domain angular separation over which the PSF is essentially unchanged. Therefore, points in the image separated by an angle larger than have significantly different PSFs due to the APEs.
Importantly, the IS and MBIR algorithms in Refs. [2] and [5] can only correct for a single PSF across the entire image field of view (FOV). Thus, in the presence of APEs, these algorithms cannot obtain diffraction-limited performance. Estimating and correcting APEs using DH data is difficult due to increased model complexity and an increased number of unknowns. Extensions of the basic IS algorithm from [2] have been demonstrated for the correction of APEs [10–12]. In particular, Tippie and Fienup demonstrated an IS algorithm for estimating APEs modeled as multiple phase screens distributed along the propagation path [11]. However, this IS algorithm still requires multishot DH data to accurately estimate the phase errors and form a focused speckle-free image.
In this paper, we propose a novel algorithm for estimating APEs and reconstructing focused speckle-free images from single-shot DH data. Our algorithm, referred to as multiplane iterative reconstruction (MIR), uses a Bayesian framework to jointly estimate both the APEs modeled as multiple phase screens distributed along the propagation path, along with the object’s reflectance, . Using MIR, we compute the joint MAP estimates of each phase screen along with the focused image. We compare the proposed algorithm to the multiplane IS (MIS) algorithm found in Ref. [11] over a range of anisoplanatic conditions using both simulated and experimental data. The results show that the MIR algorithm can more accurately estimate APEs and form focused, speckle-reduced images from single-shot DH data as compared to the leading MIS algorithm.
In what follows, we first formulate the overall concept for the MIR algorithm and review the prior models used in our Bayesian framework. Next, we formulate the execution of the MIR algorithm and our methods for quantifying performance. Results then follow with discussion and a conclusion.
2. MIR ALGORITHM
Figure 1 shows an example DH system using an off-axis image-plane recording geometry (IPRG) [4]. The system images an object with a reflectance function, , and a corresponding reflection coefficient, , through distributed-volume phase errors, . In this example, atmospheric turbulence causes anisoplanatic effects in the DH data; however, the analysis is equally applicable to more general types of APEs, as we show in the experimental results section. For mathematical simplification, it is common to treat distributed-volume phase errors, , as a finite number of discrete layers, , where each layer is a thin unit-amplitude phase screen that represents the phase perturbations along an incremental volume [8]. After encountering , the signal is imaged onto a detector array, where it is mixed with a reference field to form the digital hologram.
Figure 2 provides insight into the conventional image-processing steps used for a DH system with an off-axis IPRG. We first extract the signal, , from the spectrum of the digital hologram. Due to the discrete Fourier transform (DFT) relationship between the pupil plane and the image/detector plane, is a complex image of the signal at the pupil plane. Thus, the conventional approach to form an image is to take an inverse fast Fourier transform (FFT) of the data, . Before moving on in the analysis, the reader should note that we provide additional details about this DH IPRG model in Appendix B.
In this work, our goal is to jointly compute the MAP estimates of the reflectance, , and the phase errors, , where , from a single noisy data realization, . Here, we use vectorized notation for all 2D variables. The joint MAP estimates are given by
where and are assumed to be independent. Here, and represent conventional and conditional probability distributions, respectively. For coherent imaging applications, the log-likelihood function, , is not easily evaluated due to its computational complexity [13]. Therefore, we require a simplified approach to solve Eq. (1). In the remaining subsections, we derive the form of , develop a simplified framework using the expectation maximization (EM) algorithm, and define our distributions for and .A. Log-Likelihood Model for Coherent Imaging
For a DH system using an off-axis IPRG, we model the complex pupil image,
as a linear transform of the reflection coefficient, , plus complex-valued measurement noise, [5]. In Ref. [5], Eq. (2) was derived for a DH system using an off-axis pupil-plane recording geometry. However, as shown in Appendix B, the model also holds for the pupil image, , obtained from a DH system using an off-axis IPRG. The transform, , accounts for the propagation and measurement geometry and is dependent on .For anisoplanatic conditions, the estimation framework of Eq. (1) and data model of Eq. (2) are more complex than they are for the isoplanatic case treated in Ref. [5]. First, we must account for intermediate propagations between multiple phase screens, which adds complexity to the structure of . With this in mind, we use the split-step beam propagation method (BPM) to model as a series of Fresnel propagations between the planes of interest, as described in Appendix A [8]. Second, we must model and estimate higher dimensional phase errors, where for phase screens, we get . This makes the joint estimation of and significantly more challenging, since the number of unknowns increases relative to the measurements.
To determine the likelihood function, , we first define distributions for and ; then we relate these distributions to using Eq. (2). For a surface with reflectance, , which is rough relative to the illumination wavelength, we model the vector as
where denotes an operator that produces a diagonal matrix from its vector argument and indicates a multivariate complex normal distribution with mean, , covariance matrix, , and pseudo-covariance matrix, [5]. Here, , since the phase of is uniformly distributed. Next, to model the stochastic nature of coherent detection, we treat the measurement process as being shot-noise limited and use an additive, zero-mean, complex Gaussian model [14]. The distribution for is therefore given by where is the noise variance. Finally, using the data model from Eq. (2) and the distributions from Eqs. (3) and (4), the likelihood function for , given and , becomes where the superscript indicates the Hermitian transpose.B. Prior Models
Both the reflectance, , and the phase errors, , are modeled using Markov random fields (MRFs) with the form
where is the difference between neighboring pixel pair values, is the partition function, is the weight between pixel pairs, is the set of all pair-wise cliques falling within the same neighborhood, and is the potential function [15]. We obtain different models for the reflectance and phase errors by choosing different potential functions.1. Reflectance
For the reflectance function, , we used a Q-generalized Gaussian Markov random field (QGGMRF) potential function with the form
where is a unitless threshold value that controls the transition of the potential function from having the exponent to having the exponent , and controls the variation in [16]. The parameters of the QGGMRF potential function affect its shape, and therefore the influence neighboring pixels have an effect on one another. In general, the QGGMRF potential function is designed to preserve edges.2. Phase Errors
We assume that the layers of the distributed-volume phase error, , are statistically independent [8]. Therefore, we can write the distribution for as the product of the distributions for the individual screens, such that
Estimating phase-error functions for multiple planes can significantly increase the number of unknowns that must be estimated relative to the number of measurements, . For screens, we must estimate unknowns for , plus unknowns for , giving a total of unknowns. To reduce the number of unknowns, we model the phase errors, , using a low-resolution Gaussian Markov random field (GMRF), denoted as [5]. If is subsampled by a factor of in both dimensions to obtain , each screen will have unknowns. The value of can be a function of the propagation plane, . When the same value is used for all planes, this technique reduces the total number of unknowns to .The Gaussian potential function for is given by
where is the difference between neighboring phase samples and controls the assumed variation in . This potential function enforces a strong influence between neighboring samples, which generates a more smoothly varying output when compared to the QGGMRF model.While the phase errors are estimated as a low-resolution function, , we must up-sample the estimate to the resolution of when applying the forward-model operator, . Our interpolation scheme is given by
where is an nearest-neighbor-interpolation matrix with elements in the set {0, 1}. More complex interpolations can be used; however, nearest neighbor allows for fast computation.C. MIR Framework
Using Eqs. (5)–(10), we can write the MAP cost function from Eq. (1) as
As noted in Ref. [5], the determinant and inverse in Eq. (11) make it difficult to evaluate; however, we may use the EM algorithm to replace Eq. (11) with a more tractable surrogate function.The EM algorithm represents a special case of majorization minimization using a surrogate function. For any function , we define a surrogate function, , to be an upper-bounding function, such that , where is the current value of that determines the functional form of , and is a constant that ensures the two functions are equal at . Surrogate functions have the property that minimization of implies minimization of [15].
The EM algorithm provides a formal framework for developing a surrogate function [17]. If we choose to be the unobserved data, our surrogate function becomes
where the conditional expectation is taken over , given the data, , the current estimate of the reflectance, , and the current estimate of the phase errors, . By introducing as the unobserved data, we create a cost function that is more tractable than Eq. (11) because of the linear relationship between the data, , and the reflection coefficient, . Taking the expectation helps overcome the fact that we do not know . Evaluation of the expectation in Eq. (12), with respect to the random vector, , constitutes the E-step of the EM algorithm. This step establishes the specific functional form of . Next, we conduct the M-step to jointly minimize according to The EM algorithm consists of iteratively applying the E- and M-steps shown in Eqs. (12) and (13) and updating and . For the optimization during the M-step, we use iterative coordinate descent (ICD) to sequentially minimize the entire surrogate cost function with respect to each element, and , where the joint subscripts , indicate the th sample of the th phase screen [15]. Figure 3 shows the steps of the MIR algorithm. In Appendix D, we provide the update equations for each step. The update for each reflectance element, , has a closed form solution and the update for each phase-error element, , requires a 1D line search. Notice that the separable form of derived in Appendix C ensures that sequential optimization of and is equivalent to joint optimization of these two quantities.Since is a surrogate function for , the alternating steps of the EM algorithm converge in a stable manner to the local minima of the original MAP cost function given by Eq. (11) [13]. Furthermore, since both the MAP cost function and its surrogate are nonconvex, the algorithm is likely to converge to a local minimum that is not global. Since we cannot directly evaluate the MAP cost function to compare different local minima, we cannot easily select among different local minima based on the values of their MAP cost function. Thus, to ensure we find a good local minimum, in the sense that it produces a focused and useful image, we pay close attention to our initial conditions, and we employ an iterative initialization process first developed in Ref. [13]. For a set number of outer-loop iterations, , we allow the basic MIR algorithm to run for iterations. Each time the MIR algorithm is called, is reinitialized, but the previous value of is used. After the outer loop runs times, we again run the basic MIR algorithm. For this final reconstruction run, we stop the algorithm when a metric, , falls below a threshold value of , where
3. RESULTS
In this section, we demonstrate the utility of our proposed MIR algorithm using both simulated and experimental data. In both cases, we compare the performance of MIR to the MIS algorithm from [11] over a range of anisoplanatic conditions. Additionally, we include reconstructions from the MBIR algorithm [5], which assumes IPEs. For the reader’s convenience, we provide an overview of the MIS algorithm in Appendix E.
A. Simulation
To generate synthetic data, we simulated the scenario shown in Fig. 1 in which atmospheric turbulence causes APEs. We used the split-step BPM MATLAB code from [8] to simulate the propagation of light from the object to the pupil plane. For atmospheric turbulence, we used three Kolmogorov phase screens of equal strength, distributed over the propagation path. We quantified the strength of the distributed-volume phase errors using two parameters, and . The first parameter, , is the ratio of the pupil diameter, , to the spherical-wave coherence length, [9]. High values of correspond to strong phase errors, and low values correspond to weak phase errors. The second parameter, , is the isoplanatic angle, , normalized by the diffraction-limited viewing angle, . As previously noted, provides a gauge for the image-domain angular separation over which the PSF is essentially unchanged [9]. Thus, high values of indicate neighboring points in the image have similar PSFs. As , points that are just barely resolved will have different PSFs, making reconstruction difficult. After propagating the fields to the pupil plane, we simulated the detection process for a DH system using an off-axis IPRG. The signal-to-noise ratio (SNR) was fixed for all data sets, and in Appendix F, we provide additional details about the simulations.
To numerically quantify algorithm performance in simulation, we used two metrics, normalized root mean square error (NRMSE) and peak Strehl ratio. NRMSE measures the distortion between the reconstructed images, , and the simulation input, . It is given by
where is the least-squares fit for any multiplicative scalar offset between and .The peak Strehl ratio, , measures the quality of our phase-error estimate, . It is important to note that we should not simply compare to the simulated phase errors, , directly. The imaging system is underdetermined, and there are infinitely many solutions that will produce the same effects. It is possible to obtain an estimate, where , that can still conjugate the effects of the distributed-volume phase errors and form a focused image. Therefore, to assess the quality of the reconstructed phase errors, we measure how well conjugates the effects of during backpropagation. Specifically, we backpropagate the pupil function, , with uniform phase, from to the object plane at , through a set of phase screens, . Our estimate, , should ideally cancel out the effects of during backpropagation. We define the peak Strehl ratio as
where indicates the backpropagation through a vacuum and indicates that we take the maximum value of the argument. In Eq. (17), the numerator is the maximum intensity in the object plane that occurs when we backpropagate the function , using to conjugate the effects of at each partial-propagation plane. The denominator represents the maximum intensity that occurs in the object plane when there are no turbulence effects (i.e., backpropagation through a vacuum). Thus, quantifies how close we are to the diffraction limit, ignoring the effects of tilt. Note that this Strehl-ratio metric only measures the phase-error conjugation for a single propagation angle; however, we observed that other angles have similar degrees of correction. Furthermore, we found this Strehl-ratio metric to be more useful at assessing the quality of compared to image-quality metrics, such as those used in Ref. [11].Using the methods described above, we generated 10 independent and identically distributed (i.i.d.) data realizations for 10 different values of , and three different values of . Each i.i.d. data realization had unique turbulence, speckle, and noise realizations. Figures 4–6 show example reconstructions. The images displayed were those with the median value of from the 10 i.i.d. realizations at each value of and . We also show the original uncompensated image in the left column. Figures 7 and 8 show and NRMSE, averaged over the 10 i.i.d. realizations, as a function of , for each . We include the Strehl ratio for the uncompensated phase errors, computed using Eq. (17) with , and the NRMSE for the original uncompensated image obtained according to .
The results show that the polynomial-based MIS algorithm is able to correct aberrations in the less challenging conditions when is low and is high. For a fixed value of , the MIS algorithm’s performance tapers off as decreases. For fixed values of , MIS performance falls off quickly as increases. Similar to MIS, for a fixed value of , the MIR performance tapers off as decreases. However, the MIR algorithm is able to achieve much higher peak Strehl ratios and lower NRMSE than the MIS algorithm in all scenarios. The significantly lower NRMSE produced by the MIR algorithm can be attributed to both a better estimate of and a reduction in speckle variation, as seen in Fig. 4. Unlike the MIS algorithm, MIR performance does not fall off as rapidly as increases for fixed values of . Lastly, we see that in most cases, the MBIR algorithm, which again assumes IPEs, performs as well or better than the MIS algorithm in terms of average Strehl ratio and outperforms the MIS algorithm in terms of average NRMSE. When visually compared to the MIR algorithm, the MBIR reconstructions look similar in the weaker turbulence cases and differ more drastically in the stronger turbulence cases. Upon close inspection of the weaker turbulence cases, we see more irregularities in the MBIR reconstructions than in the MIR reconstructions. This is due to the large difference in average Strehl ratios of the two algorithms, as shown in Fig. 7. The MBIR algorithm is simply not able to produce as accurate estimates of the APEs as the MIR algorithm. This large difference is not as apparent in the images, since some of the energy that gets spread out by the turbulence is removed during the MBIR regularization process.
To assess algorithm run time, we measured as a function of time for two different cases. In Case 1, the turbulence was relatively weak, with and . For Case 2, the turbulence was strong, with and . In both cases, we used MATLAB on a PC with a 2.2 GHz Intel Xeon processor. Figure 9 plots the results for both the MIS and MIR algorithms, and Fig. 10 shows the intermediate estimates, , as a function of time for Case 1. The plots show that in general, the MIR algorithm achieves higher values of much faster than the MIS algorithm. Additionally, we see that, for Case 1, the MIR estimate changes very little in the last 90% of the total run time; thus, we can shorten the run time and still get similar results. However, as noted in Appendix F, for simplicity we fixed the number MIR and MBIR outer-loop iterations at and the inner-loop iterations at for all cases. During the last outer loop, we let the algorithms run until the stopping criteria were met. On average, MIS ran 60 iterations in 67 min, MBIR ran 2565 iterations in 11 min, and MIR ran 2674 iterations in 28 min.
B. Experiment
To demonstrate the utility of the MIR algorithm with real data, we also conducted a laboratory experiment, as outlined in Appendix F. A laser with wavelength was used to illuminate an Air Force resolution chart 2.5 m away. To simulate APEs over the propagation path, we introduced multiple phase screens made from plastic CD jewel cases. For the experiment, we started with three screens located at , then incrementally added additional screens at distances 2.1 m, 2.0 m, and 1.9 m. However, to simplify reconstruction, we only estimated four screens located at for all data. Figure 11 shows the images reconstructed from the experimental data for three, four, five, and six plastic phase screens. The left column shows the original images reconstructed with no phase-error mitigation. Due to the nonuniformity of the plastic phase screens used, we see anisotropic (not to be confused with anisoplanatic) phase-error effects in the images. Here, the multiplane phase errors encountered in the experiment had stronger effects in the x axis than in the y axis. In Fig. 11, we also see that the phase errors were stronger in some areas of the image than others, making for a more challenging reconstruction problem.
Since the underlying statistics of the plastic phase screens used in the experiment were unknown, we cannot easily quantify performance in terms of atmospheric turbulence parameters such as and . However, we can compare against the simulated results to get a general idea of the experimental conditions. Doing so, we see that the overall turbulence strength most closely resembles our simulated data when . Furthermore, including the reconstructions from the MBIR algorithm provides a visual representation of the degree of anisoplanatic effects, since it assumes IPEs.
Figure 11 shows that the MBIR algorithm is not able to fully correct the multiplane phase errors encountered in the experiment. Alternatively, the MIS algorithm is able to improve upon the original image in all cases, although the sharpness is reduced as the number of screens increases. Similar to the simulations, the MIR algorithm is able to produce sharper and higher-contrast images with less speckle. Note that while the experimental phase errors were clearly anisoplanatic and anisotropic, the conditions were not as challenging as those produced in the simulation.
4. CONCLUSION
In summary, we have developed a novel algorithm for jointly estimating APEs and an object’s reflectance from single-shot DH data using a Bayesian framework. Our MIR algorithm models the APEs using multiple phase screens distributed along the propagation path and computes the MAP estimates of each phase screen along with the focused reflectance image.
To test our algorithm, we designed a simulation and experiment to produce data representative of a DH system using an off-axis IPRG. We generated data for a range of anisoplanatic conditions and compared the MIR algorithm to the polynomial-based MIS algorithm from [11]. To assess the results, we presented both qualitative and quantitative comparisons.
In all cases, we observed that the polynomial-based MIS algorithm works moderately well for the less severe APEs (e.g., when the value of is high and is low). However, the performance of the MIS algorithm falls off quickly as the APEs worsen (i.e., with decreases in and increases of ). Additionally, we observed that the MIR algorithm produces high-quality images and phase-error estimates over a wide range of anisoplanatic conditions. We therefore believe the MIR algorithm provides a robust alternative to the MIS algorithm for estimating APEs and for producing focused, speckle-reduced images from single-shot DH data.
APPENDIX A: PROPAGATION MODEL
In Eq. (2), the linear transform, , represents the propagation of light from the object to the observer through a number of intermediate planes containing phase screens. Using the split-step BPM found in Ref. [8], we model as a series of partial propagations between consecutive planes. This technique uses the angular spectrum method (ASM) to model each partial propagation followed by the application of a phase screen. In this section we first review the ASM, then we use the split-step BPM to provide an explicit definition for in Eq. (A9).
1. ASM
To model the vacuum propagation of a monochromatic field, , at some arbitrary distance, , we use the ASM [18]. Here, we use tildes to indicate continuous functions. The ASM computes the convolutional form of Fresnel diffraction integral by multiplying the input field in the spatial-frequency domain with the free-space transfer function. Note that the free-space transfer function is given by
where is the wavelength and and are the spatial frequency coordinates [18]. Thus, at a distance away, the diffraction pattern from is given by where and indicate the continuous-space Fourier transform (CSFT) operation and its inverse, respectively.In practice, we leverage the discrete version of Eq. (A2) found in Ref. [8]. Here, we write that discrete model, using vector notation, as
where is the vectorized input field, is the vectorized output field, and is a 2D DFT matrix. The diagonal matrix, , is the matrix form of the free-space transfer function and is given by Here, is an operator that forms a 1D vector from a 2D argument, is the sample spacing, is the number of samples in each dimension of the 2D transfer function, and are indices in the spatial-frequency domain. For simplicity, we have dropped any complex scalars and assumed a square sampling grid where .We can only use Eq. (A3) when the input- and output-field sample spacing and grid lengths are the same (i.e., ). Since this is restrictive, we use a modified version of the ASM found in Ref. [8] to control the ratio , and thus adjust the output grid size, . The modified propagation model is given by
where and Here, and are diagonal matrices that apply quadratic phase factors, is the ratio of the sample spacing in each grid, and are indices in the spatial domain. We add the subscript 0,1 to indicate that is the free-space-transfer function between planes and .2. Split-Step BPM Model
The split-step BPM model combines multiple partial propagations between phase screens using the modified ASM approach described by Eq. (A5). It accounts for differences between the object-plane sample spacing, , and the pupil-plane sample spacing, . If we convert the spit-step BPM operators found in Ref. [8] to vector notation, we can write the field passing through the pupil in the th propagation plane as
where and Here, is the sample spacing in the th plane, is the number of samples in each dimension of the 2D functions, are indices in the spatial domain, and are indices in the spatial-frequency domain. Additionally, in Eq. (A9), is the phase screen in the th plane, is the partial-propagation distance between plane and plane , and . We add the subscript to indicate that is the free-space-transfer function from planes to . As previously noted, is a diagonalization of the pupil-aperture vector, , and is a 2D DFT. Note that we have dropped any complex scalars and assumed a square sampling grid for simplicity (i.e., ).In Eq. (A9), we see that in order to model the propagation using , we must have knowledge of the geometry. Specifically, we must know the total distance, , and the sample spacings, and . In practice, these values can be determined from knowledge of the object range, sensor FOV, and pupil diameter. The partial-propagation distances can be an arbitrary set, depending on how many screens are to be modeled and their desired locations. Following [8], the values of through are set automatically from the values of and .
APPENDIX B: MEASUREMENT MODEL
In this section, we derive a discrete input–output relationship for a DH system using an off-axis IPRG. We start with the complex signal image, , a continuous quantity located at the focal-plane array (FPA) of the DH sensor. For the off-axis IPRG, we perform spatial-heterodyne detection by mixing the signal image with a reference field from a local oscillator (LO) that is split off from the master-oscillator laser. For a point source LO, located at point in the pupil plane, the reference field at the FPA is given by
where is the amplitude, assumed to be spatially uniform. For brevity, we ignore the quadratic phase factor, , applied to the reference, , by the imaging lens. Since this phase is applied to both the reference and the signal, it cancels out of Eq. (B2). In Eq. (B1), and are the spatial frequencies associated with the modulation from the off-axis reference, and is the focal length of the imaging lens. The resulting hologram intensity is given by where indicates the complex conjugate. Substituting Eq. (B1) into Eq. (B2) gives We model the discrete and noisy digital hologram, measured from the intensity in Eq. (B3), as where and are the sample spacing on the FPA and and are the sample indices in and dimensions. Here, we have simplified the model by ignoring the blurring effects of the finite-sized pixels. Furthermore, we assume that the DH system is shot-noise limited with noise , driven by the power of the reference field [14]. Thus, we model as additive, zero-mean, white Gaussian noise. Substituting Eq. (B3) into Eq. (B4) gives where is the digitized reference amplitude, and and are the modulation frequencies associated with the digitized reference field.Following the example shown in Fig. 2, we extract the signal by first taking a DFT of the digital hologram, . The resulting hologram spectrum is given by
Here, is the 2D DFT operator, is a discrete 2D delta function, denotes a discrete convolution, and are the spatial-frequency shifts corresponding to the modulation from the off-axis reference, , and is the DFT of the complex-signal image, . Due to the Fourier-transform relationship between the image and pupil planes, is a discrete representation of the signal field in the pupil plane [4].From Eq. (B6), we see that the modulation from the off-axis reference shifts the spectrum of the two cross terms away from the spectrum of the first two terms, which are centered at D.C. in . Thus, we can spatially extract a grid around the term of interest (fourth term) and divide by to get our signal data, such that
where is the extracted and scaled noise. Finally, we can represent Eq. (B7) using vector notation as where , , and are in the set and are obtained by vectorizing , , and , respectively. To obtain our final measurement model, given by Eq. (2), we substitute the discrete version of the propagated field, obtained using Eq. (A8), into Eq. (B8) for .APPENDIX C: SURROGATE FUNCTION
In this section, we derive the EM surrogate for the MAP cost function. We start by writing Eq. (12) as
where we have used Bayes’ theorem inside the expectation and the fact that . Next, we substitute in the likelihood and prior models from Eqs. (5)–(10). This gives where is a constant with respect to and .In Ref. [5], we showed that the conditional posterior distribution of is complex Gaussian with mean
and covariance Using Eqs. (C3) and (C4), we evaluate the expectations in Eq. (C2) to get the final form of our EM surrogate function. This gives us where is the th element of the posterior mean and is the th diagonal element of the posterior covariance.APPENDIX D: MIR ALGORITHM
In this section, we present details for executing each step of the MIR algorithm. The discussion describes how we evaluate the surrogate function in the E-step and how we optimize it in the M-step. We also present initialization and stopping procedures.
1. E-Step: Surrogate Function Formation
The E-step forms Eq. (C5) so that it can be minimized during the M-step. Specifically, this consists of computing the posterior mean and covariance according to Eqs. (C3) and (C4). We simplify this process by approximating the posterior covariance as
which requires that . The matrix can be decomposed into several matrices that are all unitary, with the exception of the , which applies the pupil-aperture function. Therefore, when and/or , the approximation of Eq. (D1) is exact. Typically, this is not the case, since is circular, but we approximate it as such to dramatically simplify computations. In practice, we have found this to work well. Alternatively, a square subset of the DH data can be used for which and Eq. (D1) is exact [6].2. M-Step: Surrogate Function Optimization
During the M-step, we minimize the surrogate function according to Eq. (13). Specifically, we employ iterative coordinate descent (ICD) to sequentially minimize Eq. (C5) with respect to each element, and . Here, the joint subscripts , indicate the th sample of the phase screen [15].
To update each reflectance element, , we ignore all terms in Eq. (C5) that do not contain . The resulting 1D cost function can be minimized with a line search over ; however, a closed-form solution can be obtained using the symmetric bound method of majorization described in Ref. [15]. Using this technique, we replace the QGGMRF potential function from Eq. (7) with a quadratic surrogate function. The resulting cost function is given by
where are the modified weights between neighboring pixel pairs given by [15] To obtain a closed-form solution for , we differentiate the right-hand side of Eq. (D2) with respect to , set it equal to zero, and multiply both sides by . The solution, , can then be found by rooting the polynomial , where Thus, we update each sample using the closed-form solution for the roots of a third-order polynomial. If multiple real roots exist, we chose the one with the lowest cost according to Eq. (D2).While the reflectance can be updated using a closed-form solution, the phase-error update requires a line search. To compactly express the cost function for the update of the phase-error element, , it is helpful to write the forward-model operator as
where represents all the matrices on the left-hand side of in Eq. (A9), and represents all the matrices on the right-hand side of . Using this notation, we can write the cost function for the th phase screen as If we consider just the terms in Eq. (D6), which depend on the element , we obtain the 1D cost function for that element given by where In Eq. (D7), is an index over neighboring low-resolution phase samples, is the set of vector indices of the high-resolution phase, , corresponding to the low-resolution element , and indicates the phase of the complex-valued scalar, . Intuitively, is the difference between the fields on either side of the screen , summed over all high-resolution points that correspond to the single low-resolution point being updated.The cost function of Eq. (D7) is the superposition of a cosine term, which has an infinite number of minima with equal cost, and a quadratic term. The global minima for the surrogate will be near the minimum of the quadratic term, bracketed on either side by the period of the cosine. Thus, we minimize Eq. (D7) using a 1D line search over , where minimizes the quadratic term in Eq. (D7). After repeating this for all , where is the set of samples in the lower-resolution phase-error grid, we obtain the full-resolution estimate of the unwrapped phase screen, , according to Eq. (10).
3. Initialization and Stopping Criteria
Figure 12 summarizes the steps of the basic MIR algorithm. The parameters and are set automatically, although the regularization of can be adjusted using a unitless parameter, . Note that while we initialize using a simple back-projection of the noisy data, we use an iterative process to initialize , based on the heuristic developed in Refs. [5,13].
Figure 13 details the steps of this iterative initialization process, referred to as the MIR initialization. The initial estimate of the phase-error vector is simply . Then, for a set number of outer-loop iterations, , we allow the basic MIR algorithm to run for iterations. Each time the MIR algorithm is called, is reinitialized, but the previous value of is used. After the outer loop runs times, we again run the basic MIR algorithm. For this final reconstruction run, we stop the algorithm when , in accordance with Eq. (14). It is important to note that while the MIR algorithm has many parameters, we set most of them automatically. In the methods section, we provide settings that we have found to robustly provide good results over a wide range of anisoplanatic conditions.
APPENDIX E: MULTIPLANE IS ALGORITHM
To investigate the performance gain of the MIR algorithm over conventional techniques, we considered two different MIS algorithms that differ in their representation of the phase errors, . A 15th-order polynomial is used for each phase screen in Ref. [11], while a point-by-point representation is used in Ref. [12]. Unfortunately, neither of these algorithms was designed for single-shot data, and only the polynomial-based algorithm produced usable results. Even with the bootstrapping initialization approach, we found that the point-by-point MIS algorithm produces poor results by over-sharpening single-shot images to a few points. Conversely, the polynomial-based MIS algorithm has fewer unknowns and is naturally constrained by the 15th-order representation. We found it to produce more useful results for single-shot data. Therefore, we chose the polynomial-based algorithm in Ref. [11] as our reference for comparisons.
The polynomial-based MIS algorithm computes the phase estimate according to
where is a matrix with columns corresponding to the first Zernike polynomials and indicates a matrix direct sum. The number of matrices direct-summed is equal to the number of phase-error planes modeled. In Eq. (E1), the value is an estimate of the polynomial coefficients for all three phase screens found according to Here, indicates the L1 norm, is an image-sharpness parameter, is a scalar that balances the two optimization terms, indicates a dependence of on , and is a vector of weights that applies a high-pass filter to the magnitude squared of the back-projected image spectrum, . The second term in Eq. (E2) is a regularization term that prevents object demagnification [11,12,20].Following the algorithm in Ref. [11], we use conjugate gradient (CG) to optimize Eq. (E2) in an iterative manner by first estimating only the coefficients corresponding to a third-order polynomial for all the screens, thenfourth, and so on, continuing up to the 15th order. We generated by first creating an inverse rectangular function with width (i.e., , where and are the spectral indices in the set ; then we applied a raised-cosine transition from the inner section of zeros to the outer section of ones. Through trial and error, we found to be optimal.
APPENDIX F: METHODS
In this section, we provide details about the methods used to produce data. We begin with a full description of the wave-optics simulation used to produce synthetic data and then describe our laboratory setup used to produce experimental data.
1. Simulation Details
Figure 1 of the paper shows the reflectance function, , used as the simulation input. We generated a reflection coefficient according to . Next, we propagated a distance, , through three phase screens, located at [8]. We generated the three phase screens using a Kolmogorov model for the refractive-index power spectral density (PSD) [8].
The strength of each individual screen was parameterized using the plane-wave coherence length, , also known as the Fried parameter [9]. Each of the three screens had equal values of . We quantified the overall strength of the distributed-volume phase errors with the parameter . Note that the spherical-wave coherence length is more appropriate than the plane-wave coherence length for distributed-volume phase errors. The value of specifies the degrees of freedom in the phase-error function at the pupil plane and allows us to quantify the cumulative effects of the distributive-volume phase errors. High values of correspond to strong phase errors and low values correspond to weak phase errors. For our experiment, we controlled by varying the value of for the three screens.
We quantified the anisoplanatic conditions with the parameter . For our simulations, we controlled by adjusting . Specifically, we adjusted the grid spacing in the object plane, , and set the corresponding pupil-plane sample grid spacing according to . The intermediate grid spacings, and , were set according to the split-step BPM MATLAB code from [8].
To simulate detection by a DH system using an off-axis IPRG, we padded the propagated field to obtain an array size of . Next, we applied a binary circular pupil-aperture function, , which had a circle of ones, 256 pixels wide, centered in an array of zeros. After, we first applied a thin-lens phase function that collimated the propagated light, then applied an FFT to form an image in the focal plane. Next, we mixed the image with a reference field given by Eq. (B1) and detected the resultant power. Figure 2(a) shows an example of the resulting digital hologram. The reference-beam power was set at approximately 80% of the well depth per detector element [i.e., 80% of ]. We also modeled Gaussian read noise with a standard deviation of 40 pe and digitized the output to 12 bits.
After detection, we took an FFT of the digital hologram data and isolated the signal of interest, which was a complex image of the signal in the pupil plane. Figure 2(b) shows the magnitude squared of the digital hologram spectrum. The signal of interest, , is highlighted by a white dotted line and is also shown in Fig. 2(c) after it is removed from the spectrum. The SNR was set to 100 for all data sets, where
and is the sample-variance operator.For simulated reconstructions using the MIR algorithm, we allowed the outer initialization loop to run times, with EM iterations each time. Once the iterative initialization process was complete, we set a stopping criterion of and let the EM algorithm run to completion. During initialization, we used , , , , and , where indicates a Gaussian kernel with standard deviation, . For the final MIR reconstruction, we set , , and . Additionally, for the phase-error prior model parameters, we used , , and set .
2. Laboratory Experiment Details
Figures 14 and 15 describe our experimental setup. For our master oscillator (MO) laser, we used a LightWave nonplanar ring oscillator with 700 mW of continuous-wave power and approximately 1 km of coherence length, with a wavelength of 1064 nm. Past a Faraday isolator and two alignment mirrors, we used a half-wave plate and polarized beam-splitting cube to create two optical legs. In the first optical leg, we coupled the beam from the MO laser into a single-mode, polarization-maintaining fiber to create an off-axis LO. This off-axis LO created a quasi-uniform reference beam with the appropriate tilt for digital-holographic detection in the off-axis IPRG. The mean strength of this reference beam was set to approximately 75% of the camera’s pixel-well depth. In the second optical leg, we focused the beam from the MO laser through a pinhole and then used an adjustable collimating lens to flood-illuminate the object after reflection from a steering mirror.
After setting up the imaging system with an object distance of approximately 2.5 m and an image distance of approximately 10 cm, we simulated the effects of anisoplanatic aberrations by placing six different phase screens in front of the imaging lens with a spacing of 10 cm in between. Each experimental phase screen contained two crisscrossed plastic plates (from CD jewel cases) placed back-to-back. For digital-holographic detection, we used a single lens to image the scattered light from the object onto a Guppy PRO f-125 FireWire camera. We then saved the recorded digital hologram in a , 16-bit TIFF file. Note that our recorded digital hologram had an SNR of approximately 10 according to Eq. (F1). Also note that we positioned the object, a chrome-on-glass Air Force resolution chart backed by Labsphere Spectralon, at an angle to minimize specular reflections and ensured that it was uniformly illuminated.
For all experimental reconstructions, MIS, MIR, and MBIR, we chose to digitally estimate four phase screens located at . Note that we observed better results when we included a phase screen at the aperture plane located at . We allowed the outer initialization loop to run times, with EM iterations each time. Once the iterative initialization process was complete, we set a stopping criterion of and let the EM algorithm run to completion. During initialization, we used , , , , and . For the final MIR reconstructions, we set , , , and . Additionally, for the phase-error prior model parameters, we used , and set . For the MBIR reconstructions, we set the parameters equal to the MIR parameters, except that only the phase screen was reconstructed.
Acknowledgment
The authors would like to thank B. T. Plimmer and D. E. Thornton for their assistance in conducting the laboratory experiments and in reviewing this paper.
REFERENCES
1. T.-C. Poon and J.-P. Liu, Introduction to Modern Digital Holography: With MATLAB (Cambridge University, 2014).
2. S. T. Thurman and J. R. Fienup, “Phase-error correction in digital holography,” J. Opt. Soc. Am. A 25, 983–994 (2008). [CrossRef]
3. J. C. Marron, R. L. Kendrick, N. Seldomridge, T. D. Grow, and T. A. Höft, “Atmospheric turbulence correction using digital holographic detection: experimental results,” Opt. Express 17, 11638–11651 (2009). [CrossRef]
4. M. F. Spencer, R. A. Raynor, M. T. Banet, and D. K. Marker, “Deep-turbulence wavefront sensing using digital-holographic detection in the off-axis image plane recording geometry,” Opt. Eng. 56, 031213 (2017). [CrossRef]
5. C. Pellizzari, M. F. Spencer, and C. A. Bouman, “Phase-error estimation and image reconstruction from digital-holography data using a Bayesian framework,” J. Opt. Soc. Am. A 34, 1659–1669 (2017). [CrossRef]
6. C. Pellizzari, M. T. Banet, M. F. Spencer, and C. A. Bouman, “Demonstration of single-shot digital holography using a Bayesian framework,” J. Opt. Soc. Am. A 35, 103–107 (2018). [CrossRef]
7. J. W. Goodman, Speckle Phenomena in Optics: Theory and Applications (Roberts and Company, 2006).
8. J. D. Schmidt, Numerical Simulation of Optical Wave Propagation with Examples in MATLAB (SPIE, 2010).
9. L. Andrews and R. Phillips, Laser Beam Propagation through Random Media, 2nd ed. (SPIE, 2005).
10. S. T. Thurman and J. R. Fienup, “Correction of anisoplanatic phase errors in digital holography,” J. Opt. Soc. Am. A 25, 995–999 (2008). [CrossRef]
11. A. E. Tippie and J. R. Fienup, “Phase-error correction for multiple planes using a sharpness metric,” Opt. Lett. 34, 701–703 (2009). [CrossRef]
12. A. E. Tippie and J. R. Fienup, “Multiple-plane anisoplanatic phase correction in a laboratory digital holography experiment,” Opt. Lett. 35, 3291–3293 (2010). [CrossRef]
13. C. J. Pellizzari, R. Trahan, H. Zhou, S. Williams, S. E. Williams, B. Nemati, M. Shao, and C. A. Bouman, “Synthetic aperature ladar: a model-based approach,” IEEE Trans. Comput. Imaging 3, 901–916 (2017). [CrossRef]
14. V. V. Protopopov, Laser Heterodyning, Vol. 149 of Springer Series in Optical Sciences (Springer, 2009).
15. C. A. Bouman, Model Based Image Processing, 2013, https://engineering.purdue.edu/~bouman/publications/pdf/MBIP-book.pdf.
16. J. B. Thibault, K. Sauer, C. Bouman, and J. Hsieh, “A three-dimensional statistical approach to improved image quality for multi-slice helical CT,” Med. Phys. 34, 4526–4544 (2007). [CrossRef]
17. A. P. Dempster, N. M. Laird, and D. B. Rubin, “Maximum likelihood from incomplete data via the EM algorithm,” J. R. Stat. Soc. Ser. B 39, 1–38 (1977).
18. J. W. Goodman, Introduction to Fourier Optics, 3rd ed. (Roberts and Company, 2005).
19. C. Forbes, M. Evans, N. Hastings, and B. Peacock, Statistical Distributions (Wiley, 2011).
20. A. E. Tippie, “Aberration correction in digital holography,” Ph.D. dissertation (University of Rochester, 2012).