Anne-Flore Perrin
email: anne-flore.perrin@irisa.fr

Lu Zhang

Olivier Le Meur
email: olemeur@irisa.fr

How well current saliency prediction models perform on UAVs videos?

Keywords: Benchmark, Salience, Dynamic saliency models, Unmanned Aerial Vehicles (UAV), Videos

It is exciting to witness the fast development of Unmanned Aerial Vehicle (UAV) imaging which opens the door to many new applications. In view of developing rich and efficient services, we wondered which strategy should be adopted to predict salience in UAV videos. To that end, we introduce here a benchmark of off-the-shelf state-of-theart models for saliency prediction. This benchmark studies comprehensively two challenging aspects related to salience, namely the peculiar characteristics of UAV contents and the temporal dimension of videos. This paper enables to identify the strengths and weaknesses of current static, dynamic, supervised and unsupervised models for drone videos. Eventually, we highlight several strategies for the development of visual attention in UAV videos.

Introduction

Unmanned Aerial Vehicles (UAVs) are propitious for a broad range of applications. Drone racing tournaments, new autonomous delivery services, wildfire detection, and private house surveillance is a sample of these highly promising services. However, we reach the point that service improvement lies in technology-, imagery-and context-based solutions. Indeed, UAV imagery presents specificities, detailed later on, that are worth exploiting, especially when using videos.

Services mentioned above heavily rely on object detection, segmentation and even more on saliency detection. Salience expresses the extent of importance of an area, approximating perceptual processes of visual attention. It is represented as the prediction of the fixation probability density of a multimedia content [START_REF] Kummerer | Saliency benchmarking made easy: Separating models, maps and metrics[END_REF]. Its prediction is essential in diverse fields, from content-aware re-targeting and compression to advertisement analyses. As UAV applications may further benefit from visual attention theory, we investigated whether state-of-the-art saliency models are suitable for predicting visual attention in UAV videos.

Several works already applied saliency prediction models in UAV services. An automatic salient object detection is implemented in [START_REF] Sokalski | Automatic salient object detection in uav imagery[END_REF] through mean-shift segmentation and edge-detection operators. In [START_REF] Trinh | Efficient uav video event summarization[END_REF], the combination of frame alignment, forward/backward difference and blob detection creates a UAV video event summarization, tracking salient objects. Deep architectures are used for real-time autonomous indoor navigation [START_REF] Kim | Deep neural network for real-time autonomous indoor navigation[END_REF][START_REF] Guo | How will humans cut through automated vehicle platoons in mixed traffic environments? a simulation study of drivers' gaze behaviors based on the dynamic areas-of-interest[END_REF] and for monitoring wildfire [START_REF] Zhao | Saliency detection and deep learning-based wildfire identification in uav imagery[END_REF], among others. Such approaches build upon basic deep-learning models, namely AlexNet [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF] and CaffeNet [START_REF] Jia | Caffe: Convolutional architecture for fast feature embedding[END_REF].

Let us stress that only very few works address the temporal dimension in traditional and UAV videos. Methods that tackle the temporal dimension comprise hand-crafted motion features [START_REF] Guo | A novel multiresolution spatiotemporal saliency detection model and its applications in image and video compression[END_REF][START_REF] Rudoy | Learning video saliency from human gaze using candidate selection[END_REF], network architecture fed with optical flow [START_REF] Bak | Spatio-temporal saliency networks for dynamic saliency prediction[END_REF], possibly in a two-layer fashion [START_REF] Bak | Spatio-temporal saliency networks for dynamic saliency prediction[END_REF][START_REF] Dutt Jain | Fusionseg: Learning to combine motion and appearance for fully automatic segmentation of generic objects in videos[END_REF], or Long Short-Term Memory (LSTM) architectures [START_REF] Bazzani | Recurrent mixture density network for spatiotemporal visual attention[END_REF][START_REF] Jiang | Deepvs: A deep learning based video saliency prediction approach[END_REF][START_REF] Wang | A deep-learning based feature hybrid framework for spatiotemporal saliency detection inside videos[END_REF][START_REF] Li | Flow guided recurrent neural encoder for video salient object detection[END_REF] to benefit from their memory functionality.

Regarding conventional imaging, elegant, elaborate and efficient solutions have been developed for saliency prediction. However, we wonder if typical schemes keep their promises when dealing with UAV video characteristics. Indeed, this new-born imagery field is distinct from typical imagery in many aspects, including the bird-point-of-view which modifies the semantic and size of objects [START_REF] Krassanakis | Monitoring human visual behavior during the observation of unmanned aerial vehicles (uavs) videos[END_REF], the loss of pictorial depth cues [START_REF] Howard | Depth perception[END_REF] i.e. the lack of horizontal line [START_REF] Foulsham | Turning the world around: Patterns in saccade direction vary with picture orientation[END_REF], and the presence of camera movements [START_REF] Krassanakis | Monitoring human visual behavior during the observation of unmanned aerial vehicles (uavs) videos[END_REF].

In this paper, we provide answers to the above questions by first introducing the saliency models under review in Section 2. In Section 3, we go through the details of the benchmark. Results are discussed in Section 4. Finally, key takeaways of this benchmark are given in Section 5.

Saliency models

Tremendously different approaches have been developed throughout the years to reach high accuracy prediction of visual attention. We propose here a taxonomy to categorize salience solutions.

Taxonomy

With the introduction of temporal dimension in videos, Borji [START_REF] Borji | Saliency prediction in the deep learning era: An empirical investigation[END_REF] proposed to classify models according to the use or not of motion features and deep learning architectures. We refine this taxonomy by considering all types of supervision and not only differentiating deep-learning-based models. Accordingly, we define five model categories:

-Static Unsupervised (SU): Itti [START_REF] Itti | A model of saliency-based visual attention for rapid scene analysis[END_REF], LeMeur [START_REF] Le Meur | Predicting visual fixations on video based on low-level visual features[END_REF], GBVS [START_REF] Harel | Graph-based visual saliency[END_REF], SUN [START_REF] Zhang | Sun: A bayesian framework for saliency using natural statistics[END_REF], Judd [START_REF] Judd | Learning to predict where humans look[END_REF], Hou [START_REF] Hou | Image signature: Highlighting sparse salient regions[END_REF], RARE2012 [START_REF] Riche | Rare2012: A multi-scale rarity-based saliency detection with its comparative statistical analysis[END_REF], BMS [START_REF] Zhang | Exploiting surroundedness for saliency detection: a boolean map approach[END_REF], -Static Deep learning (SD): Salicon [START_REF] Huang | Salicon: Reducing the semantic gap in saliency prediction by adapting deep neural networks[END_REF], DeepNet [START_REF] Pan | Shallow and deep convolutional networks for saliency prediction[END_REF], ML-Net [START_REF] Cornia | A deep multi-level network for saliency prediction[END_REF], SalGAN [START_REF] Pan | Salgan: Visual saliency prediction with generative adversarial networks[END_REF], -Dynamic Unsupervised (DU): Fang [START_REF] Fang | Video saliency incorporating spatiotemporal cues and uncertainty weighting[END_REF], OBDL [START_REF] Hossein Khatoonabadi | How many bits does it take for a stimulus to be salient?[END_REF],

-Dynamic Machine learning (DM): PQFT [START_REF] Guo | A novel multiresolution spatiotemporal saliency detection model and its applications in image and video compression[END_REF], Rudoy [START_REF] Rudoy | Learning video saliency from human gaze using candidate selection[END_REF], -Dynamic Deep learning (DD): DeepVS [START_REF] Jiang | Deepvs: A deep learning based video saliency prediction approach[END_REF], ACL-Net [START_REF] Wang | A deep-learning based feature hybrid framework for spatiotemporal saliency detection inside videos[END_REF], STSconvNet [START_REF] Bak | Spatio-temporal saliency networks for dynamic saliency prediction[END_REF], FGRNE [START_REF] Li | Flow guided recurrent neural encoder for video salient object detection[END_REF].

Deep learning models could be further classified depending on their loss function(s) or architecture attributes (e.g. Convolutional Neural Networks (CNN) design, multi-level features, encoder-decoder system, LSTM(s)).

Models under study

Models under study were selected based on an exhaustive study of the state of the art, exhibiting most used models for cross-model comparisons. Accordingly, we evaluate a total of 11 off-the-shelf models, including 8 static models and 3 dynamic models. They are briefly described hereinafter, following our taxonomy categorization.

Static unsupervised models -Itti: One of the first and most used static models is referred to as Itti's model. In [START_REF] Itti | A model of saliency-based visual attention for rapid scene analysis[END_REF], authors got inspired from primates' visual systems, especially the center-surround analytic behavior of retina and cortical lateral inhibition mechanisms. Accordingly, they investigated three modalities, namely color, intensity, and orientation, through features on variable scales. -Graph-Based Visual Saliency (GBVS) [START_REF] Harel | Graph-based visual saliency[END_REF]: GBVS is a bottom-up visual saliency model, exploiting the real-time ability of graph algorithms. It implements a feature extraction using biologically inspired filters, an activation measure by subtracting features at different scales and finally a normalization, applied based on local maxima, convolution, and non-linear weighting. -Saliency Using Natural statistics (SUN) [START_REF] Zhang | Sun: A bayesian framework for saliency using natural statistics[END_REF]: SUN mimics the visual system behavior to find potential targets by means of a Bayesian probabilistic framework. It includes computation of self-information, likelihood and location prior to reach an estimation of pointwise mutual information, which expresses the overall salience. -SIM [START_REF] Murray | Saliency estimation using a non-parametric low-level vision model[END_REF] : First, a convolution of the image with a bank of filters using a multiresolution wavelet transform yields to a scale-space decomposition. Then, a Gaussian Mixture Model (GMM), trained on eye-fixation data, simulates the inhibition mechanisms in visual cortex cells. Finally, multi-scale information is integrated using an inverse wavelet transform. -RARE2012 [START_REF] Riche | Rare2012: A multi-scale rarity-based saliency detection with its comparative statistical analysis[END_REF]: The key of RARE2012 is its multi-scale rarity mechanism applied on low-level colour and medium-level orientation features extracted beforehand. To emphasise both local contrasts and global rare regions, highly salient regions have the lowest occurrence probabilities of pixel at all scales. -Hou [START_REF] Hou | Image signature: Highlighting sparse salient regions[END_REF]: Assuming that sparsity in spatial and frequential domains discriminate the foreground from the background, the foreground is represented by the sign of the discrete cosine transform of the signal. The inverse transform of the three color channels foreground are squared, summed up and smoothed to get the saliency map.

-Boolean Map Saliency (BMS) [START_REF] Zhang | Exploiting surroundedness for saliency detection: a boolean map approach[END_REF]: First, a set of Boolean maps is created by applying uniformly distributed thresholds on the image color channels. Attention maps are derived from this set through morphological dilatation and Gaussian blurring. The final saliency map is the average attention map.

Dynamic machine learning model

-PQFT [START_REF] Guo | A novel multiresolution spatiotemporal saliency detection model and its applications in image and video compression[END_REF]: This model extends phase spectrum of Fourier transform models to a temporal approach. Chroma, luma and motion information form the quaternion of a frame. The quaternion sympletic form is converted to the frequency domain. Only the phase, which represents local information of the signal, is converted back to the spatial domain. Such information, after Gaussian filtering, produce the saliency map.

Static deep learning model

-ML-Net [START_REF] Cornia | A deep multi-level network for saliency prediction[END_REF]: This multi-level Fully Convolutional neural Network (FCN) combines low-and high-level features to predict saliency. However, the true ingenuity of this model lies in the defined loss function, which penalises more errors occurring on salient pixels.

Dynamic deep learning models

-DeepVS: Jiang et al. [START_REF] Jiang | Deepvs: A deep learning based video saliency prediction approach[END_REF] proposed an elegant architecture including two models applied sequentially. The Object-to-Motion CNN (OM-CNN) is encoding objectedness features through a complex and comprehensive network combining hierarchical spatial (coarse) and temporal saliency maps. To further compute dynamic salience of videos, a Saliency-Structured Convolutional LSTM (SS-ConvLSTM), is added. Its two successive LSTM networks leverage both short and long term correlations. -ACL-Net: Wang et al. [START_REF] Wang | A deep-learning based feature hybrid framework for spatiotemporal saliency detection inside videos[END_REF] implemented a CNN-LSTM architecture for video saliency. An attention module supervises the CNN, forcing it to learn static features and ensuring the LSTM to deal with dynamic characteristics.

3 Benchmark design

Baselines

To include comparison points for the above models, we include six baselines. First, the average saliency map over all observers and all sequences (OHM) brings out an overall representation of Human Mean (HM). Also, average saliency maps over observers for a sequence (SHM), and over all observers for all sequences but the one under study (abSHM) examine content-dependencies. Then, the Center Bias (CB) map is a centered isotropic Gaussian stretched to video frame aspect ratio [START_REF] Bylinskii | What do different evaluation metrics tell us about saliency models[END_REF]. Last, we add two chance representations, inspiring from [START_REF] Ninassi | Does where you gaze on an image affect your perception of quality? applying visual attention to image quality metric[END_REF]. We split SHMs into 16 blocks, which are redistributed in the final map following two 4x4 magic squares. Despite the lack of salience information, these shuffled maps cover a similar ratio and dynamic range than that of the ground truth.

Dataset

Specific datasets are required to conduct saliency performance analyses. Indeed, eye movements of human beings or any substitute, collected during stimuli visualization, establish the Ground Truth (GT) essential to assess the validity of predicted saliency maps.

EyeTrackUAV dataset [START_REF] Krassanakis | Monitoring human visual behavior during the observation of unmanned aerial vehicles (uavs) videos[END_REF] is the only UAV salience dataset available, to the best of our knowledge. It includes 19 different videos extracted from the UAV123 database [START_REF] Mueller | A benchmark and simulator for uav tracking[END_REF]. These video sequences were captured from a fully stabilized and manually controlled off-the-shelf professional-grade UAV (DJI S1000) flying at low-altitudes (varying between 5-25 meters). Criteria for content selection were the diversity of environment, distance and angle to the scene, size of the principal object and presence of sky. Authors of [START_REF] Krassanakis | Monitoring human visual behavior during the observation of unmanned aerial vehicles (uavs) videos[END_REF] collected highly precise binocular gaze data (1000 Hz) from 14 subjects in free viewing conditions. Stimuli were videos with a resolution of 720p, 30 fps which represent overall 26599 frames and 887 seconds.

From data, two maps have been computed for evaluation purposes. Saliency maps were inferred directly from raw data gathered by the eye tracker. Binocular gaze recordings were averaged over all observers. These maps were then filtered using a Gaussian kernel and normalized [START_REF] Le Meur | Methods for comparing scanpaths and saliency maps: strengths and weaknesses[END_REF]. Fixations were retrieved through a Dispersion-Threshold Identification (I-DT) algorithm [START_REF] Krassanakis | Eyemmv toolbox: An eye movement post-analysis tool based on a two-step spatial dispersion threshold for fixation identification[END_REF]. Fixation points, in this spatiotemporal detection, are assumed to cluster together.

Metrics

To carry out the evaluation, we use six quality metrics included in the MIT benchmark [START_REF] Bylinskii | Mit saliency benchmark[END_REF][START_REF] Le Meur | Methods for comparing scanpaths and saliency maps: strengths and weaknesses[END_REF]: Correlation Coefficient (CC) (∈ [-1, 1]), Similarity (SIM) (∈ [0, 1]) the intersection between histograms of saliency, Area Under the Curve (AUC) Judd and Borji (∈ [0, 1]), Normalized Scanpath Saliency (NSS) (∈ ]-∞, +∞[), and Kullback Leibler divergence (KL) (∈ [0, +∞[). Details of these metrics can be found in [START_REF] Le Meur | Methods for comparing scanpaths and saliency maps: strengths and weaknesses[END_REF][START_REF] Bylinskii | Mit saliency benchmark[END_REF]. Also, a very interesting comparison of metrics behavior is presented in [START_REF] Bylinskii | What do different evaluation metrics tell us about saliency models[END_REF].

Results and discussion

Analyses

Static, dynamic, non-supervised, machine-learning and deep-learning models are compared here. We discuss our results qualitatively, overall and on a frame-byframe basis.

Qualitative analysis

The qualitative verification is done on the sequence Person20, which illustrates models efficiency in a typical scenario of UAV applications. By observing Figure 1, one can note that deep learning and dynamic models detected less salient areas than unsupervised static models. This is in line with the assumed sparsity of video saliency.

(h) BMS (i) RARE2012 (j) SIM (k) Hou (l) SUN (m) GBVS (n) Itti (o) OHM (p) SHM (q) abSHM (r) CB (s) Shuffle 1 (t) Shuffle 2
Fig. 1: Saliency maps of the 11 models, the ground truth, together with fixations and source of frame 1035 in person20.

Overall analysis

Overall results are presented in Table 1. The analysis relies on mean and standard deviation of raw results while Analysis of Variance (ANOVA) and multicomparison studies have been conducted to verify the statistical significance of observed differences. We draw several observations.

Deep learning models present the best scores, especially for SIM for which differences with other models are statistically significant. This result was ex-pected as the same behavior is observed since the development of deep learning techniques in saliency prediction [START_REF] Vig | Large-scale optimization of hierarchical features for saliency prediction in natural images[END_REF][START_REF] Huang | Salicon: Reducing the semantic gap in saliency prediction by adapting deep neural networks[END_REF][START_REF] Pan | Shallow and deep convolutional networks for saliency prediction[END_REF][START_REF] Cornia | A deep multi-level network for saliency prediction[END_REF][START_REF] Pan | Salgan: Visual saliency prediction with generative adversarial networks[END_REF].

Within static models, BMS model, followed closely by GBVS and Itti models, happens to be the most efficient static unsupervised model. SIM and SUN models are significantly less performing for CC, SIM, NSS and KL.

Regarding dynamic models, ACL-Net reaches the first and second position for all but AUC-Borji metric. DeepVS is also ranked high up behind the two other deep learning schemes. Contrarily, PQFT achieves statistically significantly the lowest performance of all models for all metrics and achieves similar performances as baselines. We thus recommend avoiding frequency-based solutions, even though this result must be replicated with other models and datasets for validation.

Baselines results show the potential of all models but PQFT -and SIM and SUN models with respect to KL and SIM -as they are proven statistically different from baselines. Despite being the most performing baseline, Center Bias results are quite low when compared to results on traditional contents. This raises suspicions about the suitability of this bias for UAV videos. We have indecisive results regarding overall image-dependent properties. Indeed, SHM is not statistically different from other HM baselines and is alternatively better (CC, SIM, KL) or worse (AUC-based and NSS).

From all the above observations, UAV saliency models can be built upon static and dynamic deep learning models. Moreover, performances of today's solutions are fair and could potentially reach better scores through fine-tuning or training on a UAV-dedicated database.

Sequence-based Analysis

A challenge in video quality assessment is to design metrics that are representative of the entire sequence, considering the disparity in quality scores of video frames. Saliency prediction in videos presents the same challenge. In this benchmark, we tackled this issue by reporting average and according standard deviation values of metrics over the entire sequence. To dig further in UAV video saliency, let us discuss the accuracy of predictions on a frame-by-frame basis. This in-depth investigation exhibits potential perspectives for metrics development in video saliency prediction.

In Figures 2a to 2c, a comparison of averaged scores for models sorted by categories, following the above-defined taxonomy, is provided in terms of CC scores for sequence bike3. Also, overall CC results of all models and average baseline are reported in Figure 2d for wakeboard10.

High disparity. Overall, the disparity of standard deviation and mean values spreads out over a wide range among sequence frames. Hence, giving one score per video for saliency prediction may be fallacious. Also, we claim that the lack of consideration of temporal continuity in performance estimation may prevent thorough assessments of predictions. Implementing temporal coherence metric for video saliency prediction may be a solution to reach highly accurate and constructive analyses. Deep learning models perform better. Even on a frame-by-frame basis and amongst sequences presenting statistically different results, we can observe in Figures 2a and2b the clear advantage of deep learning architectures over unsupervised and machine learning schemes, for both dynamic and static models.

Static vs dynamic deep learning models. There is no evidence of significant over-performance of static deep learning over dynamic deep learning approaches, as illustrated in Figures 2c. Actually, ANOVA has not rejected the null hypothesis when comparing the distributions of these two categories. We observe that results are content-dependent, hinting that the difference of performance between deep learning static and dynamic models lies in videos characteristics.

Event-related performance. We could not relate the varying features of EyeTrackUAV videos (overall angle and distance to the scene, environment, presence of sky, and object size) to models' efficiency. Our intuition is that model performances are event-related. We further studied our results to identify causes of abrupt changes in metrics along video frames.

It turns out that a decrease in performance may be caused by the entry of new objects of potential interest or environments variations. However, a point of interest arises in that a sudden increase in prediction accuracy often follows camera movements towards re-framing the content (change of camera angle, distance, speed or trajectory guided by an operator). Reframing reintroduces a center bias in the content. This explains the increase in prediction accuracy at the end of the scene wakeboard10 in Figure 2d. This shows the importance of understanding patterns and biases in the used imagery.

Challenges and open problems

Several different aspects of our study must be examined and discussed, from the limitations due to the used dataset, to the "normalization" of metrics for better comparison and analyses of saliency prediction strategies.

Dataset. EyeTrackUAV presents a clear lack of non-natural videos. But most importantly, this dataset is not sufficient to train a deep model, even with data augmentation. To address the need of dedicated UAV video saliency prediction, Fig. 3: Center bias for sequences bike3, boat6, building5 and car4.

one needs new larger databases with more videos, possibly with natural and unnatural images.

Likewise, the larger the population taking part in the gaze data gathering is, the more accurate should be the ground truth -on which rely performance analyses. Also, based on our results, event-related annotations, even though complicated to collect, can provide a valuable contribution.

Finally, this dataset has been designed to study bottom-up attention. However, multiple UAV applications, such as aerial surveillance, monitoring and observation with drone, require the study of top-down attention, which is taskrelated. Both bottom-up and top-down attention ground truth is necessary for the development of UAV services.

Center bias. An additional UAV feature is that the center bias [START_REF] Tseng | Quantifying center bias of observers in free viewing of dynamic natural scenes[END_REF] is less significant in EyeTrackUAV sequences than in traditional video sequences, as seen in Figure 3. It may explain, in part, the weak performance of existing saliency models. Further analyses are required to confirm this compelling finding.

Metrics. As mentioned earlier, temporal metrics must be developed to correctly assess video saliency predictions. Our main point is that the averaged prediction performance over video frames may not be representative of the prediction along scenes in videos. It might be profitable to mimic the continuum of our visual gaze deployment. Several strategies may be deployed accordingly, possibly inspiring from event summarization [START_REF] Trinh | Efficient uav video event summarization[END_REF].

With a slight shift of perspective, one may also consider to apply postprocessing to predicted saliency maps to reach comparable metric results. To reduce metric biases, Kummerer et al. [START_REF] Kummerer | Saliency benchmarking made easy: Separating models, maps and metrics[END_REF] proposed to turn a given fixation density into metric-specific saliency maps prior to metric computations.

Conclusion

UAV new services want to benefit from the substantial improvement in saliency prediction of this last decade. This benchmark reviews the ability of 11 state-ofthe-art off-the-shelf prediction schemes to identify the direction to take toward the development of UAV imagery-dedicated models.

We studied qualitatively models performance over the EyeTrackUAV dataset which includes 19 natural UAV videos together with precise eye tracking information, collected from 14 subjects. Predicted maps have been assessed using typical metrics, namely CC, SIM, AUC-Judd, AUC-Borji, NSS and KL. Overall scores of models are reported through mean and standard deviation over the entire sequence. ANOVA and multi-comparison analyses have been included to detect if models differences are statistically significant.

Several insights are provided among which are three key takeaway messages:

(1) In line with studies on traditional contents, static and dynamic deep learning models, trained on conventional contents, show the most promising results. This outcome is highly encouraging, especially with the potential of reaching higher performance through fine-tuning or training on larger UAV databases.

(2) Yet there are no significant difference between static and dynamic deep learning. This outcome has been shown to be content-dependent. Although video characteristics (angle, distance to the scene, environment type and object size) were not sufficient to explain our results, it seems that event-related annotations could help to efficiently learn saliency on such contents.

(3) The need of dedicated UAV-centric models is hereby made clear. We need to go deeper in content specificities and better encode them for future high-quality UAV-based services. Also, we found out that the center bias does not necessarily apply to UAV videos, which needs deeper exploration. A broad investigation of typical biases of attention and cognition could be carried out.

Different challenges have also been discussed, including the need to develop video-based metrics in view to further investigate video prediction performance as well as to better represent the quality of prediction along the entire sequence, or the necessity to create large datasets for UAV imagery saliency prediction.

  (a) Original (b) Fixations (c) GT (d) ACL-Net (e) DeepVS (f) MLNet (g) PQFT

Fig. 2 :

 2 Fig. 2: Temporal comparison of model categories in bike3 and results of all models in wakeboard10 for metric CC for all sequence frames. Colored envelop represents standard deviation associated with mean results of models.

Table 1 :

 1

	DD	ACL-Net [40], 2018 0.4516 0.3586 0.8199 DeepVS [19], 2018 0.3986 0.3204 0.8059	0.7717 1.9622 1.7803 0.7384 1.7904 1.9063
	SD	MLNet [6], 2016	0.4621 0.3149 0.8347 0.7866 2.1479 1.5857
	DM	PQFT [10], 2010	0.1367 0.1817 0.7054	0.5591 0.7364 2.3282
		BMS [41], 2016	0.3846 0.2482 0.8189	0.7855 1.8180 1.8053
		RARE2012 [34], 2013 0.3422 0.2566 0.7946	0.7582 1.5093 1.8240
		Hou [14], 2012	0.2811 0.2204 0.7565	0.7279 1.1529 1.9818
	SU	SIM [30], 2011	0.2213 0.1893 0.7511	0.7235 0.9899 2.1482
		SUN [42], 2008	0.2065 0.1913 0.7361	0.7059 0.9593 2.1702
		GBVS [12], 2007	0.3687 0.2431 0.8125	0.7854 1.5915 1.8115
		Itti [17], 1998	0.3687 0.2431 0.8125	0.7855 1.5915 1.8115
		OHM	0.1889 0.1883 0.7208	0.6934 0.6881 2.3526
		SHM	0.1917 0.1944 0.7051	0.6790 0.6880 2.2463
	Baselines	abSHM CB	0.1808 0.1868 0.7125 0.2055 0.1980 0.7493	0.6840 0.6544 2.3539 0.7172 0.7588 2.3776
		Shuffle 1	0.0282 0.1558 0.5724	0.5427 0.1315 2.4275
		Shuffle 2	-0.0085 0.1468 0.5323	0.5041 -0.0046 2.4745

Performance of saliency models in average over all videos and all frames. Best performances are in bold. (AUC-B=AUC-Borji; AUC-J=AUC-Judd) CC ↑ SIM ↑ AUC-J ↑ AUC-B ↑ NSS ↑ KL ↓

The presented work is funded by the ongoing research project ANR ASTRID DIS-SOCIE (Automated Detection of SaliencieS from Operators' Point of View and Intelligent Compression of DronE videos) referenced as ANR-17-ASTR-0009.


Christine Böhmer

Anick Abourachid

Philippe Wenger

Benjamin Fasquelle

Matthieu Furet

Christine Chevallereau

Damien Chablat

Combining precision and power to maximize performance: a case study of the woodpecker's neck

Keywords: neck, bird, tensegrity, bio-inspired, robotics

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The neck in birds has been characterized as a highly complex system that positions the head during all kinds of behavior [START_REF] Böhmer | What is a long neck? The effects of scaling relationships between skeletal dimensions and body size in birds[END_REF][START_REF] Zweers | Motor organization of the avian head-neck system[END_REF]. This involves a variety of highly demanding tasks such as feeding, manipulation, preening, sexual display, nest building, and combat behavior [START_REF] Böhmer | Correlated evolution of neck length and leg length in birds[END_REF]. This versatility is an opportunity for bioengineering, but prior to designing a technological model, the biological system must be understood; in particular regarding the form-function relationships. One behavioral pattern in birds that appears to be relatively constrained is pecking in woodpeckers. There is evidence that the drilling trajectory is essentially planar [START_REF] Spring | Climbing and pecking adaptations in some North American woodpeckers[END_REF]. From a morphofunctional point of view, woodpeckers are highly interesting because they are very lightweight (mean body mass = 300 g), but are capable to dig into dense tree trunks [START_REF] Puverel | This is my spot: characteristics of trees bearing Black Woodpecker cavities International Woodpecker Conference[END_REF]. They strike their beak against a tree repeatedly with high accelerations, high velocity and large force transmissions [START_REF] May | Woodpecker drilling behavior. An endorsement of the rotational theory of impact brain injury[END_REF][START_REF] Wang | Why do woodpeckers resist head impact injury: a biomechanical investigation[END_REF]. We took this bird species as a first example for a biologically inspired robot arm that combines precision and power to maximize performance. Adding biological knowledge into the design process requires a simplification of the complexity of the biological entity [START_REF] Whitesides | Bioinspiration: something for everyone[END_REF]. In this context, we (1) analyzed the neck musculature which supplies force for movement; (2) established a planar robotic model using several stacked tensegrity crossed bar mechanisms and (3) integrated all data into an actuated model of the bird neck.

Methods

Dissection

The neck musculature of one generalist bird (Corvus monedula) and three species of woodpeckers was comparatively investigated by quantitative dissection: the green woodpecker (Picus viridis), the great spotted woodpecker (Dendrocopos major) and the black woodpecker (Dryocopus martius). They differ in pecking performance [START_REF] Jenni | Das Skelettmuskelsystem des Halses von Buntspecht und Mittelspecht Dendrocopos major und medius[END_REF]) which allows identification of muscular adaptations specific for high power. After removing the skin and fascia, each muscle was identified and systematically detached from both the origin and insertion. Muscles were assigned to main functional categories indicating the direction and vertebral region of action. The craniocervical muscular system was illustrated in a schematic diagram.

Kinematic and dynamic modelling

The woodpecker neck is modelled in the sagittal plane with a serial stack of Snelson's X-shape tensegrity mechanisms [START_REF] Wenger | Kinetostatic analysis and solution classification of a planar tensegrity mechanism 7th International Workshop on Computational Kinematics[END_REF], which is a class 2 tensegrity mechanism composed of an antiparallelogram made of 4 rigid links. The sides of the mechanism are attached with two springs that ensure a stable rest position, and cables run through the springs. By modifying the cable tension with motors, different positions can be reached. The kinematic model links the output of the modelled stacked modules (i.e., the position and orientation of the tip of the beak) to the input of the robot (the orientation of each X-shape mechanism). The dynamic modelling takes into account the influence of the rigid links masses and inertia, spring stiffnesses and cable tensions and links the motion of the model (position, velocity, acceleration) to the forces needed to move it [START_REF] Furet | Optimal design of tensegrity mechanisms used in a bird neck model[END_REF].

Results and discussion

Muscular adaptations

The overall muscular system of the neck in the analyzed woodpeckers corresponds to an earlier anatomical description of Dendrocopus provided by [START_REF] Jenni | Das Skelettmuskelsystem des Halses von Buntspecht und Mittelspecht Dendrocopos major und medius[END_REF]. Our quantitative analysis revealed that the adaptation to pecking is particularly evident in the craniocervical and ventral muscle system. As opposed to generalist birds, the muscles of the craniocervical system in woodpeckers show adaptations specific for precision of force control. The highly specialized tendinous architecture of the ventral muscle system in woodpeckers enables forceful and rapid movements. Compared to the other species of woodpeckers, the tendons are most prominently developed in the black woodpecker. In summary, the anatomical arrangement of muscle and tendon in the neck of woodpeckers facilitates efficient pecking by combining precision and power.

Neck movement

Videos of alive woodpeckers do not allow us to precisely detect the neck movement because of the feathers. As a consequence, a film of the corpse of a green woodpecker without its feathers and moved by hand was recorded and used to capture a neck movement, with the aim of reproducing it with our robot in simulation.

To convert the video into a robot movement, the frames are analyzed independently. For each of them, the middle line of the neck is detected, then the mechanisms of the robots are placed one by one from the bottom of the neck. The obtained positions give us the trajectory.

Simulation

The manipulator is made of 11 identical mechanisms that is the number of vertebrae in the green woodpecker. The dimensions were chosen to have a manipulator with the same proportions than the woodpecker neck. Simulations are performed with Matlab Simulink and Simscape Multibody. The use of the dynamic model in the control law allows obtaining an accurate tracking with a full-actuated robot. In this case, the forces applied have similar shapes for each mechanism, which is promising for under-actuation. Simple under-actuation strategies also succeed at tracking the movement with a good precision, and the neck musculature of the woodpeckers will be used to create related under-actuation strategies.

Conclusions

The aim of the present work was to extract the anatomical parameters that are essential to perform precise and powerful movements. The conceptualization of the extremely complex musculature of the neck in birds revealed such a specialization. Eventually, the prototype of the woodpecker's neck will serve as a starting point for future invention and application such as industrial softrobots that execute their task in collaboration with humans in an efficient and safe way.

Acknowledgements

We thank the Chene Association and the Goupil Association for providing bird specimens for dissection. This study was financially supported by the Agence National de la Recherche (ANR): Project ID #ANR-16-CE33-0025 (AVINECK), Project coordinator: A. Abourachid.


Agniva Sengupta

Alexandre Krupa

Eric Marchand

Tracking of Non-Rigid Objects using RGB-D Camera

A method to accurately track deformable objects using a RGB-D camera with the help of a coarse object model is presented in this paper. The deformation model is based on corotational FEM formulation. The physical model of the object does not need to be exact, nor do we require the precise physical properties for accurately tracking the object. The position of the vertices of the surface mesh of the tracked object is deformed using a set of virtual forces. A point-to-plane distance based geometric error between the pointcloud and the mesh is minimized with respect to these virtual forces. The point of application of force is determined by analysis of the error obtained from rigid tracking, which is done in parallel with the non-rigid tracking. This architecture also enables the overall system to be realtime. The proposed approach is evaluated on a synthetic data with ground-truth for deformation at every frame, as well as on real data.

I. INTRODUCTION

Interaction with soft objects using machines remain a challenging problem in the field of robotics. To precisely interact with deforming objects using any robotic manipulator, it is necessary to track the surface of the object with high accuracy. Non-rigid object tracking also finds implementation in the field of augmented/mixed reality systems. This is an open area of research and the available literature that aims to tackle this problem is not very extensive.

Real-time tracking of the surface of non-rigid objects using RGB-D camera has become possible only very recently. However, most of the proposed approaches for non-rigid object tracking either requires a detailed model or considers tracking and reconstruction as intertwined problems. The tracking and reconstruction algorithm usually generates impressive visual results, but the accuracy of the final reconstruction does not necessarily imply accurate frame-to-frame tracking of the non-rigid object. Moreover, for model-based tracking, a detailed model of the physical properties of the object may not always be available in practice.

In this paper, we propose a method to handle the problem of accurate tracking of the surface of non-rigid objects undergoing deformation. A commodity-level RGB-D camera is used for sensing. It is assumed that we know the visualsurface model of the object (which can be a CAD model), but this model does not have to be precise.

A. Related Work

We propose to track the object using its 3D model, which consists of a mathematical representation of the surface of the object in 3D. To track the complete surface of deforming objects using depth (or monocular) camera when its 3D model is available, there are three components which are generally used: a) an algorithm to track each face of the model w.r.t the observed 3D depth, possibly along with the corresponding image intensity values of the pixels, b) a mechanism to impart coherence to the motion of these individual faces, which is done either with the help of a geometric shape preserving function (also called a regularizer), such as [START_REF] Alexa | As-rigid-as-possible shape interpolation[END_REF], [START_REF] Chui | A new algorithm for non-rigid point matching[END_REF] or [START_REF] Sumner | Embedded deformation for shape manipulation[END_REF], or with some assumptions about the underlying physical properties of the material of the object, such as [START_REF] Royer | Real-time tracking of deformable target in 3d ultrasound images[END_REF] or [START_REF] Petit | Tracking fractures of deformable objects in real-time with an rgb-d sensor[END_REF], and c) some strategy to combine the tracking of the surface with the regularizer or the physical model of the object (e.g: 'deterministic annealing' in case of [START_REF] Chui | A new algorithm for non-rigid point matching[END_REF]). Some of the popular methods for geometrically regularizing the deformations includes as-rigid-as-possible (ARAP) regularizer [START_REF] Sorkine | As-rigid-as-possible surface modeling[END_REF], embedded deformation graphs [START_REF] Li | Global correspondence optimization for non-rigid registration of depth scans[END_REF] and thin plate spline model [START_REF] Lee | Intensity-based visual servoing for non-rigid motion compensation of soft tissue structures due to physiological motion using 4d ultrasound[END_REF]. The physics based framework for deformation modelling includes mass-spring-damper systems [START_REF] Royer | Real-time target tracking of soft tissues in 3d ultrasound images based on robust visual information and mechanical simulation[END_REF], finite element model (FEM) [START_REF] Petit | Real-time tracking of 3d elastic objects with an rgb-d sensor[END_REF] and kinematic chain [START_REF] Schmidt | Dart: Dense articulated real-time tracking[END_REF]. For model based tracking of non-rigid objects, the representation of the model is usually done by a set of triangular faces representing the surface of the object. However, tetrahedral, volumetric model is preferred [START_REF] Paloc | Online multiresolution volumetric mass spring model for real time soft tissue deformation[END_REF] for tracking using mass-spring-damper systems, while the FEM based systems available in the literature typically preferred to use a combination of surface model and volumetric model. Haouchine et. al. [START_REF] Haouchine | Impact of soft tissue heterogeneity on augmented reality for liver surgery[END_REF] proposed a linear tetrahedral corotational FEM based model for tracking large deformations. [START_REF] Petit | Capturing deformations of interacting non-rigid objects using rgb-d data[END_REF] uses co-rotation FEM based model to track deformation by estimating the direction and magnitude of elastic force acting on the object using depth information. As summarized in [START_REF] Nadon | Multi-modal sensing and robotic manipulation of non-rigid objects: A survey[END_REF], FEM based approaches suffer from excessive dependency on the availability of the accurate physical properties of the object being tracked. The approach we propose in this paper removes this dependency using a novel, closed loop minimization technique.

The recent advancements in non-rigid object tracking should also include a survey of tracking and reconstruction algorithms. [START_REF] Newcombe | Dynamicfusion: Reconstruction and tracking of non-rigid scenes in real-time[END_REF] uses ARAP to regularize a warp field from the current state of deformation to a canonical model. The reconstructed model is visually coherent and accurate. [START_REF] Innmann | Volumedeform: Real-time volumetric non-rigid reconstruction[END_REF] improves upon this work and proposes to enhance the accuracy by introducing SIFT based sparse correspondence. [START_REF] Slavcheva | Killingfusion: Nonrigid 3d reconstruction without correspondences[END_REF] achieves the same objective, but uses Approximately Killing Vector Field (AKVF) to regularize the deformation. However, these approaches do not evaluate the per-frame tracking accuracy. The use of Truncated Signed Distance Field (TSDF) makes it easier to reconstruct the model, despite minor inaccuracies in tracking.

B. Contribution

We propose a method for tracking the entire visible surface of a deforming, non-rigid object in 3D, such that it:

• requires only a very approximate estimate of the physical properties (Young's modulus and Poisson ratio) of the object. • performs accurate frame-to-frame tracking of deformation of the non-rigid object.

• has been validated on a simulated data with ground truth, as well as real data. The method that we present in this paper is closely related to [START_REF] Petit | Capturing deformations of interacting non-rigid objects using rgb-d data[END_REF]. We use co-rotational FEM to simulate the physics behind the deformation of non-rigid objects and we use a tetrahedral mesh as the physical model of the object. But, unlike the approaches in the literature, we extend the method to perform a closed-loop optimization based on a point-toplane geometric error. This additional step of minimization effectively reduces the dependency of the system on the accuracy of the physical properties of matter used for the finite element modelling. Moreover, we dissociate the rigid and non-rigid tracking into two parallel modules. The rigid tracking is purely based on the surface model, while the FEM is used strictly for tracking the non-rigid deformations. This parallel system is capable of real-time performance. Our approach has been both quantitatively and qualitatively evaluated using a combination of ground-truthed, simulated data and real data.

II. METHOD

We will now describe the details of the methodology that we use for tracking deformable objects. We begin by describing the preliminary notations and the elastic deformation model that we use for interpreting the physics behind the objects we track. This is followed by a brief description of the rigid registration module (section II-C). The description of matching the deformable object between consequent frames is discussed next (section II-D). We follow this up with the mechanism for computing the Jacobian that links the variation in the geometric error with the variation of force applied on a particular vertex. We use this Jacobian for minimizing the error using an Iteratively Reweighted Least Square (IRLS) formulation. We finish this section by summarizing the steps for minimizing the non-rigid error.

A. Notation

We use two types of 3D model for the non-rigid tracking. The first is a visual-surface model, denoted by V, and the other one is a tetrahedral, volumetric, internal model, denoted by M. The visual-surface model is described using a set of planes which are represented using the nodes defining the boundary of the surface (usually a triangle) and the connectivity between these nodes. We can represent this surface as V = p 1 p 2 ... p M , where M is the number of vertices in V and the i-th vertex is p i = X i , Y i , Z i , the 3D coordinate of a point in the reference frame of the base model. We have a corresponding, connectivity Fig. 1: The depth camera and the base model is centered at F C , the tracked model is used for the rigid registration of the pointcloud map, given by

X V = a 1 , b 1 , c 1 ... a N , b N , c N , where a i , b i , c i ∈ [0.
.M ] denotes the indices in V and N is the number of faces in V, thereby defining connectivity between its elements. Clearly, the normal n i = (n X i , n Y i , n Z i ) and distance to origin d i for any given face can be derived from V and X V using elementary geometry. We similarly define a tetrahedral, volumetric, internal model M along with its connectivity map X M , such that every element of X M has four indices, instead of three. The visual-surface model and the internal model is maintained at the same reference frame and together, we refer to them as the base model.

In the proposed approach, a depth camera is used to track the object from the observed pointcloud. The camera is maintained at a camera-centered frame of reference F C . The intrinsic parameters of this camera are denoted by f x , f y , the focal lengths and c x , c y the optical centers of the image. The pointcloud captured by the depth camera is represented as a set of points P i = (x i , y i , z i ), which are the 3D coordinate of the points in the camera's reference frame. The visualsurface model of the object being tracked, denoted by the tuple {V, X V }, is also based on this same frame of reference, as shown in Fig. 1.

We use the homogeneous matrix

A T B = A R B A t B 0 1
to define the rigid transformation between any two arbitrary Cartesian frame F A and F B , where T ∈ SE 3 . In between frames, the rigid motion of the object with respect to its previous pose is denoted by q = ( A t B , θu), where θ and u are the angle and axis of the rotation A R B . The time derivative of q is given as v = δq, where v ∈ se(3) is the velocity screw.

C T O denotes the transformation from the camera-centered coordinate frame F C to the object in the pointcloud. For rigid tracking, we define a tracked model using the tuple

{ O V, X V }, such that: O V = O R C C V + O t C • • • O t C 3×M (1)
During non-rigid deformation of the base model ( C V) centered at F C , the transformation of ( 1) is repeated, so that the model at the object centered reference frame always stays updated according to the deformation.

B. Deformation Modelling

We use Finite Element Model (FEM) for modelling the deformation. For the object representation, we consider the visual-surface model and the internal model to be coupled with each other in a master-slave configuration. The internal model acts as the mechanical model, which is the master of the system [START_REF] Faure | Sofa: A multi-model framework for interactive physical simulation[END_REF]. We use a function J to map the position of

C M to C V such that p V i = J p M i .
We use a co-rotational FEM formulation for estimating the displacement of the nodes after being acted upon by an external force (for a detailed background on FEM, interested readers may refer to [START_REF] Cook | Finite element modeling for stress analysis[END_REF]). Beginning from the previous description of the visual and the internal model, we describe the co-rotational model for a single, tetrahedral element of M. Let the four vertices of this particular face be denoted by p 1 , p 2 , p 3 and p 4 . Let the centroid be F C , as shown in Fig. 2.

Let us assume that in the next frame, the model undergoes a rigid rotation, as well as a deformation. Let the rotated centroid at the new frame be given by C R . An arbitrary point p x in the first frame gets deformed to p in the next frame. Using the corotational model [START_REF] Felippa | A systematic approach to the element-independent corotational dynamics of finite elements[END_REF], it can be shown that:

d = R U -C + p x -p x ( 2
)
where R is the rotation matrix (also called the rotator), U is the total displacement of p x , expressed in the base frame C, and d is the deformational displacement of p x , expressed in the corotated frame C R (all variables with • are expressed in the corotated frame C R ). Now, applying similar motion to the nodes of the tetrahedron, we can represent the internal, elastic forces acting on the nodes by:

F e = R e B e ŨR (3)
Here, ŨR = ( p1 , p2 , p3 , p4 ) is the elastic forces acting on the nodes of this tetrahedral element.

To resolve the interaction between forces and their resulting displacement using the FEM, we solve a second order differential equation, given by:

Mp + D ṗ + Bp = F ext + F e (4)
Fig. 2: Demonstration of the kinematics of a single element of a corotated, tetrahedral mesh. In the corotated frame, the red outline shows the deformed surface where M and D are the model's mass and damping matrices respectively. B is the global stiffness matrix. F ext are the external forces acting on the vertices. The method for determining F ext is discussed in the subsequent sections. To solve this differential equation, we use a linear solver based on conjugate gradient descent [START_REF] Faure | Sofa: A multi-model framework for interactive physical simulation[END_REF]. To impose additional constraint on the solution of (4), it can be multiplied using a projection matrix that sets the values of certain indices to zero. We use the projective constraint to eliminate the rigid motion of the object. This is discussed in the next section It must be noted that the force F ext is just the deforming force acting on the model, and does not include the forces causing the rigid transformation of the body. In fact, the deformation model proposed here is completely independent of rigid motion, the effects of gravity and interaction with other contact surfaces. As described in section II-C, the rigid motion is tracked separately. This separation of rigid and non-rigid tracking method, along with the minimization described section II-D, makes the overall system independent from the inaccuracies of the physical parameters.

C. Rigid Registration

The rigid registration is a joint minimization of two error terms: depth based geometric error and keypoint based feature tracking [START_REF] Trinh | A modular framework for model-based visual tracking using edge, texture and depth features[END_REF]. The two error terms are explained below:

1) Depth based geometric error: Assuming that we know the accurate C T n-1 O at the (n -1)-the frame and n-1 T n gives the initial estimate of transformation between previous and current frame, the error is given as:

e D ( n-1 q n ) = n-1 R n P + n-1 t n • n k -d k (5)
where n k and d k are the normal and distance to origin respectively, for the k-th planar face that corresponds with the point P in the pointcloud. The Jacobian that links the variation of the error e D ( n-1 q n ) with the time varying n-1 qn , is given by:

J D i = n k n k × P i (6)
2) Feature based minimization: The Harris corner features are used for tracking with keypoints. This is based on the classical KLT algorithm [START_REF] Baker | Lucas-kanade 20 years on: A unifying framework[END_REF]. Let u = (x, y, 1) be the homogeneous 2D coordinate of a feature point in the (n-1)th image and u * = (x * , y * , 1) be the matched coordinate for the same point in the n-th frame. We define the error term as:

e K ( n-1 q n ) = x( n-1 q n ) -x * y( n-1 q n ) -y * (7)
where (x( n-1 q n ), y( n-1 q n ), 1) = n-1 H n p, given that:

n-1 H n = n-1 R n + n-1 t n d n (8)
where n and d are interpreted in the same way as [START_REF] Petit | Tracking fractures of deformable objects in real-time with an rgb-d sensor[END_REF]. The corresponding Jacobian that relates the variation of e K ( n-1 q n ) with the time variation of n-1 q n is given by:

J K = -1 Z 0 x Z xy -(1 + x 2 ) y 0 -1 Z y Z (1 + y 2 ) -xy -x (9)
Given these two errors and their corresponding Jacobian, the combined error e( n-1 q n ) is obtained by stacking the vectors e D ( n-1 q n ) and e K ( n-1 q n ), while the combined Jacobian J is obtained by stacking J D on J K . The combined error is minimized with the update given by:

v = -λ WJ + We ( 10
)
where W is the weight matrix for outlier removal using the Tukey m-estimator [START_REF] Beaton | The fitting of power series, meaning polynomials, illustrated on band-spectroscopic data[END_REF], and v ∈ se(3). This becomes an iteratively re-weighted least squares problem.

D. Non Rigid Tracking

Before starting the non-rigid tracking using a minimization technique, it is necessary to determine the possible points of application of forces that needs to be tracked. It is sufficient to track only a small subset of nodes, depending upon a set of criteria. These criteria, as detailed below, is obtained by clustering the error from the rigid tracking step. Only regions surrounding the clusters with higher value of errors are considered to be relevant for non-rigid tracking.

At the end of the rigid registration process, we can easily construct a map that links the points in the 3D pointcloud to their corresponding geometric error from [START_REF] Petit | Tracking fractures of deformable objects in real-time with an rgb-d sensor[END_REF]. Given that e D i denotes the error for the i-th point P i in the entire pointcloud P, it is possible to derive a new pointcloud P * using a linear thresholding operation, such that P * = {P i ∈ P|e D i ≥ θ D }, where θ D is a threshold that depends on the geometry of the object being tracked. We subject the pointcloud P * to a clustering step using Eucledian distances [START_REF] Hartigan | Algorithm as 136: A k-means clustering algorithm[END_REF]. Very small clusters and clusters spanning more than half of the size of the entire pointcloud are discarded. Let us assume that we obtain j clusters from this operation, denoted by K 1 , K 2 , • • •, K j , and their corresponding centroids are represented by k 1 , k 2 , •••, k j , where k = (x, y, z, 1) is the homogeneous 3D coordinate of the point w.r.t F C . These points are inverse transformed with the last estimate of C T O obtained after minimization with [START_REF] Petit | Real-time tracking of 3d elastic objects with an rgb-d sensor[END_REF]

, such that k * i = C T O -1 k i .
For all these j centroids, we determine their nearest neighbors in C M using the k-nearest neighbor algorithm [START_REF] Fukunage | A branch and bound algorithm for computing k-nearest neighbors[END_REF]. Let the nearest neighbors of k

* 1 , k * 2 , • • •, k * j in C M be denoted by C M k1 , C M k2 , • • •, C M kj .
For the following discussion, let us consider an arbitrary k i -th cluster alone, which we refer to as C M i for the sake of simplicity. At this stage, we must modify (5) to accommodate a slightly different variant of the same error function. We redefine it as e N O V i = n r • P s -d r assuming the r-th plane corresponds to the s-th point in the pointcloud, and the normal n r and the distance to origin d r is derived from

O V i .
Jacobian Computation: The Jacobian that relates the variation of e N O V i with the variation of the applied force is computed numerically by perturbing the node C M i by a very small, constant force ∆F J successively along the three axes, as shown in Fig. 3. After applying the force on the node, a simulation of FEM is done using the minimization described in section II-B -(4). After the conjugate gradient Fig. 3: The node closest to the centroid of a cluster is perturbed by a small force along the three axes in both positive and negative direction, producing six deformed configurations per node solver optimizes (4), the system will attain a static condition.

There will be six such configurations of the mesh, obtained by the perturbation along the positive and the negative direction of the three axes. We denote these new configurations of

C M i as C M X+ Ji , C M X- Ji , C M Y + Ji , C M Y - Ji , C M Z+ Ji and C M Z-
Ji . This is immediately propagated down to the visualsurface model via the map C V Ji = J ( C M Ji ). We take these modified, visual mesh and transform it back to the registered pointcloud using [START_REF] Alexa | As-rigid-as-possible shape interpolation[END_REF]. The relation of the variation of the external force with the variation in error can be expressed as ėN ( O V Ji ) = J i Ḟext The term J i is obtained numerically by finite difference computation. The final Jacobian used is:

J i = J X i J Y i J Z i .
The temporary deformation of the model for Jacobian computation is discarded, once the Jacobian has been determined.

Minimization: So far, j nodes have been selected for application of external force to deform the model. We have n points in the pointcloud, and their correspondence with the model is known. At the first iteration, the initial estimate of the vertically stacked force vector F n-1 ext 3j×1 is set to zero. These Jacobian matrices are stacked up horizontally, such that

J n×3j = J 1 J 2 • • • J j The update is computed as: ∆F = -λ W n×n J + We N ( O V) n×1
The force vector is updated by

F n ext = F n-1 ext + ∆F.
This force is applied on the base model and the final node displacements are determined, once again, using (4).

III. IMPLEMENTATION

There are a few practical aspects of implementing this algorithm which needs to be discussed. It is necessary to associate every 3D point P i with one of the planes of the visual-surface model O V. At every frame, we know the value of C T O obtained at the previous frame. We project the 3D points O V for all visible planes into the image obtained in the current frame. The visibility is checked using the classical ray-casting algorithm [START_REF] Roth | Ray casting for modeling solids[END_REF]. The ray-casting is also necessary for imposing the projective constraints. The nodes invisible to the camera are considered to be immobile. We do not consider the volumetric, internal mesh to be an input to the system. The mesh is rather generated apriori, using Dirichlet tessellation, followed by Delaunay tetrahedralization of the input mesh using the Bowyer-Watson algorithm [START_REF] Si | Tetgen, a delaunay-based quality tetrahedral mesh generator[END_REF].

The approach proposed here does not require a visual segmentation for separating the region of interest from the background. This is done by initializing the pose of the object at the first frame using pre-trained markers on the object. It is done using the ViSP library [START_REF] Marchand | Visp for visual servoing: a generic software platform with a wide class of robot control skills[END_REF], by matching the keypoints detected in the very first image with those extracted in the training images using an approach similar to [START_REF] Choi | Real-time 3d model-based tracking using edge and keypoint features for robotic manipulation[END_REF].

The algorithm is implemented using a parallel framework, where two different processes are involved with the rigid and non-rigid tracking of the object. The rigid tracking module is capable of processing each frame at < 100 ms, thereby ensuring real-time interaction. The un-optimized, non-rigid tracking code runs at 800ms to 1.3s (approximately) per frame, depending on the size of the object and its proximity to the camera. The results reported here have been achieved using processes running on a single core of a i7-6600U CPU with 16GB of RAM. An Intel Skylake GT2 GPU has been utilized, but only for the ray-casting. The FEM solvers are implemented using the SOFA library [START_REF] Faure | Sofa: A multi-model framework for interactive physical simulation[END_REF], without using CUDA. The rigid tracking process runs at > 10 fps, while the non-rigid tracking module handles every 10-th/12-th frame coming from the sensor.

IV. RESULTS

The results are validated on two sets of simulated data, as well as on multiple real objects 1 .

Simulation: For the simulated data, the deformation of the objects were generated using simulation of co-rotational FEM. Two objects are considered, a cube and a rectangular board. They are made to undergo simple but large deformations, as shown in Fig. 4. We made the dataset available at: github.com/lagadic/nr-dataset. The visual models produced as an output of the simulation were subjected to texturing, shading and rendering using the Blender2 software, followed by the generation of the pointcloud using a RGB-D camera simulator. In the simulation, both the objects were modelled using Young's modulus (YM) of 50000 Pa and Poisson's ratio of 0.3. Rayleigh mass was assumed to be 0.1 and Rayleigh stiffness to be 0.3.

The simulated dataset is used not only to validate our approach, but also to compare the proposed method to a partial re-implementation of [START_REF] Petit | Capturing deformations of interacting non-rigid objects using rgb-d data[END_REF], where a co-rotational FEM is simulated for non-rigid object tracking. The external, elastic forces, f ext , acting on the object is given by f ext,i = k ext,i (x i -y i ), where k ext is the stiffness corresponding to the external elastic force and y i is a point in the pointcloud which has been matched to x i , a vertex on the mesh. We reimplemented this method. Moreover, to reduce ambiguity, the Since the ground-truth for the deformation is known, it is possible to compute the Hausdorff distance metric [START_REF] Henrikson | Completeness and total boundedness of the hausdorff metric[END_REF] between the output of the tracking and the ground-truth, to quantify the error between the ideal and the actual output. Fig. 6 delineates the Hausdorff distance of the output for various values of Young's modulus.

As given in Table I, the mean error in terms of Hausdorff distance using PA and SOA are 0.305 and 1.622 units for the cube dataset and 0.834 and 1.100 units for the board dataset respectively (the length of the largest diagonal is 69.28 and 55.19 units for the cube and the board models respectively). Hence, the tracking accuracy of the proposed approach is 81.18% better than SOA in the cube data and 24.20% better than SOA in the board data. The tracking accuracy of the proposed approach in terms of Hausdorff distance varies only by an average of 5.006% when the Young's modulus is increased by a factor of 10 8 . It can therefore be claimed that the system is significantly robust to error in the estimation of the physical properties used to model the deforming object. It must be noted that the results shown for SOA contains only the error in estimating the magnitude of deformation of the objects. As stated before, the direction of application of force is provided to SOA from the ground-truth in our re-implementation (while PA makes no such exception), which explains the impressive accuracy of SOA when using Young's modulus equal to the ground- Real Data: The real data has been captured using the Intel RealSense SR300 RGB-D camera. It is a commodity level depth sensor that produces RGB-D images of reasonably good resolution between 20 to 150 cm distance. Three objects have been tracked for the validation: a) a pizza, b) a cuboidal soft toy and c) a rectangular sponge. A Poisson's ratio of 0.3 was assigned to all the objects uniformly. All

V. DISCUSSION

The approach proposed in this paper accurately tracked all the deforming objects that it had been experimented upon. The results are promising and the validation tests of the approach had been successful. This approach has shown significant improvement over equivalent methods seen in the literature. Using the proposed approach, FEM based deformation tracking no longer needs to be fine tuned for the precise physical parameters of the objects being tracked. The overall algorithm, especially the Jacobian computation, remains highly parallelizable. There remains significant scope for performance optimization, which could be an important future work.

VI. CONCLUSION

We present an algorithm to track deformable objects using RGB-D camera. The FEM based deformation tracking methods available in the literature suffers from the drawback of being excessively dependant on the accurate knowledge However, the approach proposed in this paper removes this dependency, while accurately tracking the deforming object. This will be beneficial when tracking deformation in objects with unknown physical properties. The algorithm is capable of real-time interaction using a parallelized architecture.

The approach has been validated on simulated objects with ground-truth, as well on real objects of unknown physical properties. The simulated dataset with ground-truth and all real data used in this paper has been made publicly available for future research and comparison.

ACKNOWLEDGEMENT

We gratefully acknowledge the support from The Research Council of Norway through participation in iProcess-255596 project and from Inria.

12×1,

  B e 12×12 is the stiffness matrix and R e 12×12 is the block diagonal matrix of four R rotation matrices stacked diagonally. F e 12×1

Fig. 4 :

 4 Fig. 4: The undeformed model of the board and the cube. The arrows in the image shows the approximate direction of application of force

8 ( 8 (Fig. 5 :

 885 Fig. 5: Comparison between the proposed approach (PA) and state-of-the-art (SOA) for the cube and the board sequence. The images show the tracked model, positioned in the pointcloud with the current estimate of C T O , thereby causing the model to be partially occluded by the pointcloud. The small, red cubes are the model vertices. In the experiments, SOA fails to track accurately when YM is lower than the ground-truth value (5 × 10 4 Pa) and tends to overcompensate the deformation when YM is higher. PA maintains consistently accurate tracking, regardless of the changes in the value of YM

Fig. 6 :

 6 Fig. 6: Comparison of Hausdorff distance between output and ground-truth for PA and SOA for multiple values of Young's modulus Young's Modulus: 5 Pa 500 Pa

  the objects were tracked using multiple values of Young's modulus which were set coarsely and empirically. The pizza data was tracked with an Young's modulus of 5 × 10 3 Pa, 5 × 10 4 Pa and 5 × 10 5 Pa, the toy was tracked using 5 × 10 4 Pa, 5 × 10 5 Pa and 5 × 10 6 , and the sponge was tracked using 8 × 10 5 Pa, 8 × 10 6 Pa and 8 × 10 7 Pa. The tracking results are shown in Fig. 7. Given these results, the average, per pixel point-to-plane error for multiple values of Young's modulus are given in Fig. 8.

Fig. 7 :

 7 Fig. 7: The tracking out for the pizza (top), toy (middle) and sponge (bottom). Input denotes the pointcloud obtained from the RGB-D camera, projected on the image plane, the second row shows the tracked model placed inside the pointcloud with the current estimate of C T O and the third row depicts the base model separately

Fig. 8 :

 8 Fig. 8: Summary of point-to-plane error obtained from experiments on real data

TABLE I :

 I

Table for comparison of Hausdorff distance between output and ground-truth for PA and SOA for multiple values of Young's modulus

All results can be viewed at: youtu.be/2aqrGtLcrqU

https://www.blender.org/


Mia Baise

Phillip M Maffettone

Fabien Trousselet

Nicholas P Funnell

Franc ¸ois-Xavier Coudert

Andrew L Goodwin

Negative Hydration Expansion in ZrW 2 O 8 : Microscopic Mechanism, Spaghetti Dynamics, and Negative Thermal Expansion

We use a combination of X-ray diffraction, total scattering and quantum mechanical calculations to determine the mechanism responsible for hydration-driven contraction in ZrW2O8. Inclusion of H2O molecules within the ZrW2O8 network drives the concerted formation of new W-O bonds to give one-dimensional (-W-O-)n strings. The topology of the ZrW2O8 network is such that there is no unique choice for the string trajectories: the same local changes in coordination can propagate with a large number of different periodicities. Consequently, ZrW2O8•H2O is heavily disordered, with each configuration of strings forming a dense aperiodic 'spaghetti'. This new connectivity contracts the unit cell via large shifts in the Zr and W atom positions. Fluctuations of the undistorted parent structure towards this spaghetti phase emerge as the key NTE phonon modes in ZrW2O8 itself. The large relative density of NTE phonon modes in ZrW2O8 actually reflects the degeneracy of volumecontracting spaghetti excitations, itself a function of the particular topology of this remarkable material.

Zirconium tungstate, ZrW 2 O 8 , is an important material because it harbours a number of useful and counterintuitive mechanical responses to external stimuli. First, on heating, its molar volume decreases such that the material is 2% denser at 1050 K than at 0.3 K [1][2][3][4]. This negative thermal expansion (NTE) effect is about as strong as the usual positive thermal expansion (PTE) of conventional materials, and so can be exploited in athermal composites [5][6][7]. Second, under hydrostatic pressure, ZrW 2 O 8 amorphises unusually easily, transforming into a disordered phase that is denser by more than 20% [8][9][10]. This switching between crystalline and amorphous states is related to the mechanism exploited for data storage in phase-change chalcogenides [11]. And, third, hydration of ZrW 2 O 8 gives the monohydrate ZrW 2 O 8 •H 2 O, which has a molar volume ∼10% smaller than ZrW 2 O 8 itself [12,13]. In this sense the material-which has no appreciable pore network-is nonetheless a kind of 'inverse sponge' that shrinks instead of swells when it takes up water. We call this effect "negative hydration expansion" (NHE). NHE is known to occur in anisotropic 'breathing' frameworks such as MIL-53 [14,15], but is rare in cubic systems [16,17]; pore sizes are typically 5-20 Å in all cases. The phenomenon is largely unexplored-and particularly so in ZrW 2 O 8 -yet NHE may eventually have applications in counteracting hydration-driven swelling of construction materials (e.g. concretes and clays).

While the mechanism responsible for NHE in ZrW 2 O 8 remains unclear [12,13], both NTE and pressure-induced amorphisation (PIA) have been studied extensively and are understood to arise from a peculiar phonon spectrum. The material supports an extensive family of low-energy modes with large and negative Grüneisen parameters [18][19][20][21][22][23] (i.e. their frequencies soften with pressure)-the thermodynamic requirement for NTE in cubic systems [6]. Under hydrostatic pressure the energies of these modes fall to zero [8,24]. Because the modes are distributed throughout the Brillouin zone (BZ), the resulting structural transition is associated with a family of modulation periodicities, and hence amorphisation [10].

Despite the many experimental and computational studies of ZrW 2 O 8 , the microscopic origin of its phonon anomalies remains contentious. It is the sheer number of NTE vibrational modes that is difficult to explain: each one appears to involve a remarkably different combination of polyhedral translations, rotations, and distortions [22]. Early indications of the importance of 'rigid unit modes' (RUMs) [25][26][27] have increasingly lost experimental [28,29] and computational [22,23] support, with consensus only emerging regarding the role of polyhedral distortions and the 'skipping rope' or 'tension' mechanism proposed in Ref. space actually reflect the degeneracy of volume-contracting spaghetti excitations, itself a function of network topology.

The starting point for our study was to prepare and characterise polycrystalline samples of both ZrW 2 O 8 and its monohydrate, following the hydrothermal methodologies of Ref. [START_REF] Banek | Autohydration of Nanosized Cubic Zirconium Tungstate[END_REF]. The corresponding X-ray powder diffraction patterns [Fig. 1] reflect all hydration-driven structural responses noted in earlier studies [12,13]: the crystal symmetry increases on hydration from P 2 1 3 to P a 3, the lattice parameter a decreases from 9.153(3) to 8.876(4) Å, and there is a marked increase strain-driven peak broadening. We found hydration/dehydration cycles (∼36 h) to be reversible in all respects, including the appearance and disappearance of strain effects. Rietveld refinement of the parent ZrW 2 O 8 diffraction pattern using the published structure [3] gave an excellent fit, but we encountered the same difficulty noted in [12] of obtaining a structural model for the monohydrate: the scattering density at the Zr and W sites indicated severe positional disorder and the O atom positions were impossible to determine.

Given the difficulty of directly interpretating the X-ray diffraction pattern of ZrW 2 O 8 •H 2 O, we used quantum mechanical calculations to determine an energetically-sensible candidate structure for this phase. The structure of ZrW 2 O 8 contains a single symmetry-related set of cavities, with each cavity of the correct size to host one H 2 O molecule [Fig. 2(a)]. There is one cavity per formula unit, and four formula units per unit cell. The cavities are positioned on the 4a Wyckoff sites and coincide with three-fold rotation axes of the P 2 1 3 space group. Clearly this point symmetry is incompatible with that of the H 2 O molecule. So there is a choice of three possible orientations for each of the four water molecules to be placed in a single unit cell. This situation contrasts that of NH 3 binding in ZrW 2 O 8 •NH 3 , for which NH 3 orientations are ordered and P 2 1 3 symmetry is maintained [START_REF] Cao | [END_REF]-the molecular point symmetry of NH 3 (3m) being a supergroup of the 4a site symmetry (3) [START_REF]molar volume of ZrW2O8•NH3 is only[END_REF]. Instead the arrangement of H 2 O molecules that results in the least severe symmetry lowering gives P 2 1 2 1 2 1 symmetry. We will come to show that ZrW 2 O 8 •H 2 O is in fact heavily disordered, with P a 3 symmetry emerging as a configurational average over states related to this P 2 1 2 1 2 1 arrangement we now consider in detail.

Starting with this highest-symmetry H 2 O arrangement, we carried out density functional theory (DFT) relaxations at two different levels of theory (LDA and PBESOL0; see [START_REF]See Supplemental Material for details concerning methods, the P 212121 DFT structural model of ZrW2O8•H2O, and configurational entropy calculations[END_REF]). Energy minimisation was sensitive to the precise starting positions of the additional H 2 O molecules within the framework cavities, yet we were able to identify a series of closely-related low-energy structures with molar volumes smaller by 6-8% than that of ZrW 2 O 8 . The corresponding enthalpies of water adsorption were in the range 60-90 kJ mol -1 (PBESOL0 functional), indicating strong physisorption [START_REF] Bolis | Calorimetry and Thermal Methods in Catalysis[END_REF]. Hence our calculations are entirely consistent with a NHE effect of the same magnitude as that observed experimentally.

The structure of our lowest-energy P 2 1 2 1 2 1 model for ZrW distinct W atoms [W1 in Fig. 2(b)] along a direction that is no longer parallel to a triad axis of the original unit cell. This new coordination induces a series of coupled rotations and translations of both types of W-centred polyhedra that collectively results in a set of new W-O-W connections. These connections draw the framework in on itself: the terminal W1-O and W2-O bonds of the original structure are forced within bonding distance of nearby W2 and W1 atoms, respectively. Not only does this connectivity reduce the average W. . .W separation-and in turn the unit-cell volume-but it results in substantial displacement of both Zr and W atoms away from their original sites. Indeed, provided we can account for the P a 3 crystal symmetry, these displacements are fully consistent with experiment: Fig. 1(b) shows the excellent match between the experimental diffraction pattern and that derived from our DFT coordinates. Moreover, our model accounts for the NTE/PTE switch on hydration: zone-centre phonon calculations show a dramatic reduction in the number of NTE modes, with the corresponding coefficient of thermal expansion now positive [START_REF]See Supplemental Material for details concerning methods, the P 212121 DFT structural model of ZrW2O8•H2O, and configurational entropy calculations[END_REF].

Given the existence of this low-energy P 2 1 2 1 2 1 structure, why would the experimental system retain cubic symmetry? We suggest the key is the existence of a diverging number of different ways of propagating the local changes in coordination identified by our DFT calculations. This degeneracy, we claim, leads to a strongly disordered state. We proceed to describe this state and in turn test our claims using experimental X-ray pair distribution function (PDF) measurements.

Consider the new W-O-W bonds formed during hydration [Fig. 2(b)]. While there is only one W1 atom accessible to each terminal W2-O bond, there are three equally-spaced W2 atoms which any terminal W1-O bond might approach. These three choices of W2 atom are related by the same triad axis that is broken by the symmetry of an H 2 O molecule. So, in selecting one of three possible states, the particular orientation of each H 2 O molecule effectively determines the trajectory of a corresponding W1-O→W2 linkage [Fig. 2(c)]. These new W-O-W bonds form one-dimensional strings that in our DFT structure reflect the same P 2 1 2 1 2 1 crystal symmetry of the H 2 O molecule arrangements [Fig. 2(d)]. But more disordered H 2 O arrangements are associated with equally disordered arrangements of W-O-W strings (we will call collections of strings a 'spaghetti', noting the relationship to the 'loop' models [START_REF] Jaubert | [END_REF] of e.g. Coulomb phases [37,38]). Such states are not random. Note that in the DFT structure each W1-O bond approaches only one W2 atom, and each W2 atom is approached by only one W1-O bond. Preserving these local constraints enforces correlations between the new W-O-W bonds formed, and in turn between the orientations of successive H 2 O molecules. It can be shown [START_REF]See Supplemental Material for details concerning methods, the P 212121 DFT structural model of ZrW2O8•H2O, and configurational entropy calculations[END_REF] that the number of configurations satisfying these local rules grows exponentially with system size to give a macroscopic number of viable hydration states; the corresponding configurational entropy is S config R ln(9/8). Hence there is an entropic driving force for ZrW 2 O 8 •H 2 O to adopt spaghetti states, with the bulk crystal symmetry reflecting the (cubic) configurational average.

There is a robust analogy between this model and icelike states on the pyrochlore lattice, which also preserve cubic symmetry [37,39,40].

In order to test the relevance of this spaghetti model to the NHE mechanism of ZrW 2 O 8 •H 2 O, we proceeded to measure X-ray total scattering data for both ZrW 2 O 8 and its hydrate. Such data are directly sensitive to the local structure changes induced by hydration. We then developed a supercell model of a disordered spaghetti state, informed by the DFT structure we have already described, from which an X-ray-weighted PDF might be calculated and compared against experiment.

Our total scattering measurements were performed using a PANalytical Empyrean X-ray diffractometer fitted with an Ag anode (Q max = 20 Å-1 ) and a GaliPIX3D detector; the corresponding PDFs [41], normalised using GUDRUNX [42], are shown in Fig. 3(a). We were able to model the PDF of ZrW 2 O 8 in PDFGUI [43] using the crystal structure of Ref. 3. By contrast, the P 2 1 2 1 2 1 model for ZrW 2 O 8 •H 2 O is capable of fitting convincingly only the very lowest-r region of the corresponding PDF, and cannot account for the positional disorder evident at distances r 4 Å.

Because it was not computationally feasible to explore a disordered spaghetti configuration at the DFT level, we used harmonic lattice-dynamical methods, informed by the bonding geometries observed in our P 2 1 2 1 2 1 DFT structure, to generate a 5 × 5 × 5 approximant [44] [Fig. 3(b)]. First we determined an appropriate disordered spaghetti of selfavoiding W-O-W pathways using a custom Monte Carlo algorithm, guided by the connectivity rules of Fig. 2(c). This spaghetti was used to assign both connectivity and H 2 O orientations within a supercell of the P 2 1 3 ZrW 2 O 8 structure. Atomic positions and unit cell dimensions were then relaxed using the GULP code [START_REF]See Supplemental Material for details concerning methods, the P 212121 DFT structural model of ZrW2O8•H2O, and configurational entropy calculations[END_REF]45], with harmonic potentials constructed from the equilibrium bond distances and angles of our earlier DFT calculation. The cell volume decreased during relaxation, reflecting activation of the NHE mechanism as described above. We consider this GULP-relaxed configuration to be representative of the degree of structural disorder implied by the spaghetti model. The connectivity of the relaxed structure is shown in Fig. 3(b), and the PDF calculated from the full model compared against experiment in Fig. 3(a). The quality of this fit-to-data shows that the displacive disorder introduced by spaghetti W-O-W connectivity can indeed account for the broad oscillations in the experimental PDF.

A necessary consequence of the macroscopic configurational entropy of spaghetti formation on the ZrW 2 O 8 lattice is that the corresponding hydration-driven distortions are dense throughout the BZ. Hence the structure of ZrW 2 O 8 •H 2 O reflects a superposition of spaghetti states with arbitrarily large periodicities. Anisotropy vanishes in this limit such that the corresponding global symmetry is cubic, with an average structure given by the projection of all possible distortions onto a single unit cell [Fig. 1(b)]. The large displacements involved account for the substantial crystal strain observed experimentally, and also for the difficulty of interpreting even high-quality neutron diffraction data [12]. In turn, proton migration and topological defects in spaghetti state formation account for the introduction of inversion symmetry. quently we conclude that NHE is driven by a degenerate family of 'spaghetti'-like lattice distortions modulated with periodicities distributed continuously throughout the BZ.

Having established this mechanism, we sought to determine whether the same distortions that drive NHE are related to the phonons responsible for NTE in ZrW 2 O 8 . Published data [23] allow access to the frequencies ω(k, ν), eigenvectors e(k, ν), and Grüneisen parameters

γ(k, ν) = - V ω(k, ν) ∂ ∂V ω(k, ν) (1)
for each phonon branch ν and wave-vector k. NTE modes are those for which γ < 0. By comparing the atomic coordinates of each ZrW 2 O 8 framework atom in the anhydrous and hydrated DFT structures, we formed a normalised, massweighted displacement vector u that describes the volumereducing framework deformations induced by hydration. The projections c(k, ν) = u • e(k, ν) partition these deformations amongst the normal modes; for the (k = 0) DFT structure the c(k, ν) necessarily vanish for k = 0. Remarkably, we find that 88% of the displacements described by u can be accounted for in terms of two of the lowest-energy k = 0 NTE modes. The same result is found for zone-centre phonons generated by our quantum mechanical calculations [START_REF]See Supplemental Material for details concerning methods, the P 212121 DFT structural model of ZrW2O8•H2O, and configurational entropy calculations[END_REF]. Extending our analysis to wave-vectors k = 0 requires access to supercell configurations. So, using the same GULP approach described above, we generated a family of ten different modes (mostly acoustic) account for ∼45% of the hydrationdriven displacement, and on including all modes with negative Grüneisen parameters this total projection rises to 87% of the total framework response [Fig. 4] [START_REF]See Supplemental Material for details concerning methods, the P 212121 DFT structural model of ZrW2O8•H2O, and configurational entropy calculations[END_REF].

So the most important NTE phonon branches are strongly implicated in the NHE response of ZrW 2 O 8 , but the question remains of whether the spaghetti displacements invoked for NHE are necessarily involved in the NTE dynamics. The c(k, ν) play a crucial role in answering this question: whereas previous studies have focussed either on individual collective modes [22] or on the correlated motion of isolated structural fragments [29], the c(k, ν) allow us to reconstruct combination NTE modes with intuitive real-space realisations. For example, the P 2 1 2 1 2 1 DFT configuration of Fig. 2(b) gives a framework displacement vector u that-as we have said-projects almost entirely onto just two zonecentre modes ν 1 , ν 2 (with ω 1 = ω 2 = 44.1 cm -1 and γ 1 = γ 2 = -4.14). The combination mode c(0, ν 1 )ν 1 + c(0, ν 2 )ν 2 then describes one specific collective vibration of ZrW 2 O 8 contributing to NTE. Crucially, the associated displacement vector e = c(0, ν 1 )e(0, ν 1 ) + c(0, ν 2 )e(0, ν 2 ) does indeed involve fluctuations towards the very same spaghetti states assumed by ZrW 2 O 8 •H 2 O; an animation of this combination mode is provided as SI [START_REF]See Supplemental Material for details concerning methods, the P 212121 DFT structural model of ZrW2O8•H2O, and configurational entropy calculations[END_REF]. The same result also holds for the more complex family of NTE modes associated with the 3 × 3 × 3 spaghetti configurations described above [START_REF]See Supplemental Material for details concerning methods, the P 212121 DFT structural model of ZrW2O8•H2O, and configurational entropy calculations[END_REF].

The significance of this result is that the key real-space mechanism for NTE in ZrW 2 O 8 is the correlated translation, rotation and distortion of ZrO 6 /WO 4 polyhedra towards the volume-reducing spaghetti states adopted in ZrW 2 O 8 •H 2 O. The configurational degeneracy of these states-a property of the peculiar ZrW 2 O 8 topology-is such that these fluctuations are described by a family of phonons necessarily spread throughout the BZ, with individual spaghetti distortions described by the population of a family of NTE phonons. Consequently a 'spaghetti dynamics' mechanism for NTE in ZrW 2 O 8 finally explains (i) why NTE modes are distributed throughout the BZ, (ii) why it is so difficult to interpret in-dividual NTE phonons in isolation, and (iii) how the local polyhedral rotations/translations identified in previous studies [27,29] actually propagate throughout the crystal lattice. Moreover, in determining this mechanism, we have established an unexpected link between the correlated structural disorder associated with NHE, on the one hand, and the property of NTE, on the other hand, that is of significant conceptual importance in the broader context of understanding disorder-property relationships in functional materials [40].

FIG. 1 :

 1 FIG. 1: X-ray diffraction patterns (λ = 1.54056 Å) of (a) anhydrous and (b) hydrated ZrW2O8: data in black, fits in red, difference (data-fit) in blue, and reflection positions as tick marks. The fit in (a) is from a Rietveld refinement using the coordinates of Ref. 3; * = minor impurity. The fit in (b) used our DFT coordinates, with P a 3 symmetry applied (inset: the corresponding structural model represents the configurational average of a disordered spaghetti).

  FIG. 2: (a) The crystal structure of ZrW2O8: Zr in green, W in blue, O in red, and the four cavities of each unit cell in surface representation. (b) DFT geometries of the key structural fragment before (top) and after (bottom) hydration. Note the formation of additional W-O-W bonds. W2. . .W1 pairs are indicated by a grey arrow; the location of these pairs within the structure is shown in (a). (c) On hydration, neighbouring W2. . .W1 pairs connect to form W-O-Wchains. Note the symmetry-equivalent choices for these new connections (dashed lines). The red cylinders denote the connections made in the DFT (P 212121) structure. (d) The final P 212121 structure, with W-O-W strings shown in grey and H2O molecules shown in space-filling representation.

FIG. 3 :

 3 FIG. 3: (a) X-ray PDFs G(r) of anhydrous (top) and hydrated (bottom) ZrW2O8: data as points and fits as lines. For ZrW2O8 a satisfactory fit is obtained using the model of Ref. 3 (red lines). For the hydrate, the DFT P 212121 model fits at low r but cannot reproduce the broad oscillations at r > 4 Å (thin black line). Instead, the PDF calculated from our disordered spaghetti model (red line) captures all features. (b) The relaxed GULP model for a ZrW2O8•H2O spaghetti configuration as described in the text: only W-O-W linkages are shown, using grey cylinders [cf. Fig. 2(d)].

FIG. 4 :

 4 FIG. 4: Hydration-driven framework displacements in ZrW2O8•H2O project strongly onto NTE phonon eigenvectors in ZrW2O8. The average magnitude of projection coefficients c (vertical bars) is shown as a histogram of γ/ω; nearly 90% of the cumulative projection (blue shaded region) is accounted for by NTE modes.

  3. Intriguingly, hydration switches off NTE in ZrW 2 O 8 : ZrW 2 O 8 •H 2 O is a conventional PTE material [13, 30]. One implication is that the mechanism responsible for NHE saturates the same set of displacements from which NTE arises. Consequently we sought to determine the structure of ZrW 2 O 8 •H 2 O, establish a mechanism for NHE, and explore the relationship between NHE and NTE in ZrW 2 O 8 itself. In this Letter, we show that hydration of ZrW 2 O 8 proceeds via concerted formation of one-dimensional strings of new W-O bonds. The topology of the ZrW 2 O 8 network is such that there is no unique choice for the trajectory of these

	strings: the same local changes in coordination can propagate
	with a macroscopically large number of different periodici-
	ties. Each configuration of strings forms a dense 'spaghetti'
	that describes an additional connectivity responsible for con-
	tracting the unit cell via surprisingly large shifts in the Zr and
	W atom positions. So ZrW 2 O 8 •H 2 O is heavily disordered, al-
	beit in a highly correlated fashion. Remarkably, fluctuations
	of the anhydrous parent towards this spaghetti phase emerge
	as the key NTE phonon modes in ZrW 2 O 8 . Consequently the
	density and distribution of NTE modes throughout reciprocal

  3 × 3 × 3 spaghetti configurations. For each ZrW 2 O 8 •H 2 O configuration we extracted the corresponding framework displacement vectors u and projected these onto the phonon eigenvectors e(k, ν) of ZrW 2 O 8 for k ∈ n1

	3 , n2 3 , n3 3

*

, n i ∈ (0, ±1). We find that the projections onto the very lowest energy

The authors gratefully acknowledge financial support from the E.R.C. (Grant 279705), the E.P.S.R.C. (U.K.), the Leverhulme Trust (Grant No. RPG-2015-292), a CNRS/Oxford collaboration grant, and from the Marshall Foundation to P.M.M. This work benefitted from access to HPC platforms provided by a GENCI grant (A0030807069). ALG thanks J. S. O. Evans (Durham), J. P. Attfield (Edinburgh), B. Slater (UCL), and T. Ledsam (Oxford) for valuable discussions.