## 1.

## Introduction

When imaging at long ranges through the atmosphere, the acquired images are highly susceptible to degradations from atmospheric turbulence.^{1}2.^{–}^{3} Fluctuations in the index of refraction along the optical path length, driven by temperature and pressure variations, give rise to spatially and temporally varying blur and warping. This is a well-researched area with respect to astronomical imaging.^{3} In astronomical imaging, with narrow fields of view, the degradation caused by the atmosphere can usually be modeled as isoplanatic. That is, the atmospheric effects are uniform across the image. This gives rise to warping that is a global image shift, and the blurring can be modeled with a spatially invariant point spread function (PSF). Wide field-of-view imaging, at long ranges through the atmosphere, generally leads to anisoplanatic imaging conditions. Here, the atmospheric PSF varies significantly across the field of view of the imaging sensor.

Adaptive optics has proven to be effective in treating the isoplanatic problem.^{4} However, for the anisoplanatic problem, mitigation methods are generally based on acquiring and digitally processing a sequence of short-exposure (SE) images.^{5} Short integration time means that warping during integration is minimized, reducing a major source of blurring. However, under anisoplanatic imaging conditions, there is still temporally and spatially varying warping and blurring to contend with. Note that with long-exposure (LE) images, turbulence-induced image warping gets temporally integrated. While this has the advantage of “averaging out” the geometric warping, to reveal the correct scene geometry, it also leads to high levels of blurring that may be difficult to treat effectively with image restoration.

One important class of turbulence mitigation algorithms is bispectral speckle imaging.^{6}7.8.9.10.11.12.13.14.15.^{–}^{16} This method seeks to recover the ideal image in the Fourier domain, by estimating the magnitude and phase spectrum separately. The magnitude spectrum is obtained with an inverse filter, or pseudoinverse filter, based on the LE optical transfer function (OTF). The phase is estimated using properties of the bispectrum.^{7}^{,}^{9}^{,}^{10}^{,}^{16}

Another class of turbulence mitigation algorithms uses some form of dewarping, fusion, and then blind deconvolution.^{8}^{,}^{15}16.17.18.19.20.^{–}^{21} Other related methods can also be found in the literature.^{22}23.24.25.26.27.^{–}^{28} With most of these methods, a motion-compensated temporal average of video frames is computed first. The motion compensation, prior to temporal averaging, reduces the motion blurring that might otherwise be seen in LE imaging. In the case of a static scene, the true geometry can be revealed with a prototype image obtained with a sufficiently long temporal average (or LE). Input frames can be registered to the prototype to provide turbulent motion compensation. As we shall show, even global frame registration can be of benefit. Performing such registration with a dynamic scene, containing moving objects, presents additional challenges.^{24}^{,}^{26}^{,}^{29} The current paper limits its scope to static scenes. Fusion is often done next by simple temporal averaging. This reduces noise and averages the spatially and temporally varying speckle PSFs in the individual frames. The result is an image that appears to be blurred with a spatially invariant PSF (with less blurring than an LE PSF). A blind image restoration process is then used to jointly estimate the spatially invariant PSF and true image. Note that using blind deconvolution has its challenges. First, it can be very computationally demanding. Also, unless a significant amount of *a priori* knowledge is incorporated, the recovered PSF and image may not be accurate.^{30}

Here, we present a block-matching and Wiener filtering (BMWF) approach to atmospheric turbulence mitigation for long-range imaging of extended scenes. We seek to leverage the rich theoretical work on atmospheric turbulence to aid in the design of a practical image restoration algorithm. We evaluate the proposed method, along with some benchmark methods, using simulated and real-image sequences. The simulated data are generated with a simulation tool developed by one of the current authors.^{31} These data provide objective truth and allow for a quantitative error analysis. The proposed turbulence mitigation method takes a sequence of SE frames of a static scene and outputs a single restored image. The images are globally registered to the temporal average and then reaveraged. This forms our prototype with the approximately correct geometry. A block-matching algorithm (BMA) is used to align the individual input frames to the prototype. We discuss how atmospheric statistics can help in setting the tuning parameters of the BMA. The BMA method here also uses a prefilter on the individual frames, so they better match the power spectrum of the prototype image for improved registration. The BMA registered frames are then averaged to generate a fused image. The final step is deconvolving the fused image using a Wiener filter.

An important aspect of the proposed method lies in how we model the degradation PSF. We use a parametric model that takes into account the level of geometric correction achieved during image registration. This is unlike any method we are aware of in the literature. By matching the PSF to the level of registration in this way, the Wiener filter is able to fully exploit the reduced blurring achieved by registration. We also describe a method for estimating the atmospheric coherence diameter (or Fried parameter) from the same estimated motion vectors used for restoration. We provide a detailed performance analysis that illustrates how the key tuning parameters impact the BMWF system performance. The proposed BMWF method is relatively simple computationally, yet it has excellent performance in comparison with state-of-the-art benchmark methods in our study.

The remainder of this paper is organized as follows. In Sec. 2, we present our observation model. This includes key statistics and the OTF models. The proposed BMWF turbulence mitigation approach is described in Sec. 3. The efficacy of the BMWF turbulence mitigation, in comparison with some benchmark methods, is demonstrated in Sec. 4. Finally, we offer conclusions in Sec. 5.

## 2.

## Optical Turbulence Modeling

## 2.1.

### Atmospheric Turbulence Statistics

One of the most important statistics that can be derived from the widely used Kolmogorov turbulence model is the atmospheric coherence diameter (or Fried parameter).^{3}^{,}^{32} This is given as

## (1)

$${r}_{0}={[0.423{\left(\frac{2\pi}{\lambda}\right)}^{2}{\int}_{z=0}^{z=L}{C}_{n}^{2}(z){\left(\frac{z}{L}\right)}^{5/3}\mathrm{d}z]}^{-3/5},$$Another very salient statistic is the tilt variance for a point source. This is the angle of arrival variance of a point source due to turbulence. An expression for the one-axis tilt variance, for the spherical wave case, is given as^{33}

## (2)

$${\sigma}_{\phi}^{2}=2.91{D}^{-1/3}{\int}_{z=0}^{L}{C}_{n}^{2}(z){\left(\frac{z}{L}\right)}^{5/3}\mathrm{d}z,$$## 2.2.

### Optical Transfer Functions

When imaging in atmospheric turbulence, the overall camera OTF can be modeled to include the atmospheric OTF and the diffraction OTF. This is given as

where $\rho =\sqrt{{u}^{2}+{v}^{2}}$ and $u$ and $v$ are the spatial frequencies in units of cycles per unit distance. The atmospheric OTF model typically used is given as^{3}

## (5)

$${H}_{\mathrm{atm},\alpha}(\rho )=\mathrm{exp}\{-3.44{\left(\frac{\lambda l\rho}{{r}_{0}}\right)}^{5/3}[1-\alpha {\left(\frac{\lambda l\rho}{D}\right)}^{1/3}\left]\right\}.$$^{34}

## (6)

$${H}_{\mathrm{dif}}(\rho )=\{\begin{array}{ll}\frac{2}{\pi}\left[{\mathrm{cos}}^{-1}\right(\frac{\rho}{2{\rho}_{c}})-\frac{\rho}{2{\rho}_{c}}\sqrt{1-{\left(\frac{\rho}{2{\rho}_{c}}\right)}^{2}}]& \rho \le {\rho}_{c}\\ 0& \text{otherwise}\end{array},$$## (7)

$${H}_{\mathrm{LE}}(\rho )={H}_{0}(\rho )={H}_{\mathrm{atm},0}(\rho ){H}_{\mathrm{dif}}(\rho ).$$## (8)

$${H}_{\mathrm{SE}}(\rho )={H}_{1}(\rho )={H}_{\mathrm{atm},1}(\rho ){H}_{\mathrm{dif}}(\rho ).$$^{3}An alternative SE OTF is given by Charnotskii.

^{35}

With the two main transfer functions defined, we now highlight a very interesting and important relationship between them that comes from the original development of Eq. (5). That is, it can be shown that

where ${G}_{\alpha}(\rho )$ is a Gaussian, given as## (10)

$${G}_{\alpha}(\rho )=\mathrm{exp}\{-(1-\alpha )3.44\frac{{(\lambda l)}^{2}}{{r}_{0}^{5/3}{D}^{1/3}}{\rho}^{2}\}.$$## (11)

$${H}_{\alpha}(\rho )={H}_{\mathrm{atm},1}(\rho ){G}_{\alpha}(\rho ){H}_{\mathrm{dif}}(\rho )={H}_{\mathrm{SE}}(\rho ){G}_{\alpha}(\rho ).$$Based on the Fourier transform properties of a Gaussian, the spatial-domain function resulting from the inverse Fourier transform ${g}_{\alpha}(r)={\mathrm{FT}}^{-1}[{G}_{\alpha}(\rho )]$ is also Gaussian. This is given as

## (13)

$${g}_{\alpha}(r)=\frac{2{\sigma}_{g}^{2}(\alpha )}{\pi}\text{\hspace{0.17em}}\mathrm{exp}\{-\frac{{r}^{2}}{2{\sigma}_{g}^{2}(\alpha )}\},$$## (14)

$${\sigma}_{g}^{2}(\alpha )=(1-\alpha ){\left(\frac{.4175\lambda l}{{r}_{0}^{5/6}{D}^{1/6}}\right)}^{2}.$$Thus, Eq. (12) shows that the parametric atmospheric PSF, ${h}_{\alpha}(r)$, is the SE PSF convolved by a Gaussian motion blur impulse response. When $\alpha =0$ (or equivalently $\beta =1$), the Gaussian motion blur standard deviation is the theoretical tilt standard deviation in Eq. (3). That is, ${\sigma}_{g}(\alpha )={\sigma}_{r}$, and we get the LE PSF

In the frequency domain, we have When the motion blur standard deviation is zero (i.e., $\alpha =1$ or equivalently $\beta =0$), Eq. (12) gives the SE PSF. For $0<\alpha <1$, Eq. (12) gives us the SE PSF convolved with a Gaussian motion blur somewhere between full tilt compensation and no tilt compensation. In the literature, it is typically the LE OTF ($\alpha =0$ or $\beta =1$) that is used for image restoration. However, if some level of registration is applied to the SE images (even if only global image registration), prior to fusion and deconvolution, we show that better results can be achieved by tuning $\alpha $ to the level tilt motion compensation. This gives us a powerful way to match the deconvolution step with the tilt correction processing step.Examples of atmospheric PSFs, ${h}_{\alpha}(r)$, from Eq. (12) are shown in Fig. 1. The optical system parameters corresponding to these plots are from the simulated data used in Sec. 4. The specific parameters are listed in Table 1. Figure 1(a) is for ${C}_{n}^{2}=0.25\times {10}^{-15}\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{m}}^{-2/3}$ (${r}_{0}=0.1097\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$) and Fig. 1(b) is for ${C}_{n}^{2}=1.00\times {10}^{-15}\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{m}}^{-2/3}$ (${r}_{0}=0.0478\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$). Note how the choice of $\alpha $ has a very significant impact on the width of the PSF for both levels of turbulence. As $\alpha $ is reduced, the Gaussian blurring component smoothes out and widens the PSF. We seek to match the $\alpha $ (or equivalently $\beta $), used in our PSF model, to the level of tilt correction provided by the registration stage of processing.

## Table 1

Optical parameters for simulated data.

Parameter | Value |
---|---|

Aperture | $D=0.2034\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$ |

Focal length | $l=1.2\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$ |

$f$-number | $f/\#=5.9\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$ |

Wavelength | $\lambda =0.525\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ |

Spatial cutoff frequency | ${\rho}_{c}=322.8410\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{cyc}/\mathrm{mm}$ |

Object distance | $L=7\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{km}$ |

Nyquist pixel spacing (focal plane) | ${\delta}_{N}=1.5488\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ |

## 2.3.

### Observation Model

Given the information in the preceding subsections, we are now ready to define the model we use to relate observed SE frames to a truth image. Note that the model is similar to that described earlier by Fraser et al.^{17} The observed frames are expressed in terms of a spatially varying blur operator and a spatially varying geometric warping operator. In particular, observed frame $k$ is given as

## (22)

$$\u27e8{f}_{k}(x,y)\u27e9={g}_{0}(x,y)*{h}_{\mathrm{SE}}(x,y)*z(x,y)={h}_{\mathrm{LE}}(x,y)*z(x,y).$$Now, consider the case where perfect tilt correction is applied to the SE frames. Let this tilt correction operator be expressed as ${\tilde{s}}_{k}^{-1}(x,y)[\xb7]$. Applying this to Eq. (19) and comparing this to Eq. (21), we get

## (23)

$$\u27e8{\tilde{s}}_{k}^{-1}(x,y)[{f}_{k}(x,y)]\u27e9=\u27e8{\tilde{h}}_{k}(x,y)[z(x,y)]\u27e9={h}_{SE}(x,y)*z(x,y).$$^{36}

^{,}

^{37}Thus, we define a partial tilt correction operator as ${\tilde{s}}_{k,\alpha}^{-1}(x,y)[\xb7]$. Applying this to the SE frames, and applying an ensemble mean, yields

## (24)

$$f(x,y)=\u27e8{\tilde{s}}_{k,\alpha}^{-1}(x,y)[{f}_{k}(x,y)]\u27e9={g}_{\alpha}(x,y)*{h}_{\mathrm{SE}}(x,y)*z(x,y)={h}_{\alpha}(x,y)*z(x,y).$$## 3.

## Turbulence Mitigation Approach

## 3.1.

### Block-Matching and Wiener Filtering Turbulence Mitigation

A block diagram representing the proposed BMWF turbulence mitigation algorithm is provided in Fig. 2. The input is a set of $N$ SE frames ${f}_{k}(x,y)$, for $k=1,2,\dots ,N$. We assume that these frames are sampled such that they are free from aliasing. Treating turbulence and aliasing simultaneously has been explored in the literature^{25}^{,}^{26}^{,}^{38}^{,}^{39} but it is not addressed here.

The input frames are buffered and averaged. Next, robust global translational registration is used to align the $N$ frames to the average. A least-squares gradient-based registration algorithm is used. This method is based on Lucas and Kanade^{40} but includes the robust multiscale processing described by Hardie et al.^{41}^{,}^{42} The frames are reaveraged after this global alignment to produce the prototype image with the desired geometry. This step also gives us the opportunity to compensate for any camera platform motion. For ground-based systems, translations may be sufficient. For airborne applications, affine registration at this stage may be appropriate.^{41}

Next, a BMA algorithm^{43} is used to estimate the local motion vectors for each pixel within each frame. The images are then interpolated, based on the motion vectors, to match the geometry of the prototype. Note that there is a mismatch between the level of blurring in the raw frames and the protoype being matched. One of the features of our method is that we prefilter the raw frames, using the Gaussian tilt blur, ${g}_{0.5}(x,y)$ from Eq. (13). Note that this is only done for the purposes of BMA, and we revert to the raw frames for subsequent processing.

As discussed in Sec. 2.3, the registration will not be ideal, and the accuracy of the registration is quantified by the parameter $\alpha $ (or equivalently $\beta $). The BMA registered frames are then expressed as ${\tilde{s}}_{k,\alpha}^{-1}(x,y)[{f}_{k}(x,y)]$, as shown in Fig. 2. Let us define the BMA block size as $B\times B$ pixels, and let the search window be $S\times S$ pixels in size (as defined by the position of the block centers). We use an exhaustive search within the search window, using the full rectangular blocks, and employ the mean absolute difference metric. We present results using a whole pixel search and subpixel search. The subpixel results are obtained by upsampling the images with bicubic interpolation. Because of its widespread use in image compression, much work has been done regarding performance analysis, speed enhancements, and hardware implementations of BMA.^{43} We leverage that work by incorporating the BMA here. The key parameters for BMA are $B$ and $S$. If knowledge of the atmospheric coherence diameter, ${r}_{0}$, is available, we can predict the amount of motion using Eq. (3). Exploiting this, we employ a search window that includes $\pm 2$ standard deviations, giving $S=2\lfloor 2{\sigma}_{r}/{\delta}_{N}\rfloor +1$ pixels, where ${\delta}_{N}$ is the pixel spacing measured on the focal plane.

With regard to block size, the larger these are, the less sensitive the BMA is to noise and warping. However, with increased size, there tends to be an increased underestimation of the true local motion from atmospheric tilt.^{36}^{,}^{37} Thus, a balance is required. The exact amount of underestimation will depend on the block size, the particular ${C}_{n}^{2}(z)$ profile, and optical parameters.^{36}^{,}^{37} Notwithstanding this, we have found that a fixed block size of $B=15\text{\hspace{0.17em}\hspace{0.17em}}\text{pixels}$ is effective for the range of turbulence conditions used in the simulated imagery. Furthermore, our results show that the corresponding residual RMS tilt factor is approximately a constant $\beta =0.1$ in the simulated imagery.

The next step of the BMWF method is to simply average the registered frames, as shown in Fig. 2. This gives rise to the result in Eq. (24). This step is important for two main reasons. First, it reduces noise and reduces the impact of any BMA errors. Second, by averaging the spatially varying blurring, it allows us to accurately model the resulting blur as spatially invariant,^{17} as shown in Eq. (24). This justifies the use of a spatially invariant deconvolution step. The deconvolution step is implemented here using a Wiener filter. The frequency response of the Wiener^{44} filter is given as

^{45}The images are padded symmetrically to minimize ringing artifacts associated with the circular convolution that results from FFT products.

Examples of the atmospheric OTF, ${H}_{\alpha}(\rho )$, from Eq. (11), are shown in Fig. 3. The optical system parameters corresponding to these plots are listed in Table 1. Figure 3(a) is for ${C}_{n}^{2}=0.25\times {10}^{-15}\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{m}}^{-2/3}$ (${r}_{0}=0.1097\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$), and Fig. 3(b) is for ${C}_{n}^{2}=1.00\times {10}^{-15}\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{m}}^{-2/3}$ (${r}_{0}=0.0478\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$). Also shown in Fig. 3 are the degradation OTFs multiplied by the Wiener filter transfer function in Eq. (25) for $\mathrm{\Gamma}=0.0002$ (the value used for the simulated and real data with 200 frames). Clearly, as $\alpha $ approaches 1 (equivalently $\beta $ approaches 0), the degradation OTF is more favorable to high spatial frequencies. The signal will be above the noise floor out to a higher spatial frequency. Consequently, the Wiener filter is able to provide gain out to a higher spatial frequency, without being overwhelmed with noise. This greatly extends the effective restored system OTF. When the degradation OTF value gets below the noise floor, governed by $\mathrm{\Gamma}$, the Wiener filter in Eq. (25) succumbs, as shown in Fig. 3. Note that with the illustrated NSR, the effective bandwidth of the sensor is nearly doubled, going from $\alpha =0$ (no registration or tilt correction) to $\alpha =1$ (full tilt correction). Matching the degradation model to the level of registration is essential to exploiting the full benefits of the registration.

Another important thing to note from Fig. 3 is that the degradation OTF has no zeros up to the optical cutoff frequency. Thus, the blur degradation is theoretically invertible (barring noise and quantization). Given this, and the fact that we often can expect a low NSR because many frames are averaged, the computationally simple Wiener filter tends to give excellent performance. More complex deblurring methods may be warranted when additional ill-posed blurring functions and/or very high levels of noise are present. However, we observed negligible performance gains using more complex regularization-based image restoration methods in our experiments.

## 3.2.

### Estimating the Atmospheric Coherence Diameter

To define the degradation OTF in Eq. (11) and the corresponding Wiener filter in Eq. (25), we require the parameter ${r}_{0}$. This can be measured using a scintillometer. However, in most practical imaging scenarios, it will not be known. Estimating this parameter from observed imagery is an active area of research.^{30}^{,}^{36}^{,}^{37} In some applications, it may be possible to set this parameter manually, based on subjective evaluation of the restoration results.

Here, we propose a method for estimating ${r}_{0}$ from the BMA motion vectors used in the BMWF algorithm. Based on Eq. (3), it is clear that ${r}_{0}$ is directly related to the warping motion in the observed imagery. The BMA motion vectors can give us an estimate of the warping motion. However, note that it is important that we exclude any camera platform motion or within scene motion when doing this. If the residual RMS tilt is $\beta {\sigma}_{r}$ from Eq. (16), the reduction in tilt due to registration is $(1-\beta ){\sigma}_{r}$. Let the BMA single-axis RMS motion, converted from pixels to distance on the focal plane, be denoted as ${\sigma}_{\mathrm{BMA}}$. Now we can estimate ${\sigma}_{r}$ as ${\widehat{\sigma}}_{r}={\sigma}_{\mathrm{BMA}}/(1-\beta )$. Using this and Eq. (3), we obtain an estimate of the atmospheric coherence diameter as

## 4.

## Experimental Results

In this section, we present a number of experimental results. These include results for both simulated and real SE image sequences. We also provide a comparison with state-of-the art benchmark methods. One benchmark is the bispectral speckle imaging method.^{7}^{,}^{9}^{,}^{16} Our implementation uses apodization with $16\times 16\text{\hspace{0.17em}\hspace{0.17em}}\text{pixel}$ tiles^{9} and incorporates local registration alignment with each tile, as this gave the best bispectrum method performance. Another benchmark is the method of Zhu and Milanfar.^{21} Our results for this method come from publicly available MATLAB^{®} code, provided courtesy of the authors. They use a B-spline nonrigid image registration (NRIR) of SE frames. This is followed by temporal regression to produce what the authors refer to as a near-diffraction-limited (NDL) image. Zhu and Milanfar^{21} suggest that blind deconvolution be applied to the NDL image. However, blind deconvolution code is not provided by those authors. Here, we have exact knowledge of the diffraction PSF, and, therefore, we apply a parameter-optimized Wiener filter to deconvolve diffraction-only blur from the NDL image. We also compute the temporal average of the NRIR frames as an alternative (bypassing the NDL regression operation), as an additional comparison.

## 4.1.

### Simulated Data

The simulated data are generated using an anisoplanatic simulation tool described by Hardie et al.^{31} The optical parameters used are listed in Table 1, and the simulation parameters are listed in Table 2. Five different levels of turbulence are simulated, and some statistical parameters for these scenarios are listed in Table 3. Additive-independent Gaussian noise, with a standard deviation of ${\sigma}_{\eta}=2$ digital units, is added to each simulated frame. The metric we use to evaluate the simulated data results is peak signal-to-noise-ratio (PSNR). Our first set of results is for $N=200$ temporally independent frames; then, we show results for $N=30$ frames.

## Table 2

Simulation parameters used in generating simulated frames.31

Parameter | Value |
---|---|

Path length | $L=7\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{km}$ |

Propagation step | $\mathrm{\Delta}z=700\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$ |

Cropped screen samples | $N=256$ |

Propagation screen width | $X=0.9699\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$ |

Propagation sample spacing | $\mathrm{\Delta}x=0.0038\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$ |

Number of phase screens | $N=10$ (9 nonzero) |

Phase screen type | Modified von Kármán with subharmonics |

Inner scale | ${l}_{0}=0.01\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$ |

Outer scale | ${L}_{0}=300\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$ |

Image size (object plane) | $2.3218\times 2.3218\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$ |

Image size (pixels) | $257\times 257\text{\hspace{0.17em}\hspace{0.17em}}\text{pixels}$ |

Pixel skip | 4 pixels ($65\times 65$ PSF array) |

Input image dynamic range | 256 digital units |

Gaussian noise variance | 4 digital units |

## Table 3

Theoretical statistical parameters for the different levels of simulated atmospheric turbulence and related restoration parameters.

Cn2×1015 (m−2/3) | |||||
---|---|---|---|---|---|

Parameter | 0.10 | 0.25 | 0.50 | 1.00 | 1.50 |

Isoplanatic angle ${\theta}_{0}$ (Pixels) | 6.62 | 3.82 | 2.52 | 1.66 | 1.30 |

BMA search size $S$ (Pixels) | 3.00 | 5.00 | 7.00 | 11.00 | 13.00 |

BMA block size $B$ (Pixels) | 15.00 | 15.00 | 15.00 | 15.00 | 15.00 |

Theoretical RMS G-tilt ${\sigma}_{r}$ (Pixels) | 0.8831 | 1.3962 | 1.9746 | 2.7925 | 3.4201 |

Estimated RMS G-tilt ${\widehat{\sigma}}_{r}$ (Pixels) | 0.8269 | 1.3781 | 1.9601 | 2.9022 | 3.5321 |

Theoretical ${r}_{0}(m)$ | 0.1902 | 0.1097 | 0.0724 | 0.0478 | 0.0375 |

Estimated ${\widehat{r}}_{0}(m)$ | 0.2058 | 0.1115 | 0.0730 | 0.0456 | 0.0360 |

## 4.1.1.

#### Results using 200 simulated frames

The PSNR results using 200 temporally independent frames are reported in Table 4. For the BMA algorithm, the search window size is set to $S=2\lfloor 2{\sigma}_{r}/{\delta}_{N}\rfloor +1\text{\hspace{0.17em}\hspace{0.17em}}\text{pixels}$. The block size is a constant $B=15\text{\hspace{0.17em}\hspace{0.17em}}\text{pixels}$. We report results for both whole pixel BMA and subpixel BMA in Table 4. We use the theoretical ${r}_{0}$ for each sequence for our OTF model. For the Wiener filter, we also report two sets of results. One where the optimum $\mathrm{\Gamma}$ and $\beta $ are searched for and used and another where fixed parameters are employed. The fixed NSR is $\mathrm{\Gamma}={\sigma}_{\eta}^{2}/(100N)$, where $N$ is the number of frames. The fixed residual tilt factors are: $\beta =1.0$ ($\alpha =0.00$) for the Wiener filter applied to the raw frame average (i.e., LE PSF), $\beta =0.5$ ($\alpha =0.75$) for the Wiener filter applied to the global registration average, and $\beta =0.1$ ($\alpha =0.99$) for the Wiener filter applied to BMA registered average. It is interesting to see in Table 4 how the PSNR increases by incorporating different levels of registration before averaging. As might be expected, the highest PSNR values are obtained with the subpixel BMA registration. It is also clear that there is a big boost in performance by adding the Wiener filter. The best results in Table 4 are generally from the subpixel BMA + average + Wiener filter.

## Table 4

PSNR (dB) results using 200 frames of simulated data with ση=2.0.

Cn2×1015 (m−2/3) | |||||
---|---|---|---|---|---|

Method | 0.10 | 0.25 | 0.50 | 1.00 | 1.50 |

Raw frames | 22.94 | 21.24 | 20.09 | 19.11 | 18.60 |

Average | 24.10 | 22.44 | 21.18 | 20.04 | 19.44 |

Global registration + average | 25.02 | 23.53 | 22.17 | 20.83 | 20.12 |

BMA + average | 25.46 | 24.36 | 23.06 | 21.57 | 20.74 |

BMA (sub) + average | 25.57 | 24.45 | 23.13 | 21.61 | 20.77 |

NDL image^{21} | 25.25 | 24.26 | 23.14 | 21.72 | 20.73 |

NRIR^{21} + average | 25.24 | 24.11 | 22.76 | 21.13 | 20.22 |

Average + Wiener (optimized) | 31.21 | 26.87 | 24.51 | 22.61 | 21.52 |

Average + Wiener ($\mathrm{\Gamma}=0.0002$, $\beta =1.0$) | 30.49 | 26.60 | 24.27 | 22.44 | 21.31 |

Global registration + average + Wiener (optimized) | 33.94 | 30.76 | 27.68 | 24.94 | 23.55 |

Global registration + average + Wiener ($\mathrm{\Gamma}=0.0002$, $\beta =0.5$) | 33.82 | 30.76 | 27.66 | 24.92 | 23.53 |

BMA + average + Wiener (optimized) | 34.25 | 32.96 | 31.33 | 28.41 | 25.98 |

BMA + average + Wiener ($\mathrm{\Gamma}=0.0002$, $\beta =0.1$) | 34.19 | 32.91 | 31.31 | 28.41 | 25.97 |

BMA (sub) + average + Wiener (optimized) | 34.12 | 33.22 | 31.69 | 28.74 | 26.16 |

BMA (sub) + average + Wiener ($\mathrm{\Gamma}=0.0002$, $\beta =0.1$) | 34.12 | 33.22 | 31.69 | 28.72 | 26.12 |

NDL^{21} + Wiener-diffraction (optimized) | 29.51 | 26.88 | 24.51 | 22.26 | 20.99 |

NRIR^{21} + average + Wiener-diffraction (optimized) | 30.15 | 27.12 | 24.30 | 21.72 | 20.54 |

Bispectral speckle imaging^{7}^{,}^{9}^{,}^{16} | 30.34 | 27.26 | 24.89 | 22.24 | 20.81 |

The bold entries represent the highest PSNR for that level of turbulence.

An analysis of the Wiener filter performance, as a function of $\mathrm{\Gamma}$ and $\beta $, is provided in Fig. 4 for $N=200$ and ${C}_{n}^{2}=1.00\times {10}^{-15}\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{m}}^{-2/3}$. Figure 4(a) shows the PSNR surface for the Wiener filter applied to the BMA registered frame average. Note that the optimum $\beta $ is near 0.1 and the optimum $\mathrm{\Gamma}$ is near 0.0002. This suggests that the BMA registration is about 90% effective in eliminating the RMS tilt (and $\sim 10\%$ remains). A similar surface plot is provided in Fig. 4(b) for the globally registered frame average (i.e., no BMA). Here, the optimum $\beta $ is near 0.5, suggesting that the global registration is about 50% effective in eliminating the RMS tilt. Finally, Fig. 4(c) is for no registration at all. Here, the optimum $\beta $ is approaching 1.0, as would be expected for an LE image. This analysis shows that the $\beta $ parameter should be matched to the level of tilt correction provided by the registration.

An analysis of the BMA parameters is shown in Fig. 5 for $N=200$ and ${C}_{n}^{2}=1.00\times {10}^{-15}\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{m}}^{-2/3}$. In particular, this plot shows the PSNR values as a function of $B$ and $S$. Here, one can see that the optimum block size is near $B=15$, and the optimum search window size does not increase after $S=11$. It is clear that small block sizes give much lower PSNRs. This is likely due to an insufficient amount of information for accurate matching, given the atmospheric degradations. Also, one can see that larger search windows generally do not hurt performance, but they do add to the computational load.

Figure 6 shows the system PSNR as a function of the number of input frames for ${C}_{n}^{2}=1.00\times {10}^{-15}\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{m}}^{-2/3}$. Performance increases dramatically for the first 30 frames or so and then becomes more incremental. However, additional frames continue to improve performance. Note that the curve is not monotonically increasing. The drops are likely due to the introduction of frames with large shifts, relative to the truth image.

Estimation of the atmospheric coherence diameter is illustrated in Fig. 7. The continuous curve shows the relationship from Eq. (3). The five simulation turbulence levels are shown with blue circles. The red squares show the estimated parameters from Eq. (27), using the BMA motion vectors with the parameters in Table 3, and $\beta =0.1$. This result appears to show a promising level of agreement between the estimates and true ${r}_{0}$ values.

Let us now turn our attention to image results. The truth image is shown in Fig. 8. Several output images, formed using $N=200$ and ${C}_{n}^{2}=0.25\times {10}^{-15}\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{m}}^{-2/3}$, are shown in Fig. 9. Figure 9(a) shows a single raw frame. Figure 9(b) shows the temporal frame average with no registration. The subpixel BMA registered frame average is shown in Fig. 9(c). Finally, the subpixel BMA + average + Wiener filter output is shown in Fig. 9(d). Here, the fixed-parameter Wiener filter is used. Note that the temporal average in Fig. 9(b) is rather blurry, as it is effectively equivalent to the true image corrupted with the LE PSF. The BMA registered average has corrected geometry and a blur level that is comparable to the observed SE frames. We see that by matching the PSF to the BMA registered average excellent results are possible, as shown in Fig. 9(d).

A similar set of results is shown in Fig. 10. These images are the same as Fig. 9, except we have increased the turbulence to ${C}_{n}^{2}=1.00\times {10}^{-15}\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{m}}^{-2/3}$. The raw SE frame is noticeably more corrupted than in Fig. 9. Also, the blurring in Figs. 10(b) and 10(c) is far more pronounced than in the corresponding images in Fig. 9. However, despite the increased turbulence, the subpixel BMA + average + Wiener output in Fig. 10(d) maintains much of the original detail. This is a consequence of having a very high signal-to-noise ratio, by virtue of the large number of input frames, and of an effective match between the PSF model and the blur in BMA registered average.

## 4.1.2.

#### Results using 30 simulated frames

The next set of results is for the restoration methods using $N=30$ input frames. The quantitative PSNR results are shown in Table 5. The results are similar to those in Table 4, but as expected, with fewer frames the PSNR values drop somewhat. The best results in Table 5 are for the subpixel BMA + average + Wiener filter, and these results are significantly better than those of the benchmark methods.

## Table 5

PSNR (dB) results using 30 frames of simulated data with ση=2.0.

Cn2×1015 (m−2/3) | |||||
---|---|---|---|---|---|

Method | 0.10 | 0.25 | 0.50 | 1.00 | 1.50 |

Raw frames | 22.87 | 21.17 | 20.04 | 19.07 | 18.56 |

Average | 24.03 | 22.37 | 21.14 | 20.01 | 19.41 |

Global registration + average | 24.99 | 23.50 | 22.14 | 20.80 | 20.09 |

BMA + average | 25.43 | 24.34 | 23.03 | 21.54 | 20.71 |

BMA (sub) + average | 25.55 | 24.43 | 23.09 | 21.58 | 20.73 |

NDL image^{21} | 25.23 | 24.24 | 23.13 | 21.60 | 20.62 |

NRIR^{21} + average | 25.22 | 24.07 | 22.72 | 21.10 | 20.18 |

Average + Wiener (optimized) | 29.07 | 25.52 | 23.59 | 22.03 | 21.14 |

Average + Wiener ($\mathrm{\Gamma}=0.0013$, $\beta =1.0$) | 28.91 | 25.36 | 23.45 | 21.99 | 21.12 |

Global registration + average + Wiener (optimized) | 31.84 | 28.61 | 25.98 | 23.76 | 22.63 |

Global registration + average + Wiener ($\mathrm{\Gamma}=0.0013$, $\beta =0.5$) | 31.75 | 28.61 | 25.98 | 23.76 | 22.62 |

BMA + average + Wiener (optimized) | 32.28 | 30.76 | 28.94 | 26.23 | 24.31 |

BMA + average + Wiener ($\mathrm{\Gamma}=0.0013$, $\beta =0.1$) | 32.28 | 30.75 | 28.94 | 26.22 | 24.29 |

BMA (sub) + average + Wiener (optimized) | 32.55 | 31.10 | 29.24 | 26.44 | 24.42 |

BMA (sub) + average + Wiener ($\mathrm{\Gamma}=0.0013$, $\beta =0.1$) | 32.54 | 31.10 | 29.24 | 26.41 | 24.39 |

NDL^{21} + Wiener-diffraction (optimized) | 29.25 | 26.69 | 24.44 | 22.09 | 20.86 |

NRIR^{21} + average + Wiener-diffraction (optimized) | 29.72 | 26.82 | 24.14 | 21.64 | 20.48 |

Bispectral speckle imaging^{7}^{,}^{9}^{,}^{16} | 29.50 | 26.63 | 24.29 | 21.87 | 20.55 |

The bold entries represent the highest PSNR for that level of turbulence.

To allow for a subjective comparison of the proposed method and the benchmark methods, output images from several methods are shown in Fig. 11 for the $N=30$ frame input, with ${C}_{n}^{2}=1.00\times {10}^{-15}\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{m}}^{-2/3}$. Figure 11(a) shows the temporal average, followed by the Wiener filter, using the LE PSF model. Figure 11(b) shows the NRIR + NDL image from Zhu and Milanfar,^{21} followed by the Wiener filter using the diffraction-only PSF model. Figure 11(c) shows the bispectral speckle image output.^{7}^{,}^{9}^{,}^{16} Finally, Fig. 11(d) shows the BMA + average + Wiener filter output, with subpixel BMA and fixed-parameter Wiener filter. The result in Fig. 11(a) is limited because no tilt correction is used, and the LE PSF is used in the Wiener filter. The NRIR + NDL + Wiener filter image in Fig. 11(b) provides improved results, but some areas remain highly blurred and there appear to be some artifacts at the edges. The bispectrum output in Fig. 11(c) also looks better than the LE restoration in (a) but is fundamentally limited by its use of the LE PSF model in recovering the magnitude frequency spectrum. The bispectrum method also tends to suffer from tiling artifacts when treating high levels of turbulence, as can be seen in Fig. 11(c). Processing without using tiles eliminates the tiling artifacts but leads to a lower quantitative performance (hence the use of tiles here). The subpixel BMA + average + Wiener filter output in Fig. 11(d) appears to have the best overall detail, with no major artifacts. This is supported by the quantitative analysis in Table 5.

The processing time for the various algorithms and their components is provided in Table 6. Note that the proposed method has a significantly shorter run time than the benchmark methods using our MATLAB^{®} implementations. However, run times with other implementations may differ. For the bispectral imaging method, processing time can be reduced by reducing the number of tiles and eliminating tile-based registration. Furthermore, hardware acceleration can be employed to speed up the multidimensional FFTs used with this method.

## Table 6

Algorithm run times for processing 30 simulated 257×257 pixel frames to produce a single output frame. Processing was done with MATLAB® 2016a using a PC with Intel(R) Xeon(R) CPU E5-2620 v3 at 2.40 GHz, 16 GB RAM, and running Windows 10. For the BMA method S=11, B=15, and 3× subpixel BMA.

Algorithm/component | Run time (s) |
---|---|

BMA + average + Wiener | 9.60 |

BMA (sub) + average + Wiener | 203.29 |

NRIR + NDL + Wiener^{21} | 3125.26 |

Bispectrum^{7}^{,}^{9}^{,}^{16} | 613.65 |

Global registration | 1.03 |

Whole pixel BMA | 8.56 |

Subpixel BMA | 202.25 |

NRIR^{21} | 3090.25 |

Average image fusion | 0.0012 |

NDL image fusion^{21} | 35.00 |

Wiener filter | 0.0125 |

## 4.2.

### Real Data

Our final set of experimental results uses a real-image sequence acquired from a tower to a truck and an engineering resolution target at a distance of 5 km. The resolution target is made up of a sequence of vertical and horizontal three-line groups. The five large groups on the right side have bars with the following widths: 7.00, 6.24, 5.56, 4.95, and 4.91 cm. The optical parameters for this sensor are listed in Table 7. The sensor sampling is very close to Nyquist, so the Wiener filter is evaluated and implemented at the pixel pitch of the sensor (i.e., no resampling of the imagery is performed). A scintillometer is used to provide an estimate of ${r}_{0}$, as shown in Table 7. This value has been confirmed by analysis of an edge target, imaged within the larger field of view of the camera. Assuming a constant ${C}_{n}^{2}(z)$ profile, note that the isoplanatic angle, when converted to pixels, is only 0.25 pixels. This gives rise to warping that is highly uncorrelated at a small scale. This makes BMA registration somewhat less effective than we saw in the simulated data. For this reason, we have chosen to use a residual RMS tilt factor of $\beta =0.4$ (compared with $\beta =0.1$ for the simulated data). An estimate of the atmospheric coherence diameter using Eq. (27), for $\beta =0.4$, is shown in Table 7.

## Table 7

Optical parameters for the real sensor data.

Parameter | Value |
---|---|

Aperture | $D=57.150\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}$ |

Focal length | $l=926\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}$ |

$f$-number | $f/\#=16.203$ |

Wavelength | $\lambda =0.785\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ |

Spatial cutoff frequency | ${\rho}_{c}=78.620\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{cyc}/\mathrm{mm}$ |

Object distance | $L=5\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{km}$ |

Nyquist pixel spacing (focal plane) | ${\delta}_{N}=6.36\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ |

Detector pitch (focal plane) | ${\delta}_{s}=6.50\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ |

Undersampling | 1.022 |

Scintillometer path averaged ${C}_{n}^{2}(z)$ | ${C}_{n}^{2}=7.44\times {10}^{-15}\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{m}}^{-2/3}$ |

Isoplanatic angle (assuming constant ${C}_{n}^{2}(z)$) | ${\theta}_{0}=0.26$ (Pixels) |

Scintillometer Fried parameter | ${r}_{0}=28.4\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}$ |

Estimated Fried parameter (for $\beta =0.4$) | ${\widehat{r}}_{0}=26.0\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}$ |

Theoretical RMS tilt for scintillometer ${r}_{0}$ | ${\sigma}_{r}=1.46$ (Pixels) |

BMA estimated RMS tilt (for $\beta =0.4$) | ${\widehat{\sigma}}_{r}=1.61$ (Pixels) |

The image results using the real data are shown in Fig. 12. Figure 12(a) shows raw frame 1. The NRIR + NDL^{21} + Wiener filter output using $N=30$ input frames is shown in Fig. 12(b). The bispectrum output^{7}^{,}^{9}^{,}^{16} is shown in Fig. 12(c), also for $N=30$ and using the scintillometer ${r}_{0}$. The 30 frame average + Wiener filter using the LE PSF with scintillometer ${r}_{0}$ is shown in Fig. 12(d). The subpixel BMA + average + Wiener output is shown for $N=30$ and $N=200$ in Figs. 12(e) and 12(f), respectively. Here, we use the ${\widehat{r}}_{0}$ estimated from the BMA, with $\beta =0.4$ and $\mathrm{\Gamma}=0.0002$. Note that the results obtained using the scintillometer ${r}_{0}$ are very similar. The BMWF results appear to provide the best overall subjective quality and recover the resolution target lines notably better than the benchmark methods. We attribute this to a reduction in tilt blurring, by means of the BMA registration, and by the proper matching of the PSF model to the residual tilt blurring.

## 5.

## Conclusions

We have presented a block-matching and Wiener filter-based approach to optical turbulence mitigation. In addition to the restoration method, we have also presented a method for estimating the atmospheric coherence diameter from the BMA motion vectors. We demonstrate the efficacy of this method quantitatively, using simulated data from a simulation tool developed by one of the authors. Results using real data are also provided for evaluation. The proposed restoration method utilizes a parametric OTF model for atmospheric turbulence and diffraction that incorporates the level of tilt correction provided by the registration step. By matching the PSF model to the level of registration, improved results are possible, as shown in Fig. 4. For the BMA component of our algorithm, we present a few innovations. For example, we use a search window size determined by the theoretical RMS tilt, when ${r}_{0}$ is available. The BMA also uses a prefilter on the raw frames, so they better match the prototype in spatial frequency content. Compared with benchmark methods, the proposed method provides the highest PSNR restorations in our study, as shown in Tables 4 and 5.

We quantify the level of registration tilt correction by what we term the residual RMS tilt factor, $\beta $, or equivalently, the tilt variance reduction factor, $\alpha =1-{\beta}^{2}$. Recall that $\beta $ is such that the residual RMS tilt, after registration, is $\beta {\sigma}_{r}$, where ${\sigma}_{r}$ is the theoretical uncorrected RMS tilt given in Eq. (3). Given $\alpha $ and ${r}_{0}$ and the optical system parameters, the degradation OTF model is given by Eq. (11) and the Wiener filter is given by Eq. (25). We have demonstrated that $\alpha $ can have a significant impact on the degradation OTF and corresponding restored image OTF, as shown in Figs. 1 and 3, respectively.

In cases where ${r}_{0}$ is known, it is possible to infer $\beta $ from Eq. (27), using the BMA motion vectors. If one does not have knowledge of ${r}_{0}$, it can be estimated from Eq. (27) using an assumed $\beta $ and the BMA motion vectors. Thus, a complete restoration can be achieved with the assumption of only one unknown parameter, $\beta $ (or $\alpha $). Note that this parameter is linked to the size of the BMA block size $B$, along with the camera parameters, and the ${C}_{n}^{2}(z)$ profile. We achieved excellent results using a constant $B=15\text{\hspace{0.17em}\hspace{0.17em}}\text{pixels}$, assuming a corresponding $\beta =0.1$ for the simulated data and $\beta =0.4$ for the real data studied here. In practice, it may be possible to perform a search over $\beta $ and evaluate the results subjectively or by some other metric.

## Acknowledgments

The authors would like to thank Dr. Doug Droege at L-3 Communications Cincinnati Electronics for providing support for this project. We would like to also thank Matthew D. Howard at AFRL for assisting with data collection. This work has been supported in part by funding from L-3 Communications Cincinnati Electronics and under AFRL Award Nos. FA8650-10-2-7028 and FA9550-14-1-0244. It has been approved for public release (PA Approval # 88ABW-2016-4934).

## References

## Biography

**Russell C. Hardie** is a full professor in the Department of Electrical and Computer Engineering, the University of Dayton, with a joint appointment in the Department of Electro-Optics. He received the University of Dayton’s top university-wide teaching award, the 2006 Alumni Award in teaching. He received the Rudolf Kingslake Medal and Prize from SPIE in 1998 for super-resolution research. In 1999, he received the School of Engineering Award of Excellence in teaching and the first annual Professor of the Year Award in 2002 from the student chapter of the IEEE.

**Michael A. Rucci** is a research engineer at the Air Force Research Laboratory, Wright-Patterson AFB, Ohio. His current research includes day/night passive imaging, turbulence modeling and simulation, and image processing. He received his MS and BS degrees in electrical engineering from the University of Dayton in 2014 and 2012, respectively.

**Alexander J. Dapore** is a senior image processing engineer at L-3 Cincinnati Electronics. He received his BSEE and MSEE from the University of Illinois at Urbana–Champaign in 2008 and 2010, respectively. He has worked on research and development projects in many areas of digital image processing. His specific areas of interest are image restoration, image enhancement, object/threat detection and tracking, multiview computer vision, and the real-time implementation of digital image processing algorithms on image processing algorithms using general-purpose computing on graphics processing units.

**Barry K. Karch** is a principal research electronics engineer at the Multispectral Sensing and Detection Division, Sensors Directorate of the Air Force Research Laboratory, Wright-Patterson AFB, Ohio. He received his BS degree in electrical engineering in 1987, his MS degree in electro-optics and electrical engineering in 1992 and 1994, respectively, and his PhD in electrical engineering in 2015 from the University of Dayton, Dayton, Ohio. He has worked in the areas of EO/IR remote sensor system and processing development for 29 years.