DeepSUM: Deep neural network for Super-resolution of Unregistered Multitemporal images

TL;DR

DeepSUM introduces a CNN framework that jointly performs super-resolution and registration of multitemporal unregistered images, leveraging spatial and temporal correlations for enhanced image quality.

Contribution

It presents a novel end-to-end CNN architecture that integrates registration within super-resolution, specifically designed for unregistered multitemporal remote sensing images.

Findings

Won the PROBA-V super-resolution challenge

Effectively combines registration and super-resolution in a single network

Improves super-resolution quality for unregistered multitemporal images

Abstract

Recently, convolutional neural networks (CNN) have been successfully applied to many remote sensing problems. However, deep learning techniques for multi-image super-resolution from multitemporal unregistered imagery have received little attention so far. This work proposes a novel CNN-based technique that exploits both spatial and temporal correlations to combine multiple images. This novel framework integrates the spatial registration task directly inside the CNN, and allows to exploit the representation learning capabilities of the network to enhance registration accuracy. The entire super-resolution process relies on a single CNN with three main stages: shared 2D convolutions to extract high-dimensional features from the input images; a subnetwork proposing registration filters derived from the high-dimensional feature representations; 3D convolutions for slow fusion of the features from multiple images. The whole network can be trained end-to-end to recover a single high resolution image from multiple unregistered low resolution images. The method presented in this paper is the winner of the PROBA-V super-resolution challenge issued by the European Space Agency.

Tables

Table 1. TABLE I: Average mPSNR (dB) and SSIM - RegNet Performance

	Proposed without RegNet	Proposed with RegNet
NIR	47.68 / 0.98519	47.84 / 0.98578
RED	49.87 / 0.99038	50.00 / 0.99075

Table 2. TABLE II: Average mPSNR (db) and SSIM

	Bicubic	Bicubic+Mean	IBP [36]	BTV [39]	SISR	SISR+Mean	DUF [58]	DeepSUM
NIR	45.05/0.97654	45.69/0.97782	45.96/0.97960	45.93/0.97942	45.56/0.97938	46.41/0.98166	47.06/0.98417	47.84/0.98578
RED	47.61/0.98474	47.91/0.98507	48.21/0.98648	48.12/0.98606	48.20/0.98704	48.71/0.98787	49.36/0.98948	50.00/0.99075

Equations

$$I^{\text{SR}}=f(I_{[0,N-1]}^{\text{LR}},\theta ),$$

$$G_i$$

$$Z_i^{\text{IRLR}}$$

$$I_i^{\text{IRLR}}$$

$${\bar{I}}^{\text{IRLR}}$$

$$I^{\text{SR}}$$

$$L=\underset{u,v\in [0,2d]}{\min }\Vert I_{u,v}^{\text{HR}}-(I_{\text{crop}}^{\text{SR}}+b){\Vert }^2,$$

$$b=\frac{1}{(rW-d)(rH-d)}\sum _{x,y}\left(I_{u,v}^{\text{HR}}-I_{\text{crop}}^{\text{SR}}\right).$$

$$\text{mPSNR}=\underset{u,v\in [0,6]}{\max }20\log \frac{2^{16}-1}{\parallel I_{u,v}^{\text{HR}}-(I_{\text{crop}}^{\text{SR}}+b){\parallel }^2}$$

Full text

I Introduction

Super-resolution (SR) techniques reconstruct a high-resolution (HR) image from one or more low-resolution (LR) images. Remote sensing is playing a key role in mapping and monitoring the Earth and increasing the availability of high spatial resolution data is crucial for many applications such as urban mapping, military surveillance, intelligence gathering, disaster and vegetation growth monitoring. The ever increasing spatial and spectral resolution of instruments onboard of satellites generates large amounts of data which challenge compression algorithms [1, 2] to meet the available downlink bandwidth. This often results in reduced availability of HR products. Combining this issue with the high cost of hardware for smaller missions, it is clear that developing a new generation of post-processing techniques to enhance the spatial resolution is a critical objective.

The approaches to image super-resolution can be broadly framed into two main categories: single-image SR (SISR) and multi-image SR (MISR). SISR exploits spatial correlation in a single image to recover the HR version. However, the amount of information available in a single image is quite limited as some information has inevitably been lost in the LR image formation process. Certain applications provide multiple LR versions of the same scene to be combined by means of MISR techniques, where the reconstruction of high spatial-frequency details takes full advantage of the complementary information coming from different observations of the same scene.

For remote sensing problems, multiple images of the same scene can typically be acquired by a spacecraft during multiple orbits, by multiple satellites imaging the same scene at different times, or may be obtained at the same time with different sensors. In this context, developing a successful MISR model hinges on solving important problems such as image registration, invariance to absolute brightness variability, time-varying scene content (e.g., due to the time elapsed between multiple acquisitions), and unreliable data (e.g., due to cloud coverage). Deep learning methods have been proved highly successful in the SISR problem but little work has been done for the MISR problem with remote sensing data.

In this paper we present a deep learning architecture addressing MISR applied to a novel dataset provided by the European Space Agency’s Advanced Concepts Team in the context of a challenge [3]. The goal of the challenge is to super-resolve images from the PROBA-V satellite. The method presented in this paper won the challenge by achieving the highest fidelity on the reconstructed images. The unique feature of this dataset is that both LR and HR images have been acquired by the same spacecraft, as opposed to previous works where LR images are artificially down-scaled, degraded and shifted versions of an HR image. The images are not simultaneously acquired so temporal variations exist and have to be handled as well in the super-resolution process.

Our main contribution is DeepSUM, a novel CNN-based architecture to combine multiple unregistered images from the same scene exploiting both spatial and temporal correlations. Our method includes image registration inside the CNN architecture, as a subnetwork named RegNet, which dynamically computes custom filters and applies them to higher dimensional image representations. This is in contrast with the vast majority of deep-learning MISR methods in literature [4] that compensate for the motion as a preprocessing step. This approach allows the registration task to leverage the feature learning capabilities of the network in order to be more accurate and resilient to scene variations, and it also optimizes it in an end-to-end fashion for the final goal of reconstructing a single HR image. The proposed method is blind to the image degradation model as it does not require to explicitly model the blur kernel or the noise statistics, and it is robust to temporal variations in the scene as well as occlusions due to cloud coverage. The only assumption of our model is the translational nature of the shift among LR images.

A preliminary version of this work [5] addressed MISR for the PROBA-V dataset. With respect to our previous work, in this paper we add the registration as part of the network, we improve the use of the image mask information, we expand the experimental results and comparisons with alternative methods and we discuss how a variable number of LR images can be used while the network is designed to handle a fixed input size.

The remainder of this paper is organized as follows. Section II introduces related works on SISR and MISR. Section III provides details on the novel PROBA-V dataset. Sections IV and V detail the proposed framework and the training procedure. Section VI contains results and performance evaluation. Section VII draws some conclusions.

II Related work

The literature on SR techniques is extensive, both for SISR and for MISR techniques. SISR approaches can be classified into three main classes: interpolation-based methods (e.g., Lanczos kernels), optimization-based methods and learning-based methods. Optimization-based methods explicitly model prior knowledge about natural images to regularize this ill-posed inverse problem, and include low total-variation priors [6], gradient-profile prior [7, 8] and non-local similarity [9, 10, 11]. Adding prior knowledge restricts the possible solution space generating higher quality solutions. However, the performance of many optimization-based methods degrades rapidly when the upscaling factor increases, and these methods are usually computationally expensive.

Learning-based methods can be pixel-based or example-based. The latter ones are the most popular and they model the correspondence among LR and HR patches for HR patch prediction. After the early work by Freeman et al. [12] based on searching $k$-nearest neighbors LR-HR patch pairs of the input LR patch to estimate the HR patch, neighbor embedding [13, 14, 15], sparse-coding [16, 17, 18, 19, 20], anchored neighborhood regression [21], and random forest [22] methods were proposed. More recently deep convolutional neural networks (CNNs) [23, 24, 25, 26, 27, 28, 29, 30] achieved state-of-the-art results for the SISR task. The deep learning paradigm gained attention due to its natural capability of extracting high-level features from images. This is particularly important in remote sensing scenarios where images are highly detailed and their statistics can be very complex.

While most of the deep learning SISR works are related to traditional natural images, lately CNNs have been exploited for remote sensing imagery. A deep learning based method has been applied by Ma et al. [31] on remote sensing images in the frequency domain. Their CNN takes as input discrete wavelet transformed images and adopts recursive block and residual learning in global and local manners to reconstruct HR wavelet coefficients. Jiang et al. [32] proposed a generative adversarial network-based edge-enhancement network for robust satellite image SR reconstruction.

Concerning MISR, the first work was proposed by Tsai and Huang [33], who used a frequency-domain technique to combine multiple under-sampled images with sub-pixel displacements to improve the spatial resolution of Landsat TM acquisitions. Due to the drawbacks of the frequency-domain algorithms, like the difficulty to incorporate the prior information about HR images, many spatial-domain MISR techniques were proposed over the years [34]. Typical spatial-domain methods include non-uniform interpolation [35], iterative back-projection (IBP) [36], projection onto convex sets (POCS) [37, 38], regularized methods [39, 40, 41], and sparse coding [42, 43].

The iterative back projection (IBP) approach was introduced by Irani and Peleg [36]. IBP aims to improve an initial guess of the super-resolved image by back projecting the difference between simulated LR images and actual LR images to the SR image. The updates are iteratively performed attempting to invert the forward imaging process. Drawbacks come from the inability to deal with unknown or very difficult to model image degradation processes, as well as the difficulty in including image priors.

Another class of MISR methods is the projection onto convex sets (POCS) [37, 38], where restoration and interpolation problems are simultaneously solved to estimate the SR image, after accurate motion compensation. Despite allowing an easy incorporation of a priori knowledge, POCS suffers from slow convergence and high computational cost.

Regularized methods are some of the most effective multi-frame SR reconstruction approaches. In the past decades, many kinds of regularizers have been proposed to preserve edge information while removing image noise, such as Tikhonov regularizer [44, 45], Markov random field regularizer [46], total variation (TV) [47, 48, 49] and bilateral total variation (BTV) [39]. In particular, a few works have been proposed for remote sensing applications. Shen et al. [41] proposed a maximum-a-posteriori (MAP) SR method with Huber prior for MODIS images captured in different dates. Another multi-temporal SR method was proposed by Li et al.[50] for Landsat-7 PAN images. Instead, other works [51, 49] proposed SR methods for multi-angle remote sensing captures.

In recent years, sparse coding methods based on dictionary learning have successfully been applied to MISR [42, 43].

Most of the above SR methods assume a priori knowledge of the motion model, blur kernel and noise level, where both blur identification and image registration are performed as a preprocessing stage before reconstruction. However, there are many applications where knowing the image degradation process or reliably estimating it can be challenging. For this reason, many studies have been carried out on blind SR image reconstruction [52, 53]. These blind methods usually work in two stages, namely, (1) image registration from LR images, followed by (2) simultaneous estimation of both the HR image and blurring function. In order to reduce the effect of registration errors, some researchers have developed methods to simultaneously estimate the motion parameters and the reconstruction [54, 55]. Hardie et al. [54] presented a MAP framework to jointly estimate image registration parameters and the HR image. Along the same lines, Zhang et al. [55] also integrated the joint estimation of the blurring function. Moreover, Kato et al. [42] recently proposed a sparse coding method where image registration and sparse coding are performed in a unified framework reducing the image registration error.

In the last years, deep learning based methods have been proposed to solve similar MISR problems in context of video super-resolution [56, 57]. Most of these works are composed of two steps: a motion estimation and compensation procedure followed by an upsampling process, heavily relying on the prior motion estimation. Recently, Jo et al. [58] presented a novel end-to-end residual CNN to produce a SR image without explicit motion compensation. A CNN is trained to simultaneously solve motion estimation and HR image reconstruction tasks by producing a set of pixel-dependent filters and a residual correction. A similar idea was developed by Tian et al. [59]. However, little work has been done on deep learning MISR methods in the context of remote sensing, which poses specific challenges. Kawulok et al. [60] propose a MISR method that does not fully exploit the benefit of deep learning, restraining their CNN to solve a SISR problem. The fusion of the upsampled LR images is performed by the median shift-and-add method, generating a SR image that is used as initial guess for a classic regularized method. Their method is not end-to-end trainable in a supervised manner and their CNN is trained against LR images obtained by artificially degrading HR images. Inspired by the recent video super-resolution works, we aim to tackle the MISR problem on satellite images by jointly registering the input LR images and reconstructing the SR image, all within an end-to-end trainable CNN, where the two tasks are optimized jointly.

III The PROBA-V SR dataset

At present, it is difficult to find a dataset collecting both a set of real-world LR observations and the corresponding HR image for the same scene, as captured from the same platform. Many of the works found in the SR literature are based on simulated data, where LR observations for a specific scene are obtained through a degradation and down-sampling process of the HR images by assuming a sensor imaging model. This is a simplified scenario as it either assumes a non-blind problem, i.e., the degradation model can be characterized to some extent, or has the limitation that a too simple degradation model may not accurately match the real one, especially when in presence of temporal variations in the scene content.

The Advanced Concepts Team of the European Space Agency has issued a competition [61] to perform MISR for the images acquired by the PROBA-V satellite. PROBA-V is an Earth observation satellite designed to map land cover and vegetation growth across the entire globe. It was launched in 2013 into a Sun-synchronous orbit at an altitude of 820km. Its payload provides an almost global coverage with 300m LR images and 100m HR images. However, the HR images are acquired with a higher revisit time, roughly one every 5 days, instead of one per day. The dataset gathers satellite data from 74 regions located around the world from the PROBA-V mission. Images are provided as level 2A products composed of radiometrically and geometrically corrected Top-of-Atmosphere reflectance in Plate Carrée projection for the RED and NIR spectral bands. The size of the collected images is $128\times 128$ and $384\times 384$ for the LR and HR data respectively. The images have a single channel with a bit-depth of 14 bits. Each data point consists of one HR image and several LR images (ranging from a minimum of 9 to a maximum of 30) from the same scene. In total, the dataset contains 1160 scenes, 566 are from NIR spectral band and 594 are from RED band. The images of a specific scene are captured at multiple times over a maximum period of 30 days. Weather and changes in the landscape pose a limitation in the similarity of the images. Clouds, cloud shadows, ice, water, missing regions, presence of agricultural activities and, in general, human activity are the main sources of inconsistency across these images, thus posing a major challenge for any image fusion method. Moreover, each image comes with a mask, indicating which pixels in the image can be reliably used for reconstruction (e.g., they are not covered by clouds). The geometric disparity among the images can be considered as translational only. Subpixel shifts in the content of the LR images do occur and are indeed important for the MISR task.

The unique nature of this dataset (with real LR and HR images captured by the same platform at multiple times) makes for an interesting case study for SR techniques. Developing SR products from multiple, more frequent LR images could simultaneously provide enhanced resolution and higher temporal availability and is therefore an interesting application of MISR. Moreover, having real images of the same scene for both the low and high resolutions enables data-driven methods such as CNNs to learn the inversion of possibly complex degradation models and the best feature fusion strategy to handle temporal variations.

IV Proposed method

The proposed method, called DeepSUM, reconstructs a high-resolution image $I^{\text{HR}}$ given a set of $N$ LR images $I_{[0,N-1]}^{\text{LR}}$ representing the same scene:

[TABLE]

where $\theta$ represents the model parameters and $f$ represents the mapping function from LR to HR. $I_{[0,N-1]}^{\text{LR}}$ and $I^{\text{HR}}$ are represented as real-valued tensors with shape $N\times H\times W\times C$ and $1\times rH\times rW\,\times C$ respectively, where $H$ and $W$ are the height and the width of the input LR frames, $C$ is the number of channels and $r$ is the scale factor. While the LR images roughly represent the same scene as the HR image, there are several factors to be considered:

• •

the LR images are not registered with each other;

• •

the LR images and the HR image are not registered;

• •

the brightness of the HR image may be different from that of any LR image;

• •

the scene changes over multiple acquisitions;

• •

LR and HR images may be covered by different clouds and cloud shadows patterns or affected by corrupted pixels.

To tackle this problem we propose to employ a supervised deep learning approach, where a CNN learns the residual between bicubic interpolation and the ground truth. As a preprocessing step, the LR images are bicubically interpolated to the desired size and then fed into a CNN composed of three main building blocks. An overview of the network is shown in Fig. 1.

The first block, called SISRNet, is a feature extractor that can be seen as a SISR network without the output projection to a single channel. Each of the $N$ input images is processed independently by a sequence of 2D convolutional layers. The convolutional filters are shared along the temporal dimension, i.e., all the $N$ interpolated LR (ILR) images go through the same set of filters.

The second network block, called RegNet, aims at estimating a set of filters to register the $N$ higher dimensional image representations produced by the SISRNet block to each other at integer-pixel precision (notice that the network is working at the same spatial resolution as the HR image, so integer shifts correspond to sub-pixel shifts in the LR data). RegNet has been devised to align $N-1$ instances with respect to the first, taken as reference, by operating purely translational shifts. Therefore, the output is a set of $N-1$ 2D filters to be applied spatially to each feature map of the $N-1$ inputs.

Finally, the third block, called FusionNet, merges the registered image representations in the feature space in a “slow” fashion, i.e., by exploiting a sequence of 3D convolutional operations with small kernels. The output is a single super-resolved image.

In the following, we are going to describe each individual block more in detail.

IV-A SISRNet Architecture

The goal of SISRNet is to exploit spatial correlations to improve upon the initial bicubic interpolation. In doing so, the network learns to extract visual features that can be conveniently exploited by the subsequent network blocks. SISRNet has multiple 2D convolutional layers whose weights are shared among the $N$ input images, effectively processing each of them independently. Each convolutional layer is followed by Instance Normalization [62]. Instance normalization is used in place of Batch normalization [63] to make the network training as independent as possible of the contrast and brightness differences among the input images.

IV-B RegNet Architecture

RegNet is composed of two sub-blocks: a CNN, and a global dynamic convolutional layer (GDC). The CNN processes the higher dimensional image representations $Z_{[0,N-1]}^{\text{ILR}}$ generated by SISRNet block and outputs a set of $N-1$ filters $G_{[1,N-1]}$. Each filter $G_{i}$ is subsequently applied in the spatial dimensions to each of the channels of $Z_{i}^{\text{ILR}}$ by the GDC layer by means of a 2D convolution in order to register each feature map of $Z_{i}^{\text{ILR}}$ with respect to the reference one $Z_{0}^{\text{ILR}}$. The filters $G_{[1,N-1]}$ have a fixed support equal to $K\times K$ that upper bounds the maximum possible translational shift correction to $\lfloor K/2\rfloor$. Notice that there is an implicit assumption that all feature maps of an image require the same shift to be registered with the reference, so that the computed filter is shared channel-wise. The registered feature maps $Z_{[0,N-1]}^{\text{IRLR}}$ of the $N$ images are thus obtained as:

[TABLE]

being $*$ the 2D convolution operator. The same filters are also applied to the input ILR images to register them in the residual connection:

[TABLE]

The novelty of this network is twofold: firstly the filters are dynamically computed for each input image, and secondly it makes use of the features to compute the per-image optimal registration instead of performing it in image space, like most of motion estimation algorithms do. This allows to leverage the powerful feature space of the network to boost the registration performance by making it robust to scene variations. In addition, it is fully differentiable so that the whole architecture can be trained end-to-end.

More in detail, the operations performed by RegNet are depicted in Fig. 2. SISRNet outputs a tensor $Z^{\text{ILR}}$ with shape $N\times rH\times rW\,\times F$, where $F$ is the number of features, that is reshaped before being fed to RegNet. The features of the first image $Z_{0}^{\text{ILR}}$ are chosen as a reference and a new tensor of size $2(N-1)\times rH\times rW\,\times F$ is built by replicating the reference $Z_{0}^{\text{ILR}}$ $N-1$ times and interleaving each replica with the other $(N-1)$ image representations $Z_{[1,N-1]}^{\text{ILR}}$. This sequence of paired reference/unregistered features is then processed by convolutional layers to produce the filters. RegNet has a first 3D convolutional layer and a series of shared 2D convolutional layers. The first layer is the key component of registration and it processes the $2(N-1)$ image representations in pairs by using a stride equal to 2 along the temporal dimension and filters of shape $2\times 3\times 3$. This operation allows to correlate the features of each $Z_{i}^{\text{ILR}}$ with respect to the ones of the reference $Z_{0}^{\text{ILR}}$ and compute the shift. Notice that this processing in pairs is necessary to avoid any ordering ambiguity and let the network understand that the output is relative to the reference. After this 3D convolutional layer the output tensor has shape $(N-1)\times rH\times rW\,\times F$.

This tensor passes through a series of 2D convolutional layers with shared weights along the temporal dimension. The last RegNet 2D convolutional layer applies a number of kernels corresponding to the spatial size of the dynamic filters $K\times K$, obtaining a tensor with shape $(N-1)\times rH\times rW\,\times K^{2}$. Each value over the spatial dimensions can be seen as a local estimate of the desired shift based on the local image representation. Since there is a global translational shift by assumption, the values are averaged over the spatial dimensions to obtain a tensor with shape $(N-1)\times 1\times 1\times K^{2}$.

Finally, this tensor is passed through a softmax layer, so that the values over the last dimension ($K^{2}$) add up to 1. The softmax layer promotes a spiked filter with most elements set to zero [64]. The final tensor represents the $(N-1)$ dynamic filters with shape $K\times K$ to be used to register the $(N-1)$ image representations with the GDC operation, as in Fig. 3.

IV-C Mutual Inpainting

The registered and interpolated feature maps $Z_{[0,N-1]}^{\text{IRLR}}$ have regions with unreliable values due to cloud coverage, shadows, corrupted pixels and so on. A per-pixel boolean mask is assumed to be available as side information, with the purpose of mapping pixels that can be reliably used for the fusion task. An example on how to obtain such mask is to run a cloud detection algorithm on the image to segment areas with clouds. This is very important because areas occluded by clouds do not provide any useful information. In order to prevent FusionNet from combining feature maps from multiple images where some have unreliable intensities, we fill the masked areas with values from the feature maps of other images. The regions with missing or unreliable values in each feature map of each image are filled with values taken from the corresponding feature map of other images having reliable values in those regions, if any are available. In the case none of the images has feature maps with reliable values, we keep those unreliable regions as they are. Since after RegNet the masks are not aligned with the corresponding image representations, we shift the masks by an integral shift as close as possible to subpixel shift computed and operated by RegNet. This procedure is performed on both the residual image representations $Z_{[0,N-1]}^{\text{IRLR}}$ and the registered input images $I_{[0,N-1]}^{\text{IRLR}}$ right before averaging them.

IV-D FusionNet Architecture

The $N$ registered outputs $Z_{[0,N-1]}^{\text{IRLR}}$ are progressively fused by the FusionNet subnetwork. FusionNet is composed of $\lfloor N/2\rfloor$ 3D convolutional layers where convolution is performed without any padding in the temporal dimension, so that the temporal depth eventually reduces to 1. This architecture implements a “slow” fusion process in the feature space, which allows the network to learn the best space to decouple image features that are relevant to the fusion from irrelevant variations and to construct the best function to exploit spatio-temporal correlations [57]. Finally, the proposed architecture employs a global input-output residual connection. The network estimates only the high frequency details necessary to correct the bicubically-upsampled input. This is an established technique for image restoration problems using deep learning [23], including SISR. However, with respect to SISR, our proposed network is a many to one mapping, so the residual is actually added element-wise to a basic merge of $I_{[0,N-1]}^{\text{IRLR}}$ in the form of their average. Notice that registration of the input images is performed before averaging by means of the same filters produced by RegNet. Hence, the output is computed as follows:

[TABLE]

being $R$ the residual estimated by the CNN.

IV-E Loss Function

Model parameters are optimized by minimizing a loss function computed as a modified version of the Euclidean distance between the SR image and the HR target. Minimizing the Euclidean distance is optimal in terms of the mean-squared error metric. Some deep learning works on SISR attempted to use an adversarial loss [65]. While this approach produces visually pleasing results, it tends to hallucinate information, resulting in lower MSE scores and less reliable products in the context of remote sensing; hence, the adversarial approach has not been followed in the present work. As we mentioned in Sec. III, since the PROBA-V satellite does not capture LR images and HR images of a specific ground scene simultaneously, there are discrepancies coming from different weather conditions, changes in the landscape and variable absolute brightness due to the large interval between scene acquisitions. The LR images could be quite different from one another and from the corresponding HR image as well. For this reason, we must make the training objective as invariant as possible to such conditions. In particular, in order to build invariance to absolute brightness differences between ${I}^{\text{SR}}$ and ${I}^{\text{HR}}$, the modified loss function equalizes the intensities of the SR and HR images so that the average pixel brightness is the same on both images. Moreover, since ${I}^{\text{SR}}$ and ${I}^{\text{HR}}$ could be shifted, the loss embeds a shift correction. ${I}^{\text{SR}}$ is cropped at the center by $d$ pixels, i.e., as many pixels as the maximum expected shift. Then all possible patches $I_{u,v}^{\text{HR}}$ of size $(rH-d)\times(rW-d)$ for vertical and horizontal shifts $u,v$ are extracted from the target ${I}^{\text{HR}}$. All possible Euclidean distances are computed and the minimum one is taken as loss to optimize. In summary, our loss is as follows:

[TABLE]

where ${I}_{\text{crop}}^{\text{SR}}$ is the cropped version of $I^{\text{SR}}$ and $b$ represents the brightness correction:

[TABLE]

The loss is computed by utilizing only the HR image pixels that are marked as reliable by the mask provided with the dataset and the SR image pixels for which at least one out of $N$ LR images were clear. The reason for this is that a cloud in the HR image can never be predicted from terrain data in the IRLR images, so its pixels should not contribute to the loss function. Viceversa, it is also impossible to predict HR terrain if all the IRLR images have concealed regions.

V Training process

V-A Pre-training

Training the whole network end-to-end from scratch is hard due to several local minima that do not make SISRNet, RegNet and FusionNet work as expected. For example, the gradients computed during training do not sharply discriminate the RegNet task to generate registration filters from the high-resolution feature learning of SISRNet.

In order to solve this issue, it is possible to pretrain each block to handle its specific subtask, and then combine all the blocks to be fine-tuned in an end-to-end fashion.

V-A1 SISRNet pre-training

As mentioned in Sec. IV, SISRNet aims to independently super-resolve each of the $N$ input images, while providing useful higher dimensional image representations. SISRNet is pretrained by setting up a pure SISR problem (i.e., a single input image) where an additional projection layer is added at the end, in order to turn the high-dimensional feature space into a single-channel image. SISRNet with the final projection layer is trained with the same objective function of the final training, where the single image reconstruction is compared with the only HR image available for the scene. The rationale behind this is to make SISRNet exploit spatial correlations as much as possible to generate the best image features for the SISR task. Once the pretraining procedure is completed, the final layer is removed and a dataset of feature maps of the input training images is generated to pretrain RegNet.

V-A2 RegNet pre-training

The purpose of pre-training RegNet is learning to generate registration filters, i.e., filters that shift the feature maps of the $N-1$ input images with respect to the reference input. This operation would be quite challenging to learn if the whole network was trained end-to-end, so its pretraining is crucial for the overall network performance. RegNet is pre-trained by casting registration as a multi-class classification problem. Each dynamic registration filter generated by the network is viewed as a probability distribution over the possible shifts with the objective of estimating the correct shift. The number of classes is $K^{2}$ since the filter size is $K\times K$. In case of an ideal shift of an integer number of pixels, the predicted filter should be a delta function centered at the desired shift.

The input data to be used for the pretraining of RegNet are the feature maps produced by the pretrained SISRNet for the images in the training set. As described in Sec. IV-B, the input to RegNet are $N$ feature maps from images of the same scene. These feature maps are then synthetically shifted with respect to the first one by a random integer amount of pixels. The purpose is to create a balanced dataset where all possible $K^{2}$ classes (shifts) are seen by the network. The desired output is a filter with all zeros except for a one in the position corresponding to the chosen shift. A cross-entropy loss between the softmax output and the true filter is used to learn the RegNet weights.

V-B Final training

The proposed network is finally trained as a whole, end-to-end for the MISR task. FusionNet is trained from scratch while SISRNet and RegNet weights are initialized from the pretraining procedures. The concurrent optimization of all the network blocks allows SISRNet to finetune the image representations to facilitate the RegNet task that in turn finds the best registration to boost the efficiency of FusionNet.

V-C Testing phase

The network architecture presented in the previous sections has been designed to deal with a fixed number $N$ of LR images for a given scene. However, it might happen that more than $N$ images are available and exploiting them could further boost the SR reconstruction performance. Therefore, during testing, one can perform multiple forward passes by using multiple subsets of the available images. Each subset will produce a different SR estimate and, in the end, all SR estimates are averaged. Notice that the estimates should be registered to each other so it is advisable to always use the same LR image as the reference in the network (e.g., one could choose the image with fewer masked pixels). One method to produce useful subsets when more than $N$ LR images are available is to sort them by increasing number of masked pixels and then use a sliding window over $N$ images to compute SR estimates. It must be remarked that the SR estimate quality degrades with increasing number of masked pixels. Also, the estimates are clearly not independent if some images are reused multiple times, but we found consistent gains on our test set, nevertheless.

Defining the optimal function to merge SR estimates or making the network independent of the number of input images could be studied in future research.

VI Experimental results and discussions

In this section we perform an experimental evaluation of DeepSUM, comparing it with several alternative approaches. Code and pretrained models are available online111https://github.com/diegovalsesia/deepsum. We first perform an ablation study to highlight the contribution given by RegNet to the overall network performance. Then, we assess the performance of alternative approaches.

VI-A Experimental setting

In the following experiments, we employ both the NIR and RED band datasets described in Sec. III. We use 396 scenes for training and 170 for testing from the NIR band dataset and 415 for training and 176 for testing from the RED band dataset. Expanding the training set with more scenes should further improve performance as more variability can be captured by our model. Since DeepSUM is devised to work with a fixed size temporal dimension, we train the network using the minimum number of images available for each scene, i.e., $N=9$ images. When more images are available we select the 9 clearest images according to the masks. As a preprocessing step, all LR images are clipped to $2^{14}-1$ since corrupted pixels with large values occur in the LR images throughout the PROBA-V dataset.

After the bicubic interpolation, each scene is a data-cube of size $9\times 384\times 384$, from which we extract a dataset with patches of size $9\times 96\times 96$. 100 random patches are extracted from each scene, resulting in a total of 38400 samples. The patches are extracted considering the available pixel masks: a patch is accepted only if at least 9 scene images are at least 70% clear and the HR image in the same coordinates is at least 85% clear. The amount of unreliable pixels is relaxed to keep as much information as possible from the original images at the cost of training with sub-optimal patches. Separate networks are trained for RED and NIR. The proposed network is trained for around 3000 epochs with a batch size of 8 for both RED and NIR.

The Adam optimization algorithm [66] is employed for training, with momentum parameters $\beta_{1}=0.9$, $\beta_{2}=0.999$, and $\epsilon=10^{-8}$. The learning rate $\lambda$ is initialized to $5\times 10^{-6}$ for the whole network. We employ the Tensorflow framework to train the proposed network on a PC with 64-GB RAM, an Intel Xeon E5-2609 v3 CPU, and an Nvidia 1080Ti GPU. The exact number of network layers is shown in Fig. 1 and the number of filters is 64 everywhere except for the RegNet’s first layer, which has 128 filters. In order to mitigate border effects, we use reflection padding in all 2D convolutions. Each layer in the network is followed by Leaky ReLU non-linearity, except for the last layer. Each layer in SISRnet and FusionNet is followed by an Instance Norm layer. Instance normalization [62] is used in place of Batch normalization layer to make the network training as independent as possible of the contrast and brightness differences among the input images. Finally, since the network produces a residual estimate $R$, we normalize $\bar{I}^{\text{IRLR}}$ and $I^{\text{HR}}$ so that their difference gives a unit variance residual $R$, thus avoiding any scaling to be performed by the last layer of the network and improving convergence speed.

VI-B Quantitative results

The evaluation metric that we consider is a modified version of the PSNR (mPSNR), from which we derived the loss function described in Sec. IV-E.

[TABLE]

The mPSNR computation is meant only for pixels that are not concealed both in the target HR image and in the reconstructed image. Similarly to the loss function during training, this metric has been devised to cope with the high sensitivity of the PSNR to biases in brightness and with the relative translation that the reconstructed image might have with respect to the target HR image. In this case the maximum mPSNR over all possible shifts is considered for evaluation. Note that, by design of the dataset, the maximum shift in the horizontal and vertical directions is equal to 6 pixels.

We remark that this metric was also used to evaluate submissions to the ESA challenge, where the score was computed as a ratio between the mPSNR of the submission and that of the baseline approach, average over all the held-out test set.

VI-B1 Ablation study

First, we want to assess the effectiveness of the sliding window procedure described in Sec.V-C to account for more than 9 images for a given scene. Fig. 4 shows the mPSNR as function of the number of SR estimates used for computing the average. Notice that the mPSNR quickly saturates due to the lower quality of the images in the dataset (e.g., too many masked pixels). Nevertheless, averaging allow to achieve an mPSNR gain up to 0.3 dB over a single SR estimate on the NIR data and up to 0.2 dB on the RED data. All the following results have been obtained with a sliding factor equal to 5.

Then, we want to verify the effectiveness of the RegNet component of DeepSUM with respect to external registration of the images by means of cross correlation. This test should highlight the advantage of exploiting the feature space of the end-to-end trained network for the registration task. Hence, we compare two versions of our network:

• •

full network (SISRNet+RegNet+FusionNet);

• •

network without the RegNet block (SISRNet+FusionNet). We keep the registration filters but they are fixed to be a delta centered at the integer shift determined by maximum cross correlation on the ILR input images.

The full network outperforms the one without RegNet by 0.16 dB and 0.13 dB for the NIR and RED test sets, respectively, as shown in Table I. This is a significant margin and it is due to the fact that an inaccurate registration can be an important source of error for the SR reconstruction.

On the other hand, the full network, being trainable end-to-end, is able to exploit the feature space produced by SISRNet to provide a more accurate registration and help FusionNet to perform the feature merging task. We remark that the full network and the reduced network have been trained independently.

VI-B2 Comparison to State-of-the-Art

We compare the proposed MISR technique to a number of alternatives based on deep learning and model-based methods:

single image bicubic interpolation with least masked image (Bicubic); 2. 2.

averaged bicubic interpolated and registered images (Bicubic+Mean); 3. 3.

CNN-based SISR with least masked image; 4. 4.

CNN-based SISR method shared across multiple images followed by registration and averaging (SISR+Mean); 5. 5.

IBP [36]; 6. 6.

BTV [39]; 7. 7.

deep learning method based on simultaneous motion compensation and interpolation developed for video (dynamic upsampling filters (DUF) network) [58].

Table II reports the results of the comparison. It can be noticed that the proposed method outperforms all the other methods.

For all these methods, we followed the same procedure for the data preparation: bicubic interpolation and registration by phase correlation algorithm, except for DUF that computes its own registration. For MISR methods we averaged the 5 SR estimates produced by the sliding window method to ensure a fair comparison with the proposed technique.

Our IBP implementation takes as input an initial guess corresponding to our Bicubic+Mean baseline and the precomputed shifts related to the LR images using phase correlation algorithm. At each step, the LR images are estimated through the forward (HR to LR) imaging model and the error with respect to the actual LR images is back projected to the current SR image. We can observe that IBP improves over the Bicubic+Mean baseline but its performance is ultimately limited by its inability to deal with a complex and unknown degradation model. BTV implementation takes the same initial guess and precomputed shifts as in IBP with the difference that at each iteration the cost function to minimize is a L1 norm plus the bilateral regularization term. BTV shows comparable performance with respect to IBP. BTV is slightly worse due to the L1 norm data fidelity that tends to be more robust to outliers but suboptimal with respect to the mPSNR metric. The deep learning models show marked improvements over the Bicubic+Mean baseline. We consider two deep learning baselines (SISR only and SISR+Mean) that use the SISRNet architecture with the addition of a final layer projecting from the feature space to the image space, a residual connection from the (IRLR) bicubic image(s) and an increased number of parameters to roughly match the number of parameters of the full proposed architecture in order to ensure a fair comparison. The SISR+Mean result has been obtained by averaging 9 SISR images. Notice that SISR+Mean does not train the network by showing the averaged image to the loss function; it just uses the pretrained SISR network on multiple images and averages its outputs. The reason behind this choice is to provide a reference result to reader who might be interested in taking a state-of-the-art off-the-shelf SISR model, apply it to multiple images and then average the results. The comparison between SISR+Mean and the SISR only method is meant to highlight the large gain brought by exploiting both the spatial and temporal correlations, even if the LR images of a specific scene are taken under different conditions and might be wildly different from one another in terms of contrast, brightness and landscape due to temporal variations. Also, notice that SISR only is unable to improve over the simple Bicubic+Mean MISR on the NIR data. Instead, the comparison between DeepSUM and the SISR+Mean method shows the improvement brought by the introduction of FusionNet, which can exploit the slow fusion via 3D convolutions to find the best way to merge the image representations.

Another method chosen for comparison is the recent DUF network [58]. This is one of the current state-of-the-art methods for video super-resolution. DUF network processes $N$ frames in order to compute local pixel-dependent dynamic filters that are later applied on the central frame to increase its resolution and compensate motion. The network has a residual branch estimating a residual image to increase sharpness of the final SR image. The DUF network has been trained from scratch, maintaining the original structure and roughly the same number of learnable parameters with respect to our method for fair comparison. The only difference lies in using the loss function stated in Sec. IV-E instead of the one used in the original paper (Huber loss). Moreover, we always considered the first one among the 9 input LR images as central frame. The performance is worse than our proposed method and we can deduce that it highly depends on the LR input image taken to apply the dynamic local filters. We cannot know in advance which is the LR image that is closer to the HR image due to change in brightness, landscape, weather, and clouds. Involving all the LR images for HR estimation is crucial to somehow average the differences across them and try to include as much information as possible in the final SR estimate.

For completeness, we report the score achieved by DeepSUM on the unreleased test set of the PROBA-V challenge. DeepSUM achieved a score equal to 0.9474466476281652, computed as the average ratio between the mPSNR of ESA’s baseline and the mPSNR of the submitted images, over both RED and NIR data in the held-out test set.

VI-C Qualitative results

We present a set of qualitative comparisons on the RED and NIR images of our PROBA-V test set.

First of all, Figs. 5 and 8 show the multitemporal variability among the LR images and between the LR set and the HR target for the NIR and RED bands, respectively.

Figs. 6 and 9 show a visual comparison between the SR images reconstructed by the various methods for the NIR and RED bands, respectively. It can be noticed that our proposed method produces visually more detailed images, recovering finer texture and sharper edges. In order to help visualization, Figs. 7 and 10 report the absolute difference between the HR target and the SR reconstructions for the various methods after registration and compensation for absolute brightness variations (as in the mPSNR computation).

VII Conclusion

In this paper we have introduced DeepSUM, one of the first CNN architectures to deal with super-resolution from multitemporal remote sensing images. We showed that the proposed deep learning framework can successfully deal with complex degradation and temporal variation models and provide state-of-the-art performance, resulting as the best method in the PROBA-V SR challenge. Future work may focus on integrating non-local features in the network, e.g., by using graph-convolutional architectures [67], a kind of convolution that draws from ideas in graph signal processing [68, 69].