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

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³³

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

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

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

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

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}{\mathrm{m}}^{-2/3}$ ($r_0=0.1097\mathrm{m}$) and Fig. 1(b) is for $C_n^2=1.00\times {10}^{-15}{\mathrm{m}}^{-2/3}$ ($r_0=0.0478\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.

Fig. 1

Overall system PSFs with different $\alpha$ for the optical system parameters in Table 1 with (a) $C_n^2=0.25\times {10}^{-15}{\mathrm{m}}^{-2/3}$ ($r_0=0.1097\mathrm{m}$) and (b) $C_n^2=1.00\times {10}^{-15}{\mathrm{m}}^{-2/3}$ ($r_0=0.0478\mathrm{m}$). Note that as alpha increases, the PSF narrows.

Table 1

Optical parameters for simulated data.

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.¹⁷ 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

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)[·]$. Applying this to Eq. (19) and comparing this to Eq. (21), we get

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.

Fig. 2

Proposed BMWF turbulence mitigation system block diagram.

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⁴⁰ 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.⁴¹

Next, a BMA algorithm⁴³ 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.⁴³ 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{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,¹⁷ 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⁴⁴ filter is given as

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}{\mathrm{m}}^{-2/3}$ ($r_0=0.1097\mathrm{m}$), and Fig. 3(b) is for $C_n^2=1.00\times {10}^{-15}{\mathrm{m}}^{-2/3}$ ($r_0=0.0478\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.

Fig. 3

Overall system OTFs with different $\alpha$ for the optical system parameters in Table 1 with (a) $C_n^2=0.25\times {10}^{-15}{\mathrm{m}}^{-2/3}$ ($r_0=0.1097\mathrm{m}$) and (b) $C_n^2=1.00\times {10}^{-15}{\mathrm{m}}^{-2/3}$ ($r_0=0.0478\mathrm{m}$). The system OTFs multiplied by the Wiener filter frequency responses with $\mathrm{\Gamma }=0.0002$ are also shown to illustrate the effective OTF with restoration. Note that as alpha increases, the OTF widens.

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.1.

Simulated Data

The simulated data are generated using an anisoplanatic simulation tool described by Hardie et al.³¹ 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

Table 3

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

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{pixels}$. The block size is a constant $B=15\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 image21	25.25	24.26	23.14	21.72	20.73
NRIR21 + 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
NDL21 + Wiener-diffraction (optimized)	29.51	26.88	24.51	22.26	20.99
NRIR21 + average + Wiener-diffraction (optimized)	30.15	27.12	24.30	21.72	20.54
Bispectral speckle imaging7,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}{\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.

Fig. 4

Wiener filter parameter optimization using 200 simulated frames with $C_n^2=1.00\times {10}^{-15}{\mathrm{m}}^{-2/3}$. Wiener filter operating on: (a) the BMA registered frame average, (b) globally registered frame average, and (c) unregistered frame average. The highest PSNR is marked with a red “x.”

An analysis of the BMA parameters is shown in Fig. 5 for $N=200$ and $C_n^2=1.00\times {10}^{-15}{\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.

Fig. 5

BMA parameter optimization using 200 simulated frames with $C_n^2=1.00\times {10}^{-15}{\mathrm{m}}^{-2/3}$.

Figure 6 shows the system PSNR as a function of the number of input frames for $C_n^2=1.00\times {10}^{-15}{\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.

Fig. 6

Restoration performance versus the number of input frames for $C_n^2=1.00\times {10}^{-15}{\mathrm{m}}^{-2/3}$.

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.

Fig. 7

Simulation results for the atmospheric coherence diameter (Fried parameter) estimation from the BMA motion vectors using Eq. (27). Here, the $N=200$ simulated frames are used at the five $C_n^2$ levels. The same BMA parameters are used as shown in Table 3 and $\beta =0.1$.

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}{\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).

Fig. 8

Truth image used for simulated image sequence.

Fig. 9

Restoration results using $N=200$ frames with $C_n^2=0.25\times {10}^{-15}{\mathrm{m}}^{-2/3}$. (a) Raw frame 1, (b) raw frame average (no registration), (c) BMA registered frame average, and (d) BMA (sub) + average + Wiener filter output. Video 1 (MOV, 824 KB [URL: http://dx.doi.org/10.1117/1.OE.56.7.071503.1]) shows the raw frames on the left and BMA on the right.

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}{\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.

Fig. 10

Restoration results using $N=200$ frames with $C_n^2=1.00\times {10}^{-15}{\mathrm{m}}^{-2/3}$. (a) Raw frame 1, (b) raw frame average (no registration), (c) BMA registered frame average, and (b) BMA (sub) + average + Wiener filter output. Video 2 (MOV, 816 KB [URL: http://dx.doi.org/10.1117/1.OE.56.7.071503.2]) shows the raw frames on the left and BMA on the right.

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 image21	25.23	24.24	23.13	21.60	20.62
NRIR21 + 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
NDL21 + Wiener-diffraction (optimized)	29.25	26.69	24.44	22.09	20.86
NRIR21 + average + Wiener-diffraction (optimized)	29.72	26.82	24.14	21.64	20.48
Bispectral speckle imaging7,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}{\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,²¹ 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.

Fig. 11

Restoration results using 30 frames with $C_n^2=1.00\times {10}^{-15}{\mathrm{m}}^{-2/3}$. (a) Raw frame average (no registration) + Wiener filter, (b) NRIR + NDL²¹ + Wiener (diffraction), (c) bispectral speckle imaging method,^7,9,16 and (d) BMA (sub) + average + Wiener filter output.

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 + Wiener21	3125.26
Bispectrum7,9,16	613.65
Global registration	1.03
Whole pixel BMA	8.56
Subpixel BMA	202.25
NRIR21	3090.25
Average image fusion	0.0012
NDL image fusion21	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.

The image results using the real data are shown in Fig. 12. Figure 12(a) shows raw frame 1. The NRIR + NDL²¹ + 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.

Fig. 12

Restoration results using real-image sequence. (a) First raw frame, (b) 30 frame NRIR + NDL²¹ + Wiener (diffraction) output, (c) 30 frame bispectral speckle imaging output,^7,9,16 (d) 30 frame average + Wiener filter output, (e) 30 frame BMA (sub) + average + Wiener filter output, and (f) 200 frame BMA (sub) + average + Wiener filter output. Video 3 (MOV, 507 KB [URL: http://dx.doi.org/10.1117/1.OE.56.7.071503.3]) shows the raw frames on the top and BMA on the bottom.