A New Adaptive Video Super-Resolution Algorithm With Improved Robustness to Innovations

TL;DR

This paper introduces an adaptive video super-resolution algorithm that enhances robustness to outliers, outperforming traditional methods in quality while maintaining low computational costs.

Contribution

A novel cost function and two algorithms are proposed, improving robustness to innovation outliers in video super-resolution reconstruction.

Findings

Outperforms traditional LMS in robustness to outliers

Maintains computational efficiency comparable to R-LMS

Competitive with state-of-the-art SRR methods

Abstract

In this paper, a new video super-resolution reconstruction (SRR) method with improved robustness to outliers is proposed. Although the R-LMS is one of the SRR algorithms with the best reconstruction quality for its computational cost, and is naturally robust to registration inaccuracies, its performance is known to degrade severely in the presence of innovation outliers. By studying the proximal point cost function representation of the R-LMS iterative equation, a better understanding of its performance under different situations is attained. Using statistical properties of typical innovation outliers, a new cost function is then proposed and two new algorithms are derived, which present improved robustness to outliers while maintaining computational costs comparable to that of R-LMS. Monte Carlo simulation results illustrate that the proposed method outperforms the traditional and regularized versions of LMS, and is competitive with state-of-the-art SRR methods at a much smaller computational cost.

Tables

Table 1. TABLE I: Memory cost of the algorithms.

	Memory
LMS	$M^2+\left|𝐇\right|/M^2$
R-LMS	$M^2+\left(\left|𝐇\right|+\left|𝐒\right|\right)/M^2$
TSR-LMS	$2M^2+\left(\left|𝐇\right|+\left|𝐒\right|+\left|𝐌\right|+\left|𝐐\right|\right)/M^2$
LTSR-LMS	$2M^2+\left(\left|𝐇\right|+\left|𝐒\right|+\left|𝐐\right|\right)/M^2$

Table 2. TABLE II: Computational cost per iteration of the algorithms (additions and multiplications, surplus additions, and re-samplings were considered).

	Operations
LMS	$3\left|𝐇\right|+2M^2$
R-LMS	$3\left|𝐇\right|+2\left|𝐒\right|+2M^2$
TSR-LMS	$3\left|𝐇\right|+2\left|𝐒\right|+2\left|𝐐\right|+\left|𝐌\right|+2M^2$
LTSR-LMS	$3\left|𝐇\right|+2\left|𝐒\right|+2\left|𝐐\right|+3M^2$

Table 3. TABLE III: Parameter values used on the simulations with the with outlier-free sequences.

	LMS	R-LMS	TSR-LMS	LTSR-LMS
$\mu$	2	2.75	1.15	3
$\alpha$	–	5$\times {10}^{-4}$	1.5$\times {10}^{-4}$	1$\times {10}^{-4}$
${\alpha }_{\text{T}}$	–	–	82	0.02

Table 4. TABLE IV: Parameter values used on the simulations considering the presence of outliers.

	LMS	R-LMS	TSR-LMS	LTSR-LMS
$\mu$	4.7	4.2	2.2	3.4
$\alpha$	–	40$\times {10}^{-4}$	18$\times {10}^{-4}$	1$\times {10}^{-4}$
${\alpha }_{\text{T}}$	–	–	16	0.017

Table 5. TABLE V: Average processing time per frame for the videos in Example 4.

Bicubic	LMS	R-LMS	TSR-LMS	LTSR-LMS	Bayesian [15]	CNN [18]
1.56 s	0.46 s	0.53 s	0.60 s	0.62 s	49.68 s	81.72 s

Table 6. TABLE VI: Average PSNR (in decibels) for the videos in Example 4.

	Bear	Bus	Elephant	Horse	Paragliding	Sheep	Mean
Bicubic	28.43	32.66	30.56	29.91	35.58	25.02	30.36
LMS	32.78	31.06	32.51	28.73	33.24	30.46	31.46
R-LMS	33.38	33.09	33.57	30.89	35.10	30.37	32.73
TSR-LMS	34.00	34.46	34.16	32.61	36.48	30.99	33.78
LTSR-LMS	33.94	34.89	34.25	33.26	36.59	30.75	33.95
Bayesian	30.52	33.42	32.08	31.72	34.87	29.60	32.04
CNN	32.21	33.89	32.96	32.75	34.54	29.66	32.67

Equations

$$\mathbf{y}(t)=\mathbf{DHx}(t)+\mathbf{e}(t)$$

$$\mathbf{x}(t)=\mathbf{G}(t)\mathbf{x}(t-1)+\mathbf{s}(t)$$

$$\mathbf{ϵ}(t)=\mathbf{y}(t)-\mathbf{DH}\hat{\mathbf{x}}(t)$$

$${\mathcal{L}}_{\text{R-MS}}(t)=E\{\Vert \mathbf{y}(t)-\mathbf{DH}\hat{\mathbf{x}}(t){\Vert }^2|\hat{\mathbf{x}}(t)\}+\alpha \Vert \mathbf{S}\hat{\mathbf{x}}(t){\Vert }^2$$

$$\nabla {\mathcal{L}}_{\text{R-MS}}(t)=-2{\mathbf{H}}^{\text{T}}{\mathbf{D}}^{\text{T}}\{E[\mathbf{y}(t)]-\mathbf{DH}\hat{\mathbf{x}}(t)\}+2\alpha {\mathbf{S}}^{\text{T}}\mathbf{S}\hat{\mathbf{x}}(t)$$

$${\hat{\mathbf{x}}}_{k+1}(t)={\hat{\mathbf{x}}}_k(t)-\frac{\mu }{2}\nabla {\mathcal{L}}_{\text{R-MS}}(t),k=0,1,\dots ,K-1$$

$${\hat{\mathbf{x}}}_{k+1}(t)=$$

$$-\alpha \mu {\mathbf{S}}^{\text{T}}\mathbf{S}{\hat{\mathbf{x}}}_k(t)\text{,}k=0,1,\dots ,K-1$$

$${\hat{\mathbf{x}}}_0(t+1)=\mathbf{G}(t+1){\hat{\mathbf{x}}}_K(t)\text{.}$$

$$\displaystyle\hat{\mathbf{x}}_{k+1}(t)=\underset{\mathbf{z}}{\arg\min}\Big{\{}\mathcal{L}_{\tiny{\text{R-MS}}}(\hat{\mathbf{x}}_{k}(t))$$

$$\displaystyle\hskip 12.91663pt+\big{(}\mathbf{z}-\hat{\mathbf{x}}_{k}(t)\big{)}^{\top}\nabla\mathcal{L}_{\tiny{\text{R-MS}}}(\hat{\mathbf{x}}_{k}(t))+\frac{1}{\mu}\left\|\mathbf{z}-\hat{\mathbf{x}}_{k}(t)\right\|^{2}\Big{\}}$$

$${\hat{\mathbf{x}}}_{k+1}(t)=$$

$$\displaystyle+\frac{1}{\mu}\left\|\mathbf{z}-\hat{\mathbf{x}}_{k}(t)\right\|^{2}\Big{\}}$$

$${\hat{\mathbf{x}}}_1(t)=\underset{\mathbf{z}}{\mathrm{argmin}}$$

$$\displaystyle-2\mathbf{z}^{\text{\tiny T}}\mathbf{H}^{\text{\tiny T}}\mathbf{D}^{\text{\tiny T}}\boldsymbol{\epsilon}_{0}(t)+\frac{1}{\mu}\left\|\mathbf{z}-\hat{\mathbf{x}}_{0}(t)\right\|^{2}\Big{\}}$$

$$\displaystyle\hat{\mathbf{x}}_{1}(t){}\approx{}\underset{\mathbf{z}}{\arg\min}\Big{\{}$$

$$\displaystyle+\frac{1}{\mu}\left\|\mathbf{z}-\mathbf{G}(t)\mathbf{x}(t-1)\right\|^{2}\Big{\}}.$$

$$\displaystyle\hat{\mathbf{x}}_{1}(t){}\approx{}\underset{\mathbf{z}}{\arg\min}\Big{\{}$$

$$\displaystyle+\frac{1}{\mu}\left\|\mathbf{z}-\mathbf{x}(t)+\mathbf{s}(t)\right\|^{2}\Big{\}}.$$

$$\displaystyle\hat{\mathbf{x}}_{1}(t){}\approx{}\underset{\mathbf{z}}{\arg\min}\Big{\{}$$

$$\displaystyle+\frac{1}{\mu}\left\|\mathbf{z}-\mathbf{x}(t)\right\|^{2}\Big{\}}.$$

$$\Vert {\mathbf{ϵ}}_0(t)\Vert \simeq$$

$$\approx$$

$$\displaystyle\hat{\mathbf{x}}_{1}(t){}\approx{}\underset{\mathbf{z}}{\arg\min}\,\,\Big{\{}$$

$$\displaystyle+\frac{1}{\mu}\left\|\mathbf{z}-\big{(}\mathbf{x}(t)-\mathbf{s}(t)\big{)}\right\|^{2}\Big{\}}$$

$$\Vert {\mathbf{ϵ}}_0(t)\Vert \simeq$$

$$=$$

$$\displaystyle\hat{\mathbf{x}}_{1}(t){}\approx{}2\,\,\underset{\mathbf{z}}{\arg\min}\Big{\{}$$

$$\displaystyle+\frac{1}{2\mu}\left\|\mathbf{z}-\hat{\mathbf{x}}_{0}(t)\right\|^{2}\Big{\}}\,\text{.}$$

$$\displaystyle\hat{\mathbf{x}}_{k+1}(t){}={}\underset{\mathbf{z}}{\arg\min}\big{\{}\mathcal{L}_{\tiny{\text{R-MS}}}(\hat{\mathbf{x}}_{k}(t))$$

$$\displaystyle\hskip 8.61108pt+\big{(}\mathbf{z}-\hat{\mathbf{x}}_{k}(t)\big{)}^{\!\top}\nabla\mathcal{L}_{\tiny{\text{R-MS}}}(\hat{\mathbf{x}}_{k}(t))+\frac{1}{\mu}\left\|\mathbf{z}-\hat{\mathbf{x}}_{k}(t)\right\|^{2}$$

$$\displaystyle\hskip 8.61108pt+\frac{1}{\alpha_{\tiny{\text{T}}}\mu}\left\|\mathbf{z}-\mathbf{G}(t)\hat{\mathbf{x}}(t-1)\right\|^{2}\big{\}}$$

$$\displaystyle\|\hat{\mathbf{x}}_{k+1}(t)-\mathbf{G}(t)\hat{\mathbf{x}}(t-1)\|\approx\|\hat{\mathbf{x}}_{k+1}(t)-\big{(}\mathbf{x}(t)-\mathbf{s}(t)\big{)}\|$$

$$\displaystyle\big{\|}\mathbf{Q}\big{(}\hat{\mathbf{x}}_{k+1}(t)-\mathbf{G}(t)\hat{\mathbf{x}}(t-1)\big{)}\big{\|}$$

$$\begin{array}{rl}\mathbf{Q}\mathbf{x}(t)\sim & \text{details}\\ \mathbf{Q}\mathbf{s}(t)\sim & 0\end{array}$$

$$\Vert \mathbf{Q}{\hat{\mathbf{x}}}_{k+1}(t)-\mathbf{Qx}(t)+\mathbf{Qs}(t)\Vert \text{.}$$

$$\Vert \mathbf{Q}{\hat{\mathbf{x}}}_{k+1}(t)-\mathbf{Qx}(t)+\mathbf{Qs}(t)\Vert \approx \Vert \mathbf{Q}{\hat{\mathbf{x}}}_{k+1}(t)-\mathbf{Qx}(t)\Vert$$

Full text

I Introduction

Super-resolution reconstruction (SRR) is a well established approach for digital image quality improvement. SRR consists basically in combining multiple low-resolution (LR) images of the same scene or object in order to obtain one or more images of higher resolution (HR), outperforming physical limitations of image sensors. Applications for SRR are numerous and cover diverse areas, including the reconstruction of satellite images in remote sensing, surveillance videos in forensics and images from low-cost digital sensors in standard end-user systems. References [2, 3] review several important concepts and initial results on SRR.

SRR algorithms usually belong to one of two groups. Image SRR algorithms, which reconstruct one HR image from multiple observations, and video SRR algorithms, which reconstruct an entire HR image sequence. Video SRR algorithms often include a temporal regularization that constrains the norm of the changes in the solution between adjacent time instants [4, 5, 6, 7, 8, 9]. This introduces information about the correlation between adjacent frames, and tends to ensure video consistency over time, improving the quality of the reconstructed sequences [10].

A typical characteristic of SRR algorithms is the high computational cost. Recent developments in both image and video SRR have been mostly directed towards achieving increased reconstruction quality, either using more appropriate a priori information about the underlying image and the acquisition process or learning the relationship between LR and HR images from a set of training data. Examples are the non-parametric methods based on spatial kernel regression [11], non-local methods [12, 13], variational Bayesian methods [14, 15] and, more recently, deep-learning-based methods [16, 17, 18].

Although these techniques have led to considerable improvements in the quality of state of the art SRR algorithms, such improvements did not come for free. The computational cost of these algorithms is very high, which makes them unsuitable for real-time SRR applications. While deep-learning methods can be significantly faster than kernel or non-local methods, they rely on extensive training procedures with large amounts of data. Also, the training must be repeated whenever the test conditions change, or their performance may degrade significantly [19]. Moreover, registration errors plague most SRR algorithms. This motivates the use of nonlinear cost functions and more complex methodologies, which further contributes to increase the computational cost of the corresponding algorithms [20, 21].

While these traditional SRR methods achieve good reconstruction results [22], real-time video SRR applications require simple algorithms. This limitation prompted a significant interest towards developing low complexity SRR algorithms. The regularized least mean squares (R-LMS) [23, 24] is one notable example among the simpler SRR algorithms. Even though other low-complexity algorithms have been proposed for global translational image motions, such as the $\text{L}_{1}$ norm estimator in [25] and the adaptive Wiener filter of [26], their computational complexity is still not competitive with that of the R-LMS. The R-LMS presents a reasonable reconstruction quality and follows a systematic mathematical derivation. This enables a formal characterization of its behavior [27] and the specification of well defined design methodologies [28]. Furthermore, the R-LMS is also naturally robust to registration errors, which were shown to have a regularizing effect in the algorithm [29]. This makes the quality of the R-LMS algorithm in practical situations to be competitive even with that of costly and elaborated algorithms like [14], as illustrated through an example in Figure 1.

Unfortunately, the performances of these simple algorithms, including the R-LMS, tend to be heavily affected by the presence of outliers such as large innovations caused by moving objects and sudden scene changes. This can lead to reconstructed sequences of worse quality than that of the observations themselves [7]. Common strategies for obtaining robust algorithms often involve the optimization of nonlinear cost functions [4, 25, 7], and thus present a computational cost that is not comparable to that of algorithms like the R-LMS. Interpolation algorithms might seem to be a reasonable option, as their performance is not affected by motion related outliers. However, they do not offer a quality improvement comparable to SRR methods [30, 7]. Therefore, it is of interest to develop video SRR algorithms that combine good quality, robustness to outliers and a low computational cost.

This paper proposes a new adaptive video SRR algorithm with improved robustness to outliers when compared to the R-LMS algorithm. The contributions of this paper include:

1. A new interpretation of the R-LMS update equation as the proximal regularization of the associated cost function, linearized about the previous estimate, which leads to a better understanding of its quality performance and robustness in different situations. Using this representation we show that the slow convergence rate of the R-LMS algorithm (typical of stochastic gradient algorithms) establishes a trade-off between its robustness to outliers and the achievable reconstruction quality.
2. A simple model for the statistical properties of typical innovation outliers in natural image sequences is developed, which points to the desirable properties of the proposed technique.
3. A new cost function is then proposed to address the identified problems and two new adaptive algorithms are derived called Temporally Selective Regularized LMS (TSR-LMS) and Linearized Selective Regularized LMS (LTSR-LMS), which present improved robustness and similar quality at a computational cost comparable to that of the R-LMS algorithm.

The paper is organized as follows. In Section II, the image acquisition model is defined. In Section III, the R-LMS algorithm [23, 24] is derived as a stochastic gradient solution to the image estimation problem [29]. In Section IV, an intuitive interpretation of the R-LMS behavior is presented using the proximal-point cost function representation of the gradient descent iterative equation. In Section V, a new robust cost function is proposed based on a statistical model for the innovations, and two adaptive algorithms are derived. In chapter VII, computer simulations are performed to assess the performance of the algorithms. The conclusions are presented in Section VIII.

II Image Acquisition Model

Given the $N\times N$ matrix representation of an LR (observed) digital image $\mathbf{Y}(t)$, and an $M\times M$ ($M>N$) matrix representation of the original HR digital image $\mathbf{X}(t)$, the acquisition process can be modeled as [2]

[TABLE]

where vectors $\mathbf{y}(t)$ ($N^{2}\times 1$) and $\mathbf{x}(t)$ ($M^{2}\times 1$) are the lexicographic representations of the degraded and original images, respectively, at discrete time instant $t$. $\mathbf{D}$ is an $N^{2}\times M^{2}$ decimation matrix and models the sub-sampling taking place in the sensor. $\mathbf{H}$ is an $M^{2}\times M^{2}$ matrix, assumed known, that models the blurring in the acquisition process111Since $\mathbf{H}$ is assumed to be estimated independently from the SRR process, the extension of the results in this work to the case of a time-varying $\mathbf{H}(t)$ is straightforward if an online estimation algorithm for $\mathbf{H}$ is employed. Here, a possible dependence of $\mathbf{H}$ on $t$ is omitted for notation simplicity.. The $N^{2}\times 1$ vector $\mathbf{e}(t)$ models the observation (electronic) noise, whose properties are assumed to be determined from camera tests. The dynamics of the input signal is modeled by [23]

[TABLE]

where $\mathbf{G}(t)$ is the warp matrix that describes the relative displacement from $\mathbf{x}(t-1)$ to $\mathbf{x}(t)$. Vector $\mathbf{s}(t)$ models the innovations in $\mathbf{x}(t)$.

III The R-LMS SRR algorithm

Several SRR solutions are based on the minimization of the estimation error (see [2] and references therein)

[TABLE]

where $\hat{\mathbf{x}}(t)$ is the estimated HR image, and $\boldsymbol{\epsilon}(t)$ can be interpreted as the estimate of $\mathbf{e}(t)$ in (1). The LMS algorithm attempts to minimize the mean-square value of the $\text{L}_{2}$ norm of (3) conditioned on the estimate $\hat{\mathbf{x}}(t)$ [23, 27]. Thus, it minimizes the cost function $\mathbf{J}_{\text{MS}}(t)=\operatorname{E}\{\|\boldsymbol{\epsilon}(t)\|^{2}\,|\,\hat{\mathbf{x}}(t)\}$.

Since natural images are known to be intrinsically smooth, this a priori knowledge can be added to the estimation problem in the form of a regularization to the LMS algorithm by constraining the solution that minimizes $\mathbf{J}_{\text{MS}}(t)$. The R-LMS algorithm [29, 24] is the stochastic gradient version of the gradient descent search for the solution to the following optimization problem

[TABLE]

where $\mathbf{S}$ is the Laplacian operator [31, p. 182]. Note that the performance surface in (4) is defined for each time instant $t$, and the expectation is taken over the ensemble.

Following the steepest descent method, the HR image estimate is updated in the negative direction of the gradient

[TABLE]

and thus the iterative update of $\hat{\mathbf{x}}(t)$ for a fixed value of $t$ is given by

[TABLE]

where $K\in\mathbb{Z^{+}}$ is the number of iterations of the algorithm, and $\mu$ is the step size used to control the convergence speed. The factor 1/2 is just a convenient scaling.

Using the instantaneous estimate of (5) in (6) yields

[TABLE]

which is the R-LMS update equation for a fixed value of $t$. The time update of (7) is based on the signal dynamics (2), and is performed by the following expression [23]:

[TABLE]

Between two time updates, (7) is iterated for $k=0,\ldots,K-1$. The estimate $\hat{\mathbf{x}}(t)$ at a given time instant $t$ is then given by $\hat{\mathbf{x}}(t)=\hat{\mathbf{x}}_{K}(t)$.

IV R-LMS Performance in the Presence of Outliers

The R-LMS algorithm is computationally efficient when implemented with few stochastic gradient iterations (small $K$) per time instant $t$. However, the R-LMS algorithm is derived under the assumption that the solution $\mathbf{x}(t)$ is only slightly perturbed between time instants, a characteristic known as the minimum disturbance principle. This assumption is satisfied in the R-LMS algorithm close to steady-state, when the estimate $\hat{\mathbf{x}}(t)$ has already achieved a reasonable quality (i.e. $\hat{\mathbf{x}}(t)\simeq\mathbf{x}(t)$). Then, the initialization $\hat{\mathbf{x}}_{0}(t+1)$ determined by (8) for the next time instant will be already relatively close to the optimal solution, what explains the good steady-state performance of the algorithm even with a small value of $K$ in (7).

The situation is significantly different in the occurrence of innovation outliers. Experience with the R-LMS algorithm shows that the slow convergence of (7) as a function of $k$ tends to degrade the quality of the super-resolved images in the presence of innovation outliers. Visible artifacts tend to be created, and the reconstructed images may end up being of inferior quality even when compared to the originally observed LR images. This significantly compromises the quality that can be achieved in real-time super-resolution of video sequences, not just using R-LMS but most of the existing low-complexity algorithms [7]. Super-resolution algorithms designed to be robust under the influence of innovations tend to impose a high computational cost, making them unsuitable for real time applications [25, 28]. In the following we examine the R-LMS recursion under a new light, what leads to a mathematically motivated explanation for its lack of robustness to outliers. This explanation will then motivate the proposition of more robust stochastic video SRR algorithms.

An interesting interpretation of the R-LMS algorithm is possible if we view each iteration of the gradient algorithm (6) (for a fixed value of $t$) as a proximal regularization of the cost function $\mathcal{L}_{\text{R-MS}}(t)$ linearized about the estimation of the previous iteration $\hat{\mathbf{x}}_{k}(t)$. Proceeding as in [32, Section 2.2] or [33, p. 546], the gradient iteration (6) can be written as

[TABLE]

which, using the previous expressions, yields

[TABLE]

where $\operatorname{E}[\boldsymbol{\epsilon}_{k}(t)]$ is the expected value of the observation error (3) conditioned on $\hat{\mathbf{x}}(t)=\hat{\mathbf{x}}_{k}(t)$. This equivalence can be verified by differentiating the expression in the curl brackets in (9) with respect to $\mathbf{z}$, setting it equal to zero and solving for $\mathbf{z}=\hat{\mathbf{x}}_{k+1}(t)$.

Now, the presence of the squared norm within the curl brackets in (9) and  (10) means that the optimization algorithm seeks $\hat{\mathbf{x}}_{k+1}(t)$ that minimizes the perturbation $\hat{\mathbf{x}}_{k+1}(t)-\hat{\mathbf{x}}_{k}(t)$ at each iteration. Evidencing this property leads to a more detailed understanding of the dynamical behavior of the algorithm, its robustness properties and the reconstruction quality it provides. For instance, this constraint on the perturbation of the solution explains how the algorithm tends to preserve details in $\hat{\mathbf{x}}(t)$ that have been estimated in previous time iterations and that are present in $\hat{\mathbf{x}}(t-1)$. However, this same term also opposes changes from $\hat{\mathbf{x}}_{k}(t)$ to $\hat{\mathbf{x}}_{k+1}(t)$, slowing down the reduction of the observation error from $\boldsymbol{\epsilon}_{k}(t)$ to $\boldsymbol{\epsilon}_{k+1}(t)$ since changes in $\boldsymbol{\epsilon}_{k}(t)$ require changes in $\hat{\mathbf{x}}_{k}(t)$. Therefore, this algorithm cannot simultaneously achieve a fast convergence rate and preserve the super-resolved details for the practically important case of a small number of iterations per time instant (small $K$). The time sequence of reconstructed images will either converge fast and yield low temporal correlation between time estimates (leading to a solution that approaches an interpolation of $\mathbf{y}(t)$), or will converge slowly and yield a highly correlated image sequence with better quality in the absence of innovation outliers. The occurrence of outliers will result in a significant deviation from the desired signal.

To illustrate this behavior, consider for instance that the reconstructed image sequence at time instant $t-1$ is reasonably close to the real (desired) sequence, i.e. $\hat{\mathbf{x}}(t-1)\simeq\mathbf{x}(t-1)$. If we consider the video sequence to be only slightly perturbed at the next time instant such that $\|\mathbf{s}(t)\|\approx 0$ in (2), the first iteration of (10) ($k=0$ at time $t$) is given by

[TABLE]

which, using (8) with $\hat{\mathbf{x}}(t-1)\simeq\mathbf{x}(t-1)$, yields

[TABLE]

Using (2) in (12),

[TABLE]

Finally, assuming $\|\mathbf{s}(t)\|\approx 0$ (no outlier in $\mathbf{x}(t)$)

[TABLE]

Now, using $\|\mathbf{s}(t)\|\approx 0$, $\hat{\mathbf{x}}(t-1)\simeq\mathbf{x}(t-1)$, (1) and (3), the norm of $\boldsymbol{\epsilon}_{0}(t)$ can be approximated by

[TABLE]

which is small since the energy of the observation noise is much smaller than that of registration errors and outliers in most practical applications [2, 28]. The first term in the r.h.s. of (14) is due to the regularization an promotes the smoothness of the solution. Hence, $\alpha$ should be small to avoid compromising the estimation of the details of $\mathbf{x}(t)$. The last term will promote a solution close to $\mathbf{x}(t)$, especially for small values of $\mu$. Then, for reasonably small values of $\alpha$ and $\mu$, $\|\mathbf{s}(t)\|\approx 0$ and $\hat{\mathbf{x}}(t-1)\simeq\mathbf{x}(t-1)$, the first and second terms in (14) can be neglected (i.e. $|2\alpha\mathbf{z}^{\text{\tiny T}}\mathbf{S}^{\text{\tiny T}}\mathbf{S}\hat{\mathbf{x}}_{0}(t)-2\mathbf{z}^{\text{\tiny T}}\mathbf{H}^{\text{\tiny T}}\mathbf{D}^{\text{\tiny T}}\boldsymbol{\epsilon}_{0}(t)|\ll\frac{1}{\mu}\left\|\mathbf{z}-\mathbf{x}(t)\right\|^{2}$). Then, the solution of the optimization problem (14) will converge to a vector $\hat{\mathbf{x}}_{1}(t)\approx\mathbf{x}(t)$. The same reasoning can be extended to the remaining iterations for $k=2,\ldots,K-1$, which shows that, for $\|\mathbf{s}(t)\|\approx 0$, the algorithm will lead to a reconstructed image of good quality $\hat{\mathbf{x}}_{K}(t)\simeq\mathbf{x}(t)$. This explains how the R-LMS algorithm preserves the reconstructed content in time and extracts information from the different observations, attaining good reconstruction results for well behaved sequences, i.e. in the absence of large innovations.

Now, let’s consider the presence of a significant innovation outlier at time $t$, while still assuming a good reconstruction result at time $t-1$ (i.e. $\hat{\mathbf{x}}(t-1)\simeq\mathbf{x}(t-1)$). Due to the outlier at time instant $t$, $\mathbf{s}(t)$ in (2) will have a significant energy. Then, repeating (14) without the assumption $\|\mathbf{s}(t)\|\approx 0$ yields

[TABLE]

where now the observation error is given by

[TABLE]

This result clearly shows that, for large $\|\mathbf{s}(t)\|$, the solution $\hat{\mathbf{x}}_{1}(t)$ of (16) can be considerably away from $\mathbf{x}(t)$, as it should contain information introduced by $\mathbf{s}(t)$. For the desirable case of small $K$, this estimation error will hardly be significantly reduced from $\hat{\mathbf{x}}_{0}(t)$ to $\hat{\mathbf{x}}_{K}(t)$ due to the slow convergence of the gradient based recursion. This explains the poor transient performance of the algorithm in the presence of outliers.

Performance improvement in the presence of outliers could be sought by increasing the value of $\mu$ to reduce the influence of the term \frac{1}{\mu}\|\mathbf{z}-\big{(}\mathbf{x}(t)-\mathbf{s}(t)\big{)}\|^{2} in (16). However, $\mu$ cannot be made arbitrarily large for stability reasons. Hence, improvement has to come from increasing the importance of the first term in (16).

Neglecting the contribution of $\mathbf{e}(t)$ in (IV), expression (16) can be written as

[TABLE]

First, note that $\hat{\mathbf{x}}_{0}(t)=\mathbf{G}(t)\hat{\mathbf{x}}_{K}(t-1)$ does not include information on the outlier $\mathbf{s}(t)$, as it is not present in $\mathbf{x}(t-1)$. Moreover, we have already seen that the first term promotes the smoothness of the solution. Thus, increasing the value of $\alpha$ in an attempt to speed up convergence in the presence of large innovations by reducing the influence of the last term in (16) will also reduce the temporal correlation of the estimated image sequence, resulting in an overly blurred solution with lower quality in the absence of outliers. Hence, the solution $\hat{\mathbf{x}}(t)$ can hardly approach the desired solution $\mathbf{x}(t)$ in few iterations (small $K$). If, however, $K$ is made sufficiently large, the solution can adapt to track the innovations even with a large weighting for the term \frac{1}{\mu}\|\mathbf{z}-\big{(}\mathbf{x}(t)-\mathbf{s}(t)\big{)}\|^{2}. The algorithm could then achieve and maintain a good reconstruction quality both with and without the presence of outliers, but at a prohibitive computational cost.

IV-A Illustrative example

The behavior of the R-LMS algorithm is illustrated in the following example, where we consider the reconstruction of a synthetic video sequence generated through small translational displacements of a $32\times 32$ window over a larger natural image. At a specific time instant during the video (the 32nd frame), an outlier is introduced by adding a black square of size $16\times 16$ to the scene. The square remain in the scene for $3$ frames, before disappearing again. The HR sequence is convolved with a $3\times 3$ uniform blurring mask and down-sampled by a factor of $2$. Finally, a white Gaussian noise with variance $10$ is added to generate the low resolution video.

The R-LMS algorithm is applied to super-resolve the synthetic LR videos generated. The mean square error (MSE) between the original and reconstructed sequences is estimated by averaging the results from 50 realizations. To illustrate the trade-offs between the effects of different values of the step size and regularization parameter in the cost function, we reconstructed the sequence with both $\alpha=2\times 10^{-4}$ and $\alpha=100\times 10^{-4}$, for $\mu=4$. To evaluate the effect of using different numbers of iterations per time interval, we ran the algorithm with $K=2$ and $K=100$. The MSE is depicted in Figure 2-(a). From the two curves for $K=2$ one can verify that a large value for $\alpha$ (red curve) reduces the MSE in the presence of the outlier, while the greater temporal correlation induced by a small value of $\alpha$ (black curve) tends to reduce the error for small innovations and to increase it in the presence of an outlier. Comparing the blue ($K=100$) and the black ($K=2$) curves, both for $\alpha=2\times 10^{-4}$ and $\mu=4$, one verifies that the MSE can be substantially decreased by employing the R-LMS algorithm with a large $K$. The MSE is smaller than that obtained for $K=2$ both for small and for large innovations. This performance improvement is because the algorithm is allowed to converge slowly for each time interval. Figure 2-(b) shows the MSE as a function of $k$ for time instant $t=32$, when the outlier is present. These results illustrate the property that a large value of $K$ is necessary to achieve a significant MSE reduction for a fixed value of $t$.

In the light of the aforementioned limitations of the R-LMS algorithm, it is of interest to devise an algorithm that performs better both in terms of robustness and quality at a reasonable computational cost.

V Improving the Robustness to Innovations

A temporal regularization of adaptive algorithms such as the R-LMS that constrains the value of $\|\hat{\mathbf{x}}(t)-\mathbf{G}(t)\hat{\mathbf{x}}(t-1)\|$ in the SRR cost function [10, 4, 7] can be interpreted as the application of the well known least perturbation or minimum disturbance principle. This principle states that the parameters of an adaptive system should only be disturbed in a minimal fashion in the light of new data [34, p. 355]. Using this principle, the one-dimensional LMS algorithm can be shown to correspond not to an approximate solution of a gradient-based optimization problem, but to the exact solution of a constrained optimization problem [35, p. 216].

Differently from simultaneous video SRR methods, the cost function (4) of the R-LMS algorithm is defined for a single time instant. Thus, the proximal regularization described in Section IV only guarantees consistence between consecutive iterations in $k$. As the solution $\hat{\mathbf{x}}(t-1)$ at the previous time instant is only introduced during the initialization in (8), consistence between consecutive time instants is only obtained if the parameters of the R-LMS algorithm are selected such that the solution is not disturbed during all iterations $k=1,\ldots,K$ (i.e. $\hat{\mathbf{x}}_{K}(t)\simeq\hat{\mathbf{x}}_{0}(t)$). However, as illustrated in the example in section IV-A, this makes the R-LMS algorithm very sensitive to outliers. A choice of parameters leading to a superior robustness, on the other hand, compromises the estimation of the details in $\mathbf{x}(t)$.

To preserve the super-resolved details between consecutive time instants regardless of the choice of the R-LMS parameters, one might be tempted to introduce an additional temporal term to the optimization problem (10) to prevent the loss of content estimated in $\hat{\mathbf{x}}(t-1)$ when reconstructing $\hat{\mathbf{x}}(t)$. This results in the following optimization problem:

[TABLE]

where $\alpha_{\tiny{\text{T}}}$ is a weighting factor controlling the temporal disturbance. Albeit removing the dependence of its solution on the time initialization (8), the algorithm in (19) fails to achieve good results. Instead, this new regularization term makes the algorithm less robust since it prevents convergence to the desired solution $\mathbf{x}(t)$ in the presence of large innovations even for a large number of iterations (large $K$). This is clearly perceived by assuming again that $\|\mathbf{s}(t)\|$ is large and $\hat{\mathbf{x}}(t-1)\simeq\mathbf{x}(t-1)$, and examining the norm of the last term in (19) for $\mathbf{z}=\hat{\mathbf{x}}_{k+1}(t)$

[TABLE]

which will be large if $\hat{\mathbf{x}}_{k+1}(t)\simeq\mathbf{x}(t)$ not only for $k=1$, but for all iterations. Furthermore, this term would be unnecessary for small innovations. In this case the R-LMS can retain the temporal consistency even for a large number of iterations ($K$), as illustrated in the example of section IV-A for $K=100$. Hence, algorithm robustness and quality must be addressed using other approaches.

Most works in single-frame or video SRR seek robustness by considering cost functions including non-quadratic (e.g. $\text{L}_{1}$) error norms [25, 4, 7] or signal dependent regularizations [9, 36], which result in non-linear algorithms. Although these techniques achieve good reconstruction results, their increased computational cost makes real-time operation unfeasible even for the fastest algorithms. Differently from the simultaneous SRR methods, the robustness problem of the R-LMS is related with its slow convergence, since a good result is achieved for large $K$. A different approach is therefore required to adequately handle the innovations in the R-LMS algorithm.

In the following, we propose to use meaningful a priori information about the statistical nature of the innovations in deriving a new stochastic SRR method using the least perturbation principle. The proposed approach can improve the robustness of the R-LMS algorithm while retaining a reduced computational cost. By employing statistical information about $\mathbf{s}(t)$, which has been overlooked in the design of simple SRR algorithms, it becomes possible to provide robustness to the innovations while maintaining a good reconstruction quality.

V-A Constructing an Innovation-Robust Regularization

To achieve the desired effect, we propose to modify the norm being minimized in the last term of (19) through the inclusion of a weighting matrix $\mathbf{Q}$ properly designed to emphasize the image details in the regularization term. This will allow the resulting algorithm to attain a faster speed of convergence with a good quality, while at the same time reducing the influence of the innovations on the solution of the optimization problem. The new constraint is then given by

[TABLE]

where $\mathbf{Q}$ must be designed to preserve only the details of the estimated images, so that innovations will have a minimal effect upon the regularization term. Thus, it is desired that

[TABLE]

which means that the image details must lie in the column space of $\mathbf{Q}$, while the innovations lie in its nullspace. Therefore, if we assume the reconstructed image in time instant $t-1$ to be reasonably close to the real (desired) image, (i.e. $\hat{\mathbf{x}}_{K}(t-1)\simeq\mathbf{x}(t-1)$), we can write the modified restriction as

[TABLE]

If $\mathbf{Q}$ satisfies (21), $\|\mathbf{Q}\mathbf{s}(t)\|\approx 0$ even in the presence of an outlier, and (22) can be approximated by

[TABLE]

enabling the preservation of the image details even in the presence of large innovations.

The question that arises is how to design the transformation matrix $\mathbf{Q}$ to achieve the required properties. We propose to base the design of $\mathbf{Q}$ on a stochastic model for the innovations.

V-B Statistical Properties of Innovations in Natural Image Sequences

The statistical properties of natural images have been thoroughly studied in the literature. A largely employed probabilistic model for natural images is characterized by a zero-mean and highly leptokurtotic, fat-tailed distribution, with its power spectral density remarkably close to a $1/f_{r}^{\rho}$ function, where $f_{r}$ is the absolute spatial frequency and $\rho$ is close to $2$ [37]. This characterization has led, for example, to the development of sparse derivative prior models for natural images [38] that have been widely employed in image processing algorithms.

When it comes to obtaining accurate probabilistic models for the signals in the dynamic evolution of a video sequence, particularly the innovations, the task becomes more challenging. This is due to the dependence of the signal statistics on the generally unknown movement in the scene. With the motion frequently estimated from a low-resolution observed video sequence, the employed model must distinguish between errors originating from the image registration and errors caused by true changes in the scene, the latter often labeled as outliers.

The modeling of large magnitude changes in a scene has already been considered for the image matching problem. Hasler et. al. [39] proposed to consider the error patterns generated by non-coinciding regions of an aligned image pair to be similar to the error generated by comparing two random regions of the underlying scene. This relationship clearly arises in a dynamical model for a video sequence when the motion model fails to account for unpredictable changes between two adjacent images, generating an error signal that will consist of the difference between the new image and a misaligned part of the previous image. Considering the case of one dimensional signals for simplicity, the auto-correlation function of the difference between two patches of an image separated by $\Delta$ samples can be computed as

[TABLE]

where $\operatorname{E}[\,\cdot\,]$ denotes statistical expectation, $I(p)$ is a point in the one dimensional image, and $\mathbf{r}_{I}(l)$ is the image auto-correlation function. Thus, $\mathbf{r}_{\Delta}(l)$ is the auto-correlation of the simulated outlier. If the covariance between the image pixels diminishes with their distance, for a sufficiently large value of $\Delta$ the terms $\mathbf{r}_{I}(l\pm\Delta)$ will become approximately equal to the square of the mean image value. Therefore, the auto-correlation function of the simulated outlier will be similar to that of a natural image.

This interpretation can be more intuitively achieved by considering a different approach, which models the innovations assuming a scene model composed by the interactions of objects in an occlusive environment [40]. Innovations in a video sequence can be broadly described as pixels in $\mathbf{x}(t)$ that cannot be described as a linear combination of the pixels in $\mathbf{x}(t-1)$ (i.e. are statistically orthogonal). These pixels will be here divided as

[TABLE]

where $\boldsymbol{\eta}(t)$ consists of small changes in the scene, originated, for instance from specular surfaces. $\boldsymbol{\eta}(t)$ can be modeled as a low power high frequency noise. $\mathbf{d}(t)$ represents large magnitude changes (outliers) arising due to occlusions or due to objects suddenly appearing in the scene (such as image borders). $\mathbf{d}(t)$ is usually sparse and compact [41] 222Note that $\mathbf{s}(t)$ is not to be confused with registration errors due to the ill-posed nature of the motion estimation process. The latter can be shown to originate from a random linear combination of the pixels in $\mathbf{x}(t-1)$ [27]..

A region of the scene corresponding to a dis-occluded area typically reveals part of a background or object at a different depth from the camera. Hence, the nonzero pixels in $\mathbf{d}(t)$ will consist of highly correlated compact regions. Furthermore, the joint pixel statistics at these locations should actually be similar to that of natural scenes. This conclusion becomes straightforward if we consider, for instance, the Dead Leaves image formation model [40], which characterizes a natural scene by a superposition of opaque objects of random sizes and positions occluding each other. Here, a dis-occluded area would correspond to the removal of an object (or a “leaf") at random from the topmost of the z-axis. The corresponding region in the new image will therefore be composed of the next objects present on the z-axis. Since the area behind the view plane is completely filled with objects (superimposed "leaves") in this model, there is no difference between the statistical properties of a region in the foremost-top image and those of a region behind an object. This reinforces the notion of correlation obtained by considering the more generic outlier model of Hasler et al. [39].

To verify the proposed innovations model, we have determined the power spectral density (PSD) of synthetic images representing the innovations. These images were generated by pasting small pieces of the difference between two independent natural images with sizes ranging from $5\times 5$ to $15\times 15$ in random positions of a $64\times 64$ background. We have extracted the small pieces from 20 different natural images, so that they emulate small regions appearing in the occluded regions of a video sequence. The PSD is computed by averaging 200 realizations of a Monte Carlo (MC) simulation. Figure 3 shows the obtained result. It can be clearly seen that the energy is concentrated in the lower frequencies of the spectrum, resulting in a highly correlated signal.

V-B1 Choosing the Operator $\mathbf{Q}$

Natural scene innovations tend to be highly correlated in space. Thus, their energy tends to be primarily concentrated in the low spatial frequencies. Hence, the operator $\mathbf{Q}$ should in general emphasize the high frequency components to accomplish the design objectives in (21). Unfortunately, the specific scenes to be processed are not known in advance, what hinders the accurate determination of the statistical properties of the innovations, and thus of the optimal operator $\mathbf{Q}$. A simple solution with reduced computational complexity is to use a basic high-pass filter with small support, such as a differentiator or a Laplacian. For simplicity, the Laplacian filter mask will be employed during the remaining of this work. Thus, we shall use $\mathbf{Q}=\mathbf{S}$, leaving the search for an optimal operator for a future work.

VI The Temporally Selective Regularized LMS (TSR-LMS) Algorithm

To derive the new algorithm, we propose a new cost function that minimizes the perturbation only on the details of the reconstructed image, while at the same time observing the objectives of the R-LMS algorithm. Differently from (19), the new cost function allows for more flexibility for the component of the solution in the subspace corresponding to the outlier while retaining its quality. Such strategy leads to an increased algorithm robustness.

We propose to solve the following optimization problem:

[TABLE]

Calculating the gradient of the cost function with respect to $\mathbf{z}$, setting it equal to 0, solving for $\mathbf{z}=\hat{\mathbf{x}}_{k+1}(t)$ and approximating the statistical expectations by their instantaneous values yields the iterative equation for the TSR-LMS algorithm:

[TABLE]

where the time update is based on the signal dynamics (2) and performed by $\hat{\mathbf{x}}_{0}(t+1)=\mathbf{G}(t+1)\hat{\mathbf{x}}_{K}(t)$ [23].

Algorithm (26) generalizes R-LMS and the least perturbation approach (19). It collapses to these solutions if $\alpha_{\tiny{\text{T}}}\rightarrow\infty$ or $\mathbf{Q}={\mathbf{I}}$, respectively. It should perform well both with and without outliers, at the cost of little extra computational effort. Though the matrix inversion can be made a priori and the resulting inverse might be sparse, its storage is still rather costly. If $\mathbf{Q}$ is chosen to be block circulant (BC) (such as a Laplacian), then \big{(}{\mathbf{I}}+\frac{1}{\alpha_{\tiny{\text{T}}}}\mathbf{Q}^{\text{\tiny T}}\mathbf{Q}\big{)}^{-1} is block circulant [42] and can be computed as a convolution, leading to important memory savings.

Although (26) may resemble the Gradient Projection Method (GPM) [43], this is not generally true, as \mathbf{M}=\big{(}{\mathbf{I}}+\frac{1}{\alpha_{\tiny{\text{T}}}}\mathbf{Q}^{\text{\tiny T}}\mathbf{Q}\big{)}^{-1} is not necessarily a projection matrix (i.e. $\mathbf{M}^{2}\neq\mathbf{M}$). Hence, the convergence and stability properties of these algorithms are not the same in general.

VI-A The Linearized Temporally Selective Regularized LMS (LTSR-LMS)

Whereas algorithm (26) should present a good cost-benefit ratio, the aforementioned limitations motivates the pursuit of another algorithm that trades a small performance loss for a simpler implementation and a more predictable performance. This section describes one possible modification.

Since the details of the solution are minimally disturbed between iterations, we can safely assume that $\mathbf{Q}\hat{\mathbf{x}}_{k+1}(t)\approx\mathbf{Q}\hat{\mathbf{x}}_{k}(t)$. Therefore, we can employ a linear approximation for the quadratic regularization introduced in the last term of (25) using a first-order Taylor series expansion of this norm with respect to the transformed variable $\mathbf{Q}\mathbf{z}$ about the point $\mathbf{Q}\mathbf{z}=\mathbf{Q}\hat{\mathbf{x}}_{k}(t)$. The resulting cost function can be written as:

[TABLE]

Note that if the algorithm initialization is selected as $\hat{\mathbf{x}}_{0}(t)=\mathbf{G}(t)\hat{\mathbf{x}}_{K}(t-1)$ [23], the linearized regularization introduced in the last term of (27) is equal to zero for the first iteration ($k=1$). Therefore, $K\geq 2$ iterations per time instant are necessary in order to have an improvement over the R-LMS algorithm. This is not the case for the algorithm proposed in (26), where an improvement can be obtained even for $K=1$.

By ignoring the constant terms in the optimization problem (27) and using (4) and (5), it can be shown that (27) corresponds to the proximal point cost function of the following Lagrangian for a single time time $t$:

[TABLE]

where $\alpha_{\tiny{\text{T}}}=1/(2\tilde{\alpha_{\tiny{\text{T}}}}\mu)$. Note that by using $\mathbf{Q}={\mathbf{I}}$ on (VI-A), the algorithm particularizes to the well known temporal regularization, commonly employed in simultaneous video SRR in order to achieve temporal consistency [10, 7]. In this case, the algorithm is not expected to be robust since the innovations are not accounted for.

Calculating the gradient of the cost function in (27) with respect to $\mathbf{z}$, setting it equal to 0, solving for $\mathbf{z}=\hat{\mathbf{x}}_{k+1}(t)$ and approximating the statistical expectations by their instantaneous values yields the iterative equation for the LTSR-LMS algorithm based on the linearized version of the proposed regularization:

[TABLE]

where the update is performed for a fixed $t$ and for $k=1,$ $\ldots,K$. Like the traditional R-LMS, the time update of (VI-A) is based on the signal dynamics (2), and performed by $\hat{\mathbf{x}}_{0}(t+1)=\mathbf{G}(t+1)\hat{\mathbf{x}}_{K}(t)$ [23].

VI-B Computational Cost of the Proposed Solution

The computational and memory costs of the (L)TSR-LMS algorithms are still comparable to those of the (R)-LMS algorithm. An important characteristic of the problem that allows a fast implementation of both the (R)-LMS and the (L)TSR-LMS methods is the spatial invariance assumption of the operators $\mathbf{M}=({\mathbf{I}}+\frac{1}{\alpha_{\text{\tiny T}}}\mathbf{Q}^{\text{\tiny T}}\mathbf{Q})^{-1}$, $\mathbf{H}$, $\mathbf{S}$ and $\mathbf{Q}$, which results in them being block-circulant or block-Toeplitz matrices. In this case, the number of nonzero elements of the matrices (denoted by $|\,\cdot\,|$) scales linearly with the number of HR image pixels (i.e. $\propto M^{2}$), and so does the number of operations of the algorithms. In this case, the computational and memory costs for the algorithms considered can be seen in Tables I and II.

VII Results

We now present four examples to illustrate the performace of the TSR-LMS and LTSR-LMS methods. Examples 1 and 2 use a controlled environment to assess differences in quality (section VII-A) and robustness (section VII-B) when compared to interpolation, LMS and R-LMS without the influence of unaccountable effects. Example 3 (section VII-C) compares the performances of the proposed methods to those of LMS, R-LMS and interpolation algorithms using real video sequences. Example 4 (section VII-D) evaluates the TSR-LMS and LTSR-LMS methods against state-of-the-art SRR algorithms in practical applications.

Example 1 evaluates the average performance of the algorithms without outliers, in a close-to-ideal environment. We used synthetic video sequences with small translational motion to enable Monte Carlo simulations and to be able to control the occurrence of modeling errors. The motion between frames was assumed known a priori, and the mean squared reconstruction error could be evaluated because we had access to the desired HR images. The simulation was also performed using a typical registration algorithm to evaluate the influence of motion estimation on the performances of algorithms (26) and (VI-A). A decline in performance was expected, as reported in [10] for the classical temporal regularization algorithms.

Example 2 evaluates the proposed algorithms in the presence of innovation outliers. A synthetic simulation emulates the case of a flying bird when an object suddenly appears in a frame or moves independently of the background, generating occlusions and leading to a high level of innovations in some specific frames of the video sequence.

Example 3 evaluates the performances of the algorithms when super-resolving real video sequences. We super-resolved a set of video sequences containing complex motion patterns and frames with large levels of innovations and registration errors. Finally, Example 4 extends Example 3 by comparing the TSR-LMS and LTSR-LMS methods with state-of-the-art SRR algorithms, namely, a Bayesian method [15] and a Convolutional Neural Network [18]. We employ recent and challenging video sequences and compare the results for robustness, quality and computational cost.

For $\mathbf{Q}={\mathbf{I}}$, the LTSR-LMS algorithm particularizes to the popular classical temporal regularization [10, 7]. We do not report this case here, as we could not obtain any improvement (quantitative or perceptual) when compared to the R-LMS ($\mathbf{Q}=0$) performance. The matrix $\mathbf{Q}$ employed in both algorithms for all examples shown here was a Laplacian filter. The SRR algorithms were compared to the bicubic interpolation algorithm proposed in [30]. Both the MSE and the structural similarity index (SSIM) were considered in the evaluation. The obtained SSIM results were qualitatively similar to the MSE results. Hence, the SSIM values are provided only for the displayed images.

The boundary condition for the convolution matrices was chosen to be circulant. This simplifies the implementation and results in the inversion of a circulant matrix in (26) [42]. We selected the boundary condition for $\mathbf{G}(t)$ in the global translational case to be circulant as well to simplify implementation. For the case of a dense motion field, the warped images were computed through the bilinear interpolation of the original image pixels.

VII-A Example 1

For a Monte Carlo simulation, each HR video sequence was created based on the translation of an $256\times 256$ window over a static image, resulting in whole-image translational movements. The window displacements consisted in a random walk process (i.i.d. unitary steps) on both horizontal and vertical directions. The still images used to generate each video sequence were natural scenes such as Lena, Cameraman, Baboon and others, and were totally distinct from each other. The resulting sequence was then blurred with a uniform unitary gain $3\times 3$ mask and decimated by a factor of $2$, resulting in LR images of dimension $N=128$. Finally, white Gaussian noise with variance $\sigma^{2}=10$ was added to the decimated images. The algorithm performances were evaluated by averaging 50 realizations, one for each input image.

The performances of the different algorithms are highly dependent on the parameters selected, as verified in section IV-A). Hence, the parameters were carefully selected to yield an honest comparison. We selected the parameters for each method to achieve the minimum steady-state MSE (i.e. for large $t$). The steady-state MSE for each set of parameters was estimated by running an exhaustive search over a small, independent set of images and averaging the squared errors for the last 5 frames. Table III shows the parameter values that resulted in the best performance for each algorithm.

We applied both standard and regularized versions of LMS and the proposed methods to super resolve the synthetic sequences, all initialized with $\hat{\mathbf{x}}(1)$ as a bicubic (spline) interpolation of the first LR image. First, we considered the motion vectors to be known a priori to avoid the influence of a registration algorithm, and used $K=2$ iterations per time instant. The super-resolved sequences were compared to the original HR sequence and the MSE was computed across all realizations.

The MSE performance is depicted in Figure 4. It can be seen that the proposed methods outperform LMS and R-LMS. Moreover, both proposed algorithms achieve the same MSE given enough time. Finally, the LTSR-LMS algorithm reaches the minimum MSE faster than the original TSR-LMS algorithm in the absence of outliers.

For a more realistic evaluation, we repeated the MC simulation considering the influence of registration errors. The Horn & Schunck registration algorithm [44, 45] was employed333The parameters were set as: lambda=1$\times$103, pyramid_levels=4, pyramid_spacing=2., with the velocity fields averaged across the entire image to compute the global displacements. The algorithm parameters were the same used in the previous simulation. The resulting MSE results are depicted in Figure 5, and an example of a reconstructed image of a resolution test chart is shown in Figure 6.

It can be verified that the proposed methods still outperform the conventional (R-)LMS algorithms, though by a smaller margin due to the effect of registration errors. It should be noted that large levels of registration errors tend to reduce the effectiveness of the TSR-LMS and LTSR-LMS methods, as we are basically preserving information (details) from previous frames that must be registered. Furthermore, the performance of the TSR-LMS algorithm showed greater sensitivity to unknown registration, as its performance degraded more when compared to the LTSR-LMS algorithm. The images reconstructed by the four evaluated algorithms were found to be perceptually similar, although a careful inspection reveals a slight improvement in the reconstruction result using the LTSR-LMS algorithm. Nevertheless, the following examples will illustrate that the proposed methods perform considerably better than the others in the presence of outliers.

VII-B Example 2

This example evaluates the robustness to innovation outliers. This was done by super-resolving synthetic video sequences containing a suddenly appearing object, which is independent of the background. The first MC simulation of Example 1 was repeated, this time with the inclusion of an $N\times N$ black square appearing in the middle of the 32nd frame and disappearing in the 35th frame of every sequence, thus emulating the behavior of a flying bird outlier on the scene.

The MSE evolution for the parameters shown in Table III is shown in Figures 7-(a) and 7-(c) (zoomed). The improvement provided by the LTSR-LMS algorithm is clearly significant, suggesting its greater robustness to outliers.

To improve the design of the proposed algorithms, we performed exhaustive searches in the parameter space of all algorithms to determine good sets of parameters for reconstructing a small independent set of images. The parameters in Table IV yielded the minimum MSE averaged between frames $30$ and $40$ for each algorithm. The MSE evolutions for these parameters are shown in Figures 7-(b) and 7-(d).

The proposed methods led to a significant performance gain compared to the other algorithms when in the presence of outliers in frames 32 and 35, with a slightly better results verified for the TSR-LMS method.

While R-LMS led to a MSE similar to that achieved by the TSR-LMS and LTSR-LMS algorithms for frames 33 and 34 (when the black square remained in the scene), its performance degraded considerably for larger $t$. Note also that the loss in steady-state performance as a result of the optimization to handle outliers was higher for LMS and R-LMS. Finally, the LTSR-LMS algorithm showed to be less sensitive to the parameter selection, performing reasonably well in both simulations. Figure 8 shows the reconstructed images for frame 32, which confirm the quantitative results. The black square introduced in the sequence is significantly better represented for the proposed methods (when it is indeed present in the HR image), and a slight improvement can be noticed in the result of the TSR-LMS algorithm when compared to that of the LTSR-LMS.

VII-C Example 3

This example illustrates the effectiveness of the proposed methods when super-resolving real video sequences. We performed a Monte Carlo simulation with 15 realizations consisting of natural video sequences, like Foreman, Harbour, News and others. In this case, the true motions of the objects and camera are unknown. We used Horn & Shunck algorithm with the same parameters shown in Table IV to estimate the dense velocity field, but now considering the displacement to be unique for each image pixel.

For a quantitative evaluation, the original videos were used as available HR image sequences. For simplicity, only the $256\times 256$ upper-right region of the original sequence was considered so that the resulting images were square. Like in Example 1, the HR sequences were blurred with an uniform unitary gain $3\times 3$ mask, decimated by a factor of $2$ and contaminated with white Gaussian noise with variance $\sigma^{2}=10$ to form the LR images. The standard LMS versions and the proposed methods were used to super-resolve the LR sequences with $K=2$ iterations per time sample and the parameters set at the values in Table IV. Hence, they were not guaranteed to be optimal, as the amount of innovations is different and the motion is not known in advance.

Figure 9 shows the MSE evolutions, which indicate a better performance for the TSR-LMS and LTSR-LMS methods, both of which performed similarly. The improvement offered by the proposed methods can be most clearly observed when a high degree of innovations is present in the scene. It can be noted that the reconstruction error exhibits a more regular behavior across the entire sequence, without being considerably influenced by the outliers, which cause significant spikes in the MSE of the LMS and R-LMS algorithms. To illustrate this scenario, Figure 10 shows the 93rd super-resolved frame of the Foreman sequence. The advantage of the proposed methods becomes apparent through a more clear reconstruction result, as opposed to a vast amount of artifacts found in the images reconstructed by the LMS and R-LMS algorithms, mainly in the regions where innovations are present. The MSE performances of the four methods are less discernible at the time intervals where the amount of innovations is less significant. Nevertheless, the difference in the perceptual quality of the reconstructed images is still noticeable. For instance, Figure 11 shows the reconstruction results of part of the 33rd frame of the Foreman sequence. Although for this frame the MSE differences between the results of the four algorithms are small, the images super-resolved by the proposed methods still offer a better perceptual quality with reduced artifacts where small localized motion occur (e.g. close to the man’s mouth).

VII-D Example 4

This example illustrates the performance of the algorithms using recent and challenging video sequences taken from the DAVIS dataset [46]. We extracted six reference HR sequences from videos with resolution of $1080\times 1920$ pixels, and generated the degraded LR sequences from them as in Example 3.

To compare performances with recent state-of-the-art algorithms, we super-resolved these sequences using interpolation, LMS, R-LMS and two more recent video SRR algorithms, namely, the adaptive Bayesian method from [15] and the Convolutional Neural Network (CNN) from [18]. For the method of [15] we fixed the solution to the blur matrix estimation of $\mathbf{H}$ at the optimal value and used the same registration algorithm of [45], which was also used for the other methods. The CNN has an embedded registration algorithm which could not be modified. The parameters used for LMS, R-LMS, TSR-LMS and LTSR-LMS methods were the same as in Table IV (which were determined based on the simulations with synthetic sequences in Example 2). The CNN [18] was implemented in Python using TensorFlow. The other methods were implemented in Matlab©. Codes for methods from [15] and [18] were provided by the respective authors. All simulations were executed on a desktop computer with an Intel Core I7 processor with 4.2Ghz and 16Gb of RAM.

We assess the performances of the different methods both quantitatively and visually. The peak signal to noise ratio (PSNR) for all algorithms and video sequences is presented in Table VI. It can be verified that the proposed methods usually led to an image quality that is at least comparable with (usually better than) that resulting from using the competing algorithms. Excerpts from reconstruction results of two sequences444More extensive results are available in the supplementary document. Illustrative video excerpts are also included as supplementary material on https://www.dropbox.com/sh/hmyp5nks2sbj6to/AADnSqRE5YpTDJheLzDlTUI3a?dl=0., shown in Figs. 12 and 13, also support this conclusion. Fig. 12 shows that the proposed methods yield a significantly better resolution than the algorithms in [15] and [18]. Fig. 13 is an example in which the competing methods performed well, specially the CNN [18] as can be noticed in the brick wall. Still, the proposed methods provided good reconstruction quality and show less influence of noise.

The methods of [15] and [18] also show robustness to innovation outliers, as can be noticed in the horse’s tail in Fig. 13 where the reconstructed images do not exhibit artifacts like the LMS and R-LMS algorithms. The proposed methods also show significantly less artifacts than the LMS and R-LMS algorithms, albeit not being as clear as the methods in [15] and [18].

Table V shows the execution times555The execution time of the CNN [18] includes image registration, since the SRR time cannot be measured separately in the provided implementation. for Example 4. Algorithms [15] and [18] were approximately two orders of magnitude slower than the remaining methods (except for bicubic interpolation). These results and the comparable reconstruction quality clearly show that the proposed algorithms are competitive in term of quality and robustness at a much lower computational cost.

The implementations of the proposed algorithms were not optimized, and thus they could not be tested in real-time. However, real-time implementation is perfectly within reach using existing devices. For instance, the LTSR-LMS algorithm needs approximately $0.137$ billion floating point operations (GFLOPS) to process each frame in the current example. Consider using a fast image registration algorithm such as [47], which has a computational complexity of $\kappa^{2}M^{2}+g_{\max}^{2}M^{2}+M^{2}$, with $\kappa$ being a small image window and $g_{\max}$ the maximum displacement amplitude. For typical values $\kappa=3$ and $g_{\max}=10$, the image registration cost is $0.228$ GFLOPS per frame. Now, for real-time performance (at 30 frames per second), the aggregate cost of image registration and SRR becomes $10.95$ GFLOPS/second, which is well within the capability of graphical processing units released almost a decade ago [48].

VIII Conclusions

This work proposed a new super-resolution reconstruction method aimed at an increased robustness to innovation outliers in real-time operation. An intuitive interpretation was proposed for the proximal-point cost function representation of the R-LMS gradient descent algorithm. A new regularization was then devised using statistical information on the innovations. This new regularization allowed for faster convergence of the solution in the subspace related to the innovations, while preserving previously estimated details. Two new algorithms were derived which present an increased robustness to outliers when compared to the R-LMS, with only a modest increase in the resulting computational cost. Computer simulations both with synthetic and real video sequences illustrated the effectiveness of the proposed methods.