N Karantzalos

B Paragios

G Le Saux

D Moser

Tuia

Student Member, IEEE, X. Zhu L Mou

Senior Member, IEEE N Paragios
email: nikos.paragios@ecp.fr

Fellow, IEEE B Le Saux

Senior Member, IEEE G Moser

Senior Member, IEEE D Tuia

M Vakalopoulou

K Karantzalos

Multitemporal Very High Resolution From Space: Outcome of the 2016 IEEE GRSS Data Fusion Contest

Keywords: Change detection, convolutional neural networks (CNN), deep learning, image analysis and data fusion, multiresolution, multisource, multimodal, random fields, tracking, video from space

published or not. The documents may come   L'archive ouverte pluridisciplinaire

I. INTRODUCTION

P UBLIC awareness about Earth observation has raised dra- matically during the latest years. Among the main driving forces, we can mention the increased availability of very high spatial resolution (VHR) data acquired by sensors on-board drones, aircrafts, and satellites. These data make it possible to address new monitoring tasks, such as object detection and classification in urban areas (which are nowadays approached at a global scale using satellite acquisitions [START_REF] Gamba | Human settlements: A global challenge for EO data processing and interpretation[END_REF]), urban area classification [START_REF] Rottensteiner | Results of the ISPRS benchmark on urban object detection and 3D building reconstruction[END_REF], and detailed mapping at decimeter scale for precision agriculture [START_REF] Yang | Using high resolution airborne and satellite imagery to assess crop growth and yield variability for precision agriculture[END_REF]. On one hand, very high spatial resolution implies new or more complex problems, typically dealing with spatial smoothness [START_REF] Schindler | An overview and comparison of smooth labeling methods for land-cover classification[END_REF], [START_REF] Moser | Land-cover mapping by Markov modeling of spatial-contextual information[END_REF], misregistrations [START_REF] Vakalopoulou | Graph-based registration, change detection, and classification in very high resolution multitemporal remote sensing data[END_REF], changes in illumination, acquisition conditions, or sensor properties [START_REF] Pacifici | The importance of physical quantities for the analysis of multitemporal and multiangular optical very high spatial resolution images[END_REF]- [START_REF] Gómez-Chova | Multimodal classification of remote sensing images: A review and future directions[END_REF], and increased complexity of the signatures of the classes [START_REF] Tuia | Classification of very high spatial resolution imagery using mathematical morphology and support vector machines[END_REF], [START_REF] Pacifici | A neural network approach using multi-scale textural metrics from very high-resolution panchromatic imagery for urban land-use classification[END_REF]. These challenges generally raise the need for more and more sophisticated methodologies. On the other hand, the advances in hardware and the increased availability of data makes it possible to use training models that are increasingly able to extract complex and abstract features. Prominent examples are deep learning methods, which train feature extraction filters directly and nonparametrically from the data [START_REF] Paisitkriangkrai | Semantic labeling of aerial and satellite imagery[END_REF], [START_REF] Volpi | Dense semantic labeling of subdecimeter resolution images with convolutional neural networks[END_REF], or feature learners, which mine the parameter spaces of existing filter banks [START_REF] Glocer | Online feature selection for pixel classification[END_REF], [START_REF] Tokarczyk | Features, color spaces, and boosting: New insights on semantic classification of remote sensing images[END_REF] or of their combinations [START_REF] Tuia | Multiclass feature learning for hyperspectral image classification: Sparse and hierarchical solutions[END_REF].

While VHR image processing is an active field of investigation, the processing of video data taken from satellite platforms is still in its infancy. On one hand, video from space (also known as "space video") is not currently widespread. A more common approximation is represented by sensors that can retarget a specific region and provide a multiangular sequence. For example, the hyperspectral CHRIS/PROBA sensor can provide angular sequences, which were used for chlorophyll [START_REF] Dorigo | Improving the robustness of cotton status characterisation by radiative transfer model inversion of multi-angular CHRIS/PROBA data[END_REF] and leaf area index retrieval [START_REF] Dorigo | Improving the robustness of cotton status characterisation by radiative transfer model inversion of multi-angular CHRIS/PROBA data[END_REF], [START_REF] Schlerf | Vegetation structure retrieval in beech and spruce forests using spectrodirectional satellite data[END_REF], and the multiangular multispectral WorldView-2 sensor was used for height retrieval [START_REF] Lemoine | Fast surface height determination using multi-angular WorldView-2 ortho ready urban scenes[END_REF], detection [START_REF] Turlapaty | A hybrid approach for building extraction from spaceborne multi-angular optical imagery[END_REF], land use classification [START_REF] Matasci | Understanding angular effects in VHR imagery and their significance for urban land-cover model portability: A study of two multi-angle in-track image sequences[END_REF], and tracking [START_REF] Meng | Object tracking using high resolution satellite imagery[END_REF]. On the other hand, with the increasing development of small and cheap commercial satellites, such as the SkySat-constellation of Skybox Imaging, wide area spaceborne remote sensing videos are becoming increasingly accessible at a reasonable cost. Indeed, while the processing of aerial videos collected by sensors onboard terrestrial vehicles or drones is not yet a major trend in remote sensing, first applications can be found in the recent literature. These applications include tracking or navigation using videos from web cameras [START_REF] Ribas | Assessing the suitability of video imaging for studying the dynamics of nearshore sandbars in tideless beaches[END_REF], terrestrial vehicles [START_REF] Leung | Localization in urban environments by matching ground level video images with an aerial image[END_REF], [START_REF] Gu | Microwave unmixing with video segmentation for inferring broadleaf and needleleaf brightness temperatures and abundances from mixed forest observations[END_REF], aircrafts [START_REF] Linger | Aerial image registration for tracking[END_REF], and satellites [START_REF] Kopsiaftis | Vehicle detection and traffic density monitoring from very high resolution satellite video data[END_REF]. Information extraction from space videos remains a major challenge, although the acquisition technology is available, as demonstrated by the Iris camera operated by Urthecast on the International Space Station (ISS). 1The Data Fusion Contest (DFC) 2016 addressed this challenge.

Since 2006, the Image Analysis and Data Fusion Technical Committee (IADF TC2 ) of the IEEE Geoscience and Remote Sensing Society (IEEE-GRSS) has organized an annual DFC, in which a dataset has been released free of charge to the international community along with a data fusion competition [START_REF] Alparone | Comparison of pansharpening algorithms: Outcome of the 2006 GRS-S data fusion contest[END_REF]- [START_REF] Vo | Processing of extremely high resolution LiDAR and RGB data: Outcome of the 2015 IEEE GRSS data fusion contest. Part B: 3D contest[END_REF]. To tackle the aforementioned challenges implied by new VHR data with very high temporal resolution, the DFC 2016 proposed a competition related to the multimodal processing of image sequences and video streams acquired from space.

The DFC 2016 released to the international community an image dataset involving VHR multiresolution and multisensor imagery and video information. The dataset was composed of a series of two VHR multispectral and panchromatic images and of a video stream acquired by the Iris camera over the city center of Vancouver, Canada (see Section II). The data were provided by the Deimos Imaging and Urthecast companies.

The competition was framed as follows: each participating team was invited to submit an original open-topic manuscript addressing an information extraction problem based on the dataset provided. All submissions were evaluated by an Award Committee, composed of the organizers of the Contest, of industrial representatives of Deimos Imaging and Urthecast, and of several past Chairs of the IADF TC. The manuscripts were ranked on the basis of scientific novelty and originality, methodological approach, experimental discussion, and quality of presentation.

In this paper, after describing the dataset of the DFC 2016 (see Section II), we will discuss first the overall scientific outcome of the contest as a whole (see Section III). Then, we will focus in more detail on the approaches proposed by the first and second place teams (see Sections IV and V, respectively). Finally, conclusions will be drawn in Section VI.

II. DATASET OF THE DFC 2016

The dataset of the DFC 2016 included VHR imagery and video from space covering an urban and harbor area in Vancouver, Canada (49 • 15'N, 123 • 6'W; see Fig. 1).

1) VHR Images: A temporal pair of VHR acquisitions (Deimos-2 standard products), each composed of a panchromatic image at 1-m point spacing and a multispectral image at 4-m point spacing, was included in the dataset. The two acquisition dates were March 31 and May 30, 2015. Deimos-2 operates from a Sun-synchronous orbit at a mean altitude of 620 km. The spacecraft design is based on an agile platform for fast and precise off-nadir imaging (up to ±30 • over nominal scenarios and up to ±45 • in emergency cases), and carries a push-broom VHR camera with five spectral channels: one panchromatic channel and four multispectral channels corresponding to red (R), green (G), blue (B), and near-infra-red (NIR) bands.

For each date, four image products were provided for the contest: the panchromatic and the multispectral images were made available at both levels 1B and 1C. A level 1B product is calibrated and radiometrically corrected but is not resampled, and geometric information is provided separately (see Fig. 2). A level 1C product is calibrated, radiometrically corrected, and manually orthorectified and resampled to a map grid, while the geometric information is contained in the GeoTIFF tags (see Fig. 3). As a consequence, the two level 1C images cover exactly the same ground area.

2) Video From Space: A full-color ultrahigh definition (UHD) video acquired by the Iris camera on-board the ISS was provided. The acquisition took place on July 2, 2015, over the harbor of Vancouver (see Fig. 4). Iris is a high-resolution camera installed on the Zvezda module of the ISS. It uses a complementary metal oxide semiconductor (CMOS) detector to capture RGB videos with a ground sample distance as fine as 1 m at three frames per second. In the Iris video of the contest, image frames were fully orthorectified and resampled to 1 m. The resulting frame format was 3840 × 2160 pixels and covered approximately 3.8 km × 2.1 km.

For the first time in the history of the DFC, the dataset included a video captured from space. This made it possible to explore new exciting applications of Earth observation, such as tracking or live monitoring. Along with the imagery that was made available at the aforementioned processing levels, the dataset allowed participants to demonstrate a broad variety of innovative methods including registration, pan-sharpening, multitemporal analysis, change detection, object detection, and image classification.

III. SUBMISSIONS AND RESULTS

A. Submissions

Twenty-four papers were submitted to the DFC 2016. They proved to be of great diversity, because they addressed heterogeneous topics and proposed various methodological solutions for each topical area. Fig. 5 summarizes both these aspects.

1) Topic-wise, a majority of papers dealt with two thematic areas: Classification / change detection and tracking. While the former had been mainstream for years for the DFC, the latter emerged as a new topical area. Indeed, it also attracted most of the participants' efforts, essentially because the data were well suited for this task. Furthermore, in addition to the two aforementioned main topics, participants also considered some other classical fusion problems (pan-sharpening and registration) as well as new topical areas such as population modeling or estimation of human activities. 2) Methodologically speaking, we observed a wide spectrum of proposals. Approaches that were quite established in remote sensing (e.g., feature extraction or sparse methods) were challenged by techniques issued from the video processing literature (e.g., optical flow and tracking) and from the deep learning and computer vision literature (e.g., convolutional neural networks (CNNs) and structured output modeling). This trend was expected because deep learning is becoming more and more prominent in the remote sensing area as well (e.g., [START_REF] Paisitkriangkrai | Semantic labeling of aerial and satellite imagery[END_REF] and [START_REF] Volpi | Dense semantic labeling of subdecimeter resolution images with convolutional neural networks[END_REF]), and allows robust models to be learned even through simple fine tuning of existing architectures drawn from other methodological areas [START_REF] Campos-Taberner | Processing of extremely high resolution LiDAR and RGB data: Outcome of the 2015 IEEE GRSS data fusion contest. Part A: 2D contest[END_REF], [START_REF] Penatti | Do deep features generalize from everyday objects to remote sensing and aerial scenes domains?[END_REF]. Another interesting observation can be made when considering the data types used by the participating teams. Three separate sources of data (multispectral VHR, panchromatic VHR, and video) were provided but no fixed topic was imposed. Therefore, each team was free to use the most relevant sources to address their problem. Fig. 6 summarizes the choices made by the teams: 14 out of 24 teams used the video, which confirms the interest and potential of this new data modality. Only three teams used all the data at their disposal, but 15 teams used more than one data source. Another relevant point is that five teams further enriched the dataset with ancillary data (e.g., building footprints from OpenStreetMap, road maps form Google Street View, digital elevation models, and areal statistics). Indeed, this choice reminded that data of great quality are available at no cost on the Internet and can be included in processing pipelines to either make them more effective or decrease their computational load.

B. Winners

The 24 submissions were evaluated by the Award Committee and four papers were awarded. They were presented during the 2016 IEEE International Geoscience and Remote Sensing Symposium (IGARSS) in Beijing, China. The papers awarded were as follows. In the following two sections, the approaches proposed by the first and second ranked teams are discussed, and their main results are summarized. More details can be found in the corresponding papers in the IGARSS 2016 proceedings [START_REF] Mou | Spatiotemporal scene interpretation of space videos via deep neural network and tracklet analysis[END_REF], [START_REF] Vakalopoulou | Simultaneous registration, segmentation and change detection from multisensor, multitemporal satellite image pairs[END_REF].

IV. SPATIOTEMPORAL SCENE INTERPRETATION OF SPACE VIDEOS VIA DEEP NEURAL NETWORK AND TRACKLET ANALYSIS

This section discusses the approaches proposed by the winners of the Contest. The work focused on a novel framework for spatiotemporal analysis of spaceborne remote sensing video data. For this purpose, a deep neural network was proposed to make use of high resolution satellite imagery in order to achieve a fine-resolution spatial scene labeling map. Moreover, a sophisticated approach was proposed to analyze activities and estimate traffic density from 150 000+ tracklets produced by a Kanade-Lucas-Tomasi (KLT) keypoint tracker. Both visual and quantitative analysis of the experimental results demonstrated the effectiveness of the proposed approach.

A. Motivation

In comparison with the widely used static satellite imagery or video data acquired from unmanned aerial vehicles, space videos have remarkable advantages. For instance, a space video gathered from the ISS or by small satellites provides an opportunity to observe the dynamic states at urban or suburban scale and usually covers large areas. However, due to the relatively coarse spatial resolution, the poor appearance of moving targets (e.g., cars) creates ambiguities in tracking and challenges scene interpretation tasks, such as activity analysis and density estimation. In this paper, the focus was the semantic scene interpretation of space videos, using spatiotemporal analysis. More specifically, the aim was to answer the following three questions.

1) What land-use categories are observed in the video and where do they appear (spatial scene labeling)? 2) Which activities can be observed in the video (temporal activity analysis)? 3) And how dense is the traffic flow (traffic density estimation)? Fig. 7. Overview of the pipeline proposed by the first ranked team. A novel deep neural network architecture tailored to pixel-wise spatial scene labeling was proposed for spatial analysis and a tracklet-based temporal parsing method was introduced for activity analysis and traffic density estimation (from [START_REF] Mou | Spatiotemporal scene interpretation of space videos via deep neural network and tracklet analysis[END_REF]).

Fig. 8. Foreground mask yielded by the frame differencing method with an estimated background model [START_REF] Loy | Time-delayed correlation analysis for multi-camera activity understanding[END_REF]. It can be seen that the camera movement makes it difficult to obtain a valid background model.

Owing to their rich temporal information, video data can be exploited for activity analysis and traffic density estimation, which are classical computer vision applications. Since videos from space have only recently become available, their use to augment the semantic interpretation of covered geographic areas has rarely been addressed in the remote sensing community so far. In a pioneering work moving from static image analysis to video parsing, Kopsiaftis and Karantzalos [START_REF] Kopsiaftis | Vehicle detection and traffic density monitoring from very high resolution satellite video data[END_REF] developed an automatic vehicle detection and traffic density estimation model for VHR satellite video data. This model was based on background estimation followed by background subtraction in every video frame. However, it could not obtain the trajectories of moving objects for further analysis of temporal activities. In addition, it was almost impossible to automatically obtain a stationary background model, since the camera was moving during the video acquisition (see Fig. 8). In this paper, the goal was to develop a sophisticated spatio-temporal scene analysis approach for space videos that was for the first time capable of performing spatial scene labeling and temporal activity analysis on multimodal image data.

In the following, the proposed approach will be summarized (see the processing chain in Fig. 7): first, a deep neural network architecture tailored to pixel-wise spatial scene labeling was proposed for spatial analysis (see Section IV-B); then, a tracklet-based temporal parsing method was introduced for activity analysis and traffic density estimation (see Section IV-C).

B. Spatial Scene Labeling via Deep Neural Network

For spatial scene labeling, the recent study [START_REF] Penatti | Do deep features generalize from everyday objects to remote sensing and aerial scenes domains?[END_REF] on the use of typical CNNs such as VGG16 [START_REF] Simonyan | Very deep convolutional networks for large-scale image recognition[END_REF] was extended to classify category-independent regions sampled from VHR images. In particular, unpooling and oversegmentation were exploited to refine the typically low resolution (LR) results produced by the CNN.

1) CNN: A typical convolutional network interleaves convolutional layers and pooling layers often followed by one or more fully connected layers. The specific architecture of the network was topologically identical to VGG16 [START_REF] Simonyan | Very deep convolutional networks for large-scale image recognition[END_REF].

The VGG16 network makes use of a stack of convolutional layers with a very small receptive field of 3 × 3, rather than using larger ones, such as 5 × 5 or 7 × 7. The reason is that 3 × 3 convolutional filters are the smallest kernels to capture patterns in different directions, such as center, up/down, and left/right, and bring multiple advantages. First, the use of stacks of small convolutional filters separated by rectified linear units (ReLUs) increases the nonlinearities inside the network. Then, the number of parameters to be learned decreases. The convolution stride in the network is fixed to 1 pixel. The spatial padding of the convolutional layer input is such that the spatial resolution is preserved after convolution, i.e., 1 pixel for the 3 × 3 convolutional layers. Spatial pooling is achieved by carrying out several max pooling layers, which follow some of the convolutional layers. Max pooling is performed over 2 × 2 pixel windows with stride 2. After several convolutional and pooling layers, the high-level reasoning in the neural network is achieved via fully connected layers. 2) Fine Tuning: Here, a fully connected layer was added and the new network was fine-tuned using the patches sampled from the multispectral image, which made the network significantly easier to train than a CNN trained from scratch. Fine tuning was performed by retraining the final fully connected layer. Fine tuning the ImageNet-trained VGG16 model for the specific task addressed here was clearly advantageous and yielded decent results, since the amount of labeled training data for the Vancouver images was insufficient. Finally, the fine-tuned network was applied to the first video frame for spatial scene labeling.

3) Refinement: Using pretrained architectures was very well suited at extracting features from data by spatially shrinking the feature maps. Pooling was necessary to allow gradually gathering information over the feature maps, and, to make the network computationally feasible. However, this produced a scene labeling map with a reduced spatial resolution. To obtain a high-resolution fine-edged scene labeling map, it was necessary to refine the LR probability maps produced by the CNN.

Several strategies have been proposed to perform this step, ranging from excluding spatial pooling in convolutional layers [START_REF] Sherrah | Fully convolutional networks for dense semantic labelling of high-resolution aerial imagery[END_REF] to learning a deconvolution network to upsample the maps [START_REF] Volpi | Dense semantic labeling of subdecimeter resolution images with convolutional neural networks[END_REF]. A simpler approach, shown in Fig. 7, was applied. The main ingredients were unpooling and region refinement.

Unpooling extended the feature maps, as opposed to pooling (see Fig. 10), and took the LR probability maps as input to produce VHR but grainy maps called coarse probability maps. Then, these coarse maps were projected to superpixels3 by as-signing, to each superpixel, the average value of its corresponding region in the coarse probability maps as

p(s i ) = 1 |s i | c∈s i p(c) ( 1
)
where s i represents the ith superpixel, |s i | is the number of pixels belonging to s i , c stands for a generic pixel in the superpixel, and p(•) denotes probability map entries. In this way, the refined fine-edged probability maps were obtained. Fig. 9 compares the original LR output (top) with the refined one obtained by the proposed approach (bottom).

C. Temporal Analysis Using Tracklets

Object tracking for temporal analysis has been well studied in computer vision for more than two decades. Video from space, however, brings new methodological challenges, such as 1) the low frame rate, which undermines the common motion smoothness assumption; and 2) the limited spatial resolution of moving objects, which prevents the use of robust and discriminative appearance-based tracking algorithms.

1) Extracting Valid Tracklets: In this paper, tracklets, instead of the trajectories obtained by the object tracking algorithms, were used to analyze the activity. A tracklet is a fragment of a trajectory obtained by a tracker with a short period. Compared to long-term object tracking, they are more reliable. In our approach, a KLT keypoint tracker [START_REF] Tomasi | Detection and tracking of point features[END_REF] was used to extract 150 000+ tracklets from the space video [see Fig. 11(a)]. If the camera was stationary, the tracklets would have been entirely caused by real object movements (positive tracklets). Unfortunately, the camera kept slightly moving throughout the entire video acquisition (see the effect of such movements over 30 s in Fig. 4), which meant that the obtained tracklets included many negative samples caused by parallax. Therefore, it became necessary to remove negative tracklets before carrying out any follow-up activity analysis. Here, a prior constraint on the tracklet denoising was enforced according to the observation that the displacements of moving object tracklets were usually larger than the displacements of still object tracklets (e.g., buildings) caused by the camera movement. In this paper, the threshold value was fixed to 11 pixels, which corresponded to a minimum car speed of 20 km/h, given the spatial resolution of the video 2) Activity Analysis and Traffic Density Estimation: For the remaining tracklets, the angle of motion θ = arctan(ζ/ε) was calculated according to the displacement vector (ε, ζ) of the tracklet, and 2-D information was transformed into the (ε, ζ, θ) space, in which k-means was used for clustering. The output could be interpreted as activity information of every cluster. The traffic density was estimated on the superpixels superimposed on the scene. For every superpixel, the number of positive tracklets was calculated and the density was estimated for the whole scene at every frame. Finally, the normalized density estimation map was obtained by dividing the number of positive tracklets in each superpixel by the maximum number among all superpixels.

D. Results and Discussion

In order to quantitatively evaluate the performance of the proposed approach, a ground truth for spatial scene classification 4 was built by manually labeling the first frame of the video to four land-use scene categories, according to a detailed visual analysis and some prior information. The spatial scene labeling result is illustrated in the left panel of Fig. 12. Cyan, red, yellow, and blue represent urban, vegetative cover, harbor, and river, respectively. 400 labeled patches (100 in each category) with a size of 224 × 224 were randomly selected from the multispectral image as the training set. On the held out test data, the developed approach can achieve an average accuracy of 90.5%, overall accuracy of 96.5%, and Kappa of 0.9353. 4 The manually labeled ground truth for spatial scene labeling and the corresponding code for visualization can be found at http://www.sipeo. bgu.tum.de/downloads/gt4dfc16video.rar.

Regarding activity analysis, the standard quantitative measure of precision is based on the number of true positives (TP) and false positives (FP) and is computed as precision = TP/(TP + FP). For quantitative evaluation, the whole visible area of the video was used and tracklet ground truth was manually determined for 30 s. The precision of the temporal activity result reached 95.3%. Fig. 13 provides a close zoom of the activity analysis in which two moving cars were detected and their moving directions and speeds could be estimated. In addition, the traffic density estimation result is presented in the bottom right image of Fig. 12. The estimated density highly correlates with the actual traffic situations, e.g., high values in main roads and intersections.

It is also worth noting that, here, the focus was on addressing two relatively independent problems, i.e., spatial scene labeling and temporal activity analysis. Accordingly, an "open-loop" solution in which the scene labeling was accomplished without making use of the temporal analysis results was designed. Indeed, a "closed-loop" system, in which the spatial scene labeling is further refined using the outcome of the temporal analysis, will represent an interesting future extension of this work.

V. SIMULTANEOUS REGISTRATION, SEMANTIC SEGMENTATION, AND CHANGE DETECTION FROM MULTISENSOR, MULTITEMPORAL SATELLITE IMAGE PAIRS

This section details the method proposed by the second ranked team of the Contest.

A. Motivation

The current generation of space-borne and airborne imaging sensors are generating nearly continuous streams of massive, multitemporal, high-resolution remote sensing data. However, in order to efficiently exploit these datasets, their accurate coregistration is the first indispensable processing step along with any further analysis and change detection procedures. In the framework of the DFC 2016, a methodology was proposed that was able to ingest information from multisensor and multitemporal satellite image pairs, by addressing concurrently the tasks of registration, semantic segmentation, and change detection. These three tasks are usually addressed separately.

A number of studies were focused on jointly tackling the problem of registration and semantic segmentation for mainly video sequences or medical images [START_REF] Chen | Multiphase joint segmentationregistration and object tracking for layered images[END_REF]- [START_REF] Mayer | A supervised framework for the registration and segmentation of white matter fiber tracts[END_REF]. Similar research efforts were focused on jointly addressing the tasks of segmentation and tracking in image video sequences [START_REF] Tsai | Motion coherent tracking using multi-label MRF optimization[END_REF]- [START_REF] Karantzalos | Implicit Free-Form-Deformations for Multi-frame Segmentation and Tracking[END_REF]. However, such formulations cannot exploit sparse multitemporal datasets with changes in-between the various acquisition dates. To tackle such problems, a framework was designed for jointly register images, detect changes, and producing semantic segmentation maps at the same time.

In particular, the formulation of [START_REF] Vakalopoulou | Graph-based registration, change detection, and classification in very high resolution multitemporal remote sensing data[END_REF] and [START_REF] Vakalopoulou | Simultaneous registration and change detection in multitemporal, very high resolution remote sensing data[END_REF], was extended by adding a graph tackling the semantic segmentation problem. The developed method jointly integrated the following three types of energies into an energy minimization framework: Fig. 12. Final results. From left to right: spatial scene labeling (cyan, red, yellow, and blue represent urban, vegetative cover, harbour, and river, respectively.), temporal activity analysis (each color corresponds to a cluster), and traffic density estimation (from [START_REF] Mou | Spatiotemporal scene interpretation of space videos via deep neural network and tracklet analysis[END_REF]). Fig. 13. Zoomed area. From left to right and top to bottom: frame #1, frame #61, a red-cyan color composite used to illustrate the pixel-wise difference, all tracklets, positive tracklets, and final activity analysis including their moving speeds (from [START_REF] Mou | Spatiotemporal scene interpretation of space videos via deep neural network and tracklet analysis[END_REF]).

1) data-driven costs issued from classification scores5 (various classification algorithms can be employed, e.g., [START_REF] Lagrange | Benchmarking classification of earth-observation data: From learning explicit features to convolutional networks[END_REF]); 2) registration metrics (e.g. similarity metrics); 3) change detection scores. These energies were efficiently coupled with local geometric constraints in the context of a higher order graph. Reduction methods were used to map this graph into a pairwise one, which was optimized using efficient linear programming. The experimental results were very promising: less than 2 pixels in terms of mean displacement error for the registration; above 77% in most cases regarding the completeness and correctness rates of the semantic segmentation; and change detection accuracy higher than 70%.

B. Methodology

1) Graph Formulation: Without loss of generality, three distinct graphs with exactly the same topology/geometry and of lower resolution than the image dimensions were considered. They corresponded to a lower resolution grid than the input image and every original pixel was assigned to the closest graph node to compute the related energy contributions. Following the notations of [START_REF] Vakalopoulou | Graph-based registration, change detection, and classification in very high resolution multitemporal remote sensing data[END_REF] and [START_REF] Vakalopoulou | Simultaneous registration and change detection in multitemporal, very high resolution remote sensing data[END_REF], the first graph, G reg , involved nodes where the labels corresponded to deformation vectors from the registration process, i.e., a mapping between the source and the target images. The second graph, G ch , referred to nodes with binary labels expressing changes in the temporal domain. The third graph, G seg , which was the one introduced in this work, referred to the labels representing the semantic segmentation of the image. In terms of connectivity, each graph was endowed with a first-order neighborhood system accounting for local consistency on the corresponding label space. Furthermore, two intergraph connections were considered, one associating the corresponding nodes in the detection and registration graphs, and the other involving hypercliques that linked the corresponding nodes of all graphs.

This graph structure was superimposed on the image in such a way that each node of the graph depended on a subset of pixels in its vicinity. The objective function was defined on the entire image domain, and every pixel contributed to the graph with a weight depending on its distance from a node in the support domain. The graph structure defined the support domain and the computational complexity of the implementation, thus, a finer grid generated a higher computational load. In particular, the dimensions of the graph were related to the image size, thus representing a tradeoff between accuracy and computational complexity.

2) Energy Formulation: By integrating the three graphs in a single framework, the proposed energy was a function E reg,ch,seg (l reg , l ch , l seg ), in which the labels of each node p of the coupling graph were l p = [l ch p , l seg p , l reg p ]. The term l ch p ∈ {0, 1} represents the label for change detection, l seg p ∈ {0, 1} represents the label for binary semantic segmentation, and l reg p ∈ Δ is the registration label. Δ = {d 1 , . . . , d n } corresponds to all possible displacements. The label space could be described as L = {0, 1} × {0, 1} × Δ. In the next section, we detail the key ideas of all the terms used to calculate the energy function.

C. Energy Terms

The energy was composed of a series of terms related to the three tasks of change detection, registration, and semantic segmentation. The energy terms could be either singletons (i.e., functions of the score of one specific node), pairwise (i.e., depending on the spatial neighborhood for one single task), or coupled (i.e., using the results of two or more tasks on a single node simultaneously).

1) Singleton and Pairwise Terms: Both the registration and change detection terms depended on the two considered images simultaneously, and therefore, only had pairwise energy contributions. Regarding the semantic segmentation process, the goal was to assign the correct semantic segmentation label to each node of the target image. In particular:

1) the registration (V pq ,reg (l reg p , l reg q )) and change detection (V pq ,ch (l ch p , l ch q )) pairwise terms followed the same formulation as in [START_REF] Vakalopoulou | Graph-based registration, change detection, and classification in very high resolution multitemporal remote sensing data[END_REF] and [START_REF] Vakalopoulou | Simultaneous registration and change detection in multitemporal, very high resolution remote sensing data[END_REF] and penalized neighboring nodes with different registration or change detection labels; 2) the semantic segmentation graph contained a term with the classification score for each node (i.e., a singleton term) (V seg (l seg p )) and a pairwise term (V pq ,seg (l seg p , l seg q )), which penalized different semantic segmentation labels in neighboring nodes. The two terms followed the same formulation as in [START_REF] Vakalopoulou | Simultaneous registration, segmentation and change detection from multisensor, multitemporal satellite image pairs[END_REF].

2) Coupled Intergraph Energy Terms: As stated previously, two intergraph connections were considered. They translate into two coupled terms addressing two and three tasks simultaneously.

1) The coupling between the registration and change detection terms (V reg,ch (l reg p , l ch p )) was achieved through the interconnection between the two graphs. In the absence of change, the cost was calculated using a similarity function, while in the presence of change, a fixed cost was used.

2) The coupling between the three terms was performed using one potential term (V reg,ch,seg (l reg p , l ch p , l seg p )) that penalized different semantic segmentation labels in the source and target images for all possible displacements in cases of no-change and the same ones in cases of change.

3) Final Global Energy: The global energy function was a linear combination of the aforementioned energy terms, weighted by six nonnegative weight parameters w 1 , w 2 , . . . , w 6 . These parameters should be optimized prior to the minimization of the energy with respect to the labels l ch p , l seg p , and l reg p of all nodes p. Details on the specific formulation of the energy can be found in [START_REF] Vakalopoulou | Simultaneous registration, segmentation and change detection from multisensor, multitemporal satellite image pairs[END_REF].

D. Experimental Results and Evaluation

The evaluation of the developed framework was performed on all the data described in Section II. From now on, we will refer to the first Deimos-2 acquisition (March 2015) as [D1], to the second Deimos-2 image (May 2015) as [D2], and to the Iris video as [V]. The level 1C image products were used, and were further radiometrically corrected, and then, pan-sharpened based on the standard high-pass filter method, resulting into an overlapping image pair of approximately 12 760 × 11000 pixels. The area of the overlap between the image pairs and the Iris video frames was approximately 4 720 × 2 680 pixels. In order to employ an additional image/map, which could serve as a reference/target map while contributing to the automation of the subsequent training procedure, an image mosaic [G] and the corresponding land map 6 were downloaded from Google Maps APIs. 7 6 By "land map" we intend the raster simplified view that can be used in Google Maps as an alternative to the satellite image. 7 All the considered raw unregistered data and ground truth data, and several experimental results can be viewed at: http://users.ntua. gr/karank/Demos/DemoContest16.html for all node p ∈ G do 8:

for all pixel j in the vicinity of p do 9:

z pj ← η(j) 10:

end for 11:

end for 12:

for all label l p ∈ {0, 1} × {0, 1} × Δ do 13:

for all pixel j do 14:

sum j = w 1 V reg,ch (l p ) + w 3 V seg (l p )+ w 2 V reg,ch,seg (l p ) 15: U p (l p ) + = z pj • sum j 16:
end for 17:

end for 18:

for all label l p do 19:

for all label l q do 20:

B pq (l p , l q ) = w 4 V pq ,reg (l p , l q ) +w 5 V pq ,ch (l p , l q ) + w 6 V pq ,seg (l p , l q ) 21:

end for 22:

end for 23:

l ← FastPD-Optimization(U, B) 24:

for all node p do 25:

for all pixels j in the vicinity p do 26:

l j + = z pj • l p 27: end for 28:

end for 29:

Transform the source image 30: end while 1) Implementation Details: The parameters of the developed framework were not constant across different datasets. Grid search over a small area of the dataset was used for fine tuning all the parameters. Then, this set of parameter values was applied to the entire dataset. This was of course a limitation of the developed method, however, with satellite or aerial data of similar spatial and spectral resolution, significant differences were not expected. Details on techniques and implementation can be found in [START_REF] Vakalopoulou | Simultaneous registration, segmentation and change detection from multisensor, multitemporal satellite image pairs[END_REF], while for the optimization procedure, FastPD was employed [START_REF] Komodakis | Approximate labeling via graph cuts based on linear programming[END_REF]. For an overview of the method, see the pseudocode in Algorithm 1.

2) Classification Scores: In order to estimate the classification scores, a patch-based deep learning framework was used. It reported high classification accuracy rates for a number of cases [START_REF] Lagrange | Benchmarking classification of earth-observation data: From learning explicit features to convolutional networks[END_REF], [START_REF] Mnih | Machine learning for aerial image labeling[END_REF], [START_REF] Vakalopoulou | Building detection in very high resolution multispectral data with deep learning features[END_REF]. These deep architectures require a large number of training data. To this end, the Roads, Buildings, Sea, and Vegetation classes were automatically extracted from the IEEE JOURNAL OF SELECTED TOPICS IN APPLIED EARTH OBSERVATIONS AND REMOTE SENSING Fig. 14. Chessboard visualization with the unregistered inputs (left) and the registered results (right) multisensor data. From [START_REF] Vakalopoulou | Simultaneous registration, segmentation and change detection from multisensor, multitemporal satellite image pairs[END_REF]. Google land map (based on the different colors). To do so, all datasets were registered to the Google image mosaic in order to relate every pixel to one of the Google land map colors. Furthermore, polygons for the Vessels, Building shadow, Vegetation shadow, Soil, and Clouds classes were manually annotated. For the Iris video sequence, Vegetation shadows was merged with Vegetation. Spectral analysis on the derived numerous polygons and probabilities was employed (as in [START_REF] Mnih | Machine learning for aerial image labeling[END_REF]) in order to create the final training and validation sets for the eight terrain classes, i.e., Roads, Buildings, Building Shadows, Soil, Sea, Ship/vessels, Vegetation, and Vegetation shadows.

The training for the experiment considering the D1 and D2 images was performed on the large (12 760 × 11000 pixels) overlapping region with eight classes. On the contrary, seven classes were employed for V. Based on the derived polygons, numerous patches of size 21 × 21 including all the spectral bands were created by centering each patch on the annotated pixels [START_REF] Vakalopoulou | Building detection in very high resolution multispectral data with deep learning features[END_REF]. Approximately 200 000 randomly selected patches per class were used for the Deimos-2 data, while 50 000 were used for the Iris video. A CNN with ten layers was trained: the network was made of two convolutional layers, each followed by tangent and max pooling layers, and ended with two fully connected, a tangent, and a linear layer. The model was trained with a learning rate equal to 1 for 40 epochs; afterwards, the learning rate was halved every 2 epochs. The same CNN architecture and configuration settings were used for the classification of the Iris video frames.

3) Registration Results: For the validation of the registration results (see Fig. 14), several ground control points were manually collected in all resulting image pairs. In Table I, the These results are compared with those obtained by the algorithm in [START_REF] Karantzalos | Efficient and automated multi-modal satellite data registration through MRFs and linear programming[END_REF] for both tasks of registering the different sources to the Google mosaic and registering each frame of the video sequence to the first frame. One can observe that both methods resulted in subpixel accuracy, indicating that the coupling of registration priors with semantic segmentation and change detection priors does not harm the quality of the registration.

4) Semantic Segmentation Results:

The framework was validated for the detection of three different classes out of the original eight, namely: Buildings, Ship/vessels, and Vegetation. For the semantic segmentation and change detection results, the quantitative evaluation of the framework was performed using the completeness, correctness, and overall quality criteria at the object level. In particular, after an intensive manual photointerpretation procedure, ground-truth/reference data were created for the three aforementioned classes in four different smaller regions. The resulting true positives, false negatives, and false positives were calculated on the validation dataset after the application of the developed framework. The overall accuracy (OA) and average accuracy (AA) were also calculated. After the optimization and based on the polygons of the Sea class, which were derived automatically from Google Maps, all objects segmented as Buildings in the Sea and as Ship/vessels in the land areas were neglected.

As one can observe in Fig. 15, although the original classification scores constrained significantly the result, the developed framework could improve the semantic segmentation results in several image regions as compared with the classification accuracies obtained by the CNN. 8 In particular, the quantitative evaluation (see Table II) indicated that the detection completeness rates were above 78% (apart from the Buildings class in the D1-D2 pair), and the detection correctness rates were above 72% in all cases. The highest rates were for the Vegetation class, indicating that the NIR Deimos-2 band significantly contributed to class separation. Most semantic segmentation errors were due to false alarms near the port, pier, and ship wake on the 8 http://users.ntua.gr/karank/Demos/DemoContest16.html sea. Moreover, the Buildings and Roads classes were confused in certain cases. 5) Change Detection Results: Similar qualitative errors were observed on the change detection results (see Fig. 16) obtained from both image pairs. Quantitative results (see Table III) suggested lower completeness and correctness rates as compared to the semantic segmentation task, as expected. These results were mainly due to a number of false positives in the dense urban regions where the relief displacements were significant due to the tallest buildings and skyscrapers.

Additionally, experiments with the unsupervised change detection method in [START_REF] Vakalopoulou | Graph-based registration, change detection, and classification in very high resolution multitemporal remote sensing data[END_REF] were also performed to compare the performance with and without semantic segmentation labels. The proposed method not only provided additional information (i.e., semantic segmentation maps and from-to change trajectory maps) but also resulted in higher accuracy rates. In particular, the developed system led to an improvement of at least 6% in the mean overall quality compared to the method in [START_REF] Vakalopoulou | Graph-based registration, change detection, and classification in very high resolution multitemporal remote sensing data[END_REF], indicating that the use of the semantic segmentation labels allows the number of false alarms to be reduced and more accurate change maps to be produced.

Finally, the results of the proposed method were compared with those of [START_REF] Vakalopoulou | Graph-based registration, change detection, and classification in very high resolution multitemporal remote sensing data[END_REF] and [START_REF] Karantzalos | Efficient and automated multi-modal satellite data registration through MRFs and linear programming[END_REF] in terms of computational performance. For the proposed higher order formulation, reduction methods were employed to reformalize as a pairwise model, which is indeed quite computationally demanding. In particular, for a 1000 × 1000 pixel image, the deformable registration component alone (see [START_REF] Karantzalos | Efficient and automated multi-modal satellite data registration through MRFs and linear programming[END_REF]) required 2-3 min. Then, the registration and the unsupervised change detection (like in [START_REF] Vakalopoulou | Graph-based registration, change detection, and classification in very high resolution multitemporal remote sensing data[END_REF]) converged after 6 min. The proposed framework converged after approximately 12 min with the same configuration setting.

VI. CONCLUSION

In this paper, we discussed the scientific results of the 2016 IEEE GRSS DFC organized by the IEEE GRSS Image Analysis and Data Fusion Technical Committee. We described the dataset and described the overall scientific outcomes of the competition, by first presenting its overall results, and then, focusing on the strategies proposed by the first and second place teams. These teams made use of all the data available and developed methodologies rooted in the latest advances in computer vision and machine learning: The winners proposed a CNN and a tracking algorithm, while the second place team developed a random field model performing multiple tasks simultaneously. Both methodologies were effective and tailored to the peculiarities of the new generation remote sensing data that were released for the contest, including in particular video from space (introduced for the first time in this annual competition).

Moreover, the results of the DFC 2016 also exceeded these two individual methods: The organizers appreciated the originality of all the submissions and the hard work of all the participating teams. The participants' efforts in stretching the current limits of the discipline were remarkable, both methodologically-by integrating advanced video processing and computer vision techniques-and application wise. It also was highly appreciated that new challenging fields of study, such as population density estimation (proposed by the third place team [START_REF] Kelbe | Sensor-agnostic photogrammetric image registration with applications to population modeling[END_REF]) and human activity modeling (from the winning team), were addressed within the contest.

The data will remain downloadable for free from the IEEE GRSS website. 9 Ground references were made available by the two winning teams (see Sections IV and V for the URLs). We do hope that these data will serve to push the boundaries of remote 9 http://www.grss-ieee.org/community/technical-committees/data-fusion, under the 'Past Contests' tab.

sensing data fusion even further and make video form space a new valuable data source for remote sensing applications.

Fig. 1 .

 1 Fig. 1. Ground coverage of the Deimos-2 multispectral (full frame) and panchromatic images (yellow frame) and of the Iris video (red frame).

Fig. 2 .

 2 Fig.2. Level 1B (i.e., calibrated and radiometrically corrected, but not orthorectified) data for both dates and both the panchromatic (Pan) and multispectral (MS) modes.

Fig. 3 .

 3 Fig.3. Level 1C (i.e., calibrated, radiometrically corrected and manually orthorectified, and resampled to a map grid) data for both dates and both the panchromatic (Pan) and multispectral (MS) modes.

Fig. 4 .

 4 Fig. 4. Two frames of the UHD video acquired by the Iris camera on-board the International Space Station. (a) Iris video frame at t = 1 s, (b) Iris video frame at t = 30 s.

Fig. 5 .

 5 Fig.5. Problems studied (each paper can cover more than one; left panel) and main methodological approaches adopted in the submitted papers (right panel).

Fig. 6 .

 6 Fig.6. Use of the data sources in the 24 submitted manuscripts.

Fig. 9 .

 9 Fig. 9. Spatial analysis: Example of region refinement for resolution improvement. The first and second row show the coarse and refined probability maps of different scenes, respectively (from left to right: urban, vegetation, harbor, and river).

Fig. 10 .

 10 Fig. 10. Illustration of pooling (left) and unpooling (right).

Fig. 11 .

 11 Fig. 11. Visualization of (a) all tracklets produced by the KLT keypoint tracker and (b) positive tracklets (from [39]). and the duration time of the tracker. Fig. 11(b) shows the result of tracklet denoising.2) Activity Analysis and Traffic Density Estimation: For the remaining tracklets, the angle of motion θ = arctan(ζ/ε) was calculated according to the displacement vector (ε, ζ) of the tracklet, and 2-D information was transformed into the (ε, ζ, θ) space, in which k-means was used for clustering. The output could be interpreted as activity information of every cluster. The traffic density was estimated on the superpixels superimposed on the scene. For every superpixel, the number of positive tracklets was calculated and the density was estimated for the whole scene at every frame. Finally, the normalized density estimation map was obtained by dividing the number of positive tracklets in each superpixel by the maximum number among all superpixels.

Algorithm 1 :

 1 General formulation of the algorithm proposed by the second ranked team. Require: param: parameters of the algorithm Ensure: Registration, Semantic Segmentation, and Detection labels. 1: i ← param.mingridlevel 2: while i++ < param.maxgridlevel do 3) ← projection function (see [40]) 7:

Fig. 15 .

 15 Fig. 15. Semantic segmentation of the multitemporal multisensor data: (left) Deimos-2 March 2015 and Deimos-2 May 2015 images (D1-D2, from [40]), (right) Iris video sequence (first frame) and Deimos-2 May 2015 image (V-D2).

Fig. 16 .

 16 Fig. 16. Change Detection from multitemporal multisensor data (from [40]): (left) Deimos-2 March 2015 and Deimos-2 May 2015 images (D1-D2), (right) Iris video sequence (first frame) and Deimos-2 May 2015 image (V-D2).

1) First Place: Lichao Mou and Xiaoxiang Zhu from the

	German Aerospace Center (DLR), for the paper entitled
	Spatiotemporal scene interpretation of space videos via
	deep neural network and tracklet analysis [39].
	2) Second Place: Maria Vakalopoulou, Christos Platias,
	Maria Papadomanolaki, Nikos Paragios, and Konstanti-
	nos Karantzalos from the National Technical University of
	Athens, Greece, and the Ecole Centrale de Paris, France,
	for their paper entitled Simultaneous registration, segmen-
	tation and change detection from multi-sensor, multitem-
	poral satellite image pairs [40].
	3) Third Place: Dave Kelbe, Devin White, Andrew Hardin,
	Jessica Moehl, and Melanie Phillips, from the Oak Ridge
	National Laboratory, USA, for their paper entitled Sensor-
	agnostic photogrammetric image registration with appli-

cations to population modeling

[START_REF] Kelbe | Sensor-agnostic photogrammetric image registration with applications to population modeling[END_REF]

. 4) Fourth Place: Zuming Huang, Guangliang Cheng, Hongzhen Wang, Haichang Li, Limin Shi, and Chunhong Pan from the Institute of Automation of the Chinese Academy of Sciences, China, for their paper entitled Building extraction from multi-source remote sensing images via deep deconvolution neural networks

[START_REF] Huang | Building extraction from multi-source remote sensing images via deep deconvolution neural networks[END_REF]

.

TABLE I QUANTITATIVE

 I EVALUATION OF THE REGISTRATION RESULTS (MEAN DISPLACEMENT ERRORS IN PIXELS; FROM [40])

			Method of [60]		Proposed
		D1 to G	D2 to G	V to G	V Frame to Frame	D1 to D2	V to D2
	Dx	1.09	1.22	0.93	0.84	1.12	1.04
	Dy	1.62	1.49	1.73	0.92	1.59	1.61
	DS	1.95	1.93	1.97	1.24	1.94	1.92

TABLE II QUANTITATIVE

 II EVALUATION OF THE SEMANTIC SEGMENTATION RESULTS (PER CLASS RESULTS FROM [40])

		Deimos-2 March [D1]-Deimos-2 May [D2]
		Completeness	Correctness	Overall Quality
	Ship/vessel	81.4%	78.0%	66.2%
	Vegetation	83.9%	88.3%	75.6%
	Buildings	68.9%	77.4%	57.4%
	OA	85.3%
	AA	81.2%
		Iris July [V]-Deimos-2 May [D2]
		Completeness	Correctness	Overall Quality
	Ship/vessel	79.0%	77.9%	65.6%
	Vegetation	82.5%	86.2%	72.8%
	Buildings	78.8%	72.2%	60.5%
	OA	80.2%
	AA	78.8%
	mean displacement errors for both axis (Dx and Dy) and the dis-
	tance (DS) in pixels are presented. It should be mentioned that
	the registration process did not manage to address the largest
	relief displacements of the tallest buildings/skyscrapers of this
	part of Vancouver, and these errors hindered the results of both
	semantic segmentation and change detection. All other building
	rooftops, roads, and terrain classes were registered with subpixel
	accuracy. A quantitative comparison is also provided in Table I:

TABLE III QUANTITATIVE

 III

		Completeness	Correctness	Overall Quality
	Ship/vessel	68.6%	66.7%	66.2%
	Vegetation	88.2%	82.3%	75.6%
	Buildings	69.2%	67.4%	57.4%
	Overall (proposed)	75.3%	72.1%	66.4%
	OA, method of [6]	70.3%	67.1%	52.3%
		Iris July [V]-Deimos-2 May [D2]
	Ship/vessel	70.6%	69.5%	53.8%
	Vegetation	81.1%	79.6%	67.2%
	Buildings	71.3%	65.6%	51.9%
	Overall (proposed)	74.3%	71.6%	57.6%
	OA, method of [6]	70.4%	66.3%	51.8%

EVALUATION OF THE CHANGE DETECTION RESULTS (PER CLASS RESULTS FROM [40]) Deimos-2 March [D1]-Deimos-2 May [D2]

https://www.urthecast.com/data#iris

http://www.grss-ieee.org/community/technical-committees/data-fusion/

A superpixel can be defined as a set of locally connected similar pixels that preserve detailed edge structures for a fine segmentation. The code used for generating superpixels is available at http://coewww.rutgers. edu/riul/research/code/EDISON/[START_REF] Comaniciu | Mean shift: A robust approach toward feature space analysis[END_REF]. In this paper, spatial and range bandwidth parameters were set to 7 and 6.5, respectively.

Throughout Section V, we will use the word classification for the scores issued from this independent classification and semantic segmentation for those obtained by the proposed energy minimization method.

ACKNOWLEDGMENT

The authors would like to express their greatest appreciation to Deimos Imaging and Urthecast, for acquiring and providing the data used in the competition and for indispensable contribution to the organization of the Contest, and the IEEE GRSS for continuously supporting the annual Data Fusion Contest through funding and resources.

The work of M. Vakalopoulou was supported by the "ELKE" Ph.D. Scholarship from the National Technical University of Athens. The work of L. Mou and X. Zhu was supported in part by the European Research Council (ERC) under the European Union Horizon 2020 research and innovation programme under Grant ERC-2016-StG-714087 and in part by Helmholtz Association under the framework of the Young Investigators Group "SiPEO" (VH-NG-1018, www.sipeo.bgu.tum.de). The work of D. Tuia was supported by the Swiss National Science Foundation for financial support through Grant PP00P2-150593. L. Mou.)