Ahmet Iscen

Yannis Avrithis

Giorgos Tolias

Teddy Furon

Ondřej Chum

Hybrid Diffusion: Spectral-Temporal Graph Filtering for Manifold Ranking

State of the art image retrieval performance is achieved with CNN features and manifold ranking using a k-NN similarity graph that is pre-computed off-line. The two most successful existing approaches are temporal filtering, where manifold ranking amounts to solving a sparse linear system online, and spectral filtering, where eigen-decomposition of the adjacency matrix is performed off-line and then manifold ranking amounts to dot-product search online. The former suffers from expensive queries and the latter from significant space overhead. Here we introduce a novel, theoretically well-founded hybrid filtering approach allowing full control of the space-time trade-off between these two extremes. Experimentally, we verify that our hybrid method delivers results on par with the state of the art, with lower memory demands compared to spectral filtering approaches and faster compared to temporal filtering.

Introduction

Most image retrieval methods obtain their initial ranking of the database images by computing similarity between the query descriptor and descriptors of the database images. Descriptors based on local features [START_REF] Sivic | Video Google: A text retrieval approach to object matching in videos[END_REF][START_REF] Philbin | Object retrieval with large vocabularies and fast spatial matching[END_REF] have been largely replaced by more efficient CNN-based image descriptors [START_REF] Gordo | End-to-end learning of deep visual representations for image retrieval[END_REF][START_REF] Radenović | Fine-tuning CNN image retrieval with no human annotation[END_REF]. Regardless of the initial ranking, the retrieval performance is commonly boosted by considering the manifold structure of the database descriptors, rather than just independent distances of query to database images. Examples are query expansion [START_REF] Chum | Total recall: Automatic query expansion with a generative feature model for object retrieval[END_REF][START_REF] Arandjelovic | Three things everyone should know to improve object retrieval[END_REF] and diffusion [START_REF] Donoser | Diffusion processes for retrieval revisited[END_REF][START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF][START_REF] Iscen | Fast spectral ranking for similarity search[END_REF]. Query expansion uses the results of initial ranking to issue a novel, enriched, query [START_REF] Chum | Total recall: Automatic query expansion with a generative feature model for object retrieval[END_REF] on-line only. Diffusion on the other hand, is based on the k-NN graph of the dataset that is constructed off-line, so that, assuming novel queries are part of the dataset, their results are essentially pre-computed. Diffusion can then be seen as infinite-order query expansion [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF].

The significance of the performance boost achieved by diffusion has been recently demonstrated at the "Large-Scale Landmark Recognition"1 challenge in conjunction with CVPR 2018. The vast majority of top-ranked teams have used query expansion or diffusion as the last step of their method.

Recently, efficient diffusion methods have been introduced to the image retrieval community. Iscen et al . [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF] apply diffusion to obtain the final ranking, in particular by solving a large and sparse system of linear equations. Even though an efficient conjugate gradient (CG) [START_REF] Hackbusch | Iterative solution of large sparse systems of equations[END_REF] solver is used, query times on large-scale datasets are in a range of several seconds. A significant speed-up is achieved by truncating the system of linear equations. Such an approximation, however, brings a slight degradation in the retrieval performance. Their method can be interpreted as graph filtering in the temporal domain.

In the recent work of Iscen et al . [START_REF] Iscen | Fast spectral ranking for similarity search[END_REF], more computation is shifted to the offline phase to accelerate the query. The solution of the linear system is estimated by low-rank approximation of the k-NN graph Laplacian. Since the eigenvectors of the Laplacian represent a Fourier basis of the graph, this is interpreted as graph filtering in the spectral domain. The price to pay is increased space complexity to store the embeddings of the database descriptors. For comparable performance, a 5k-10k dimensional vector is needed per image.

In this paper, we introduce a hybrid method that combines spectral filtering [START_REF] Iscen | Fast spectral ranking for similarity search[END_REF] and temporal filtering [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF]. This hybrid method offers a trade-off between speed (i.e., the number of iterations of CG) and the additional memory required (i.e., the dimensionality of the embedding). The two approaches [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF][START_REF] Iscen | Fast spectral ranking for similarity search[END_REF] are extreme cases of our hybrid method. We show that the proposed method pairs or outperforms the previous methods while either requiring less memory or being significantly faster -only three to five iterations of CG are necessary for embeddings of 100 to 500 dimensions.

While both temporal and spectral filtering approaches were known in other scientific fields before being successfully applied to image retrieval, to our knowledge the proposed method is novel and can be applied to other domains.

The rest of the paper is organized as follows. Related work is reviewed in Section 2. Previous work on temporal and spectral filtering is detailed in Sections 3.1 and 3.2 respectively, since the paper builds on this work. The proposed method is described in Section 4 and its behavior is analyzed in Section 5. Experimental results are provided in Section 6. Conclusions are drawn in Section 7.

Related work

Query expansion (QE) has been a standard way to improve recall of image retrieval since the work of Chum et al . [START_REF] Chum | Total recall: Automatic query expansion with a generative feature model for object retrieval[END_REF]. A variety of approaches exploit local feature matching and perform various types of verification. Such matching ranges from selective kernel matching [START_REF] Tolias | Visual query expansion with or without geometry: refining local descriptors by feature aggregation[END_REF] to geometric consensus [START_REF] Chum | Total recall: Automatic query expansion with a generative feature model for object retrieval[END_REF][START_REF] Chum | Total recall II: Query expansion revisited[END_REF][START_REF] Joly | Logo retrieval with a contrario visual query expansion[END_REF] with RANSAC-like techniques. The verified images are then used to refine the global or local image representation of a novel query.

Another family of QE methods are more generic and simply assume a global image descriptor [14-18, 7, 6]. A simple and popular one is average-QE [START_REF] Chum | Total recall: Automatic query expansion with a generative feature model for object retrieval[END_REF], recently extended to α-QE [START_REF] Radenović | Fine-tuning CNN image retrieval with no human annotation[END_REF]. At some small extra cost, recall is significantly boosted. This additional cost is restricted to the on-line query phase. This is in contrast to another family of approaches that considers an off-line pre-processing of the database. Given the nearest neighbors list for database images, QE is performed by adapting the local similarity measure [START_REF] Jégou | A contextual dissimilarity measure for accurate and efficient image search[END_REF], using reciprocity con-straints [START_REF] Delvinioti | Image retrieval with reciprocal and shared nearest neighbors[END_REF][START_REF] Qin | Hello neighbor: Accurate object retrieval with k-reciprocal nearest neighbors[END_REF] or graph-based similarity propagation [START_REF] Zhang | Query specific fusion for image retrieval[END_REF][START_REF] Donoser | Diffusion processes for retrieval revisited[END_REF][START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF]. The graph-based approaches, also known as diffusion, are shown to achieve great performance [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF] and to be a good way for feature fusion [START_REF] Zhang | Query specific fusion for image retrieval[END_REF]. Such on-line re-ranking is typically orders of magnitude more costly than simple average-QE.

The advent of CNN-based features, especially global image descriptors, made QE even more attractive. Average-QE or α-QE are easily applicable and very effective with a variety of CNN-based descriptors [START_REF] Radenović | Fine-tuning CNN image retrieval with no human annotation[END_REF][START_REF] Gordo | End-to-end learning of deep visual representations for image retrieval[END_REF][START_REF] Tolias | Particular object retrieval with integral maxpooling of cnn activations[END_REF][START_REF] Kalantidis | Cross-dimensional weighting for aggregated deep convolutional features[END_REF]. State-of-the-art performance is achieved with diffusion on global or regional descriptors [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF]. The latter is possible due to the small number of regions that are adequate to represent small objects, in contrast to thousands in the case of local features.

Diffusion based on tensor products can be attractive in terms of performance [START_REF] Bai | Regularized diffusion process on bidirectional context for object retrieval[END_REF][START_REF] Bai | Ensemble diffusion for retrieval[END_REF]. However, in this work, we focus on the page-rank like diffusion [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF][START_REF] Zhou | Ranking on data manifolds[END_REF] due to its reasonable query times. An iterative on-line solution was commonly preferred [START_REF] Donoser | Diffusion processes for retrieval revisited[END_REF] until the work of Iscen et al . [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF], who solve a linear system to speed up the process. Additional off-line pre-processing and the construction and storage of additional embeddings reduce diffusion to inner product search in the spectral ranking of Iscen et al . [START_REF] Iscen | Fast spectral ranking for similarity search[END_REF]. This work lies exactly in between these two worlds and offers a trade-off exchanging memory for speed and vice versa.

Fast random walk with restart [START_REF] Tong | Fast random walk with restart and its applications[END_REF] is very relevant in the sense that it follows the same diffusion model as [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF][START_REF] Zhou | Ranking on data manifolds[END_REF] and is a hybrid method like ours. It first disconnects the graph into distinct components through clustering and then obtains a low-rank spectral approximation of the residual error. Apart from the additional complexity, parameters etc. of the off-line clustering process and the storage of both eigenvectors and a large inverse matrix, its online phase is also complex, involving the Woodbury matrix identity and several dense matrixvector multiplications. Compared to that, we first obtain a very low-rank spectral approximation of the original graph, and then solve a sparse linear system of the residual error. Thanks to orthogonality properties, the online phase is nearly as simple as the original one and significantly faster.

Problem formulation and background

The methods we consider are based on a nearest neighbor graph of a dataset of n items, represented by n × n adjacency matrix W . The graph is undirected and weighted according to similarity: W is sparse, symmetric, nonnegative and zero-diagonal. We symmetrically normalize W as W : = D -1/2 W D -1/2 , where D : = diag(W 1) is the diagonal degree matrix, containing the row-wise sum of W on its diagonal. The eigenvalues of W lie in [-1, 1] [START_REF] Chung | Spectral graph theory[END_REF].

At query time, we are given a sparse n × 1 observation vector y, which is constructed by searching for the nearest neighbors of a query item in the dataset and setting its nonzero entries to the corresponding similarities. The problem is to obtain an n × 1 ranking vector x such that retrieved items of the dataset are ranked by decreasing order of the elements of x. Vector x should be close to y but at the same time similar items are encouraged to have similar ranks in x, essentially by exploring the graph to retrieve more items.

Temporal filtering

Given a parameter α ∈ [0, 1), define the n × n regularized Laplacian function by

L α (A) : = (I n -αA)/(1 -α) (1)
for n × n real symmetric matrix A, where I n is the n × n identity matrix. Iscen et al . [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF] then define x as the unique solution of the linear system

L α (W)x = y, (2)
which they obtain approximately in practice by a few iterations of the conjugate gradient (CG) method, since L α (W) is positive-definite. At large scale, they truncate W by keeping only the rows and columns corresponding to a fixed number of neighbors of the query, and re-normalize. Then, (2) only performs re-ranking within this neighbor set. We call this method temporal filtering because if x, y are seen as signals, then x is the result of applying a linear graph filter on y, and the iterations of CG are a set of recurrence relations that define the filter. While W is computed and stored off-line, ( 2) is solved online (at query time), and this is expensive.

Spectral filtering

Linear system (2) can be written as x = h α (W)y, where the transfer function h α is defined by

h α (A) : = L α (A) -1 = (1 -α)(I n -αA) -1 (3)
for n × n real symmetric matrix A. Given the eigenvalue decomposition W = U ΛU of the symmetric matrix W, Iscen et al . [START_REF] Iscen | Fast spectral ranking for similarity search[END_REF] observe that h α (W) = U h α (Λ)U , so that (2) can be written as

x = U h α (Λ)U y, (4)
which they approximate by keeping only the largest r eigenvalues and the corresponding eigenvectors of W. This defines a low-rank approximation h α (W) ≈ U 1 h α (Λ 1 )U 1 instead. This method is referred to as fast spectral ranking (FSR) in [START_REF] Iscen | Fast spectral ranking for similarity search[END_REF]. Crucially, Λ is a diagonal matrix, hence h α applies element-wise, as a scalar function

h α (x) : = (1 -α)/(1 -αx) for x ∈ [-1, 1].
At the expense of off-line computing and storing the n × r matrix U 1 and the r eigenvalues in Λ 1 , filtering is now computed as a sequence of low-rank matrix-vector multiplications, and the query is accelerated by orders of magnitude compared to [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF]. However, the space overhead is significant. To deal with approximation errors when r is small, a heuristic is introduced that gradually falls back to the similarity in the original space for items that are poorly represented, referred to as FSRw [START_REF] Iscen | Fast spectral ranking for similarity search[END_REF].

We call this method spectral filtering because U represents the Fourier basis of the graph and the sequence of matrix-vector multiplications from right to left in the right-hand side of (4) represents the Fourier transform of y by U , filtering in the frequency domain by h α (Λ), and finally the inverse Fourier transform to obtain x in the time domain by U .

Hybrid spectral-temporal filtering

Temporal filtering ( 2) is performed once for every new query represented by y, but W represents the dataset and is fixed. Could CG be accelerated if we had some very limited additional information on W?

On the other extreme, spectral filtering (4) needs a large number of eigenvectors and eigenvalues of W to provide a high quality approximation, but always leaves some error. Could we reduce this space requirement by allocating some additional query time to recover the approximation error?

The answer is positive to both questions and in fact these are the two extreme cases of hybrid spectral-temporal filtering, which we formulate next.

Derivation

We begin with the eigenvalue decomposition W = U ΛU , which we partition as Λ = diag(Λ 1 , Λ 2 ) and U = (U 1 U 2 ). Matrices Λ 1 and Λ 2 are diagonal r × r and (n -r) × (n -r), respectively. Matrices U 1 and U 2 are n × r and n × (nr), respectively, and have the following orthogonality properties, all due to the orthogonality of U itself:

U 1 U 1 = I r , U 2 U 2 = I n-r , U 1 U 2 = O, U 1 U 1 + U 2 U 2 = I n . (5)
Then, W is decomposed as

W = U 1 Λ 1 U 1 + U 2 Λ 2 U 2 . (6)
Similarly, h α (W) is decomposed as

h α (W) = U h α (Λ)U (7) = U 1 h α (Λ 1 )U 1 + U 2 h α (Λ 2 )U 2 , (8)
which is due to the fact that diagonal matrix h α (Λ) is obtained by element-wise application, hence decomposed as

h α (Λ) = diag(h α (Λ 1 ), h α (Λ 2 )).
Here the first term is exactly the low-rank approximation that is used by spectral filtering, and the second is the approximation error

e α (W) : = U 2 h α (Λ 2 )U 2 (9) = (1 -α) U 2 (I n-r -αΛ 2 ) -1 U 2 + U 1 U 1 -U 1 U 1 (10) = (1 -α) U 2 (I n-r -αΛ 2 )U 2 + U 1 U 1 -1 -U 1 U 1 (11) = (1 -α) I n -αU 2 Λ 2 U 2 -1 -U 1 U 1 (12) = h α (U 2 Λ 2 U 2 ) -(1 -α)U 1 U 1 . ( 13
)
We have used the definition (3) of h α in [START_REF] Hackbusch | Iterative solution of large sparse systems of equations[END_REF] and [START_REF] Joly | Logo retrieval with a contrario visual query expansion[END_REF]. Equation ( 12) is due to the orthogonality properties [START_REF] Chum | Total recall: Automatic query expansion with a generative feature model for object retrieval[END_REF]. Equation [START_REF] Tolias | Visual query expansion with or without geometry: refining local descriptors by feature aggregation[END_REF] follows from the fact that for any invertible matrices A, B of conformable sizes,

(U 1 AU 1 + U 2 BU 2 ) -1 = U 1 A -1 U 1 + U 2 B -1 U 2 , (14)
which can be verified by direct multiplication, and is also due to orthogonality. Now, combining (8), ( 13) and ( 6), we have proved the following.

Theorem 1. Assuming the definition (3) of transfer function h α and the eigenvalue decomposition (6) of the symmetrically normalized adjacency matrix W, h α (W) is decomposed as

h α (W) = U 1 g α (Λ 1 )U 1 + h α (W -U 1 Λ 1 U 1 ), ( 15
)
where

g α (A) : = h α (A) -h α (O) = (1 -α) (I n -αA) -1 -I n ( 16
)
for n × n real symmetric matrix A.

For x ∈ [-1, 1] in particular, g α (x) : = h α (x) -h α (0) = (1 -α)αx/(1 -αx).
Observe that Λ 2 , U 2 do not appear in [START_REF] Jégou | A contextual dissimilarity measure for accurate and efficient image search[END_REF] and indeed it is only the largest r eigenvalues Λ 1 and corresponding eigenvectors U 1 of W that we need to compute. The above derivation is generalized from h α to a much larger class of functions in appendix A.

Algorithm

Why is decomposition (32) of h α (W) important? Because given an observation vector y at query time, we can express the ranking vector x as

x = x s + x t , (17)
where the first, spectral, term x s is obtained by spectral filtering

x s = U 1 g α (Λ 1 )U 1 y, (18)
as in [START_REF] Iscen | Fast spectral ranking for similarity search[END_REF], where g α applies element-wise, while the second, temporal, term x t is obtained by temporal filtering, that is, solving the linear system

L α (W -U 1 Λ 1 U 1 )x t = y, (19)
which we do by a few iterations of CG as in [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF]. The latter is possible because

L α (W -U 1 Λ 1 U 1 ) is still positive-definite, like L α (W).
It's also possible without an explicit dense representation of U 1 Λ 1 U 1 because CG, like all Krylov subspace methods, only needs black-box access to the matrix A of the linear system, that is, a mapping z → Az for z ∈ R n . For system [START_REF] Tolias | Particular object retrieval with integral maxpooling of cnn activations[END_REF] in particular, according to the definition (1) of L α , we use the mapping

z → z -α Wz -U 1 Λ 1 U 1 z /(1 -α), (20)
where product Wz is efficient because W is sparse as in [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF], while U 1 Λ 1 U 1 z is efficient if computed right-to-left because U 1 is an n × r matrix with r n and Λ 1 is diagonal as in [START_REF] Iscen | Fast spectral ranking for similarity search[END_REF].

Discussion

What is there to gain from spectral-temporal decomposition (17) of x?

First, since the temporal term [START_REF] Tolias | Particular object retrieval with integral maxpooling of cnn activations[END_REF] can recover the spectral approximation error, the rank r of U 1 , Λ 1 in the spectral term [START_REF] Zhang | Query specific fusion for image retrieval[END_REF] can be chosen as small as we like. In the extreme case r = 0, the spectral term vanishes and we recover temporal filtering [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF]. This allows efficient computation of only a few eigenvectors/values, even with the Lanczos algorithm rather than the approximate method of [START_REF] Iscen | Fast spectral ranking for similarity search[END_REF]. Most importantly, it reduces significantly the space complexity, at the expense of query time. Like spectral filtering, it is possible to sparsify U 1 to compress the dataset embeddings and accelerate the queries online. In fact, we shall see that sparsification is much more efficient in our case.

Second, the matrix W -U 1 Λ 1 U 1 is effectively like W with the r largest eigenvalues removed. This improves significantly the condition number of matrix L α (W -U 1 Λ 1 U 1 ) in the temporal term [START_REF] Tolias | Particular object retrieval with integral maxpooling of cnn activations[END_REF] compared to L α (W) in the linear system (2) of temporal filtering [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF], on which the convergence rate depends. In the extreme case r = n, the temporal term vanishes and we recover spectral filtering [START_REF] Iscen | Fast spectral ranking for similarity search[END_REF]. In turn, even with small r, this reduces significantly the number of iterations needed for a given accuracy, at the expense of computing and storing U 1 , Λ 1 off-line as in [START_REF] Iscen | Fast spectral ranking for similarity search[END_REF]. The improvement is a function of α and the spectrum of W, and is quantified in section 5.

In summary, for a given desired accuracy, we can choose the rank r of the spectral term and a corresponding number of iterations of the temporal term, determining a trade-off between the space needed for the eigenvectors (and the off-line cost to obtain them) and the (online) query time. Such choice is not possible with either spectral or temporal filtering alone: at large scale, the former may need too much space and the latter may be too slow.

Analysis

How "easier" for CG is W -

U 1 Λ 1 U 1 in (19) compared to W in (2)?
In solving a linear system Ax = b where A is an n × n symmetric positivedefinite matrix, CG generates a unique sequence of iterates x i , i = 0, 1, . . . , such that the A-norm e i A of the error e i : = x * -x i is minimized over the Krylov subspace K i : = b, Ab, . . . , A i b at each iteration i, where x * : = A -1 b is the exact solution and the A-norm is defined by x A : = √

x Ax for x ∈ R n . A well-known result on the rate of convergence of CG that assumes minimal knowledge of the eigenvalues of A states that the A-norm of the error at iteration i, relative to the A-norm of the initial error e 0 : = x * , is upper-bounded by [START_REF] Trefethen | Numerical linear algebra[END_REF]

e i A e 0 A ≤ φ i (A) : = 2 κ(A) -1 κ(A) + 1 i , (21)
where κ(A) : = A A -1 = λ 1 (A)/λ n (A) is the 2-norm condition number of A, and λ j (A) for j = 1, . . . , n are the eigenvalues of A in descending order.

In our case, matrix L α (W) of linear system (2) has condition number

κ(L α (W)) = 1 -αλ n (W) 1 -αλ 1 (W) = 1 -αλ n (W) 1 -α . ( 22
)
The first equality holds because for each eigenvalue λ of W there is a corresponding eigenvalue (1 -αλ)/(1 -α) of L α (W), which is a decreasing function.

The second holds because λ 1 (W) = 1 [START_REF] Chung | Spectral graph theory[END_REF]. Now, let W r : = W -U 1 Λ 1 U 1 for r = 0, 1, . . . , n -1, where Λ 1 , U 1 represent the largest r eigenvalues and the corresponding eigenvectors of W respectively. Clearly, W r has the same eigenvalues as W except for the largest r, which are replaced by zero. That is, λ 1 (W r ) = λ r+1 (W) and λ n (W r ) = λ n (W). The latter is due to the fact that λ n (W) ≤ -1/(n -1) ≤ 0 [START_REF] Chung | Spectral graph theory[END_REF], so the new zero eigenvalues do not affect the smallest one. Then,

κ(L α (W r )) = 1 -αλ n (W) 1 -αλ r+1 (W) ≤ κ(L α (W)). ( 23
)
This last expression generalizes [START_REF] Bai | Ensemble diffusion for retrieval[END_REF]. Indeed, W = W 0 . Then, our hybrid spectraltemporal filtering involves CG on L α (W r ) for r ≥ 0, compared to the baseline temporal filtering for r = 0. The inequality in ( 23) is due to the fact that |λ j (W)| ≤ 1 for j = 1, . . . , n [START_REF] Chung | Spectral graph theory[END_REF]. Removing the largest r eigenvalues of W clearly improves (decreases) the condition number of L α (W r ) relative to L α (W). The improvement is dramatic given that α is close to 1 in practice. For α = 0.99 and λ r+1 (W) = 0.7 for instance, κ(L α (W r ))/κ(L α (W)) = 0.0326. More generally, given the eigenvalues λ r+1 (W) and λ n (W), the improvement can be estimated by measuring the upper bound φ i (L α (W r )) for different i and r. A concrete example is shown in Figure 1, where we measure the eigenvalues of the adjacency matrix W of a real dataset, remove the largest r for 0 ≤ r ≤ 300 and plot the upper bound φ i (L α (W r )) of the relative error as a function of rank r and iteration i given by ( 21) and [START_REF] Zhou | Ranking on data manifolds[END_REF]. Clearly, as more eigenvalues are removed, less CG iterations are needed to achieve the same relative error; the approximation error represented by the temporal term decreases and at the same time the linear system becomes easier to solve. Of course, iterations become more expensive as r increases; precise timings are given in section 6.

Experiments

In this section we evaluate our hybrid method on popular image retrieval benchmarks. We provide comparisons to baseline methods, analyze the trade-off between runtime complexity, memory footprint and search accuracy, and compare with the state of the art.

Experimental setup

Datasets. We use the revisited retrieval benchmark [START_REF] Radenović | Revisiting oxford and paris: Large-scale image retrieval benchmarking[END_REF] of the popular Oxford buildings [START_REF] Philbin | Object retrieval with large vocabularies and fast spatial matching[END_REF] and Paris [START_REF] Philbin | Lost in quantization: Improving particular object retrieval in large scale image databases[END_REF] datasets, referred to as ROxford and RParis, respectively. Unless otherwise specified, we evaluate using the Medium setup and always report mean Average Precision (mAP). Large-scale experiments are conducted on ROxford +R1M and RParis +R1M by adding the new 1M challenging distractor set [START_REF] Radenović | Revisiting oxford and paris: Large-scale image retrieval benchmarking[END_REF].

Image representation. We use GeM descriptors [START_REF] Radenović | Fine-tuning CNN image retrieval with no human annotation[END_REF] to represent images. We extract GeM at 3 different image scales, aggregate the 3 descriptors, and perform whitening, exactly as in [START_REF] Radenović | Fine-tuning CNN image retrieval with no human annotation[END_REF]. Finally, each image is represented by a single vector with d = 2048 dimensions, since ResNet-101 architecture is used. Baseline methods. We consider the two baseline methods described in Section 3, namely temporal and spectral filtering. Temporal filtering corresponds to solving a linear system with CG [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF] and is evaluated for different numbers of CG iterations. It is used with truncation at large scale to speed up the search [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF] and is denoted by Temporal †. Spectral filtering corresponds to FSR and its FSRw variant [START_REF] Iscen | Fast spectral ranking for similarity search[END_REF]. Both FSR variants are parametrized by the rank r of the approximation, which is equal to the dimensionality of the spectral embedding. Implementation details. Temporal ranking is performed with the implementation2 provided by Iscen et al . [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF]. The adjacency matrix is constructed by using top k = 50 reciprocal neighbors. Pairwise similarity between descriptors v and z is estimated by (v z) 3 + . Parameter α is set to 0.99, while the observation vector y includes the top 5 neighbors. The eigendecomposition is performed on the largest connected component, as in [START_REF] Iscen | Fast spectral ranking for similarity search[END_REF]. Its size is 933,412 and 934,809 for ROxford +R1M and RParis +R1M, respectively. Timings are measured with Matlab implementation on a 4-core Intel Xeon 2.60GHz CPU with 200 GB of RAM. We only report timings for the diffusion part of the ranking and exclude the initial nearest neighbor search used to construct the observation vector.

Comparisons

Comparison with baseline methods. We first compare performance, query time and required memory for temporal, spectral, and hybrid ranking. With respect to the memory, all methods store the initial descriptors, i.e. one 2048dimensional vector per image. Temporal ranking additionally stores the sparse regularized Laplacian. Spectral ranking stores for each vector an additional embedding of dimensionality equal to rank r, which is a parameter of the method. Our hybrid method stores both the Laplacian and the embedding, but with significantly lower rank r. We evaluate on ROxford +R1M and RParis +R1M with global image descriptors, which is a challenging large scale problem, i.e. large adjacency matrix, where prior methods fail or have serious drawbacks. Results are shown in Figure 2. Temporal ranking with CG is reaching saturation near 20 iterations as in [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF]. Spectral ranking is evaluated for a decomposition of rank r whose computation and required memory are reasonable and feasible on a single machine. Finally, the proposed hybrid method is evaluated for rank r = 400, which is a good compromise of speed and memory, as shown below.

Temporal [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF] Temporal † [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF] Spectral (FSRw) [START_REF] Iscen | Fast spectral ranking for similarity search[END_REF] Spectral ranking variants (FSR, FSRw) are not performing well despite requiring about 250% additional memory compared to nearest neighbor search in the original space. Compared to hybrid, more than one order of magnitude higher rank r is required for problems of this scale. Temporal ranking achieves good performance but at much more iterations and higher query times. Our hybrid solution provides a very reasonable space-time trade-off.

Runtime and memory trade-off. We report the trade-off between number of iterations and rank r, representing additional memory, in more detail in Figure 3. It is shown that the number of iterations to achieve the top mAP decreases as the rank increases. We achieve the optimal trade-off at r = 400 where we only need 5 or less iterations. Note that, the rank not only affects the memory and the query time of the spectral part in a linear manner, but the query time of the temporal part too [START_REF] Kalantidis | Cross-dimensional weighting for aggregated deep convolutional features[END_REF].

Sparsification of spectral embeddings is exploited in prior work [START_REF] Iscen | Fast spectral ranking for similarity search[END_REF]. We sparsify the embeddings of our hybrid method by setting the smallest values of U 1 to zero until a desired level of sparsity is reached. We denote this method by SHybrid. This sparse variant provides memory savings and an additional speedup due to the computations with sparse matrices. Figure 4 shows that performance loss remains small even for extreme sparsity e.g. 99%, while the results in Figure 2 show that it offers a significant speedup in the global descriptor setup. Performance-memory-speed comparison with the baselines is shown in Table 1. Our hybrid approach enjoys query times lower than those of temporal with truncation or spectral with FSRw, while at the same time achieves higher performance and requires less memory than the spectral-only approach. We summarize our achievement in terms of mAP, required memory, and query time in Figure 5. Temporal ranking achieves high performance at the cost of high query time and its truncated counterpart saves query time but sacrifices performance. Spectral ranking is not effective at this scale, while our hybrid solution achieves high performance at low query times.

Comparison with the state of the art. We present an extensive comparison with existing methods in the literature for global descriptors at small and large scale (1M distractors). We choose r = 400 and 5 iterations for our hybrid method, 20 iterations for temporal ranking, r = 2k and r = 5k for spectral ranking at small and large scale, respectively. Temporal ranking is also performed with truncation on a shortlist size of 75k images at large scale. The comparison is presented in Table 2. Our hybrid approach performs the best or second best right after the temporal one, while providing much smaller query times at a small amount of additional required memory.

Conclusions

In this work we have tested the two most successful manifold ranking methods of temporal filtering [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF] and spectral filtering [START_REF] Iscen | Fast spectral ranking for similarity search[END_REF] on the very challenging new benchmark of Oxford and Paris datasets [START_REF] Radenović | Revisiting oxford and paris: Large-scale image retrieval benchmarking[END_REF]. It is the first time that such methods are evaluated at the scale of one million images. Spectral filtering, with both its FSR and FSRw variants, fails at this scale, despite the significant space required for additional vector embeddings. It is possible that a higher rank would work, but it wouldn't be practical in terms of space. In terms of query time, temporal filtering is only practical with its truncated variant at this scale. It works pretty well in terms of performance, but the query time is still high. Our new hybrid filtering method allows for the first time to strike a reasonable balance between the two extremes. Without truncation, it outperforms temporal filtering while being significantly faster, and its memory overhead is one order of magnitude less than that of spectral filtering. Unlike spectral filtering, it is possible to extremely sparsify the dataset embeddings with only negligible drop in performance. This, together with its very low rank, makes our hybrid method even faster than spectral, despite being iterative. More importantly, while previous methods were long known in other fields before being applied to image retrieval, to our knowledge our hybrid method is novel and can apply e.g. to any field where graph signal processing applies and beyond. Our theoretical analysis shows exactly why our method works and quantifies its space-timeaccuracy trade-off using simple ideas from numerical linear algebra.

A General derivation

The derivation of our algorithm in section 4.1 applies only to the particular function (filter) h α (3). Here, as in [START_REF] Iscen | Fast spectral ranking for similarity search[END_REF], we generalize to a much larger class of functions, that is, any function h that has a series expansion

h(A) = ∞ i=0 c i A i . ( 24
)
We begin with the same eigenvalue decomposition (6) of W and, assuming that h(W) converges absolutely, its corresponding decomposition

h(W) = U 1 h(Λ 1 )U 1 + U 2 h(Λ 2 )U 2 , (25)
similar to [START_REF] Iscen | Efficient diffusion on region manifolds: Recovering small objects with compact cnn representations[END_REF], where U 1 , U 2 have the same orthogonality properties [START_REF] Chum | Total recall: Automatic query expansion with a generative feature model for object retrieval[END_REF]. Again, the first term is exactly the low-rank approximation that is used by spectral filtering, and the second is the approximation error

e α (W) : = U 2 h(Λ 2 )U 2 (26) = ∞ i=0 c i U 2 Λ i 2 U 2 (27) = ∞ i=0 c i U 2 Λ 2 U 2 i -c 0 U 1 U 1 (28) = h(U 2 Λ 2 U 2 ) -h(0)U 1 U 1 . (29)
Again, we have used the series expansion [START_REF] Tong | Fast random walk with restart and its applications[END_REF] of h in ( 27) and (29). Now, equation ( 28) is due to the fact that

(U 2 Λ 2 U 2 ) i = U 2 Λ i 2 U 2 (30)
for i ≥ 1, as can be verified by induction, while for i = 0,

U 2 Λ 0 2 U 2 = U 2 U 2 = I n -U 1 U 1 = (U 2 Λ 2 U 2 ) 0 -U 1 U 1 . (31)
In both (30) and (31) we have used the orthogonality properties [START_REF] Chum | Total recall: Automatic query expansion with a generative feature model for object retrieval[END_REF]. Finally, combining [START_REF] Chung | Spectral graph theory[END_REF], ( 29) and ( 6), we have proved the following.

Theorem 2. Assuming the series expansion [START_REF] Tong | Fast random walk with restart and its applications[END_REF] of transfer function h and the eigenvalue decomposition (6) of the symmetrically normalized adjacency matrix W, and given that h(W) converges absolutely, it is decomposed as

h(W) = U 1 g(Λ 1 )U 1 + h(W -U 1 Λ 1 U 1 ), ( 32
)
where g(A) : = h(A) -h(O) (33)

for n × n real symmetric matrix A. For h = h α and for x ∈ [-1, 1] in particular, g α (x) : = h α (x) -h α (0) = (1 -α)αx/(1 -αx).

This general derivation explains where the general definition of function g (33) is coming from in [START_REF] Delvinioti | Image retrieval with reciprocal and shared nearest neighbors[END_REF] corresponding to our treatment of h α in section 4.1.

Fig. 1 .

 1 Fig. 1. (a) In order, eigenvalues of adjacency matrix W of Oxford5k dataset of n = 5, 063 images with global GeM features by ResNet101 and k = 50 neighbors per point (see section 6). Eigenvalues on the left of vertical red line at j = 300 are the largest 300 ones, candidate for removal. (b) Contour plot of upper bound φi(Lα(Wr)) of CG's relative error as a function of rank r and iteration i for α = 0.99, illustrating the space(r)-time(i) trade-off for constant relative error.

Fig.Fig.

  Fig. mAP vs. CG iteration i and mAP vs. time for temporal, spectral, and hybrid filtering. Sparsified hybrid is used with sparsity 99%.

Fig.Fig. 5 .

 5 Fig. mAP vs. level of sparsification on our hybrid method for r = 400. Dashed horizontal line indicates no sparsification.

Table 1 .

 1 Performance, memory and query time comparison on ROxford +R1M with GeM descriptors for temporal (20 iterations), truncated temporal (20 iterations, 75k images in shortlist), spectral (r = 5k), and hybrid ranking (r = 400, 5 iterations). Hybrid ranking is sparsified by setting the 99% smallest values to 0. Reported memory excludes the initial descriptors requiring 8.2 GB. U1 is stored with double precision.

	SHybrid

Table 2 .

 2 mAP comparison with existing methods in the literature on small and large scale datasets, using Medium and Hard setup of the revisited benchmark.

		ROxford	ROxf +R1M	RParis	RParis +R1M
		Medium Hard Medium Hard Medium Hard Medium Hard
	NN-search	64.7	38.5	45.2	19.9	77.2	56.3	52.3	24.7
	α-QE [27]	67.2	40.8	49.0	24.2	80.7	61.8	58.0	31.0
	Temporal [8]	72.3	44.1	62.4	35.2	88.4	77.8	84.5	71.1
	Temporal † [8]	-	-	59.0	31.4	-	-	79.7	65.2
	FSR [9]	70.4	41.8	22.6	6.0	88.6	77.9	67.1	41.5
	FSRw [9]	65.9	34.9	42.4	20.2	88.6	77.9	67.1	41.3
	Hybrid (ours)	70.5	40.3	62.6	34.4	88.7	78.1	82.0	66.8

https://landmarkscvprw18.github.io/

https://github.com/ahmetius/diffusion-retrieval/


Pierre Allain

Laurent Guillo

Christine Guillemot

Fast light field inpainting

come

Introduction

As the capture of 4D light fields from real scenes is gaining in popularity, the need for efficient editing tools is expected to rise as well. However, the very large volume of data which they represent, as well as the need to maintain texture and structure consistency across views, raise challenging complexity issues in these processing tasks.

In this paper, we focus on light field inpainting for applications such as object removal. Although the problem of image inpainting has already been extensively studied, only a few methods in the literature address the specific case of 4D light fields. One can find some work on the related topic of multiview image inpainting, however the methods as in [START_REF] Baek | Multiview image completion with space structure propagation[END_REF], [START_REF] Thonat | Multi-view inpainting for image-based scene editing and rendering[END_REF] are designed for captures with large baselines and do not generalize well for dense light fields.

It is only recently that methods have been proposed for dense light fields inpainting. The authors in [START_REF] Chen | Light field image editing by 4D patch synthesis[END_REF] extend 2D patch-based methods to 4D patches and ensure consistency by minimizing a 4D patch bi-directional similarity measure. This method progresses patch per patch in a greedy fashion and suffers from a high computational complexity. In addition, the global consistency of the entire light field is not guaranteed. In [START_REF] Zhang | Plenopatch: Patch-based plenoptic image manipulation[END_REF], the central view is first edited using a 2D patch-based method, and the offsets between the filled patches and their best match in the known region are propagated to the other views. But the above methods suffer from high complexity and/or from angular view coherency.

Inpainting first one view (in general the center one) and then coherently propagating the inpainted texture to the other views is an interesting approach as it allows using any 2D inpainting method on the center view to have a good initialization. This is related to the general problem of edit propagation [START_REF] Jarabo | Efficient propagation of light field edits[END_REF]. In this category of methods, the approach in [START_REF] Williem | Spatio-angular consistent editing framework for 4D light field images[END_REF] inpaints the center view using a 2D patch-based method and the other views are inpainted by searching for the best matching patches in the inpainted center view. However, because of the greedy nature of the algorithm and the fact that the views are processed separately, the inpainting results may contain angular inconsistencies. The authors in [START_REF] Frigo | Epipolar plane diffusion: An efficient approach for light field editing[END_REF] exploit smoothness priors in the epipolar plane images (EPI) to coherently propagate information across views using PDE-based diffusion, while a method based on matrix completion is proposed in [START_REF] Pendu | Light field inpainting propagation via low rank matrix completion[END_REF], assuming a low rank prior for the light field data.

One critical issue for propagating texture to all views is the estimation of disparity in the masked region. In [START_REF] Goldluecke | The variational structure of disparity and regularization of 4D light fields[END_REF], the disparity map is inpainted using a linear interpolation followed by variational regularization. However, an interpolation per EPI may lead to spatial inconsistencies in the disparity maps. In [START_REF] Frigo | Epipolar plane diffusion: An efficient approach for light field editing[END_REF], dominant orientations of known pixels are computed as weighted averages of structure tensors across all the views of a given slice. Dominant orientations of unknown pixels are found using a super-pixel guided interpolation.

In this paper, we focus on the fast propagation of the inpainted texture across the entire light field. For sake of computational efficiency, we consider a simple warping strategy which requires knowledge of scene geometry in the inpainted region, that is to say the disparity at this location, as shown in Fig. 1. We revisit the problem of disparity estimation using smoothed structure tensors. We show that, while smoothing generally increases spatio-temporal coherence of estimated disparity fields, it may also lead to misalignment between color and disparity, and in turn to incoherent disparity in the masked region. A direct structure tensor computation is used instead and it is shown that, even if the disparity field is much noisier, this noise does not affect the inpainted disparity values as it turns out to be automatically removed via the color-guided interpolation. This simple approach reduces computational complexity while improving disparity coherency in the masked region. The inpainted disparity allows us to perform a simple angular warping to propagate inpainted texture across the whole light field.

In summary, the contributions of the paper are as follows: 1) A fast disparity estimation of a masked area guided by an inpainted view. 2) The ability of applying inpainting on any view of the light field with corresponding disparity at angular position, since disparity is computed exactly on the view (no averaging).

3) A direct propagation of inpainted area over the whole light field using a 2D forward warping method. 4) An overall speed increase as compared to other methods.

The paper is organized as follows: Section 3 and 4 describe the selected disparity computation method and the method for interpolating the disparity inside the inpainted area. Section 5 presents the method used for propagating a view across the light field, and in Section 6 we discuss some experimental results.

Inpainted view+mask

Interpolated disparity

Angular warping

Light field

Fig. 1: Overview of the method. The red triangle represents inpainted object in the scene. Corresponding mask standing for inpainted area is delimited by the white line. Scene geometry (depth) behind the mask is unknown. Disparity is therefore interpolated according to subapertures color features.

Disparity estimation

Let L(x, y, u, v) denote the 4D representation of a light field, describing the radiance of a light ray parameterized by its intersection with two parallel planes [START_REF] Gortler | The lumigraph[END_REF], and where (u, v) denote the angular (view) coordinates and (x, y) the spatial (pixel) coordinates. While disparity estimation methods based on structure tensors [START_REF] Wanner | Variational light field analysis for disparity estimation and super-resolution[END_REF] usually operate on EPI smoothed with a 2D Gaussian kernel, we consider instead a direct computation of the structure tensors on the 4D light field L as

T xu = ∇ xu L ⊗ ∇ xu L * G σ (x, y, u, v), (1)
where

∇ xu L = ∂ x L ∂ u L
, denotes the (x, u) gradient of L, and where G σ (x, y, u, v) is a Gaussian kernel smoothing along the four light field dimensions. The orthogonal eigenvectors θ + and θ -with respective eigenvalues λ + and λ -(where λ + > λ -) of T xu give a robust computation of the dominant gradient orientations. We are interested in the eigenvector θ -with the smallest eigenvalue which corresponds to the isophote lines whose slopes give the amount of disparity between the views. Same applies for computing the vertical disparity with the tensor T yv .

While the use of a smoothing kernel yields disparities with enhanced spatial coherence, this also leads to mis-alignment (some cross overs) between the disparity estimates and objects contours (Fig. 2). We will see in the next section that such cross overs can produce bad estimates of disparity inside the inpainted area. The structure tensor in the approach is therefore computed with no smoothing, and in order to remove most of outliers estimates, i.e. reducing the noise, we apply a Total Variation (TV-L1) denoising [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] on the disparity estimate:

D = 1 2 θ xu- -→ x θ xu- -→ u + θ yv- -→ y θ yv- -→ v , (2)
searching to optimize:

D = arg min D ||∇ D|| 1 +λ|| D -D|| 2 , ( 3
)
where λ is the total variation regularization parameter and was set to 0.5 in our experiments. This value yields a good trade-off regarding smoothness of the regularization.

We finally obtain a disparity estimate specific to inpainted view, being spatially accurate regarding scene contours, and in a very fast way. In the next section, we will address the problem of both filling inpainted mask and dealing with remaining noise.

In order to propagate the inpainted area across all views, we need a disparity information inside the mask that does not exist in the light field scene, since inpainting a view can be seen as modifying depth geometry in the inpainted area. As proposed in [START_REF] Frigo | Epipolar plane diffusion: An efficient approach for light field editing[END_REF], we use a superpixel segmentation of the inpainted view to guide disparity inpainting, making the assumption that disparity (i.e. depth) is homogeneous in local regions having the same color. The computed superpixels S = {S s }, with {s ∈ N|s ≤ Q}, are first merged such as any of them lying even partly inside the masked area has at least K pixels outside the mask, i.e. known data. The number of superpixels is defined such as Q = 1 γ P , P being the number of pixels in each subaperture. In our experiments γ was set to 1000 pixels per superpixel, which allows to capture main structures of the image, and K to 20, so that a certain area of each superpixel exists outside the mask.

The disparity D for any pixel at position x inside the mask and belonging to a superpixel S s is estimated by interpolating known disparity values D of W ≤ K closest pixels outside the mask but from the same superpixel as:

D(x) = Σ W i w s,i (x, y i ) D(y i ) Σ W i w s,i (x, y i ) , (4)
where y i ∈ S s is a pixel position outside the mask. We also set W to 20 in our experiments. The weights are defined as:

w s (x, y i ) = exp - (d s (x, y i ) -µ s (x)) 2 σ s (x) 2 , (5)
with:

d s (x, y i ) = x -y i 2 + (I(x) -I(y i )) 2 f s (6)
accounting for the spatial distance between the two pixels (inside and outside the mask) and of the color proximity of the two pixels in the inpainted view. Knowing that the spatial distance d s (x, y i ) can greatly vary for the different pixels x within a given superpixel, we apply in Eq.( 5) a centering of weights around the minimum distance as:

µ s (x) = min( {d s,i }), σ s (x) 2 = βVar({d s,i }). (7)
Hence, pixels with known disparity and closest to the border of the mask have a higher weight. The scaling between spatial and color distance in Eq.(6 is defined by a function specific to superpixel geometry as f s = α Var({x}) Var({I(x)}) , x ∈ S s to compensate for higher spatial and color variance inside some superpixels. In our experiments, we set α = 10 -5 and β = 1 which yielded a good trade-off for interpolation smoothness and coherence with outer mask disparity.

Superpixel guided interpolation has proved to be effective for disparity reconstruction [START_REF] Frigo | Epipolar plane diffusion: An efficient approach for light field editing[END_REF]. It is nonetheless very sensitive to disparity cross overs. In Fig. 3a, we can see that the disparity interpolation uses in some cases irrelevant information because of mask and segmentation configuration. In Fig. 3b, this reconstruction error is strongly alleviated because of the disparity cross over avoidance introduced in Section 2.

As for the rest of the inpainted area, interpolation weights guided by superpixels prove to be enough for obtaining an homogeneous disparity reconstruction despite noisy inputs, in the same fashion as in [START_REF] Min | Fast global image smoothing based on weighted least squares[END_REF]. Finally, we obtain a consistent disparity estimate inside the mask which is derived from the inpainted view only and not by using information from the other views, with a few computations and approximations.

Inpainting propagation

The computed disparity in the masked area is used to propagate the inpainted area of one view to all the other views of the light field using an angular forward warping. Instead of propagating information on one angular dimension in an EPI framework, we chose to warp directly in directions (u, v). Two main advantages of such approach are as follows. First, we do not need to repeat the whole pipeline of disparity computation and inpaint propagation for each EPI stripe N +1 times (for a light field of size N × N views). And second, direct warping is very fast as compared to diffusion.

However, warping has well known drawbacks. Backward warping implies to know disparity in each destination view, i.e. computing N 2 -1 disparity estimates, which is time consuming. Forward warping on the other hand only needs disparity at the origin of propagation, but comes with ghosting and cracking effects. We propose a solution to tackle the latter issues in the scope of fast light field inpainting.

Handling ghosting artefacts: Ghosting effects, which are classical in depth image-based rendering, arise when multiple pixels with different depth and color converge to a same one in the warped image, resulting in a mean value. Ghosting artifacts are often avoided by detecting depth discontinuities, in order to separate the boundary layer (containing pixels near the contours of an object) from the main layer (containing pixels far from the contours) [START_REF] Zitnick | High-quality video view interpolation using a layered representation[END_REF]. The main layer is first projected then the boundary layer is added everywhere it is visible (i.e. where its depth value is smaller than the main layers one). Ghosting artifacts can also be avoided by estimating background and foreground contributions in the rendered view with the help of advanced matting techniques [START_REF] Hasinoff | Boundary matting for view synthesis[END_REF], or using a confidence measure for each projected pixel as in [START_REF] Jantet | Joint Projection Filling method for occlusion handling in Depth-Image-Based Rendering[END_REF].

For sake of reduced complexity, we consider instead a simple disparity mean filter to classify the projected pixels within the mask area as locally belonging to the foreground or the background. This simple mean depth classification approach gives satisfactory results in the case where the number of depth layers does not exceed two. However, in the case of more complex scenes with a higher number of layers, approaches such as z-buffering would be needed to avoid ghosting artefacts. Nonetheless, one can note that for dense light fields with small disparities, as those considered here, such case remains marginal.

Let us suppose that two pixels coming from two different objects with different depths are projected on the same pixel, when doing the texture propagation. Instead of handling a z-buffer to order the projected pixels from foreground to background, thus saving computationally expensive sorting operations, we compute the mean disparity of pixels projected on the same point. We then classify pixels having a disparity higher than the mean as foreground occluding pixels, and those having a disparity smaller than the mean as occluded pixels. In order to obtain such mean disparity, we simply warp the disparity by itself, see Fig. 4. The color is then warped with pixels classified as occluding pixels overlaying those classified as occluded pixels. This simple method also prevents aliasing of warping by still averaging occluding pixels that converge to the same point. This method is directly applied in the two angular dimensions (u, v) of light field, which also provides better continuity of propagation as compared to a sequential one dimension u, v projection.

Handling cracks: Cracking effects, oppositely, arises where a destination pixel does not receive any information due to a local divergence of disparity. This issue is dealt with the help of unsupervised inpainting as in [START_REF] Telea | An image inpainting technique based on the fast marching method[END_REF] and yields satisfying results for tested light fields.

Eventually, the computational cost of inpainting propagation for one view accounts for two forward warpings and an inpainting of cracks. The later is negligible provided cracks are small, which was the case in our tests.

Computational complexity analysis

Algorithmic complexity is presented for the differents stages of the method. First, disparity estimate obtained by Eqs.(1, 2) is computed all over the view at which the inpainting is assigned. The estimate is regularized using Eq.(3). Both algorithms have O(P ) complexity. Disparity interpolation guided by inpainted view: SLIC algorithm [START_REF] Achanta | Slic superpixels compared to state-of-the-art superpixel methods[END_REF] used for superpixel segmentation has O(P ) complexity. Fast implementations of the algorithm exist, and make in practice superpixel segmentation computational cost to be low. Merging superpixels with respect to the mask has an approximate complexity O(R log(R)), with R being the number of superpixels not containing at least K pixels outside the mask. This process therefore strongly depends on image segmentation and mask characteristics, as well as on value K. Assuming R is mostly dependant on the number of superpixels Q, we obtain O(Q log(Q)) complexity for the merging step. In our experiments, for a total number of superpixels around 400 (depending on the dataset), computation took in most complex cases only a few seconds. Disparity interpolation inside the mask defined at Eq.( 4) has complexity O(M W ), with M being the number of pixels in the mask. It is, in most cases, more costly than the last two processes.

Inpainting propagation: As presented in section 4, two forward warpings of inpainted view are performed in our method. Each of them as O(N 2 M ) complexity. The second one also involve unsupervised inpainting of cracks. This last step has O(V ) complexity, whith V being the number of cracked pixels. Their number is usally low for dense light fields, therefore the cost of this step is neglectible.

Overall: The global method then has O(P ), O(Q log(Q)), O(M W ) and O(N 2 M ) complexities. In our experiments, we found that increasing Q and W does not noticeabily change inpainting results. Therefore, our method's complexity can be approximated by O(N 2 M ). However, in case a small inpainting mask is applied, disparity estimate and superpixel segmentation become dominant, and the complexity is then O(P ).

Experimental results

Fig. 5 shows an inpainting result in comparison with method [START_REF] Williem | Spatio-angular consistent editing framework for 4D light field images[END_REF]. The experiments were performed on real light fields captured by camera Lytro 1.0. To alleviate color fluctuations between views, we first perform a histogram matching between the central and remaining views. The algorithm takes as input the mask of the area to be inpainted in the whole light field and an inpainted view (usually the central one), which here has been inpainted using the patch-based method of [START_REF] Daisy | A fast spatial patch blending algorithm for artefact reduction in pattern-based image inpainting[END_REF]. We compare our method to [START_REF] Frigo | Epipolar plane diffusion: An efficient approach for light field editing[END_REF], for which we perform 100 iterations of diffusion in the reference algorithm. Computations for the later and ours were performed using a C++ implementation on a i7-6600U without any parallelization and the execution times are presented in Tab. 1.

Our method performs the fastest inpainting for almost every dataset tested. Because the mask applied on light field Bee2 has small size, method of [START_REF] Pendu | Light field inpainting propagation via low rank matrix completion[END_REF] is faster. Indeed, the singular value decomposition used has a O(N 4 M ) complexity linearily dependent on mask's size, whereas our method has O(P ) complexity in this case.

Let us note that from a memory allocation perspective, our method consumes only a few calculation images. For comparison [START_REF] Pendu | Light field inpainting propagation via low rank matrix completion[END_REF] requires approximately 8M (N 2 + 400) floating point numbers to store in memory.

More results with videos showing the inpainting consistency across views are provided at https://www.irisa.fr/temics/demos/lightField/InpaintFast/ index.html, and Fig. 6 shows additional results comparison to state of the art methods.

Experiments show very similar results regarding refocusing on focus plane of removed object and its corresponding background, while running considerably faster than other methods. 1: Execution times. Our new method performs with an average of ∼ 21 times faster than proposed in [START_REF] Frigo | Epipolar plane diffusion: An efficient approach for light field editing[END_REF], and ∼ 2.2 times than [START_REF] Pendu | Light field inpainting propagation via low rank matrix completion[END_REF]. Regarding the later, our method performs faster for larger inpainting masks. As compared to [START_REF] Williem | Spatio-angular consistent editing framework for 4D light field images[END_REF], our experiments on Lytro Illum light fields ran around ∼ 1.3 times faster, without GPU implementation.

Conclusion

In this paper we have presented a new fast approach for inpainting light fields. Our method is based on a minimum of relevant computation steps. By first obtaining a coarse disparity estimate at inpainted view coordinates in angular space, we are able to interpolate the disparity at inpainted object's location. This allows a fast propagation of inpainted view using specifically adapted warping techniques. Results quality match state of the art, while running most of the time significantly faster. As a future work, our algorithm would highly benefit of parallel programming like GPU implementation. a,e) With method proposed by [START_REF] Williem | Spatio-angular consistent editing framework for 4D light field images[END_REF]. See the incoherence of the texture in the mask area. b,f) With method proposed by [START_REF] Frigo | Epipolar plane diffusion: An efficient approach for light field editing[END_REF]. See misalignment of the edges of the central leave. c,g) With method proposed by [START_REF] Pendu | Light field inpainting propagation via low rank matrix completion[END_REF] as reference result. d,h) Using our method which result matches reference.

Fig. 2 :

 2 Fig. 2: Disparity estimate of one light field view (a), with black lines accounting for original image's segmentation of dominant shapes. (b) Using Sobel gradient and smoothing tensor with a Gaussian kernel of size 7. (c) Without gradient or tensor smoothing. (d) Without smoothing + TVL1. Despite being nosier, estimate without smoothing better fits object contours.

Fig. 3 :

 3 Fig.3: Disparity estimate of inpainted view with focus on cross over issue, white line accounting for inpainting mask boundaries. a) As in[START_REF] Frigo | Epipolar plane diffusion: An efficient approach for light field editing[END_REF] for a vertical EPI., b) Our method. Highlighted superpixel mostly use information close to the leaf due to mask configuration. Disparity retrieval is therefore extremely sensitive to estimate around object edges.

Fig. 4 :

 4 Fig. 4: Forward warping in 1D. Circles stand for image pixels. Disparity is represented in red. a) Base algorithm, b) Modified algorithm with warping of disparity by itself.

Fig. 5 :Fig. 6 :

 56 Fig. 5: Inpainting of Totoro waterfall light field. a,b) Top left view. c,d) Bottom right view. a,c) With method proposed by [6]. See the incoherence of the texture in the mask area. b-d) Using our method. Views are more coherent along angular coordinates inside the inpainting area with the proposed algorithm.

  Table

	Dataset N Resolution Mask [7]	[8]	Our method
	Bee1	9 625*434 4.4% 1'29" -	4"6
	Bee2 11 625*434 1.6% 1'17" 4"2	5"3
	Figurines 11 625*434 17.4% 3'39" 19"3	7"6
	Fruits 9 625*434 4.8% 1"38 -	4"6
	Totoro Park	7 379*379 17.2% 50" 9"2	2"5
	Totoro Waterfall	7 379*379 24% 1'2" 9"1	3"
	Tape	7 379*379 17.1% 1"12 -	2"6
	Still life 9 768*768 4.6% 4'2" 13.2" 11"9
	Butterfly 9 768*768 21.6% 6'17" 45"	19"3

This work was supported by the EU H2020 Research and Innovation Programme under grant agreement No 694122 (ERC advanced grant CLIM)


Tanguy Kerdoncuff

Thomas Galezowski

Xavier Lagrange

Mobile relay for LTE: proof of concept and performance measurements

Even with dense base station deployments, public transport users often have a low quality for mobile services. Due to the insulation of the vehicle, passengers experience little to no connectivity on their end devices and low data bit rate. We thus propose a mobile relay architecture based on two imbricated levels of LTE networks and adapted to public transport systems. This architecture is evaluated on an actual LTE test-bed, consisting in standard radio interfaces, off-theshelf terminals and the Amarisoft software suite. One of the specifics of a mobile relay is the use of a radio link between the core network and the base stations, with all the limitations of this medium. This work therefore evaluates the impact on data rate of signaling mechanisms and of tunnelling in that context.

I. INTRODUCTION

The use of wireless broadband services is rising significantly with the deployment of Long Term Evolution (LTE) networks and the generalisation of smart phones, tablet computers and other new mobile devices. People make intensive use of these devices when they are on public transport vehicles such as buses, trams, or trains. It is anticipated that, by the year 2020, it will be fairly common to have up to 50 active vehicular User Equipments (UEs) per bus and up to 300 active vehicular UEs per train [START_REF] Popovski | Scenarios, requirements and kpis for 5g mobile and wireless systems. The METIS project: Mobile and wireless communications Enablers for the Twentytwenty Information Society[END_REF].

The service quality in public transport is far from satisfactory. Vehicles are usually well shielded with coated windows, which leads to a rather high penetration loss between outdoor and in-vehicle. The penetration loss can be as high as 25 dB, and even goes up to 35 dB in Shanghai high-speed magnetic levitation train [START_REF] Liu | Position-based modeling for wireless channel on high-speed railway under a viaduct at 2.35 ghz[END_REF]. Traditionally, the UEs inside the public transport system are connected to macro base stations via wireless links, in which the penetration loss severely attenuates the signal quality and decreases the achievable data rate. Deploying mobile relays with both an outdoor antenna to communicate with the cellular network and an indoor antenna to provide a good coverage to passengers is a natural solution to provide a high Quality of Service (QoS).

Société du Grand Paris is in charge of designing and constructing the new 200km of fully automated metro lines around Paris, connecting the 3 Paris airports as well as the main suburban and innovation areas. These "Grand Paris Express" lines, interconnected with existing Paris subway, regional trains and bus lines, will carry around 2 Million passengers per day. Société du Grand Paris is seeking to provide continuous high quality telecommunication services to passengers in Grand Paris Express stations and inside trains. Thus, Société du Grand Paris is interested in developing new technologies such as mobile relays in order to serve this purpose.

Mobile relays have received a lot of attention from a theoretical point of view. In [START_REF] Sui | The potential of moving relays-a performance analysis[END_REF], the authors studied the performance of dual-hop transmission with a mobile relay and with a fixed relay by means of a theoretical analysis and proved the benefit of mobile relays when the penetration loss is above 25 dB. In [START_REF] Chen | Performance analysis of cellular networks with mobile relays under different modes[END_REF], stochastic analysis is used to compute the capacity of a cell when mobile relays are deployed. In [START_REF] Tang | Coverage performance of joint transmission for moving relay enabled cellular networks in dense urban scenarios[END_REF] joint transmission by both a mobile relay and a macro-cell base station is proposed to improve the data rate of non-vehicular UEs and the outage probability is derived. All these publications focus on the transmission issues and have a theoretical approach.

Some studies consider the mobility protocol issues. In [START_REF] Chen | Mobile relay in lte-advanced systems[END_REF], several architectures based on the Third Generation Partnership Project (3GPP) analysis [START_REF]3rd Generation Partnership Project. Evolved Universal Terrestrial Radio Access (E-UTRA); Study on mobile relay[END_REF] are described. In [START_REF] Meng-Shiuan Pan | An enhanced handover scheme for mobile relays in lte-a high-speed rail networks[END_REF], [START_REF] Huang | Mobile relay based fast handover scheme in high-speed mobile environment[END_REF], the authors optimised handover for mobile relays deployed on high speed trains by using the predictability of the movement on a railway. In [START_REF] Pietsch | Moving Relays and Mobility aspects[END_REF], the authors proposed to define "mobility areas" for mobile relays and analysed the handover latency when the mobile relay stays in the same mobility area or moves to another one. However, this proposal requires extra planning effort on the operator side and some upgrades of the standard 3GPP based protocols. To the best of our knowledge there is no paper about the performance of a system with a mobile relay with standard LTE interfaces.

3GPP investigated several possible architectures for mobile relays in [START_REF]3rd Generation Partnership Project. Evolved Universal Terrestrial Radio Access (E-UTRA); Study on mobile relay[END_REF]. Some solutions require the modification of several protocols and thus an important standardisation effort. However, the first alternative (called alt1) described in [START_REF]3rd Generation Partnership Project. Evolved Universal Terrestrial Radio Access (E-UTRA); Study on mobile relay[END_REF] is based on the fact that nodes of the Evolved Packet Core (EPC) such as eNBm, MME just need IP connectivity and that IP connectivity can be simply provided by an LTE/EPC network. The objective of this paper is to demonstrate that mobile relaying compliant with alt1 of 3GPP can be easily implemented with standard EPC nodes and to measure the performance of such a solution on a test-bed with real radio transmissions.

The paper is organised as follows. In section II we describe the architecture and the protocol stack when mobile relays are used. In section III we present our test-bed. In section IV, we measure the performance for a web-browsing service and in section V, we study the impact of a mobile-relay handover on

II. MOBILE RELAY ARCHITECTURE

An LTE/EPC mainly includes the following nodes: the Evolved Node B (eNB), the Mobile Management Entity (MME), Serving Gateway (SGW), Packet Gateway (PGW), and the Home Subscriber Server (HSS).

The wireless access for passengers is provided by the use of two LTE/EPC networks as shown in Figure 1. The first one is called the Track Evolved Packet Core (EPCt) and can be for example deployed by the transport operator. Nodes of this network are identified by suffix t (eNBt, MMEt, SGWt, etc.). We name Transport Operator terminal Equipment (TOE) a terminal of this network. There is typically one TOE in every train. The second network is the general network. It manages the UEs of passengers and is owned by the usual telecommunication operator. It has Mobile eNB (eNBm) as well as standard eNB. Note that eNBms have access to the general network through the track network: an eNBm is connected to a TOE and can thus exchange messages with Global Mobile Management Entity (MMEg) and Global Serving Gateway (SGWg).

It is worth mentioning that all these elements are standard LTE components, and that the use of this architecture requires no adaptation of protocol stacks or hardware.

There can be several general networks that use the same LTE/EPC track network in order to allow all mobile operators of a country to provide the service on board. There can be several eNBms (one per operator) but Radio Access Network (RAN) sharing can also be used to avoid a duplication of embedded hardware in trains. Note that Figure 1 shows a functional architecture. The operator of the track network can be a mobile operator and several entities can be physically in the same node (e.g. MMEt and MMEg are in the same piece of equipment)

The protocol stack when a mobile relay is used is shown in Figure 2 on the radio interface and on the S1-interface. Due to the use of two LTE/EPC networks there is an additional level of encapsulation and thus extra-headers compared to a standard architecture. The following headers are added to each packet: 8 bytes for GPRS Tunneling Protocol (GTP) [11] The total length of this extra overhead is fixed but the length of packets is variable: an important question is to determine whether this 36-byte overhead is negligible or not compared to the user data average size for a popular service as web browsing.

One important benefit of mobile relay is that mobile management is no more done on an individual terminal basis but at the relay level. When a terminal is attached to the network through the relay, the terminal is considered as being under the coverage of the mobile eNB and is fixed as long as the user stays in the train. If the train moves to another cell, a single handover command is sent to the TOE but applies for all terminals connected to the mobile eNB. An interest of performing only one handover per train, instead of one handover per passenger can be found in the messages exchanged during a handover. For example, the message UE-CapabilitiesInformation describes the LTE capabilities of an UE and can be several 100-bytes long. Sending this message only once for the TOE (instead of as many times as the number of active terminals in the train) is a valuable argument in favour of a mobile relay. However, such a TOE handover still causes an interruption in the user data transmission because soft handover is not available in LTE/EPC. An important problem is to know whether such an interruption has a major impact on the TCP behaviour.

III. AMARISOFT LTE TESTBED

Our test-bed is based on the Amarisoft solution [START_REF]Amari lte 100, software lte base station on pc[END_REF] which is a low-cost LTE base station and EPC core network running on a personal computer. Digital signal processing and the whole protocol stacks are pure software operations. We fixed the bandwidth of each radio interface to 5-MHz in order to limit the processing load of the computer. The transmission is made in frequency division duplex mode in the 2.6 GHz band. In order to avoid interference between the access and the backhaul links, two separate carriers are used. In order to perform Relay Node (RN) handovers, backhaul computer 2 is required. It runs an instance of the amarisoft lteenb program to control its eNBt. That eNB connects by ethernet to the MMEt running on the first backhaul computer.

IV. DATA TRANSFER PERFORMANCES

This section presents the performance achieved by UEs using our test-bed while loading web pages. Every year the french regulatory authority ARCEP measures the QoS of operational mobile networks. We based our tests on the 26 different web pages used in the 2016 ARCEP campaign [START_REF]Enquete sur la qualité des services de données des opérateurs mobiles en france métropolitaine[END_REF]. The objective of this experiment is to check that web browsing is possible in a relay-based architecture with acceptable quality of service in non-overloaded conditions. A secondary objective is to measure the overhead due to multiple encapsulations and the amount of signaling on a practical implementation of the protocols.

The results presented in this section focus on four websites of the aforementioned ARCEP campaign, each one of them being opened by the UE's browser. The numbers we present in this work were obtained by realising 400 such occurrences for each site through our mobile relay test-bed, and 400 more times through a direct LTE connection that consists in a single eNB and EPC and that is called the "baseline".

Figure 4a presents the time it took for the UEs to retrieve 80% of the traffic observed after opening the web page. This avoids taking into account periodic traffic generated by a page's asynchronous elements such as ads. Comparing the distribution of direct LTE to the mobile relay setup gives us an estimation of the few additional hundred milliseconds it took the UE to load most of a page's elements, and by doing so, make it readable for the end user. Mobile relays being aimed at areas where such a direct LTE coverage is physically impossible, we believe this additional delay to be acceptable.

As explained in section II (see Figure 2), there is an additional tunnelling process when mobile relays are used. In most With mobile relays, some signaling messages between the eNBm and the Global Evolved Packet Core (EPCg) are transmitted on the backhaul radio link. However, figure 4c shows that these signaling messages are negligible compared to the amount of user data: less than 0.10% in most case. The difference between the first website and the other ones, is that it maintains inactive TCP connections after the page finished loading. Closing the browser, and therefore the TCP connections at the end of each experiment reactivates the sleeping LTE link, and therefore generates additional signaling.

In our experiments, all radio links are in line-of-sight conditions. Thus, our baseline does not represent the QoS experimented by users in a transport system. These results should not be read as a decrease of the quality of service when introducing a mobile relay, since without it there can be no or very degraded service. It should, on another hand, be read as the additional traffic to take into account when using this architecture to bring cellular coverage to new, previously uncovered, areas. For example, if a standard base station can accept a 6% increase in traffic due to encapsulation and signaling, it is then acceptable to use a mobile relay.

V. RELAY NODE HANDOVER

This section studies handovers of the TOE, which happens when the train moves between cells covered by different eNBt. This experiment allowed us to confirm that only one handover operation is executed and to check whether the temporary loss of connectivity and possible link quality change impacts the on board UEs.

To ensure the UEs were in connected state (in the LTE meaning of the word) during the experiment, an HTTP down link transfer is generated from the test controller. At a given time, the source eNBt instructs the TOE to switch to a We define the handover duration from the eNBms standpoint as the duration between the sending of the LTE connection reconfiguration message on the source eNBt and the moment from which a packet can be transmitted by the relay to target eNBt. Figures 5 presents the Cumulative Distribution Function (CDF) of such handover duration. The maximal handover duration's is around 230 ms, a value that remains small compared to, for example, the initial TCP re-transmission timer (1s) [START_REF] Paxson | Computing tcp's retransmission timer[END_REF].

Another important aspect of this experiment is the confirmation that such handovers did not interfere in an unforeseen way with the train's passengers' experience : UEs are able to maintain traffic and LTE attachment without any perturbation.

Figure 6 presents the impact a TOE handover may have on an UE performing traffic. Its X axis marks time while the Y axis presents the cumulative amount of bytes the UE was able to download. The slope gives the throughput. Around t = 6s, we observe an area of non progress delimited by the Figure 7 presents the CDF of the interruptions observed in UE traffic. In each case, the duration taken into account is the longest one observed during an experiment featuring a single RN handover. In 80% of the performed experiments, the longest interruption is less than 200 ms, while the maximal observed duration is slightly over 300 ms. It is important to emphasise that this graph needs to be read as an upper bound of the handover interruption duration, as other interruptions can happen independently from the handovers, as observable in Figure 6. Still, these values are in the range of interruptions that most mobile applications already deal with by, for example, content caching on the UE.

VI. CONCLUSION

In this work, we proved that deploying mobile relays is possible in a standard LTE architecture. Using a testbed with off-the-shelf UEs, we evaluated the impact the different LTE mechanisms can have on a web service when cohabiting in a mobile relay setup. We found out that neither encapsulation nor signaling on a radio link was prohibitive. We also studied mobility, in particular the mobility of the eNBm to which the associated UEs are blind. We observed that, on one hand, such handovers were not disruptive, and on another hand, that they occurred under a reasonable amount of time.

In this test, only two terminals were used. A perspective is to study the behaviour of such a system when a train is full of passengers. Some additional simulations would be needed to test a full charged vehicle in a close to real field conditions.

Fig. 1 :

 1 Fig. 1: Architecture of a network with a mobile relay

Fig. 2 :

 2 Fig. 2: Protocol stack with a mobile relay

Fig. 3 :

 3 Fig. 3: Testbed Architecture

Fig. 4 :

 4 Fig. 4: HTTP Transfer analysis

Fig. 5 :

 5 Fig. 5: Cumulative Distribution Function (CDF) of RN Handover duration from the eNBm standpoint

Fig. 6 :Fig. 7 :

 67 Fig. 6: Evolution of transmission and RN Handover impact on UE traffic (the vertical lines show the beginning and the end of the handover)