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[·]$ 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.

 

Fig. 1. Example DH system using an off-axis IPRG. A laser source is used to flood-illuminate an object that has a reflectance function, $r$, and corresponding reflection coefficient, $g$. In this example, the return signal is corrupted by atmospheric turbulence, which causes APEs, $\varphi$. We model $\varphi$ as a layered atmosphere consisting of multiple phase screens, $\varphi =[{\varphi }_1,{\varphi }_2,\dots ,{\varphi }_K]$, where $K=3$ in this example. In our discrete model, each propagation plane, $k$, has an associated grid sample spacing, ${\delta }_k$, and the distance between grids is given by $\mathrm{\Delta }{z}_i$. The blue and green dotted cones represent ray traces from two different points on the object. The ray traces show that light from different points on the object will encounter different phase errors, resulting in anisoplanatic conditions. The returning signal passes through the pupil-aperture function, $a$, is collimated, and is imaged onto a focal-plane detector array, where it is mixed with a reference field to form the digital hologram.

Download Full Size | PDF

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.

 

Fig. 2. Conventional image-processing steps for a DH system using an off-axis IPRG. Here, we show the magnitude squared of any complex-valued quantities. To extract the signal, we start with (a) a real-valued digital hologram, then take a 2D FFT to obtain the (b) complex-valued holographic spectrum. Note that the digital hologram, shown in (a), naturally has low contrast. Next, a section of the spectrum is isolated (c), as indicated by the dashed white line. This subset of the spectrum, $y$, represents a complex image of the signal field in the pupil plane. Finally, for basic image formation, we take an inverse FFT (IFFT) of $y$ to form the image shown in (d).

Download Full Size | PDF

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\forall 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,

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

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

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

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

 

Fig. 3. MIR algorithm for the MAP estimates of $r$ and $\varphi$. Here, $S$ is the set of all samples in $r$, $K$ is the number of partial-propagation planes, $\overline{S}$ is the set of all samples in the low-resolution phase screens, and $i$ is the iteration index.

Download Full Size | PDF

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}$.

 

Fig. 4. Example single-shot reconstructions (simulated) for $D/r_{0,sw}=5$. Images are displayed on a log-based decibel (dB) scale obtained according to $r_{\mathrm{db}}=10\log \{r/\max (r)\}$, where $\max (·)$ indicates the maximum value in the 2D array.

Download Full Size | PDF

 

Fig. 5. Example of single-shot reconstructions (simulated) for $D/r_{0,sw}=10$.

Download Full Size | PDF

 

Fig. 6. Example of single-shot reconstructions (simulated) for $D/r_{0,sw}=15$.

Download Full Size | PDF

 

Fig. 7. Average peak Strehl ratio, $S_p$, as a function of ${\tilde{\theta }}_0$ for $D/r_{0,sw}=5$ (top), $D/r_{0,sw}=10$ (middle), and $D/r_{0,sw}=15$ (bottom). The error bars show the standard deviation at each point.

Download Full Size | PDF

 

Fig. 8. Average NRMSE as a function of ${\tilde{\theta }}_0$ for $D/r_{0,sw}=5$ (top), $D/r_{0,sw}=10$ (middle), and $D/r_{0,sw}=15$ (bottom). The error bars show the standard deviation at each point.

Download Full Size | PDF

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.

 

Fig. 9. Peak Strehl ratio as a function of run time for both weak (blue) and strong (red) APEs.

Download Full Size | PDF

 

Fig. 10. Example of the intermediate samples of $\widehat{r}$ taken at various times during the algorithm run time corresponding to Case 1 in Fig. 9. The top images are the MIS estimates at $t=0,341,706,2086,3134$, and 4501 s. The bottom images are the MIR estimates at $t=0,9.9,21.4,33.2,163$, and 1324 s (ordered from left to right first, then top to bottom).

Download Full Size | PDF

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\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]\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]\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.

 

Fig. 11. Experimental results for three, four, five, and six plastic phase screens. The original reconstruction with no phase-error compensation is shown in the left column, the MIS reconstructions are in the second column, and the MIR reconstructions are in the third column. Additionally, we show reconstructions from the isoplanatic MBIR algorithm.

Download Full Size | PDF

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.