Phaseerror estimation and image reconstruction from digitalholography data using a Bayesian framework
Abstract
The estimation of phase errors from digitalholography data is critical for applications such as imaging or wavefront sensing. Conventional techniques require multiple i.i.d. data and perform poorly in the presence of high noise or large phase errors. In this paper we propose a method to estimate isoplanatic phase errors from a single data realization. We develop a modelbased iterative reconstruction algorithm which computes the maximum a posteriori estimate of the phase and the specklefree object reflectance. Using simulated data, we show that the algorithm is robust against high noise and strong phase errors.
josaa
1 Introduction
Digital holography (DH) uses coherentlaser illumination and detection to gain significant improvements in sensitivity for remotesensing applications. Detection involves measuring the modulation of a strong reference field by a potentiallyweak signal field. This modulation allows for the detection of signals with energies equivalent to a single photon or less [1].
In practice, DH systems are sensitive to phase errors imparted by atmospheric perturbations or flaws in the optical system. For imaging applications, we must estimate and remove these phase errors to form a focused image [2, 3, 4, 5, 6, 7, 8, 9, 10]. Similarly, for wavefront sensing applications, estimation of the phase errors represents the desired sensor output [11, 12, 13, 14]. In either case, accurately estimating the phase errors is a necessary and critical task.
Stateoftheart techniques for estimating phase errors from DH data involve maximizing an image sharpness metric with respect to the phase errors [3, 4, 5, 7, 8, 9, 10]. Typically, multiple data realizations are obtained for which the imagetoimage speckle variation is decorrelated but the phase errors are identical. For applications using pupilplane detection, the data realizations are corrected using an initial estimate of the phase errors then inverted using a Fast Fourier Transform (FFT). The inversion produces estimates of the complexvalued reflection coefficient, , which are incoherently averaged to reduce speckle variations. The sharpness of this speckleaveraged image is then maximized to obtain the final phaseerror estimate.
Image sharpening (IS) algorithms are designed to use multiple data realizations. However, Thurman and Fienup demonstrated that in favorable conditions, the algorithms can still be used when only one realization is available [3]. For a high signaltonoise ratio (SNR), with several hundred detected photons per detector element, their algorithm was able to estimate atmospheric phase functions with residual rootmeansquare (RMS) errors as low as .
While IS algorithms have been successfully demonstrated for estimating both isoplanatic and anisoplanatic phase errors, there remains room for improvement. First, IS algorithms use relatively simple inversion techniques to reconstruct the magnitude squared of the reflection coefficient, , rather than the realvalued reflectance, r. The reflectance, given by , where indicates the expected value, is a smoother quantity with higher spatial correlation between elements in the image. We are accustomed to seeing the reflectance in conventional images and it is of greater interest for many imaging applications. Conversely, reconstructing , leads to images with highspatialfrequency variations known as speckle. By reconstructing , we can leverage its higher spatial correlation to better constrain the estimation process and potentially produce more accurate estimates of the phase errors [15].
Another limitation of IS algorithms is that the process of estimating the phase errors is not tightly coupled to image reconstruction. Information about the object can be further leveraged during the phaseestimation process. For example, this information can help rule out contributions to the phase caused by noise. By jointly estimating both the phase errors and the reflectance, we can incorporate additional information into our estimates which helps reduce uncertainty [15].
Lastly, in many practical applications it may not be possible to obtain multiple data realizations. Particularly, in cases where there is not enough relative movement between the object and receiver to decorrelate the speckle variations between images, the phase errors are rapidly changing, or when timing requirements prohibit multiple images from being obtained. Without multiple realizations, IS algorithms can struggle to produce accurate phase error estimates [3] and an alternate method is needed.
In this paper, we propose an improved method of DH phase recovery based on a framework of joint image reconstruction and phaseerror estimation. Our approach builds on the work of [15] by reconstructing the spatiallycorrelated reflectance, , rather than the reflection coefficient, , which allows for better estimates of the phase errors from less data. We focus on reconstructing from a single data realization under isoplanatic atmospheric conditions for the offaxis pupil plane recording geometry (PPRG) [13]. Our major contributions include:

We jointly compute the maximum a posteriori (MAP) estimates of the image reflectance, , and phase errors, from a single data realization. Joint estimation reduces uncertainty in both quantities. Additionally, estimating rather than helps to further constrain the problem and produce images without speckle.

We derive the forward model for a DH system using the offaxis PPRG. The model ensures our estimates are consistent with the physics and statistics of the remote sensing scenario.

We develop a technique to compute the 2D phase errors that occur in DH. Our approach allows for the estimation of phase errors on a lower resolution grid to help reduce the number of unknowns.

We model the phase errors as a random variable using a Gaussian Markov Random Field (GMRF) prior model. This approach allows us to compute the MAP estimate of the phase errors which constrains the estimate to overcome high noise and strong turbulence.

We compare the proposed algorithm to the image sharpening approach in [3] over a range of SNRs and atmospheric turbulence strengths.
Overall, our experimental results using simulated, singleframe, isoplanatic data demonstrate that the proposed MBIR algorithm can reconstruct substantially higherquality phase and image reconstructions compared to IS algorithms.
2 Estimation Framework
Figure 1 shows an example scenario for a DH system using the offaxis PPRG. Our goal is to compute the MAP estimates of the reflectance, , and the phase errors, , from the noisy data, . Note that we are using vectorized notation for these 2D variables. The joint estimates are given by
(1)  
where represents the jointlyfeasible set and the quantities and are assumed to be independent.
To evaluate the cost function for Eq. (1), we must derive the model for . The reader is referred to App. A for a more detailed derivation. In short, we can represent the data using an additive noise model given by
(2) 
where is the measurement noise, is the linear forward model operator, and is the complex field in the object plane. For an object with reflectance function , is defined as
(3) 
where is a diagonal matrix that applies the objectplane, quadratic phase factor from the Fresnel propagation integral [16]. Note that , and therefore for all . By representing the data as a function of , rather than , we avoid having to explicitly define in our model.
In Eq. (2), the matrix accounts for the propagation and measurement geometry. If we ignore the blurring effects caused by sampling the signal with finitesized pixels, can be decomposed as
(4) 
Here denotes an operator that produces a diagonal matrix from its vector argument, is the entrancepupil transmission function, and is the phase error function previously defined. For our purposes, is a binary circular function defining the entrance pupil of the imaging system, as shown in Fig 1. Finally, we choose the reconstruction parameters such that is a twodimensional discrete Fourier transform (DFT) matrix scaled so that .
For surfaces which are rough relative to the illumination wavelength, the vector can be modeled as a conditionally complex Gaussian random variable. Given the underlying reflectance, , of the scene, the conditional distribution is given by
(5) 
where for any random vector , is the multivariate complex normal distribution with mean, , and covariance matrix, , such that
(6) 
Here, superscript indicates the Hermitian transpose. The vector, , then has an identical distribution given by
(7)  
The equality in Eq. (7) results from the commutative property of diagonal matrices and from the fact that .
It is common for coherent detection systems to use a reference beam which is much stronger than the return signal. For such cases, shot noise, driven by the power of the reference beam, is the dominate source of measurement noise and can be modeled as additive, zeromean, complex Gaussian white noise [17, 18]. The distribution of is therefore given by
(8) 
where is the noise variance.
From Eqs. 28, the likelihood function of the data, given the reflectance and the phase errors, is distributed according to [19]
(9) 
Therefore, the MAP cost function is given by
(10)  
Unfortunately, the determinate and inverse make direct optimization of Eq. (10) an extremely computationally expensive task.
Equation (10) is similar to the MAP cost function in [15] except that the forward model has changed and we include the distribution of in the cost function. Despite these differences, we use a similar optimization approach which leverages the expectation maximization (EM) algorithm and allows us to replace Eq. (10) with a moretractable surrogate function.
In order to use the EM algorithm, we need to first introduce the concept of surrogate functions. For any function , we define a surrogate function, , to be an upperbounding function, such that , where is the current value of which 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 . That is
(11) 
Surrogate functions are useful because in some cases it is possible to find a simple surrogate function, , that we can iteratively minimize in place of a more difficult to compute function, . In fact, the EM algorithm provides a formal framework for constructing such a surrogate function. More specifically, for our problem the surrogate function of the EM algorithm takes the following form:
(12)  
where and are the current estimates of and , respectively. Here we choose to be the missing data. Evaluation of the expectation in Eq. (12) constitutes the Estep of the EM algorithm. Fig. 2 shows the alternating minimization approach used for implementing the Mstep of the EM algorithm. The proposed algorithm shares the same convergence and stability properties as the algorithm in [15]. More specifically, the EM algorithm converges in a stable manner to a local minima.
Since the cost function is nonconvex, the result of the EM iterations will depend on the initial conditions and may not be the global minimum. While it is not generally possible to compute a guaranteed global minima of the cost function in Section 4.C, we present a multistart algorithm that we empirically found to robustly compute good local minima.
3 Prior Models
In this section, we present the distributions used to model both the reflectance, , and the phase errors, . In both cases, we choose to use Markov random fields (MRFs) with the form
(13) 
where is the partition function, is the weight between neighboring pixel pairs ( and , or and ), is the difference between those pixel values, is the set of all pairwise cliques falling within the same neighborhood, and is the potential function [20]. We obtain different models for the reflectance and phase errors by choosing different potential functions.
3.1 Reflectance
For the reflectance function, , we used a QGeneralized Gaussian Markov Random Field (QGGMRF) potential function with the form
(14) 
where is a unitless threshold value which controls the transition of the potential function from having the exponent to having the exponent [21]. The variable, , controls the variation in . The parameters of the QGGMRF potential function affect its shape, and therefore the influence neighboring pixels have on one another.
While more sophisticated priors exist for coherent imaging [22], the use of a QGGMRF prior is relatively simple and computationally inexpensive, yet effective [15]. The low computational burden is important for DH applications where we desire nearrealtime, phaseerror estimation, such as with wavefront sensing.
3.2 Phase Errors
In practice, the phaseerror function may not change much from sample to sample, or we may wish to correct the phase at a defined realtime frame rate using a deformable mirror with fewer actuators than sensor measurements. Therefore, to reduce the number of unknowns, we allow the phaseerror function, , to be modeled on a grid that has lower resolution than the measured data. We denote the subsampled function as .
Using allows us to estimate only a single phase value for multiple data samples. This technique reduces the number of unknowns when estimating and from to , where is the factor of subsampling used in both dimensions. To scale to the resolution of for use in the forward model, we use a nearestneighborinterpolation scheme given by
(15) 
Here, is an interpolation matrix with elements in the set .
To model the distribution of unwrapped phase error values in , we use a Gaussian potential function given by
(16) 
where controls the variation in . The Gaussian potential function enforces a strong influence between neighboring samples which generates a moresmoothly varying output when compared to QGGMRF. Since the unwrapped phase function will not have edges or other sharp discontinuities, the GMRF is more appropriate than QGGMRF in this case. An added benefit is that it is less computationally expensive than the QGGMRF model which allows for fast evaluation.
4 Algorithm
In this section, we present the details for how we execute the proposed EM algorithm shown in Fig. 2. The discussion describes how we evaluate the surrogate function in the Estep and how we optimize it in the Mstep.
4.1 Estep  Surrogate Function Evaluation
Using the models given by Eqs. (79) and (1316), along with the conditional posterior distribution of given and , we can evaluate the surrogate function of Eq. (12). Its final form is given by
(17)  
Here, indicates that the matrix is dependent on , is the element of the posterior mean given by
(18) 
and is the diagonal element of the posterior covariance matrix given by
(19) 
We simplify the evaluation of Eq. (17) by approximating the posterior covariance as
(20) 
which requires that
(21) 
For a square entrancepupil function filling the detector’s field of view, and the approximations of Eqs. (20) and (21) become exact as both and are unitary matrices. When represents a circular aperture, then the relationship of Eq. (21) is approximate. However, since it dramatically simplifies computations, we will use this approximation and the resulting approximation of Eq. (20) in all our simulations when computing the function for the EM update. In practice, we have found this to work well.
4.2 Mstep  Optimization of the Surrogate Function
The goal of the Mstep is to minimize the surrogate function according to
(22) 
The alternating optimization shown in Fig. 2 is used to minimize with respect to and during each iteration of the EM algorithm.
We use Iterative Coordinate Descent (ICD) to update the reflectance [20]. In particular, ICD works by sequentially minimizing the entire cost function with respect to a single pixel, . By considering just the terms in Eq. (55) which depend on a pixel , the cost function is given by
(23) 
where the sum over indicates a sum over all neighbors of the pixel . We carry out the minimization of Eq. (23) with a 1D line search over .
To minimize Eq. (17) with respect to the phaseerrors, we consider just the terms which depend on . The phaseerror cost function becomes
(24) 
We use ICD to sequentially minimize the cost with respect to each element of the subsampled phaseerror function, . For element , corresponding to a group of data samples with indices in the set , the cost function becomes
(25) 
where
(26) 
and is an index over neighboring phase samples on the lower resolution grid. We minimize Eq. (25) using a 1D line search over , where minimizes the prior term. By minimizing Eq. (25), we obtain an estimate of the unwrapped phase errors. Finally, we obtain the fullresolution estimate of the unwrapped phase errors, , according to Eq. (15).
4.3 Initialization and Stopping Criteria
Figure 3 summarizes the steps of the EM algorithm. To determine when to stop the algorithm, we can use either a set number of iterations, , or a metric such as
(27) 
where is the iteration index and we stop the algorithm when falls below a threshold value of .
The EM algorithm in Fig. 3 is initialized according to
(28)  
where indicates the elementwise magnitude square of a vector, is a unitless parameter introduced to tune the amount of regularization in , and computes the sample variance of a vector’s elements [23].
We use a heuristic similar to that used in [15] to iteratively initialize the phaseerror estimate. Fig. 4 details the steps of this iterative process. The initial estimate of the phaseerror vector is simply . Next, for a set number of outerloop iterations, , we allow the EM algorithm to run for iterations. At the beginning of each outerloop iteration, we reinitialize according to Eq. (28).
After the outer loop runs times, we again reinitialize according to Eq. (28) and run the EM algorithm until it reaches the stopping threshold . A Gaussian prior model for works best in the outer loop for the initialization of . Specifically, we use , , , , and , where indicates a Gaussian kernel with standard deviation, . Once the initialization process is complete, we can use different priormodel parameters for the actual reconstruction.
5 Methods
To compare the proposed algorithm to an IS algorithm from [3], we generated simulated data using an approach similar to that of [3]. Figure 5 (a) shows the reflectance function that we multiplied by an incident power, , to get the object intensity, . Figure 5 (b) shows the magnitude squared of the corresponding reflection coefficient generated according to .
In accordance with the geometry shown in Fig. 1, we used a fast Fourier transform (FFT) to propagate the field, , from the object plane to the entrancepupil plane of the DH system. The field was padded to prior to propagation. Next, we applied a phase error, and the entrance pupil transmission function, , shown in Figure 5 (c). Figure. 5 (d) is an example of the intensity of the field passing through the entrance pupil which has a diameter, , equal to the grid length.
To generate the atmospheric phase errors, we used an FFTbased technique described in [24]. Using a Kolmogorov model for the refractiveindex powerspectral density (PSD), we generated random coefficients for the spatialfrequency components of the atmospheric phase, then we used an inverse FFT to transform to the spatial domain. To set the turbulence strength, we parameterized the Kolmogorov PSD by the coherence length, , also known as the Fried parameter [25, 24]. The value of specifies the degree to which the phase of a planewave passing through the turbulent medium is correlated. Two points in the phase function which are separated by a distance greater than will typically be uncorrelated. In this paper, we report turbulence strength using the ratio which is related to the degrees of freedom in the atmospheric phase. We simulated data for values of and . Figure 6 shows examples of the wrapped phase errors for each case.
After we added phase errors to the propagated field and applied the aperture function shown in Fig. 5, we mixed the signal with a modulated reference beam and detected the resultant power. Following [3], the referencebeam power was set at approximately 80% of the well depth per detector element, i.e., 80% of . We modeled Gaussian read noise with a standard deviation of and digitized the output to . After detection, we demodulated the signal to remove the spatialfrequency offset from the reference beam, lowpass filtered to isolate the signal of interest, and decimated to obtain a data array.^{1}^{1}1It is typical for this process to be carried out by taking an FFT, windowing a small region around the desired signal spectrum, and taking an inverse FFT. The resulting data was represented by Eq. (2) after vectorization.
We generated data over a range of SNRs which we define as
(29) 
where is the samplevariance operator used in Eq. (28). For opticallycoherent systems, SNR is well approximated by the average number of detected signal photons per pixel [17, 26]. At each turbulence strength, and at each SNR, we generated 18 i.i.d. realizations of the data. We then processed each i.i.d. realization independently and computed the average performance over the 18 independent cases.
To measure the distortion between the reconstructed images, , and the simulation input, , we used normalized root mean square error (NRMSE) given by
(30) 
where
(31) 
is the leastsquares fit for any multiplicative offset between and .
To measure distortion between the reconstructed phase error, , and the actual phase error, , we calculated the Strehl ratio according to
(32) 
where is the FFT operator and indicates that we take the maximum value of the argument. The function, , is a binary function that represents the aperture in the observation plane. It takes on the value inside the white dotted circle shown in Fig. 5 and outside.
Digital phase correction of DH data is analogous to using a pistononly, segmented deformable mirror (DM) in adaptive optics. Therefore, using the DM fitting error from [27] and Maréchal’s approximation, we can compute the theoretical limit of the Strehl ratio as
(33) 
where is the spacing between correction points in .
Using NRMSE and Strehl ratio, we compared performance of the proposed algorithm to the IS approach presented in [3] using a pointbypoint estimate of and the sharpness metric. The algorithm computes the phaseerror estimate according to
(34) 
where indicates the application of an exponent to each element. Following the process described in [3], we used 20 iterations of conjugate gradient to optimize Eq. (34) and the algorithm was initialized with a order Zernike polynomial found through an iterative process.
For the proposed MBIR algorithm, we allowed the outer initialization loop to run times, with EM iterations each time. We kept constant for all reconstructions. Once the iterative initialization process was complete, we set a stopping criteria of and let the EM algorithm run to completion. We used , , , , and as the QGGMRF prior parameters for image reconstruction. Additionally, we used , and for the phase error prior parameters. Using gives a total number of unknowns of . While varying with turbulence strength may produce more optimal results, we kept it fixed for simplicity. We found these parameters to work well over a wide range of conditions.
6 Experimental Results
Figure 7 shows example reconstructions for a subset of the results. Each block of images shows the reconstructions corresponding to the median Strehl ratio of the 18 i.i.d. data sets. Note that we only show five of the 20 SNR levels for each turbulence strength. The top row of each image block shows the original blurry images, the middle shows the IS reconstructions, and the bottom shows the MBIR reconstruction. To compress the large dynamic range that occurs in coherent imagery, we present the results using a logbased dB scale given by , where is the normalized reflectance function. The residual phase errors, wrapped to , are shown below each image block. To aid in viewing the higherorder residual phase errors in Fig. 7, we removed the constant and linear phase components, piston, tip, and tilt. The values were estimated according to
(35) 
where is an operator that wraps its argument to , is the residual phase error, is a constant, and are the linear components in the and dimensions, respectively.
The examples in Figure 7 show that the IS algorithm is able to correct most of the phase errors for and some of the errors at . For and , the IS images are blurred beyond recognition. In contrast, the proposed MBIR algorithm produced focused images for all but the lowest SNR reconstructions. In addition, we see how the MBIR reconstruction has significantly less speckle variation. The MBIR results also show that in many cases, the residualphase errors contain largescale patches separated by wrapping cuts. However, as Fig. 7 shows, the patches are approximately modulo equivalent and still produce focused images.
Figure 8 shows the resulting Strehl ratios and NRMSE for each algorithm as a function of SNR and . The curves in Fig. 8 show the average results for all 18 i.i.d. realizations along with error bars which span the standard deviation. We also plot the Strehl ratio limits given by Eq. (33) where we used for the IS algorithm and for the MBIR algorithm.
The results show that the IS algorithm works in a limited range of conditions for singleshot data. For , IS peaks at a Strehl ratio of when the SNR exceeds . For lower SNRs or stronger turbulence, the IS algorithm’s performance tapers off quickly. Conversely, the proposed MBIR algorithm is able to obtain Strehl ratios much higher than IS, even for low SNRs and strong turbulence. In the case, MBIR reaches the IS Strehl limit of 0.8 with about less SNR.
7 Conclusion
In this paper, we presented an inverse, modelbased approach to estimating phase errors and reconstructing images from digital holography (DH) data. We designed the algorithm for cases where only a single data realization is available and the phase errors are isoplanatic. Rather than estimating the spatially uncorrelated reflection coefficient, , we estimate the highlycorrelated reflectance function, . This allows us to better constrain the underdetermined system and produce specklereduced images and accurate phaseerror estimates with less data. Using first principals, we derived a discrete forward model for use in the MAP cost function. To obtain a moretractable surrogate function, we used the EM algorithm. Additionally, we introduced a GMRF prior model for the phase error function modeled on a lower resolution grid and presented optimization schemes for the joint estimates.
We compared the proposed algorithm to a leading image sharpness approach over a range of conditions. The results showed that the MBIR algorithm produced phase estimates with higher Strehl ratios and lower image distortion than the IS technique. For cases of high noise and large phase errors, IS was not effective, while MBIR was still able to produce accurate estimates. In conclusion, we showed that the proposed algorithm is an effective alternative to IS algorithms for estimating phase errors from single shot DH data and reconstructing images.
Appendix A: Forward Model Derivation
In this section we derive the linear forward model expressed by Eq. (2). Following the scenario depicted in Fig. 1, we derive a continuous model for the data, then a discrete representation.
Assume that we use a monochromatic plane wave to illuminate an object which has a complexvalued reflection coefficient, . We denote continuous functions using tildes. The return field passing through the entrance pupil can be represented using the Fresnel diffraction integral [16] given by
(36)  
where is a complex constant, represents the circular transmission function for the system’s entrance pupil, is the atmospheric phaseerror function assumed to be concentrated in a single layer in the entrancepupil plane (i.e. an isoplanatic phase error), is the wave number, and is the wavelength.
To simplify notation, we define the function
(37) 
and its corresponding continuous space Fourier transform (CSFT) given by
(38) 
We can also combine the pupilplane quadratic phase function of Eq. (36) with the atmospheric phase errors to get a single unknown phase term given by
(39) 
Given the definitions in Eqs. (37) through (39), the return signal of Eq. (36) can be written as
(40) 
Optical heterodyne detection is performed by mixing the returned signal with a local reference beam given by
(41) 
where is the amplitude and , are factors which control the spatialfrequency modulation. We combine the signal and reference onto the detector which measures the intensity given by
(42)  
where indicates the complex conjugate. We assume the system is shotnoise limited with noise , driven by the power of the reference beam [18]. Thus, we model as additive, zeromean, white Gaussian noise [17]. After detection we demodulate the detector output to remove the spatialfrequency offset of the reference, lowpass filter, and decimate the data. This isolates the last term of Eq. (42) which gives us our signal of interest ^{2}^{2}2In practice, the output of the detection process is only proportional to the righthand side of Eq. (43) by some unknown constant.,
(43) 
Equation (43) is continuous in both the and coordinate systems. However, we wish to represent the signal as discrete measurements, , generated from a discretespace signal, . We start by representing the discrete field in the object plane as
(44) 
where and are the spatialsampling periods in the object plane. Furthermore, if we sample the signal with a focalplane array and ignore the blurring effects of finitesized pixels, then the discrete measurements can be represented as
(45) 
where and are the spatialsampling periods in the measurement plane. Combining Eqs. (43)(45), we get
(46) 
where , , and are discrete versions of , , and , respectively, and
(47) 
Equation (46) can represented more compactly using matrixvector notation as
(48) 
where is the vectorized measurement noise and is the vectorized field in the object plane. The matrix can be decomposed as
(49) 
Here denotes an operator that produces a diagonal matrix from its vector argument, is the vectorized entrancepupil transmission function, and is the vectorized phase error function. The matrix is defined in Eq. (47).
Since the sum in Eq. (46) represents the forward propagation of the the field , we have scaled by so that . This ensures we conserve energy when propagating between the object and entrancepupil planes. Furthermore, we choose our reconstruction parameters such that and , where and are the grid sizes in the and dimensions, respectively. Thus, is exactly a Discrete Fourier Transform (DFT) kernel and can be efficiently implemented using a Fast Fourier Transform (FFT).
Appendix B: Derivation of the EM Surrogate Function
In this section, we derive the EM surrogate for the MAP cost function. We start by writing Eq. (12) as
(50)  
where we have used Bayes’ theorem inside the expectation and the fact that . Next, we substitute in the forward and prior models specified in sections 2 and 3. This gives
(51)  
where indicates the matrix is dependent on and the variable is a constant with respect to and .
To evaluate the expectations in Eq. (51), we must specify the conditional posterior distribution of . Using Bayes’ theorem,
(52)  
where is the partition function which absorbs any exponential terms that are constant with respect to . By completing the square, we can show that the posterior distribution is a complex Gaussian with mean
(53) 
and covariance
(54) 
Acknowledgment
We would like to thank Dr. Samuel T. Thurman for his helpful discussions and support of this work. In addition, we would also like to thank the Maui High Performance Computing Center (MHPCC) for their assistance and resources which allowed us to test the proposed algorithm over a wide range of conditions.
References
 [1] P. J. Winzer and W. R. Leeb, “Coherent LIDAR at low signal powers: basic considerations on optical heterodyning,” Journal of modern Optics 45, 1549–1555 (1998).
 [2] J. C. Marron and K. S. Schroeder, “Holographic laser radar,” Optics Letters 18, 385–387 (1993).
 [3] S. T. Thurman and J. R. Fienup, “Phaseerror correction in digital holography,” JOSA A 25, 983–994 (2008).
 [4] S. T. Thurman and J. R. Fienup, “Correction of anisoplanatic phase errors in digital holography,” JOSA A 25, 995–999 (2008).
 [5] 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,” Optics express 17, 11638–11651 (2009).
 [6] J. Marron, R. Kendrick, S. Thurman, N. Seldomridge, T. Grow, C. Embry, and A. Bratcher, “Extendedrange digital holographic imaging,” in “SPIE Defense, Security, and Sensing,” (International Society for Optics and Photonics, 2010), pp. 76841J–76841J.
 [7] A. E. Tippie and J. R. Fienup, “Phaseerror correction for multiple planes using a sharpness metric,” Optics letters 34, 701–703 (2009).
 [8] A. E. Tippie and J. R. Fienup, “Multipleplane anisoplanatic phase correction in a laboratory digital holography experiment,” Optics letters 35, 3291–3293 (2010).
 [9] A. E. Tippie, “Aberration correction in digital holography,” Ph.D. thesis, University of Rochester (2012).
 [10] J. R. Fienup, “Phase error correction in digital holographic imaging,” in “Digital Holography and ThreeDimensional Imaging,” (Optical Society of America, 2014), pp. DM1B–1.
 [11] R. A. Muller and A. Buffington, “Realtime correction of atmospherically degraded telescope images through image sharpening,” JOSA 64, 1200–1210 (1974).
 [12] M. F. Spencer, I. V. Dragulin, D. S. Cargill, and M. J. Steinbock, “Digital holography wavefront sensing in the presence of strong atmospheric turbulence and thermal blooming,” in “SPIE Optical Engineering+ Applications,” (International Society for Optics and Photonics, 2015), pp. 961705–961705.
 [13] M. T. Banet, M. F. Spencer, R. A. Raynor, and D. K. Marker, “Digital holography wavefront sensing in the pupilplane recording geometry for distributedvolume atmospheric aberrations,” in “SPIE Optical Engineering+ Applications,” (International Society for Optics and Photonics, 2016), pp. 998208–998208.
 [14] M. F. Spencer, R. A. Raynor, M. T. Banet, and D. K. Marker, “Deepturbulence wavefront sensing using digitalholographic detection in the offaxis image plane recording geometry,” Optical Engineering 56, 031213–031213 (2017).
 [15] C. Pellizzari, R. Trahan III, H. Zhou, S. Williams, S. Williams, B. Nemati, M. Shao, and C. A. Bouman, “Synthetic aperture LADAR: A model based approach,” Submitted to IEEE Transactions on Computational Imaging (TBD).
 [16] J. W. Goodman, Introduction to Fourier Optics (Roberts and Company, Englewood Colorado, 2005).
 [17] R. L. Lucke and L. J. Rickard, “Photonlimited synthetic aperture imaging for planet surface studies,” Applied Optics 41, 5,084–5,095 (2002).
 [18] V. V. Protopopov, Laser heterodyning, vol. 149 (Springer, 2009).
 [19] S. M. Kay, Fundamentals of Statistical Signal Processing, Detection Theory, vol. II (Prentice Hall Signal Processing Series, 1993).
 [20] C. A. Bouman, Model Based Image Processing (Unpublished, 2015).
 [21] J. B. Thibault, K. Sauer, C. Bouman, and J. Hsieh, “A threedimensional statistical approach to improved image quality for multislice helical ct.” Medical Physics 34, 4526–4544 (2007).
 [22] C. Pellizzari, R. Trahan III, H. Zhou, S. Williams, S. Williams, B. Nemati, M. Shao, and C. A. Bouman, “Optically coherent image formation and denoising using plug and play inversion framework,” JOSA (2017).
 [23] C. Forbes, M. Evans, N. Hastings, and B. Peacock, Statistical distributions (John Wiley & Sons, 2011).
 [24] J. D. Schmidt, “Numerical simulation of optical wave propagation with examples in matlab,” in “Numerical simulation of optical wave propagation with examples in MATLAB,” (SPIE Bellingham, WA, 2010).
 [25] L. C. Andrews and R. L. Phillips, Laser beam propagation through random media, vol. 52 (SPIE press Bellingham, WA, 2005).
 [26] H. A. Haus, Electromagnetic noise and quantum optical measurements (Springer Science & Business Media, 2012).
 [27] J. W. Hardy, Adaptive optics for astronomical telescopes (Oxford University Press on Demand, 1998).