## 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, $g$, 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, $\varphi $, and the real-valued reflectance, $r$, from a single data realization—a process referred to here as single-shot DH [5]. We define the reflectance as $r=E[{|g|}^{2}]$, where $E[\xb7]$ indicates the expected value [7]. In general, $r$ is smoother and has higher spatial correlation when compared to the reflection coefficient, $g$. Reconstructing $r$, rather than $g$, 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, ${\theta}_{0}$ [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 ${\theta}_{0}$ 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, $r$. 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, $r$, and a corresponding reflection coefficient, $g$, through distributed-volume phase errors, $\varphi $. 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, $\varphi $, as a finite number of discrete layers, $[{\varphi}_{1},{\varphi}_{2},\dots ,{\varphi}_{K}]$, where each layer is a thin unit-amplitude phase screen that represents the phase perturbations along an incremental volume [8]. After encountering $\varphi $, 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, $y$, from the spectrum of the digital hologram. Due to the discrete Fourier transform (DFT) relationship between the pupil plane and the image/detector plane, $y$ 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, $y$. 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, $r\in {\mathbb{R}}^{M}$, and the phase errors, $\varphi =[{\varphi}_{1},{\varphi}_{2},\dots ,{\varphi}_{K}]$, where ${\varphi}_{k}\in {\mathbb{R}}^{M}\text{\hspace{0.17em}}\forall \text{\hspace{0.17em}}k$, from a single noisy data realization, $y\in {\mathbb{C}}^{M}$. Here, we use vectorized notation for all 2D variables. The joint MAP estimates are given by

#### 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, $g\in {\mathbb{C}}^{N}$, plus complex-valued measurement noise, $w\in {\mathbb{C}}^{M}$ [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, $y$, obtained from a DH system using an off-axis IPRG. The transform, ${A}_{\varphi}\in {\mathbb{C}}^{M\times M}$, accounts for the propagation and measurement geometry and is dependent on $\varphi $.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 ${A}_{\varphi}$. With this in mind, we use the split-step beam propagation method (BPM) to model ${A}_{\varphi}$ 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 $K$ phase screens, we get $\varphi \in {\mathbb{R}}^{KM}$. This makes the joint estimation of $r$ and $\varphi $ significantly more challenging, since the number of unknowns increases relative to the measurements.

To determine the likelihood function, $p(y|r,\varphi )$, we first define distributions for $g$ and $w$; then we relate these distributions to $y$ using Eq. (2). For a surface with reflectance, $r$, which is rough relative to the illumination wavelength, we model the vector $g$ as

where $\mathcal{D}(\xb7)$ denotes an operator that produces a diagonal matrix from its vector argument and $CN(\mu ,C,\mathrm{\Gamma})$ indicates a multivariate complex normal distribution with mean, $\mu $, covariance matrix, $C$, and pseudo-covariance matrix, $\mathrm{\Gamma}$ [5]. Here, $\mathrm{\Gamma}=0$, since the phase of $g$ 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 $w$ is therefore given by where ${\sigma}_{w}^{2}$ is the noise variance. Finally, using the data model from Eq. (2) and the distributions from Eqs. (3) and (4), the likelihood function for $y$, given $r$ and $\varphi $, becomes where the superscript $H$ indicates the Hermitian transpose.#### B. Prior Models

Both the reflectance, $r$, and the phase errors, ${\varphi}_{k}$, are modeled using Markov random fields (MRFs) with the form

### 1. Reflectance

For the reflectance function, $r$, we used a Q-generalized Gaussian Markov random field (QGGMRF) potential function with the form

### 2. Phase Errors

We assume that the layers of the distributed-volume phase error, $\varphi $, are statistically independent [8]. Therefore, we can write the distribution for $\varphi $ 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, $M$. For $K$ screens, we must estimate $KM$ unknowns for $\varphi $, plus $M$ unknowns for $r$, giving a total of $M(K+1)$ unknowns. To reduce the number of unknowns, we model the phase errors, ${\varphi}_{k}$, using a low-resolution Gaussian Markov random field (GMRF), denoted as ${\overline{\varphi}}_{k}$ [5]. If ${\varphi}_{k}$ is subsampled by a factor of ${n}_{b}$ in both dimensions to obtain ${\overline{\varphi}}_{k}$, each screen will have $M/{n}_{b}^{2}$ unknowns. The value of ${n}_{b}$ can be a function of the propagation plane, $k$. When the same value is used for all planes, this technique reduces the total number of unknowns to $M(K/{n}_{b}^{2}+1)$.The Gaussian potential function for ${\overline{\varphi}}_{k}$ is given by

While the phase errors are estimated as a low-resolution function, $\overline{\varphi}$, we must up-sample the estimate to the resolution of $\varphi $ when applying the forward-model operator, ${A}_{\varphi}$. Our interpolation scheme is given by

where $P$ is an $M\times M/{n}_{b}^{2}$ 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

The EM algorithm represents a special case of majorization minimization using a surrogate function. For any function $c(x)$, we define a surrogate function, $Q(x;{x}^{\prime})$, to be an upper-bounding function, such that $c(x)\le Q(x;{x}^{\prime})+\kappa $, where ${x}^{\prime}$ is the current value of $x$ that determines the functional form of $Q$, and $\kappa $ is a constant that ensures the two functions are equal at ${x}^{\prime}$. Surrogate functions have the property that minimization of $Q(x;{x}^{\prime})$ implies minimization of $c(x)$ [15].

The EM algorithm provides a formal framework for developing a surrogate function [17]. If we choose $g$ to be the unobserved data, our surrogate function becomes

Since $Q(r,\overline{\varphi};{r}^{\prime},{\overline{\varphi}}^{\prime})$ is a surrogate function for $c(r,\overline{\varphi})$, 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, ${N}_{L}$, we allow the basic MIR algorithm to run for ${N}_{K}$ iterations. Each time the MIR algorithm is called, $r$ is reinitialized, but the previous value of $\varphi $ is used. After the outer loop runs ${N}_{L}$ times, we again run the basic MIR algorithm. For this final reconstruction run, we stop the algorithm when a metric, $\epsilon $, falls below a threshold value of ${\epsilon}_{T}$, 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, $D/{r}_{0,sw}$ and ${\tilde{\theta}}_{0}={\theta}_{0}/(\lambda /D)$. The first parameter, $D/{r}_{0,sw}$, is the ratio of the pupil diameter, $D$, to the spherical-wave coherence length, ${r}_{0,sw}$ [9]. High values of $D/{r}_{0,sw}$ correspond to strong phase errors, and low values correspond to weak phase errors. The second parameter, ${\tilde{\theta}}_{0}$, is the isoplanatic angle, ${\theta}_{0}$, normalized by the diffraction-limited viewing angle, $\lambda /D$. As previously noted, ${\theta}_{0}$ provides a gauge for the image-domain angular separation over which the PSF is essentially unchanged [9]. Thus, high values of ${\tilde{\theta}}_{0}$ indicate neighboring points in the image have similar PSFs. As ${\tilde{\theta}}_{0}\to 1$, 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, $\widehat{r}$, and the simulation input, $r$. It is given by

The peak Strehl ratio, ${S}_{p}$, measures the quality of our phase-error estimate, $\widehat{\varphi}$. It is important to note that we should not simply compare $\widehat{\varphi}$ to the simulated phase errors, $\varphi $, 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 $\widehat{\varphi}\ne \varphi $, 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 $\widehat{\varphi}$ conjugates the effects of $\varphi $ during backpropagation. Specifically, we backpropagate the pupil function, $a$, with uniform phase, from $k=3$ to the object plane at $k=0$, through a set of phase screens, ${\varphi}_{c}=\varphi -\widehat{\varphi}$. Our estimate, $\widehat{\varphi}$, should ideally cancel out the effects of $\varphi $ during backpropagation. We define the peak Strehl ratio as

Using the methods described above, we generated 10 independent and identically distributed (i.i.d.) data realizations for 10 different values of ${\tilde{\theta}}_{0}\in [100,1]$, and three different values of $D/{r}_{0,sw}\in \{5,10,15\}$. 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 ${S}_{p}$ from the 10 i.i.d. realizations at each value of ${\tilde{\theta}}_{0}$ and $D/{r}_{0,sw}$. We also show the original uncompensated image in the left column. Figures 7 and 8 show ${S}_{p}$ and NRMSE, averaged over the 10 i.i.d. realizations, as a function of ${\tilde{\theta}}_{0}$, for each $D/{r}_{0,sw}$. We include the Strehl ratio for the uncompensated phase errors, computed using Eq. (17) with $\widehat{\varphi}=0$, and the NRMSE for the original uncompensated image obtained according to $\widehat{r}={|{A}_{\varphi =0}^{H}y|}^{\circ 2}$.

The results show that the polynomial-based MIS algorithm is able to correct aberrations in the less challenging conditions when $D/{r}_{0,sw}$ is low and ${\tilde{\theta}}_{0}$ is high. For a fixed value of $D/{r}_{0,sw}$, the MIS algorithm’s performance tapers off as ${\tilde{\theta}}_{0}$ decreases. For fixed values of ${\tilde{\theta}}_{0}$, MIS performance falls off quickly as $D/{r}_{0,sw}$ increases. Similar to MIS, for a fixed value of $D/{r}_{0,sw}$, the MIR performance tapers off as ${\tilde{\theta}}_{0}$ 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 $\widehat{\varphi}$ and a reduction in speckle variation, as seen in Fig. 4. Unlike the MIS algorithm, MIR performance does not fall off as rapidly as $D/{r}_{0,sw}$ increases for fixed values of ${\tilde{\theta}}_{0}$. 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 ${S}_{p}$ as a function of time for two different cases. In Case 1, the turbulence was relatively weak, with $D/{r}_{0,sw}=5$ and ${\tilde{\theta}}_{0}=100$. For Case 2, the turbulence was strong, with $D/{r}_{0,sw}=15$ and ${\tilde{\theta}}_{0}=2.8$. 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, $\widehat{r}$, as a function of time for Case 1. The plots show that in general, the MIR algorithm achieves higher values of ${S}_{p}$ 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 ${N}_{L}=50$ and the inner-loop iterations at ${N}_{K}=20$ 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 $\lambda =1064\text{\hspace{0.17em}}\mathrm{nm}$ 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 $z=[2.2,2.3,2.4]\text{\hspace{0.17em}}\mathrm{m}$, 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 $z=[1.9,2.1,2.3,2.5]\text{\hspace{0.17em}}\mathrm{m}$ 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 $D/{r}_{0,sw}$ and ${\tilde{\theta}}_{0}$. 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 $D/{r}_{0,sw}=5$. 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 ${\tilde{\theta}}_{0}$ is high and $D/{r}_{0,sw}$ is low). However, the performance of the MIS algorithm falls off quickly as the APEs worsen (i.e., with decreases in ${\tilde{\theta}}_{0}$ and increases of $D/{r}_{0,sw}$). 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, ${A}_{\varphi}$, 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 ${A}_{\varphi}$ 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 ${A}_{\varphi}$ in Eq. (A9).

## 1. ASM

To model the vacuum propagation of a monochromatic field, ${\tilde{u}}_{0}({x}_{0},{y}_{0})$, at some arbitrary distance, $\mathrm{\Delta}z$, 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

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 ${u}_{0}\in {\mathbb{C}}^{M}$ is the vectorized input field, ${u}_{1}\in {\mathbb{C}}^{M}$ is the vectorized output field, and $F\in {\mathbb{C}}^{M\times M}$ is a 2D DFT matrix. The diagonal matrix, $H$, is the matrix form of the free-space transfer function and is given byWe can only use Eq. (A3) when the input- and output-field sample spacing and grid lengths are the same (i.e., ${\delta}_{1}={\delta}_{0}$). Since this is restrictive, we use a modified version of the ASM found in Ref. [8] to control the ratio ${\beta}_{1}={\delta}_{1}/{\delta}_{0}$, and thus adjust the output grid size, ${\delta}_{1}N$. The modified propagation model is given by

where## 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, ${\delta}_{0}$, and the pupil-plane sample spacing, ${\delta}_{K}$. 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 $K$th propagation plane as

whereIn Eq. (A9), we see that in order to model the propagation using ${A}_{\varphi}$, we must have knowledge of the geometry. Specifically, we must know the total distance, $\mathrm{\Delta}z$, and the sample spacings, ${\delta}_{0}$ and ${\delta}_{K}$. 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 ${\delta}_{1}$ through ${\delta}_{K-1}$ are set automatically from the values of ${\delta}_{0}$ and ${\delta}_{K}$.

## 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, ${\tilde{u}}_{S}(x,y)$, 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 $({x}_{p}={r}_{x},{y}_{p}={r}_{y})$ in the pupil plane, the reference field at the FPA is given by

where ${\tilde{R}}_{0}$ is the amplitude, assumed to be spatially uniform. For brevity, we ignore the quadratic phase factor, ${e}^{j\frac{{\pi}^{2}}{\lambda f}({x}^{2}+{y}^{2})}$, applied to the reference, ${u}_{R}(x,y)$, 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), ${f}_{{r}_{x}}={r}_{x}/(\lambda f)$ and ${f}_{{r}_{y}}={r}_{y}/(\lambda f)$ are the spatial frequencies associated with the modulation from the off-axis reference, and $f$ is the focal length of the imaging lens. The resulting hologram intensity is given byFollowing the example shown in Fig. 2, we extract the signal by first taking a DFT of the digital hologram, $i(m,n)$. The resulting hologram spectrum is given by

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 $I(p,q)$. Thus, we can spatially extract a grid around the term of interest (fourth term) and divide by ${R}_{0}$ to get our signal data, such that

where $w(p,q)$ is the extracted and scaled noise. Finally, we can represent Eq. (B7) using vector notation as where $y$, $u$, and $w$ are in the set ${\mathbb{C}}^{M}$ and are obtained by vectorizing $y(p,q)$, ${U}_{s}(p,q)$, and $w(p,q)$, 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 $U$.## APPENDIX C: SURROGATE FUNCTION

In this section, we derive the EM surrogate for the MAP cost function. We start by writing Eq. (12) as

In Ref. [5], we showed that the conditional posterior distribution of $g$ is complex Gaussian with mean

and 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

## 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, ${r}_{s}$ and ${\overline{\varphi}}_{k,s}$. Here, the joint subscripts $k$, $s$ indicate the $s$th sample of the ${k}^{th}$ phase screen [15].

To update each reflectance element, ${r}_{s}$, we ignore all terms in Eq. (C5) that do not contain ${r}_{s}$. The resulting 1D cost function can be minimized with a line search over ${\mathbb{R}}^{+}$; 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

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, ${\overline{\varphi}}_{k,s}$, it is helpful to write the forward-model operator as

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 ${\overline{\varphi}}_{k,s}\in [{\overline{\varphi}}_{k,s}^{**}-\pi ,{\overline{\varphi}}_{k,s}^{**}+\pi ]$, where ${\overline{\varphi}}_{k,s}^{**}$ minimizes the quadratic term in Eq. (D7). After repeating this for all $s\in \overline{S}$, where $\overline{S}$ is the set of samples in the lower-resolution phase-error grid, we obtain the full-resolution estimate of the unwrapped phase screen, ${\widehat{\varphi}}_{k}$, according to Eq. (10).

## 3. Initialization and Stopping Criteria

Figure 12 summarizes the steps of the basic MIR algorithm. The parameters ${\sigma}_{r}$ and ${\sigma}_{w}$ are set automatically, although the regularization of $r$ can be adjusted using a unitless parameter, $\gamma $. Note that while we initialize $r$ using a simple back-projection of the noisy data, we use an iterative process to initialize $\varphi $, 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 ${\varphi}^{0}=\mathbf{0}$. Then, for a set number of outer-loop iterations, ${N}_{L}$, we allow the basic MIR algorithm to run for ${N}_{K}$ iterations. Each time the MIR algorithm is called, $r$ is reinitialized, but the previous value of $\varphi $ is used. After the outer loop runs ${N}_{L}$ times, we again run the basic MIR algorithm. For this final reconstruction run, we stop the algorithm when $\epsilon <{\epsilon}_{T}$, 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, $\varphi $. 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 $Z\in {\mathbb{R}}^{M\times {N}_{z}}$ is a matrix with columns corresponding to the first ${N}_{z}$ Zernike polynomials and $\oplus $ indicates a matrix direct sum. The number of $Z$ matrices direct-summed is equal to the number of phase-error planes modeled. In Eq. (E1), the value $\widehat{c}\in {\mathbb{R}}^{3{N}_{z}}$ is an estimate of the polynomial coefficients for all three phase screens found according toFollowing 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 $W$ by first creating an inverse rectangular function with width $0.8N$ (i.e., $(1-\mathrm{rect}[p/0.8N])(1-\mathrm{rect}[q/0.8N])$, where $p$ and $q$ are the spectral indices in the set $[-N/2,N/2-1)$; 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 $\alpha =1\times {10}^{-6}$ 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 $256\times 256$ reflectance function, $r$, used as the simulation input. We generated a reflection coefficient according to $g\sim CN(0,\mathcal{D}(r))$. Next, we propagated $g$ a distance, $\mathrm{\Delta}z=80\text{\hspace{0.17em}}\mathrm{m}$, through three phase screens, located at $z=[26.4,52.8,80]\text{\hspace{0.17em}}\mathrm{m}$ [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, ${r}_{0,pw}$, also known as the Fried parameter [9]. Each of the three screens had equal values of ${r}_{0,pw}$. We quantified the overall strength of the distributed-volume phase errors with the parameter $D/{r}_{0,sw}$. Note that the spherical-wave coherence length is more appropriate than the plane-wave coherence length for distributed-volume phase errors. The value of $D/{r}_{0,sw}$ 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 $D/{r}_{0,sw}$ correspond to strong phase errors and low values correspond to weak phase errors. For our experiment, we controlled $D/{r}_{0,sw}$ by varying the value of ${r}_{0,pw}\in [3\times {10}^{-4},10]\text{\hspace{0.17em}}\mathrm{m}$ for the three screens.

We quantified the anisoplanatic conditions with the parameter ${\tilde{\theta}}_{0}={\theta}_{0}/(\lambda /D)$. For our simulations, we controlled ${\tilde{\theta}}_{0}$ by adjusting $(\lambda /D)$. Specifically, we adjusted the grid spacing in the object plane, ${\delta}_{0}\in [1\times {10}^{-4},1.2\times {10}^{-2}]\text{\hspace{0.17em}}\mathrm{m}$, and set the corresponding pupil-plane sample grid spacing according to ${\delta}_{3}=\lambda \mathrm{\Delta}z{\delta}_{0}/256\text{\hspace{0.17em}}\mathrm{m}$. The intermediate grid spacings, ${\delta}_{1}$ and ${\delta}_{2}$, 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 $512\times 512$. Next, we applied a $512\times 512$ binary circular pupil-aperture function, $a$, 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 $5\times {10}^{4}\text{\hspace{0.17em}}\text{photoelectrons}(\mathrm{pe})$]. 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 $256\times 256$ complex image of the signal in the pupil plane. Figure 2(b) shows the magnitude squared of the digital hologram spectrum. The $256\times 256$ signal of interest, $y$, 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 ${s}^{2}(\xb7)$ is the sample-variance operator.For simulated reconstructions using the MIR algorithm, we allowed the outer initialization loop to run ${N}_{L}=50$ times, with ${N}_{K}=20$ EM iterations each time. Once the iterative initialization process was complete, we set a stopping criterion of ${\u03f5}_{T}=1\times {10}^{-4}$ and let the EM algorithm run to completion. During initialization, we used $q=2$, $p=2$, $T=1$, $\gamma =2$, and $b=G(0.8)$, where $G(\sigma )$ indicates a $3\times 3$ Gaussian kernel with standard deviation, $\sigma $. For the final MIR reconstruction, we set $p=1.1$, $T=0.1$, and $b=G(0.1)$. Additionally, for the phase-error prior model parameters, we used $b=G(0.1)$, ${\sigma}_{\overline{\varphi}}=[0.05,0.10,0.15]$, and set ${n}_{b}=4$.

## 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 $512\times 512$, 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 $z=[1.9,2.1,2.3,2.5]\text{\hspace{0.17em}}\mathrm{m}$. Note that we observed better results when we included a phase screen at the aperture plane located at $z=2.5\text{\hspace{0.17em}}\mathrm{m}$. We allowed the outer initialization loop to run ${N}_{L}=100$ times, with ${N}_{K}=50$ EM iterations each time. Once the iterative initialization process was complete, we set a stopping criterion of ${\u03f5}_{T}=1\times {10}^{-4}$ and let the EM algorithm run to completion. During initialization, we used $q=2$, $p=2$, $T=1$, $\gamma =2$, and $b=G(0.8)$. For the final MIR reconstructions, we set $p=1.1$, $T=0.05$, $\gamma =2.5$, and $b=G(0.1)$. Additionally, for the phase-error prior model parameters, we used ${\sigma}_{\overline{\varphi}}=[0.15,0.30,0.45,0.60]$, and set ${n}_{b}=[10,10,5,5]$. For the MBIR reconstructions, we set the parameters equal to the MIR parameters, except that only the $z=2.5\text{\hspace{0.17em}}\mathrm{m}$ 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).