2.1.

Algorithm Description

Figure 1 shows a block diagram summarizing the FIF SR method, originally from Karch and Hardie,¹⁷ and adapted here for joint SR and TM. In our implementation, short-exposure LR observed frames are registered using one of a variety of registration methods appropriate to the application. Next, single-frame interpolation is used to upsample the individual input images by a factor of $M$ to the Nyquist sampling grid, based on the diffraction-limited optical cut-off frequency. The interpolated images are formed in alignment to a common reference frame or frame average.^9,16 We shall consider different interpolation kernels here in our analysis. The next block is where the multiframe fusion takes place. Each interpolated pixel in each frame gets a weight based on subpixel alignment for each interpolated pixel. Finally, a Wiener filter is applied to the fused image to produce a single restored image. The Wiener filter makes use of a PSF model that incorporates atmospheric parameters, optical system parameters, level of tilt reduction in the registration step, and the interpolation blurring.

Fig. 1

Block diagram of the FIF SR method. Observed frames are registered, interpolated individually, and then fused based on a subpixel weighting. A Wiener filter provides restoration based on an OTF model that incorporates registration accuracy and subpixel weighting.

Note that the FIF SR algorithm may be viewed as a type of nonuniform interpolation SR method.² The FIF SR method and other nonuniform interpolation SR methods simplify the SR problem by separating it into registration and nonuniform interpolation followed by OTF restoration. Applying the Wiener filter for OTF restoration after the nonuniform interpolation is justified by the assumption that the warping and burring operators commute in the degradation observation model. In the case of translational motion, the assumption is fully valid.² For some other types of motion, such as affine, the warping and burring operators have been shown by Hardie et al.²¹ to approximately commute.

Let us now consider the heart of the FIF SR method, which is the fusion step. Assume there are $K$ input frames and one fused output frame. We define the interpolated input pixel $i$ in the frame $k$ as $f_k(i)$, and this corresponds to a sample of $f_k(x,y)$ in Fig. 1. Let the fused pixel $i$ on the Nyquist grid be denoted as $g(i)$. This represents a sample of $g(x,y)$ in Fig. 1. With this notation, the fused image pixels are given by

Fig. 2

Nyquist interpolation grid (red squares) and the noninterpolated pixel positions from a single LR frame (blue circles). A larger distance implies a larger interpolation error for that frame at that pixel, and consequently, a lower weight using Eq. (2).

Note that Eq. (2) gives a Gaussian weighting, with more weight given to frames that have a lower distance to an observed pixel. This is because larger distances tend to give larger interpolation errors.²² The weighting function is plotted in Figs. 3(a) and 3(b) for $\beta =0.25$ and 0.10, respectively. One can see that increasing $\beta$ makes the fusion less selective, giving increased weight to frames with larger interpolation distances. In the limiting case of $\beta =\infty$, equal weight is given to all frames and the fusion becomes a simple temporal average. Using a larger $\beta$ provides more temporal averaging in the fusion process that can help to reveal the proper scene geometry in moderate to heavy turbulence, and make the atmospheric blurring more spatially invariant in the fused image.^9,16 This also helps to attenuate temporal noise. On the other hand, a small $\beta$ gives a high level of selectivity, where significant weight is only given to interpolated pixels that are close to observed samples in that frame. A small $\beta$ would be expected to give the best SR results with minimal interpolation error, provided a sufficient diversity of samples is available with a high signal-to-noise ratio.¹⁷

Fig. 3

Interpolation weighting from Eq. (2) for (a) $\beta =0.25$ and (b) $\beta =0.10$. Binning weighting from Eq. (3) is shown in (c) for $M=4$.

The tunability provided by the parameter $\beta$ sets it apart from most other nonuniform interpolation SR methods² and makes it well suited to perform SR in the presence of turbulence. For example, it is interesting to compare the FIF SR approach of fusing interpolated frames with methods that use binning. Binning methods register the input LR frames and populate an HR grid by putting the LR pixels into discrete bins on the HR grid.^2,21 LR pixels are assigned to the nearest HR bin in a simple quantization process or by interpolation to the nearest HR bin. With such binning methods, care must be taken to address the very common scenario of empty bins.²¹ By fusing interpolated frames and using Gaussian weighting, the FIF SR method does not have this issue and will never have any empty output pixels. Using a different weighting function, however, it is possible for the FIF SR framework to perform fusion that is equivalent to binning. Specifically, this binning operation is achieved with the weighting function,

Further insight may be gained with regard to the parameter $\beta$ by considering what we term the “averaging power” of the fusion. By this, we mean the variance reduction factor for independent and identically distributed (i.i.d.) temporal samples, relative to a standard average. Consider the fusion of $K$ i.i.d. temporal samples with variance ${\sigma }^2$. The output variance would depend on the fusion weights and is given by

Fig. 4

Averaging power factor from Eq. (6) as a function of $\beta$.

It should be noted that the original FIF SR paper¹⁷ considers other weighting components, in addition to that in Eq. (2). One additional weighting term is designed to exploit color correlation. Another term is included to address in-scene motion to minimize motion blur and distortion.¹⁷ However, since our focus here is on the simultaneous SR and TM, we limit the scenarios under consideration to single-band imagery (i.e., no color), and no in-scene motion. Treating color and in-scene motion are certainly important problems that we hope to address in future work, in the context of joint TM and SR.

2.2.

Interpolation Impulse Response

In this subsection, we derive an OTF model for the fused interpolation blur. First, note that each sample in $g(x,y)$ in Fig. 1 may be expressed as a weighted sum of observed pixels. The specific weighting will depend on the exact motion and interpolation kernel used as well as the weighting from Eq. (2). For anything other than rigid translational motion, the weights will vary spatially.^21,23 However, if we assume a large number of frames with motion significantly greater than one LR pixel, we believe it is reasonable to assume a uniform subpixel distribution of motion. This allows us to model the interpolation blur with a spatially invariant impulse response. This is very similar to how the turbulence blur can be treated as spatially invariant after fusing a large number of frames.^9,16 Computing a spatially invariant interpolation impulse response for the fused image $g(x,y)$ is powerful in that it can be incorporated into the overall degradation model, as shown in Fig. 1. This allows us to mitigate some of the nonideal aspects of the interpolation step. It also adds to our overall understanding of the FIF SR method and further informs the selection of $\beta$.

The choice of interpolation kernel impacts the resulting interpolation blur model. Consider three common 1-D interpolation kernels²⁴ that are illustrated in Fig. 5. The zero-order hold (ZOH), or nearest neighbor, interpolation kernel is given by

Fig. 5

Continuous interpolation functions as a function of LR pixel spacings for ZOH, linear, and cubic interpolation.

Next, let the LR sampling function (i.e., integers corresponding to LR sample positions) be expressed as

Using a 1-D version of the weighting function from Eq. (2), we have $w(s)=e^{-s^2/{\beta }^2}$. Putting this into Eq. (14) gives

Examples of the frame-averaged interpolation impulse response from Eq. (15) are plotted in Fig. 6 for ZOH, linear, and cubic interpolation kernels. One can see that with $\beta =\infty$, the impulse response is simply the interpolation kernel, as shown in Eq. (16). As $\beta$ gets smaller, the width of the impulse response also gets smaller. The interpolation OTFs are shown in Fig. 7 for the same kernels and values of $\beta$ as those in Fig. 6. Note that in these plots, the folding frequency for SR with upsampling by a factor of $M$ is $0.5M$ cycles/(LR spacing). From Fig. 7, it is clear that the interpolation blur for all but very small $\beta$ can be quite significant. The good news is that much of the frequency content, while attenuated, is available for restoration, provided that the interpolation OTF is included in overall degradation model. To allow for SR restoration of the higher frequencies, it is important to use the smallest $\beta$ possible, while meeting the other algorithm requirements.

Fig. 6

Frame-averaged interpolation impulse response from Eq. (15) with multiple $\beta$ for (a) ZOH interpolation, (b) linear interpolation, and (c) cubic interpolation.

Fig. 7

Frame-averaged interpolation frequency response cross-section from Eq. (18) with multiple $\beta$ for (a) ZOH interpolation, (b) linear interpolation, and (c) cubic interpolation.

3.1.

Overall OTF Model

In this subsection, we present the overall OTF model used for the FIF SR restoration. The model is spatially invariant and is meant to capture the degradation as seen in the fusion image $g(x,y)$ in Fig. 1. The spatial invariance is justified because $g(x,y)$ represents a weighted sum of many short-exposure images.^9,16 The overall model and its components are expressed as

In Fried’s derivation, $\alpha =0$ in Eq. (21) gives the long exposure OTF. The average tilt-corrected short-exposure OTF under near field conditions is given by Eq. (21) when $\alpha =1$. The near field condition is defined by Fried to be $D\gg \sqrt{L\lambda }$, where $L$ is the optical path length. As shown by Tofsted’s analysis,^27,28 most sensors operate in regimes, where Fried’s near field condition should be invoked for the short exposure case (even for a far-field optical condition). Since we are applying this OTF to partially tilt-corrected imagery resulting from the potentially imperfect registration step in Fig. 1, we follow the approach in Hardie et al.⁹ and treat the parameter $\alpha$ as a continuous tilt reduction factor. At the extremes, a value of $\alpha =0$ is used for no registration (i.e., the long exposure OTF), and a value of $\alpha =1$ would be used for ideal registration (i.e., the short-exposure OTF). This parameter may be estimated based on the type of registration used.⁹

Referring to the block diagram in Fig. 1, the ideal image, $z(x,y)$, relates to the fusion image as

An example OTF from Eq. (19), showing the various components, is provided in Fig. 8 for $\alpha =0.5$ and $\beta =0.3$ using cubic interpolation. The optical parameters for this figure are listed in Table 1. The parameters in Table 1 are also used for the simulated data presented in Sec. 4.1. Also shown in Fig. 8 is the native sensor folding frequency for this particular example. Note that any signal energy above the folding frequency will be aliased during sampling. It is interesting to note from Fig. 8 that the frame-averaged interpolation OTF, $|H_{\beta }(u,0)|$ (shown in red), is comparable to the isolated atmospheric OTF from Eq. (21) (shown in green). Thus, it is clear that with $\beta$ values on this order, mitigating the impact of the interpolation is important for restoration, especially for SR, where frequency content out to the diffraction-limited cutoff frequency is sought.

Fig. 8

Example OTF from Eq. (19) showing the various components for $\alpha =0.5$ and $\beta =0.3$. The optical parameters for this figure are listed in Table 1.

Table 1

Optical parameters used for the OTF plot in Fig. 8 and for the simulation results in Sec. 4.1.

As can be seen in Eqs. (19)–(21), the parametric form of the OTF used for FIF SR restoration has several parameters. In practice, it is generally reasonable to assume that the optical parameters, $\lambda$, $D$, and $l$, would be known a priori. The remaining parameters are $\alpha$, $\beta$, and $r_0$. We have observed generally good results with $\beta =0.25$ for a wide range of datasets. The tilt reduction parameter, $\alpha$, typically ranges from 0 for no registration to about 0.5 for optical flow registration. Imperfect values of $\alpha$ can usually be reasonably well compensated for by altering the employed value of $r_0$. Thus, the bottom line is that OTF estimation and resulting restoration performance are mainly sensitive to $r_0$. Fixing the other parameters and searching over the $r_0$ space can often lead to very useful results. Subjective evaluation or other no-reference metrics³⁰ can be used for selection. Other scene-based methods for estimating $r_0$ can be found in the literature.^31,32

3.2.

Aliasing and Turbulence

As noted earlier, the Nyquist criterion dictates that the sampling frequency be greater than two times the highest spatial frequency in the image to guarantee no aliasing. The sampling frequency is given by ${\rho }_s=1/p$, where $p$ is the detector pitch, and the highest spatial frequency is limited by the diffraction-limited optical cutoff frequency ${\rho }_c$ from Eq. (20). Thus, the diffraction-limited sampling status of an imaging system can be characterized by the diffraction-limited sampling factor parameter,³³

While the metric in Eq. (25) is very useful, it does not take atmospheric turbulence into account. As can be seen from Eq. (21), the atmospheric optical turbulence acts like a nonideal low pass filter with no absolute cutoff frequency. Nevertheless, the OTF signal energy above the folding frequency diminishes in Eq. (21) as $r_0$ goes down, reducing the potential aliasing for a given sampling frequency. In order to capture this effect, we propose a modified $Q$ parameter, where we substitute $r_0$ for the aperture, $D$. This gives rise to what we term the turbulence-limited sampling factor,

Several overall OTFs are shown in Fig. 9 to illustrate the relationship between the focal length and $r_0$ for a fixed diffraction-limited sampling factor parameter of $Q=2/3$. For these plots, the optical parameters in Table 1 are used along with $\alpha =0.5$ and $\beta =0.3$. Focal lengths of $l=0.20$, $l=0.50$, and $l=1.00$ are shown in Figs. 9(a)–9(c), respectively. Note that with the short focal length (wide field of view), the overall OTF is less sensitive to the atmospheric effects represented by $r_0$. Each of the OTF curves in Fig. 9(a) has a significant amount of the OTF above the folding frequency. Furthermore, we see in this plot that the level of aliasing is more often limited by diffraction than turbulence (i.e., $Q>\tilde{Q}$). However, with increased focal length (reduced field of view), the effects of turbulence are magnified, as can be seen in Fig. 9(c). For the higher levels of turbulence in Fig. 9(c), the atmosphere acts as an anti-aliasing low-pass filter, essentially eliminating the possibility of aliasing. Most of the curves in this figure show that aliasing is limited more by turbulence than diffraction (i.e., $\tilde{Q}>Q$).

Fig. 9

Overall OTFs for a system with optical parameters listed in Table 1, $\alpha =0.5$, $\beta =0.3$, and (a) a focal length of $l=0.20\mathrm{m}$, (b) $l=0.50\mathrm{m}$, and (c) $l=1.00\mathrm{m}$.

The plots in Fig. 9 also illustrate that many scenarios exist where both turbulence and aliasing may be present. With short focal lengths, light to moderate turbulence is less problematic and traditional SR may be appropriate. For long focal length systems, all but light turbulence tends to effectively eliminate aliasing, making traditional TM appropriate. However, both aliasing and turbulence degradations are significant factors in heavy turbulence with short focal lengths, light turbulence with long focal lengths, and moderate focal lengths with a wide range of turbulence. This demonstrates the importance and relevance of developing algorithms capable of performing SR in the presence of turbulence, such as the one presented here.

4.1.

Simulated Data

The simulated data have been generated using the anisoplanatic optical turbulence simulation tool recently developed by Hardie et al.²⁰ The simulation tool uses numerical wave propagation and has performed well reproducing key image statistics in validation studies.²⁰ The specific simulation presented here is novel in that we are emulating an undersampled imaging system. The turbulence degraded images are first simulated at the Nyquist rate for the diffraction-limited optical system.²⁰ Next, however, we simulate detector integration and follow this with downsampling by a factor of $M=3$. After downsampling, we introduce additive Gaussian noise with a standard deviation of two digital units to imagery having a 0-255 original dynamic range. The optical system parameters used in the simulation are listed in Table 1, and the simulation parameters are listed in Table 2. We have simulated seven levels of turbulence, each with $K=200$ temporally independent frames.

Table 2

Simulation parameters used for the simulated imagery in Sec. 4.1.

Quantitative results for two different truth images are provided in Tables 3 and 4 for several algorithm variations. Table 3 is for the Kodak lighthouse image,¹⁷ and Table 4 is for the stream and bridge image.⁹ The metric we use to evaluate the simulated data results is peak signal-to-noise-ratio (PSNR). The PSNRs for the lighthouse image are also plotted in Fig. 10. For the results with parameter optimization, the parameters are found by a search to maximize the PSNR and these parameters are listed in Table 5. Note that results are shown with and without camera jitter. When camera jitter is on, we simulate camera platform motion by providing additional uniform random subpixel shifts between frames with no motion blur. The FIF SR results include the full OTF model, except where no $\beta$ restoration is indicated. The true lighthouse image and three levels of degradation are shown in Fig. 11. Various restored images are shown in Fig. 12 using the optimum parameters in Table 5. Global affine registration to the average frame is used for all of the FIF SR results reported on the simulated data.

Table 3

Lighthouse image PSNR (dB) results using simulated K=200 frames with ση=2.0 and M=3.

Table 4

Stream bridge image PSNR (dB) results using 200 simulated frames with ση=2.0 and M=3.

Fig. 10

Plot of lighthouse image PSNR (dB) results from Table 3 using 200 simulated frames with ${\sigma }_{\eta }=2.0$ and $M=3$.

Table 5

Optimum FIF SR parameters for lighthouse image results.

Fig. 11

Lighthouse image: (a) truth, (b) single-frame bicubic interpolation of a degraded image with no turbulence, (c) $C_n^2=1\times {10}^{-15}{\mathrm{m}}^{-2/3}$, and (d) $C_n^2=1\times {10}^{-14}{\mathrm{m}}^{-2/3}$. Degraded images have additive Gaussian noise with ${\sigma }_{\eta }=2$ digital units and downsampling of $M=3$.

Fig. 12

FIF SR restoration results with lighthouse image. (a) No turbulence and no jitter, (b) no turbulence with jitter, (c) $C_n^2=1\times {10}^{-15}{\mathrm{m}}^{-2/3}$ and no jitter, (d) $C_n^2=1\times {10}^{-15}{\mathrm{m}}^{-2/3}$ with jitter, (e) $C_n^2=1\times {10}^{-14}{\mathrm{m}}^{-2/3}$ and no jitter, (f) $C_n^2=1\times {10}^{-14}{\mathrm{m}}^{-2/3}$ with jitter.

The quantitative results show that the FIF SR method provides a significant boost over simple methods, such as single-frame interpolation and simple averages. The boost is seen over a wide range of turbulence levels. Another very interesting phenomenon can be seen in Fig. 10. Note that when platform jitter is turned off, increasing the $C_n^2$ turbulence level from zero actually increases the PSNR of the SR output. This may seem counterintuitive, as turbulence is generally considered to be a degradation, and not a benefit. However, what we see is that the wavefront tilt variance provided by light turbulence acts to shift the image relative to the camera, allowing for the necessary sampling diversity for SR. This phenomenon was first described by Fishbain et al.¹⁰ and Yaroslavsky et al.¹¹ As the turbulence level gets higher, we see in Fig. 10 that the degrading impact of the turbulence outweighs this sampling benefit, and the PSNR drops. Also, note that this effect is not seen when there is camera platform jitter. This is because the platform jitter provides sampling diversity more effectively than the turbulence and without the turbulence blurring. Platform jitter outperforms the corresponding no jitter scenario, and turbulence is never a benefit when platform jitter is present. As turbulence levels increase, we do see the benefit of jitter diminish as SR becomes increasingly difficult. At high turbulence levels, very little signal energy is available above the folding frequency to be recovered. Another point of interest in Tables 3 and 4, as well as Fig. 10, is that including the interpolation OTF into the Wiener filter is a clear benefit. When the optimum $\beta$ for restoration is used, this always outperforms the “no $\beta$ restoration” case. To the best of our knowledge, no simulation study and quantitative error analysis of this kind for joint TM and SR has been reported in the literature previously.

Subjective analysis of the images in Figs. 11 and 12 appears in line with the quantitative results. Figures 12(a) and 12(b) show the FIF SR output with no turbulence. Figure 12(a) has no camera jitter and Fig. 12(b) includes camera jitter. With no jitter and no turbulence in Fig. 12(a), there is a lack of sampling diversity and little is achieved in the way of true SR. Moiré patterns on the fence and jagged edges on the roof line are still quite visible here. On the other hand, Fig. 12(b) shows significant SR enhancement on the fence and roof edges. This result represents traditional multiframe SR without turbulence and is well understood.²

Perhaps, the most interesting result is that in Fig. 12(c). Here, we have light turbulence and no camera jitter. This result appears to be far better than the no-turbulence no-jitter case in Fig. 12(a). A significant amount of aliasing reduction is exhibited here, as a result of entirely turbulence-induced motion. This result is nearly as good as that with additional camera jitter in Fig. 12(d). Finally, at high turbulence levels, platform jitter makes little difference, as seen by comparing Figs. 12(e) and 12(f). It is also clear, at this high level of turbulence, that details in the fence and elsewhere are lost, but not aliased, as a result of the low-pass filtering from the turbulence.

4.2.

Real Data

In this subsection, we present results for real data using three different sensors. A summary of the data and processing parameters for each of the datasets is provided in Table 6. The Fried parameters and noise-to-signal ratios (NSRs) have been selected based on subjective evaluation of the results. In the case of the truck sequence, an edge target is used to estimate $r_0$.

Table 6

Real image data parameters for the results in Sec. 4.2.

The first dataset (coaster) is shown in Fig. 13. These data are from a midwave infrared (MWIR) sensor and show a portion of a wooden roller coaster with a significant amount of turbulence and aliasing ($Q=0.40$ and $\tilde{Q}=1.20$). Single-frame bicubic interpolation images of two regions of interest (ROIs) are shown in Figs. 13(a) and 13(c). The corresponding FIF SR outputs are shown in Figs. 13(b) and 13(d), using the parameters listed in Table 6. The impact of turbulence is quite evident in the inputs frames, such as the warped handrail visible in Fig. 13(a) near the power line pole. At the same time, aliasing artifacts are also clearly visible in the form of Moiré patterns on the groups of suspended wires. In contrast, the handrail appears to have corrected geometry in Figs. 13(b) and 13(d), as a result of the FIF SR processing. Furthermore, the individual suspended wires are clearly distinguishable and free from aliasing artifacts in the FIF SR output. We believe these data provide an excellent demonstration of successful joint TM and SR.

Fig. 13

Roller coaster infrared image sequence results. (a) ROI 1 single frame bicubic interpolation, (b) ROI 1 FIF SR output, (c) ROI 2 single-frame interpolation, and (d) ROI 2 FIF SR output. Sensor and algorithm parameters are listed in Table 5.

The next dataset (truck) shows a truck and bar target in heavy turbulence imaged with a near IR sensor. In contrast with the previous example, this sensor is nearly Nyquist sampled under diffraction-limited conditions ($Q=1.95$). With turbulence added, no significant aliasing is expected ($\tilde{Q}=3.44$). However, we still apply upsampling of $M=2$ for our restoration. Single-frame bicubic interpolation is shown in Fig. 14(a), the FIF SR output is shown with affine registration in Fig. 14(b). The result using block matching algorithm (BMA) optical flow registration⁹ is shown in Fig. 14(c). Both restorations appear to be significantly better than the single-frame interpolation, but the bar target with BMA does appear noticeably better because of the high level of local warping. This finding is consistent with previous studies using these data.⁹ One noticeable artifact is the ringing on the truck cab. This is a result of the Wiener filter operating on solar glint. Ringing is also present, to a lesser degree, in the restored images near strong edges. We are currently exploring methods for reducing these artifacts. One possibility is to employ an adaptive Wiener filter with spatially varying NSR.^34,35

Fig. 14

Truck infrared image sequence results. (a) Single-frame bicubic interpolation, (b) FIF SR using global affine registration, and (c) FIF SR using BMA registration. Sensor and algorithm parameters are listed in Table 5.

The final dataset (airborne) is a MWIR dataset of a bar target acquired from an airborne platform.³⁶ These data have been used for SR studies in prior work,^17,21 but without considering turbulence. The bar target is a series of four bar patterns. The scaling factor between bar groups is designed to be $2^{1/6}$. In contrast to the truck sequence, aliasing is a much more significant problem than turbulence in the airborne data with $Q=0.47$ and $\tilde{Q}=0.20$. Since the platform is moving between frames and the scene is approximately planar, we use global perspective registration.²¹ Single-frame bicubic interpolation is shown in Fig. 15(a) with $M=3$. The corresponding FIF SR output is shown in Fig. 15(b). The Moiré patterns are quite evident within the bicubic image in Fig. 15(a). After FIF SR processing, there appears to be an approximately $2\times$ resolution enhancement based on the resolvable bar patterns.

Fig. 15

Airborne image sequence results. (a) Single-frame bicubic interpolation and (b) FIF SR output. Sensor and algorithm parameters are listed in Table 5.