Dynamic DH-MBIR for Phase-Error Estimation from Streaming Digital-Holography Data

The goal of a directed energy system is to focus a beam of light on a distant object. In terrestrial applications, this can be difficult because atmospheric turbulence introduces random phase errors (or distortions) that tend to disperse the otherwise focused beam. A solution to this problem is to pre-distort the phase of the outgoing beam. This can be done by first performing wavefront sensing (WFS), in which the phase distortion of incoming light is measured. The outgoing light can then be pre-distorted with the conjugate phase. However, this approach requires that the WFS be performed with very low latency, so that the phase corrections can be used before the atmospheric distortion changes.

Traditionally, WFS is performed with a Shack-Hartmann wavefront sensor. However, since the Shack-Hartmann method is based on hardware, it is limited in speed and is effective only with moderate turbulence strength [1].

Alternatively, digital holography (DH) sensors can be used to estimate the wavefront from direct measurements of the light. DH sensors measure the amplitude and phase of incoming light by interfering the received light with a reference beam. The resulting complex-valued measurements can then be used to computationally estimate phase distortions by solving an inverse problem [2].

The Image Sharpening (IS) algorithm is one approach to estimating phase distortions from DH data [3]. This method is optimization-based and models the phase distortion with a Zernike expansion. However, the measured reflection from coherent light contains speckle [4], which reduces the effectiveness of the IS algorithm. A solution to this speckle problem is to take multiple shots of DH data in which the speckle is decorrelated, but this increases the latency of the WFS, which is not acceptable in directed energy systems.

The DH model-based iterative reconstruction (DH-MBIR) approach in [5] uses a probabilistic model that accounts for speckle together with prior information or regularization. This approach represents a target image using the real-valued reflectance, which is smoother than the complex-valued reflection coefficients associated with speckle; hence reflectance is more amenable to regularization. This approach yields state-of-the-art estimates of images and pupil-plane phase errors using single-shot DH data. However, the current version of DH-MBIR is designed to process a single shot of DH data, and it requires hundreds of iterations per shot, hence is not practical in a real directed energy system.

In this paper, we introduce the Dynamic DH-MBIR algorithm (DDH-MBIR) that can reconstruct a streaming sequence of single-shot DH data with only 1 iteration per shot. This results in a low-latency, high-throughput algorithm for dynamic reconstruction of phase errors due to atmospheric turbulence.

The DDH-MBIR algorithm has two key advantages:

• •

  It uses temporal correlation in both $\phi$ and $r$ to dramatically reduce the number of iterations per time-step.

• •

  It avoids poor quality solutions associated with local minima by incorporating the bootstrap step into the dynamic update process.

Our simulation results demonstrate that under realistic assumptions, the Dynamic DH-MBIR algorithm can achieve highly accurate, low latency phase error estimation in a streaming system using single-shot DH data in each frame.

A DH sensor allows us to measure the complex electromagnetic field reflected from an object using spatial-heterodyne interferometry. As in [5], we model the complex-valued measurements at time $n$ as

where $y_{n}\in\mathbb{C}^{M}$ is the rasterized vector of complex DH measurements, $g_{n}\in\mathbb{C}^{M}$ is the vector of unknown complex reflection coefficients from the illuminated object, $w_{n}\in\mathbb{C}^{M}$ is a vector of complex measurement noise, and $A_{\phi_{n}}\in\mathbb{C}^{M\times M}$ is a linear transformation that is dependent on the unknown phase error, $\phi_{n}$.

If the wavelength of the light is small compared to the resolved pixel size on the object, then the reflection coefficient, $g_{n}$, is accurately modeled as a circularly symmetric, complex Gaussian with unknown variance, $r_{n}$. In this case, $r_{n}$ is real-valued reflectance with $r_{n}=E[|g_{n}|^{2}]$ and $g_{n}\sim CN(0,r_{n})$.

We assume an isoplanatic propagation model in which the phase distortion, $\phi_{n}$, occurs close to the detector [6] (however, our proposed method also applies in the anisoplanatic case). The isoplanatic linear forward model can be decomposed as

where $D_{a}$ is a diagonal matrix that models the camera aperture, ${\cal D}(e^{j{\phi_{n}}})$ is a diagonal matrix of phase distortions, $F$ is a normalized 2D discrete Fourier transform (DFT), and $\Gamma$ is a diagonal matrix of quadratic phase factors resulting from Fresnel propagation [7].

Our goal is to estimate the unknown phase errors, $\phi_{n}$, from the measurements, $y_{n}$. To do this, we will also need to estimate the object reflectance, $r_{n}$. Importantly, the reflectance is typically a smoother quantity with higher spatial correlation relative to the reflection coefficient, $g_{n}$, hence is more amenable to accurate estimation. Thus, our goal will be to compute at each time-step the joint-MAP estimate

Direct minimization of (2) is not tractable due to the nonlinear relationship between $r_{n}$ and $y_{n}$. However, the DH-MBIR algorithm solves this problem by using the expectation-maximization (EM) algorithm to solve for the MAP estimate using the iterative application of surrogate functions [5].

The DH-MBIR algorithm has the form

$r\leftarrow 0$; $\phi\leftarrow 0$ while not converged do     $(r,\phi)\leftarrow EM(r,\phi;y)$ end while

However, the existing DH-MBIR algorithm is impractically slow for two reasons. First, it requires hundreds of iterations of the $EM()$ function. Second, since the underlying function being minimized is non-convex, the algorithm tends to become trapped in local minimum. This second problem has been largely solved using a bootstrapping technique proposed by Pellizzari [9] in which the reconstructed image $r$ is reset to a backprojected estimate every 200 iterations.

In directed energy applications, $\phi_{n}$ must be estimated at a sufficiently high sampling rate, $f_{s}$, that the atmospheric phase distortion is nearly unchanged between samples. For this, we propose the Dynamic DH-MBIR (DDH-MBIR) algorithm for estimating $r_{n}$ and $\phi_{n}$ from a streaming sequence of DH data, $y_{n}$. Figures 1 and 2 show the pseudo-code and a flow diagram for DDH-MBIR. With each iteration, the data is read, normalized, and then residual tip/tilt removed using least-squares regression. Next, the value of $r_{n}$ is initialized, and $N_{k}$ iterations of the EM algorithm are performed.

The key innovation of the DDH-MBIR algorithm is the initialization of $r_{n}$ given by

This initialization leverages the temporal correlation of both $\phi$ and $r$ to dramatically reduce the number of EM iterations per time step, $N_{k}$ (we show later that $N_{k}=1$ yields high quality estimates of the phase distortion on simulation data under practical conditions). This initialization is formed by a weighted sum of two terms where $\lambda$ and $\alpha$ are user selectable weights in the range $[0,1]$. Since $r$ does not typically change substantially between samples, the first term, $r_{n-1}$, provides a good initial condition for the estimation of $r_{n}$.

The second term uses the magnitude squared of the back-projected data. This is typically how images are formed with DH data, and in particular, this is the term that is used in the bootstrapping algorithm of [9]. Intuitively, this second term incorporates a small amount of continuous bootstrapping into the dynamic update loop, which keeps the DDH-MBIR output away from solutions associated with local minima.

Finally, the RemoveTipTilt() function removes the linear component of the phase corresponding to spatial shifts of the reconstructed image; and the Normalize() function is defined as

where $\mu$ is the mean value of $y$ and $p$ is the number of entries in $y$. This normalization allows algorithmic parameters to be chosen more robustly.

In this section, we present estimation results from synthetic data using the Dynamic DH-MBIR algorithm.

We described the Dynamic DH-MBIR algorithm for quickly estimating phase-errors due to atmospheric turbulence. Using synthetic DH data generated from frozen flow phase-errors with no boiling, we showed that this algorithm can achieve a Strehl ratio greater than $0.8$ within $150$ timeframes using only $1$ iteration of the EM algorithm per timeframe. Since the computational cost of each iteration of the Dynamic DH-MBIR algorithm is dominated by two 2D FFTs [8], this makes the Dynamic DH-MBIR algorithm practical for real-time implementation.