Anne-Laure Chauve

Jean-Philippe Pons

Jean-Yves Audibert

Renaud Keriven

Transductive Segmentation of Textured Meshes

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

Introduction

The generalization of digital cameras, the increase in computational power brought by graphical processors and the recent progress in multi-view reconstruction algorithms allow to create numerous and costless textured 3D models from digital photographs. In this work, we address the problem of segmenting a textured mesh into objects or object classes. The segmentation is an essential step of scene analysis, and can be used for semantic enrichment of architectural scenes, reverse engineering, subsequent recognition of known rigid objects... A meaningful mesh decomposition enables to retrieve objects of interest or remove undesired parts (see Fig. 4). This problem raises several interesting challenges: the variability of scenes and object types to handle; the subjectivity of a meaningful segmentation, which depends on the application; the simultaneous classification and segmentation involved by object detection. We took up theses challenges by designing a flexible and interactive framework that conducts collective classification of the mesh.

We propose an easy-to-use, graphical tool to provide labeled training seeds by drawing sketches on the mesh, select appropriate features among the set of available ones, compute the segmentation via our transductive learning algorithm based on Support Vector Machine classification and Laplacian graph regularization, visualize the resulting segmented mesh and refine the training sketches to rerun the segmentation if necessary.

Our work is related to three main research topics: mesh segmentation, collective point cloud classification and transductive image segmentation. Even if Fig. 1. Segmentation of the mesh into four classes: roof, wall, windows edges, cornice. Left: the input textured mesh with user supplied sketches. Right: the resulting segmentation using our algorithm.

we consider meshes in this paper, the type of data and the segmentation technique turn out to be quite different from the ones considered in the literature. Our issues are actually closer to a series of works on point cloud segmentation, however these are not directly applicable to our meshes whose connectivity and texture are of crucial importance. Our algorithm is rather based on a transductive segmentation method that was developed at first for 2D images.

The problem of mesh segmentation has become an important issue in various computer graphics applications, like parameterization and texture mapping, metamorphosis, 3D shape retrieval or modeling by example. Several segmentation algorithms for mesh partitioning have been compared in [START_REF] Attene | Mesh Segmentation -A Comparative Study[END_REF]. These algorithms deal with non-textured meshes of a single object, that are usually highquality meshes coming from CAD models or dense scans -thus very different from the image-based meshes of entire scenes we process. These algorithms fall into two main categories. The first one gathers geometric approaches, where the mesh is segmented into patches fitted with simple mathematical surfaces. The typical application of these methods is reverse engineering of CAD models. The second one is rather semantic-oriented: the aim is to decompose meshes of "natural" objects (e.g., a body model) into "meaningful" pieces (e.g., the head, two arms, two legs and the torso). However, these semantic approaches do not involve any learning-based or classification procedure, and do not consider the similarities between different, non connected parts. Mesh decomposition is usually a preprocessing step, and should thus be able to handle various types of input models and of target applications. In order to deal with varying application-dependent requirements and with the subjectivity of a meaningful segmentation, [START_REF] Wu | Ma: A sketch-based interactive framework for real-time mesh segmentation[END_REF] proposed an interactive framework similar to ours: the user draws sketches on the mesh that provide seeds for the algorithm. Nevertheless, as mentioned above and unlike [START_REF] Wu | Ma: A sketch-based interactive framework for real-time mesh segmentation[END_REF], we process textured meshes: this greatly reduces the required amount of user interaction (see for example Fig. 1).

A series of works [START_REF] Anguelov | Discriminative Learning of Markov Random Fields for Segmentation of 3D Scan Data[END_REF][START_REF] Triebel | Robust 3D Scan Point Classification using Associative Markov Networks[END_REF][START_REF] Triebel | Instance-based AMN Classification for Improved Object Recognition in 2D and 3D Laser Raneg Data[END_REF][START_REF] Munoz | Directional Associative Markov Network for 3-D Point Cloud Classification[END_REF] has been carried out on the problems of segmenting scan data into objects or objects classes using Markov networks. This problem differs from the one we address mainly on the type of data to process: textured meshes possess supplementary attributes like color or mesh connectivity that we exploit, while 3D point clouds are far denser, and more precise, thus targeting distinct applications and requiring different processes. However these works share several characteristics with ours. One is the use of collective classification through graph-based methods: instead of classifying each point or facet independently, the problem is thought of as a global classification task and adjacency relationships are used to enforce spatial contiguity of the labels on the graph. Our graph Laplacian transduction performs such a collective classification. Besides, in both problems, the ability to handle various type of scenes as well as a certain variability inside a given class of objects is crucial: this is achieved by taking advantage of various kinds of features. Thus, the Markov random field segmentation algorithm of [START_REF] Anguelov | Discriminative Learning of Markov Random Fields for Segmentation of 3D Scan Data[END_REF] has been tested successfully on both outdoor and indoor scenes, for real-world and synthetic scan datasets. In our framework, textured meshes possess both geometric and photometric attributes that should be chosen according to the scene type and jointly exploited. We combine these features in kernel weights of a graph Laplacian regularizer; the inherent modularity of kernel methods hence provides the desired flexibility.

Recently, several works have addressed the image segmentation problem in an interactive framework: a set of seed pixels representative of each region to be segmented is specified by the user, and the segmentation of the entire image is performed consistently with the seeds. The existing algorithms rely on computing weighted geodesic distances [START_REF] Bai | A Geodesic Framework for Fast Interactive Image and Video Segmentation and Matting[END_REF], graph cuts with discontinuity penalization [START_REF] Boykov | Interactive Graph Cuts for Optimal Boundary and Region Segmentation of Objects in N-D Images[END_REF], graph cuts with Gaussian mixture modelling of the segmented regions [START_REF] Blake | Interactive Image Segmentation Using an Adaptive GMMRF Model[END_REF], random walks on a graph and its relation to electrical resistance in a circuit [START_REF] Grady | Random Walks for Image Segmentation[END_REF][START_REF] Kim | Generative Image Segmentation Using Random Walks with Restart[END_REF] and transductive learning [START_REF] Duchenne | Segmentation by transduction[END_REF]. Due to its simplicity and its effective results obtained in the Microsoft GrabCut benchmark, the latter approach is adopted in our mesh segmentation algorithm. The segmentation results are visualized in our interactive framework, allowing for subsequent refinement of the query.

Outline. The rest of the paper is organized as follows. Section 2 presents the Laplacian regularizer on a graph to perform transduction and introduces the energy we minimize in the sequel; section 3 details our algorithm for the segmentation of textured meshes, and section 4 presents our experimental results. We conclude in section 5 with a brief discussion.

Graph Laplacian Transduction

Transductive Inference

The problem we are concerned with is a supervised classification task: given a set of labeled examples (called training set), we want to infer the labels of new input points (called test set). More specifically, we consider an input space X and an output space Y (typically Y = {0, . . . , K} for a classification problem); given a set of input-output couples {(x 1 , y 1 ), . . . , (x n , y n )}, we want to determine the label y of a new point x.

There are two different approaches for this problem, the inductive and transductive settings. Inductive inference is a two-step process: one tries first to learn a function f : X → Y to map the entire input space to the output space, and then, for every point x of the test set, predicts y = f (x) as its label. In contrast, in transductive inference, the test inputs are known beforehand. Thus, no general input-output mapping is inferred, but only the labels of the test points are predicted. This approach follows Vapnik's principle: "when solving a problem of interest, do not solve a more general problem as an intermediate step", and takes advantage of the unlabeled data to get an idea of the input space distribution. Indeed, transduction relies on this second principle, referred to as the smoothness assumption: "outputs vary a lot only on input regions having low density", or, in other words, "the decision boundary should lie in a low-density region".

Transductive inference can be compared to semi-supervised learning (SSL), where the training set is made up of both labeled and unlabeled data points. Like transduction, SSL utilizes the unlabeled points to infer the input distribution. However, in contrast to transduction, the test points are not (necessarily) known beforehand and semi-supervised algorithms can handle unseen data, thus being part of an inductive framework.

A large number of algorithms that have been proposed in the last few years for transductive inference relies on graph-based methods (see [START_REF] Belkin | Semi-Supervised Learning on Riemannian Manifolds[END_REF], [START_REF]Semi-Supervised Learning[END_REF] and references within), where nodes represent data points and edges encode similarities between them. The graph structure is used to propagate information from the known labels to the unlabeled points. We use the Laplacian graph regularizer of [START_REF] Duchenne | Segmentation by transduction[END_REF] with an unnormalized kernel (case s = 2, λ = 0): the labeling of the facets is carried out as the minimization of a quadratic cost function derived from the graph. The work [START_REF] Bengio | Label Propagation and Quadratic Criterion[END_REF] studies several cost functions for regularization on a graph and explains the links with label propagation algorithms.

Graph Laplacian Regularization and Energy Minimization

In the first place, we consider a binary classification problem, where the two classes have respectively 0 and 1 as labels (see 2.3 for generalization to the multi-class problem); 0 is the background class and 1 is the object class. Instead of directly predicting the labels of the test points (y ∈ {0, 1}), we consider a real-valued output space (y ∈ R) and assign each point to the class 1 y≥1/2 .

Graph Laplacian

The geometry of the data is represented by a graph G = (V, E) where the nodes V = {X 1 , . . . , X n } correspond to the input points coming from both the p labeled instances {X 1 , . . . , X p } of the training set, and the np unlabeled data {X p+1 , . . . , X n } of the test set, and the edges E represent similarities between them, in the form of a weight matrix W (of size n × n). The coefficients of this matrix (also called affinity or adjacency matrix) must satisfy:

-W ij = 0 if X i and X j are not "neighbours" (i.e. are not connected by an edge),

-W ij ≥ 0, -W ij = W ji .
The meaning of the neighbourhood needs to be specified and depends on the nature of the data (e.g., facet adjacency on a mesh, symmetrized k-nearest neighbour for a point cloud, etc.). We write the weights in the form

W ij = k(X i , X j )
where k : X × X → R is a symmetric positive function giving the similarity between two input points. A typical example is the Gaussian kernel

k(x, x ′ ) = exp -d(x,x ′ ) 2 2σ 2
with d a distance on the input space X and σ a positive parameter.

Let D be the diagonal matrix D ii = j W ij . The matrix L = D -W is called the unnormalized graph Laplacian. It is a discrete analogue of the Laplace-Beltrami operator on a Riemannian manifold (see ( 3)).

Energy design The goal is to find a labeling Y = (y 1 , . . . , y n ) of the graph that is consistent with both the labeled training points and the geometry of the entire data (represented by the graph structure).

We partition the vector Y and the matrices W , D and L according to their labeled and unlabeled parts :

Y = Y ℓ Y u W = W ℓℓ W ℓu W uℓ W uu D = D ℓℓ 0 0 D uu L = L ℓℓ L ℓu L uℓ L uu .
The initial labels of the training points are denoted Y 0 ℓ = (y 0 1 , . . . , y 0 p ). We want to find the optimal labeling Ŷ minimizing the following cost criterion, built up three different terms :

Ŷ = arg min Y ∈R n c p i=1 (y i -y 0 i ) 2 (α) + 1 2 n i,j=1 W ij (y i -y j ) 2 (β) + λ n i=p+1 (y i -s i ) 2 (γ) , (1)
where the scores s i in (γ) are defined hereafter.

(α) Consistency with the initial labeling:

p i=1 (y i -y 0 i ) 2 = Y ℓ -Y 0 l 2 .
The parameter c ∈ [0, +∞] expresses the confidence assigned to the training outputs (it could be different for each point, (c i ) i=1,...,p ).

(β) Consistency with the geometry of the data: This term penalizes rapid changes in Ŷ between points that are close on the graph (as given by the similarity matrix W ). It enforces the smoothness assumption along the graph.

1 2 n i,j=1 W ij (y i -y j ) 2 = 1 2 2 n i=1 y 2 i n j=1 W ij -2 n i,j=1 W ij y i y j = Y ⊤ (D -W )Y = Y ⊤ LY .
The latter term can be interpreted as a graph Laplacian regularizer in the following sense. Assume that the inputs are real vectors admitting a distribution with density p (with respect to the Lebesgue measure). Let ϕ be a smooth function such that y i = ϕ(X i ). From [START_REF] Duchenne | Segmentation by transduction[END_REF], the quantity Y T LY is a discrete approximation, up to a constant multiplicative factor, of the integral

(∇ϕ)(x) 2 p 2 (x)dx , (2)
which is equal to

ϕ(x)(∆ϕ)(x)p 2 (x)dx , (3)
for ∆ a (weighted) Laplace-Beltrami operator. From (2), we see that the regularization term Y T LY will be small only when the labels vary in low density regions of the input space, thus fulfilling the principle stated in Section 2.1.

(γ) Consistency with some additional knowledge:

n i=p+1 (y i -s i ) 2 = Y u -0 u 2 .
This terms allows to incorporate either some prior information or the output of another algorithm in the form of a score s i , measuring for each test point i the "likelihood" that it belongs to the object class. This term can be seen as an initializing score. We note S 0 u the (np)-vector (s p+1 , . . . , s n ).

The optimization problem (1) can hence be expressed in the following matrix form: Ŷ = arg min

Y ∈R n c Y ℓ -Y 0 ℓ 2 + Y ⊤ LY + λ Y u -S 0 u 2 . (4)
Sparse linear system giving the predicted labels Since we assume that the user-supplied seeds are trustworthy, we constrain the labels on the labeled data ( Ŷℓ = Y 0 ℓ ) and thus consider an infinite regularization coefficient c = +∞. Hence the minimization is carried over the labeling Y u of the test points, and the optimization problem rewrites Ŷu = arg min

Yu∈R n-p Y ⊤ LY + λ Y u -S 0 u 2 = arg min Yu∈R n-p Y 0 ℓ ⊤ L ℓℓ Y 0 ℓ + Y u ⊤ L uℓ Y 0 ℓ + Y 0 ℓ ⊤ L ℓu Y u + Y u ⊤ L uu Y u + λ Y u ⊤ Y u -Y u ⊤ S 0 u -S 0 u ⊤ Y u + S 0 u ⊤ S 0 u = arg min Yu∈R n-p 2Y ⊤ u L uℓ Y 0 ℓ -λS 0 u + Y ⊤ u (L uu + λI) Y u .
In order to minimize the cost criterion, we compute its derivative with respect to Y u . Let A = L uℓ Y 0 ℓ -λS 0 u and B = L uu +λI. The matrix B is symmetric positive definite matrix when λ > 0, since L uu is symmetric positive semi-definite:

Y u ⊤ L uu Y u = 1 2 n i,j=p+1
W ij (y iy j ) 2 ≥ 0 .

Thus the function

f : Y u → 2Y ⊤ u A + Y ⊤
u BY u is strictly convex and admits a unique minimum where its derivative equals 0:

∂f (Y u ) ∂Y u = 2A + 2BY u = 0 ⇐⇒ BY u = -A .
Hence, the optimal labeling Ŷu is the solution of the sparse linear system

(L uu + λI) Ŷu = λS 0 u -L uℓ Y 0 ℓ . (5)

Multi-class Segmentation

The extension of the previous algorithm to the multi-class case is straightforward using a one-versus-all approach. If there are d different classes, we resolve the linear system (5) for each class k against all other classes as background, thus having as initial label for a training point i, Y 0,k ℓ (i) = 1 if i is labeled k and Y 0,k ℓ (i) = 0 otherwise. We obtain d output vectors Ŷ k u , and assign the test point j to the class arg max k=1,...,d Ŷ k u (j).

Segmentation of Textured Meshes

Our Algorithm

We work on textured meshes of a 3D scene built from a set of calibrated images: we use the multi-view stereovision algorithm from [START_REF] Labatut | Efficient Multi-View Reconstruction of Large-Scale Scenes using Interest Points, Delaunay Triangulation and Graph Cuts[END_REF][START_REF] Pons | Multi-view Stereo Reconstruction and Scene Flow Estimation with a Global Image-Based Matching Score[END_REF][START_REF] Vu | Towards high-resolution large-scale multi-view stereo[END_REF] to reconstruct a 3D model of the scene and the multi-band blending algorithm from [START_REF] Allène | Seamless Image-Based Texture Atlases using Multi-band Blending[END_REF] to compute a texture atlas with minimal color discontinuities or blurring. The meshes presented in our experiments have been reconstructed from datasets provided by C. Strecha et al. [START_REF] Strecha | On Benchmarking Camera Calibration and Multi-view Stereo for High Resolution Imagery[END_REF] (castle-P19, castle-P30 and Herz-Jesu-25).

We want to classify the facets of the mesh given seed facets for each class. The seeds are provided by the user who draws sketches on the mesh. Thus we consider a graph G with a node for every facet of the mesh and an edge between any two adjacent facets. The modularity of our kernel method allows to chose various edges weights, depending on the scene type. We use a Gaussian kernel with a distance between facets constructed from one or several features characterizing the facets (see section 3.2 for more details and examples of such features).

The training points are the user-supplied seeds and we use an SVM classification on these training points as additional knowledge on the test points (see section 3.3). The linear system (5) of the transductive classification is sparse due to the sparseness of the facet adjacency relationship on the mesh. We solve it with a conjugate gradient algorithm (we use the IML++ implementation of [START_REF] Dongarra | A Sparse Matrix Library in C++ for High Performance Architectures[END_REF]). For each facet, we obtain a score for each class: the classification is given by the argmax of these scores (see section 2.3) but they carry much more information that can be exploited.

The interactive framework of our method allows to add supplementary seeds depending on misclassified facets and rerun the algorithm in order to improve the segmentation result or add new classes (see Fig. 4).

Features and Kernels

We compare two different facets through a set of attributes extracted from the mesh. We describe each facet by a feature vector of length m. These features are chosen according to the scene type and the discriminant characteristics of such scenes. A major element of this algorithm is the ability to take advantage of several different kinds of features, mixing photometric and geometric informations or any other available attribute (see Fig. 2).

Here are some examples of features that we used in the results of Section 4 or that could be used in other experiments. On the one hand, we use photometric features like the mean color of the facet or components of the mean color (for instance, without luminance to account for changes in illumination between cameras). On the other hand, we use geometric features like the normal of the facet (or solely its vertical component), the position of the center of the facet (for instance, height of the center), or the discrete curvature of the mesh. We could also also evaluate less local features that would consequently be smoother and more robust by averaging the previous features on a neighbourhood of the facet.

We do not compare globally the descriptor vectors of the facets, but feature by feature. Then we combine the component-wise distances in the kernel weight by multiplying the Gaussian kernels associated with each feature: if there are n f different features (each feature being a vector of variable length) and we note (f 1 i , . . . , f

n f i ) and (f 1 j , . . . , f n f j ) the feature vectors of facets X i and X j , W ij = k(X i , X j ) = n f k=1 exp - f k i -f k j 2 2σ 2 k = exp - n f k=1 f k i -f k j 2 2σ 2 k .
The kernel weight is thus parameterized by the n f values (σ 1 , . . . , σ n f ) which determine the trade-off between the various features.

Additional Knowledge on the Test Points

We train Support Vector Machines with the user-supplied training points in order to provide initialization information on the test points. We obtain a score s k i for each facet i and each class k with a one-versus-all SVM algorithm. We use the LibSVM implementation of the C-SV Classification (cf. [START_REF] Chang | LIBSVM: a library for support vector machines[END_REF]).

The incorporation of external knowledge in addition to the graph Laplacian regularization is essential in order to detect several objects of the same class which form several connected components. Indeed, each seed can generate at most one connected component in the segmentation. Besides, the initialization with the SVM scores accelerate the label propagation on the mesh and the convergence of the iterative resolution compared to the transductive segmentation based solely on graph Laplacian.

We can consider either the same kernel as for the graph Laplacian weights or a different one (different features, kernel type, or parameters), if we want to exploit distinct properties of the mesh.We learn the parameters of the SVM (parameters of the kernel plus parameter C of the SVM) in a cross-validation process. Hence, the tradeoff between the various features combined in the kernel is learned automatically. Moreover, we can employ these selected parameters in the Laplacian kernel for transduction (if we use the same kernel).

The Support Vector Machine approach is a classification method by itself and provides a comparison basis for our algorithm (just as the transductive segmentation without the SVM initialization). However it does not enforce any spatial contiguity on the mesh, and it produces a noisy result (see Fig. 3).

Note that we learn this term directly on the scene of interest, but we could incorporate as well a prior learned from other, previously segmented scenes, instead or in addition to this one. In this case, however, the algorithm is not anymore a purely transductive one.

Computational Complexity and Time

The computational cost of the algorithm can be split up into the cost of the SVM initialization and the cost of the transductive segmentation on the graph. Since in our framework the number of training points p is negligible compared to the total number of points n, we consider the complexity of the algorithm with respect to n only. In the SVM algorithm, the cost of the test prevails over the cost of the train (which is between O(p 2 ) and O(p 3 )), and its complexity is O(n.p) = O(n). In the energy minimization of the transduction, the main cost comes from solving the sparse linear system (5); this can be done in O(k.m) where k is the number of iterations of the conjugate gradient and m is the number of non-zero entries in the matrix, which is equal to 4n. We bound the number of iterations to limit the computational time, thus having a total computational complexity in O(n) (at the expense of precision on the convergence).

In practice, the computational time of the transduction prevails over the one of SVM. The whole segmentation process takes between fifteen seconds and five minutes on a Xeon 2.33 GHz, for meshes with 50,000 to 1,500,000 facets.

Experimental Results

Figure 1 in the introduction shows a first example of segmentation into four different classes using the mean color of a facet and the altitude of its center as features. The results are qualitatively good, and mostly agree with perceptual boundaries. Note that every window is detected while only few where initially labeled, illustrating the ability of our algorithm to detect several objects of the same class forming several connected components, thanks to the initialization of the transduction with the classification results of the SVM.

Figure 2 illustrates the importance of combining several different features in order to obtain a relevant segmentation. Indeed the use of mean color or altitude alone produces erroneous results (even if the latter seems less wrong, it is even more naive, and performs a simple thresholding on the altitude), while their combination gives a pertinent result. Hence, the set of features is chosen Fig. 2. Combining different features improves the segmentation. 3 classes: roof, wall, ground. From left to right: input textured mesh with seeds; segmentation using mean color; segmentation using altitude; segmentation using both mean color and altitude.

according to the scene of interest: in Fig. 3, the segmentation of the top row is performed using color only, while in the bottom row, the color, the curvedness and the vertical normal of the facet are combined.

Figure 3 compares the classification of the SVM and the final result after transduction: as mentioned in section 3.3, the results of the SVM alone are noisy and do not exploit the connectivity of the mesh. Our algorithm can be used to remove outliers in a scene. In the original mesh of the castle of Figure 4, obtained by multi-view reconstruction, a portion of the sky has been reconstructed running on from the roof. It can be easily remove in our interactive framework, using only two classes, one for the outliers (sky facets), one for the inliers (castle facets), and color and altitude as features. The first segmentation (middle column) is not entirely satisfying (some portions of sky remain), so we add a few labeled points and rerun the algorithm, then obtaining a good segmentation (right column).

Conclusion

We have presented an efficient procedure for the segmentation of textured meshes: we designed a sketch-based interactive framework which produces a meaningful segmentation according to the user's aims, thanks to a possible refinement of the training sketches. Our experiments demonstrate that we can take advantage of geometric and photometric features at the same time, combining a various number of appropriate features in our kernel. The graph Laplacian regularizer enforces the spatial contiguity of the labels on the mesh, producing robust decompositions of the mesh for subsequent applications while the SVM initialization allows to detect non-connected, similar objects. Future work will focus on the development of supplementary features and kernels, in order to apply this algorithm to a larger range of scenes as well as to other types of data like point clouds.

Fig. 3 .

 3 Fig. 3. Comparison of classification results using using SVM alone or SVM plus transductive segmentation. Left: textured mesh with selected facets. Middle: classification obtained with SVM; the labels are noisy and do not provide a decomposition of the mesh. Right: classification obtained with SVM initialization plus transduction.

Fig. 4 .

 4 Fig. 4. Removing outliers. Left: original textured mesh of the castle with portions of "sky" mistakenly reconstructed by stereo. Middle: removing sky facets using selected facets. Right: refinement of the previous segmentation by adding sky facets to the selection (yellow outline). Top row: selected facets (green: sky outliers, blue: castle inliers), bottom row: textured meshes.


Ehsan Aganj

Jean-Philippe Pons

Renaud Keriven

Globally optimal spatio-temporal reconstruction from cluttered videos

We propose a method for multi-view reconstruction from videos adapted to dynamic cluttered scenes under uncontrolled imaging conditions. Taking visibility into account and being based on a global optimization of a true spatio-temporal energy, it offers several desirable properties: no need for silhouettes, robustness to noise, independent from any initialization, no heuristic force, reduced flickering results, etc. Results on real-world data proves the potential of what is, to our knowledge, the only globally optimal spatio-temporal multi-view reconstruction method.

Introduction

In recent years, several methods for automatic generation of complete spatiotemporal models of dynamic scenes from multiple videos have been proposed [START_REF] Aganj | Spatio-temporal shape from silhouette using four-dimensional Delaunay meshing[END_REF][START_REF] Ahmed | Automatic generation of personalized human avatars from multi-view video[END_REF][START_REF] Ahmed | Dense correspondence finding for parametrization-free animation reconstruction from video[END_REF][START_REF] De Aguiar | Performance capture from sparse multi-view video[END_REF][START_REF] De Aguiar | Marker-less deformable mesh tracking for human shape and motion capture[END_REF][START_REF] Furukawa | Dense 3D motion capture from synchronized video streams[END_REF][START_REF] Goldlücke | Space-time isosurface evolution for temporally coherent 3D reconstruction[END_REF][START_REF] Neumann | Spatio-temporal stereo using multi-resolution subdivision surfaces[END_REF][START_REF] Pons | Multi-view stereo reconstruction and scene flow estimation with a global image-based matching score[END_REF][START_REF] Starck | Surface capture for performance based animation[END_REF][START_REF] Theobalt | Seeing people in different light-joint shape, motion, and reflectance capture[END_REF][START_REF] Varanasi | Temporal surface tracking using mesh evolution[END_REF][START_REF] Vedula | Image-based spatio-temporal modeling and view interpolation of dynamic events[END_REF][START_REF] Vedula | Shape and motion carving in 6D. CVPR[END_REF][START_REF] Vlasic | Articulated mesh animation from multi-view silhouettes[END_REF]. In particular, the most recent ones have proven effective for full-body marker-less motion capture. Many of these techniques rely on the visual hull concept [START_REF] Laurentini | The visual hull concept for silhouette-based image understanding[END_REF], among which [START_REF] Aganj | Spatio-temporal shape from silhouette using four-dimensional Delaunay meshing[END_REF][START_REF] Ahmed | Dense correspondence finding for parametrization-free animation reconstruction from video[END_REF][START_REF] Vlasic | Articulated mesh animation from multi-view silhouettes[END_REF]. Computationally efficient, they suffer from several limitations: they provide an approximate reconstruction; this one has to be a closed surface; and, above all, silhouettes have to be segmented in the videos, practically limiting the method to controlled condition capture with a known background. This latest limitation may be lifted when prior knowledge about the geometry is available: free-form deformation of a template body model [START_REF] Ahmed | Automatic generation of personalized human avatars from multi-view video[END_REF][START_REF] Theobalt | Seeing people in different light-joint shape, motion, and reflectance capture[END_REF][START_REF] Vlasic | Articulated mesh animation from multi-view silhouettes[END_REF], Laplacian deformation of a laser scan of the initial pose [START_REF] De Aguiar | Performance capture from sparse multi-view video[END_REF][START_REF] De Aguiar | Marker-less deformable mesh tracking for human shape and motion capture[END_REF], etc. Yet, these methods are unable to recover genuine geometric details such as facial expressions and clothing folds and wrinkles. An exception might be the method proposed by Furukawa et al. [START_REF] Furukawa | Dense 3D motion capture from synchronized video streams[END_REF]. Yielding visually impressive results, this method does not rely on global optimization and handles the occlusion problem via heuristics.

Our approach

In this paper, we propose a method for multi-view reconstruction from videos adapted to dynamic cluttered scenes under uncontrolled imaging conditions. Taking visibility into account and being based on a global optimization of a true spatio-temporal energy, it offers several desirable properties.

Starting from work by Labatut et al. [START_REF] Labatut | Efficient multi-view reconstruction of largescale scenes using interest points, delaunay triangulation and graph cuts[END_REF], our method might be seen as its spatio-temporal extension. It is based on modeling an evolving three-dimensional surface as a four-dimensional surface [START_REF] Aganj | Spatio-temporal shape from silhouette using four-dimensional Delaunay meshing[END_REF][START_REF] Cheung | Shape-from-silhouette across time part i: Theory and algorithms[END_REF][START_REF] Goldlücke | Space-time isosurface evolution for temporally coherent 3D reconstruction[END_REF]. More precisely, we first generate a quasi-dense 3D point cloud of the scene at each time step and merge the result into a 4D point cloud. This process is conducted in a lenient manner, thus possibly retaining many false matches. Then, we build an adaptive decomposition of the 3D+time space by computing the 4D Delaunay triangulation of this cloud. Finally, we label the Delaunay pentatopes as empty or occupied thus generating a 4D surface. Graph-cut based, this assignment is globally optimal and compatible with the visibility in input images. Optionally but not necessarily, the 3D surfaces corresponding to each time step might be obtained considering 3D slices.

Contributions

Our method has several significant advantages. First, it is not based on visual hulls:

-The videos do not have to be taken under controlled conditions. The background might be cluttered. -It can handle closed as well as open scenes. For example, it can simultaneously recover the walls and (potentially moving!) furnitures of an indoor scene and a complete reconstruction of subjects seen from all sides in the input images.

Second, it is based on a global optimum:

-It is robust and does not depend on some initialization.

-It exploits visibility information to guide the position of the surface. As a result, it avoids the minimum cut solution from being an empty surface. Hence it exonerates from the usual techniques proposed so far to solve this problem (ballooning term, silhouette information, photo-flux, etc.). Moreover, this visibility information is not enforced as a hard constraint but integrated in the very optimization framework, hence yielding robustness to outliers.

Finally, and mainly, compared to the independent frame-by-frame computations of [START_REF] Labatut | Efficient multi-view reconstruction of largescale scenes using interest points, delaunay triangulation and graph cuts[END_REF], it profits from the 4D representation:

-Regularization is handled both in space and time, yielding more robustness to noise both in geometry and in visibility reasoning.

-Points extracted at one given time step transfer information to the surrounding time steps. As a result, more details are obtained at each time step. -Flicking artifacts in synthesized views are reduced, as consecutive 3D slices have similar geometry and connectivity by construction. -The temporally continuous representation, which is defined at any time, optionally enables interpolation of objects shape between consecutive frames.

The output of our method might be use for different purposes: as a 4D compact representation; as a list of consecutive 3D meshes; as an initialization for variational spatio-temporal stereovision methods [START_REF] Goldlücke | Space-time isosurface evolution for temporally coherent 3D reconstruction[END_REF].

The remainder of this paper is organized as follows. Section 2 gives some background on the different techniques needed in our approach. In Section 3, we describe in detail the different steps of our algorithm. Section 4 discusses numerical experiments that demonstrate the potential of our approach for reconstructing cluttered scenes from real-world data.

Background

Delaunay triangulation

Most of the following definitions are taken from [START_REF] Boissonnat | Provably good sampling and meshing of surfaces[END_REF]. We also refer the interested reader to some computational geometry textbooks [START_REF] Boissonnat | Algorithmic Geometry[END_REF][START_REF] De Berg | Computational Geometry, Algorithms and Applications[END_REF]. In the sequel, we call k-simplex the convex hull of k + 1 affinely independent points. For example, a 0-simplex is a point, a 1-simplex is a line segment, a 2-simplex is a triangle and a 3-simplex is a tetrahedron. In this paper, we will also consider 4-simplices: they are known as pentachorons or pentatopes. Let E = {p 1 , . . . , p n } be set of points in R d . The Voronoi region, or Voronoi cell, denoted by V (p i ), associated to a point p i is the region of space that is closer from p i than from all other points in E:

V (p i ) = {p ∈ R d : ∀j, p -p i ≤ p -p j } .
The Voronoi diagram of E, denoted by Vor(E), is the partition of space induced by the Voronoi cells V (p i ). See Figure 1(a) for a two-dimensional example of a Voronoi diagram. The Delaunay triangulation Del(E) of E is defined as the geometric dual of the Voronoi diagram: there is an edge between two points p i and p j in the Delaunay triangulation if and only if their Voronoi cells V (p i ) and V (p j ) have a non-empty intersection. It yields a triangulation of E, that is to say a partition of the convex hull of E into d-dimensional simplices (i.e. into triangles in 2D, into tetrahedra in 3D, into pentatopes in 4D and so on). See Figure 1(b) for a twodimensional example of a Delaunay triangulation. The fundamental property of the Delaunay triangulation is called the empty circle (resp. empty sphere in 3D, resp. empty hypersphere in higher dimensions) property: in 2D (resp. in 3D, resp. in 4D), a triangle (resp. tetrahedron, resp. pentatope) belongs to the Delaunay triangulation if and only if its circumcircle (resp. circumsphere, resp. circumscribed hypersphere) does not contain any other points of E in its interior.

Energy minimization by graph cuts

Given a finite directed graph G = (V, E) with nodes V and edges E with nonnegative weights (capacities), and two special vertices, the source s and the sink t, an s-t-cut C = (S, T ) is a partition of V into two disjoint sets S and T such that s ∈ S and t ∈ T . The cost of the cut is the sum of the capacity of all the edges going from S to T : c(S, T ) = (p,q)∈S×T |p→q∈E w pq . The minimum s-t-cut problem consists in finding a cut C with the smallest cost: the Ford-Fulkerson theorem [START_REF] Ford | Flows in Networks[END_REF] states that this problem is equivalent to computing the maximum flow from the source s to the sink t and many classical algorithms exist to efficiently solve this problem. Such a cut can be viewed as a binary labeling of the nodes: by building an appropriate graph, many segmentation problems in computer vision can be solved very efficiently [START_REF] Greig | Exact maximum a posteriori estimation for binary images[END_REF]. More generally, global minimization of a whole class of energy is achievable by graph cuts [START_REF] Kolmogorov | What energy functions can be minimizedvia graph cuts?[END_REF].

Kirsanov and Gortler [START_REF] Kirsanov | A discrete global minimization algorithm for continuous variational problems[END_REF] first proposed to use graph cuts on complexes to globally optimize surface functionals and also developed the idea of using random sparse complexes for their flexibility over regular subdivisions: this differs from the graphs commonly used in computer vision, which are often regular grids in the input images or in the bounding volume of the scene. Our approach similarly relies on a sparse complex-based graph: this graph however directly derives from an adaptive space decomposition efficiently provided by the Delaunay triangulation. Moreover the specifies of our graph construction are quite different and tailored to the multi-view reconstruction problem.

Method

Our spatio-temporal reconstruction algorithm consists in four steps: (i) a quasidense 3D point cloud is generated for each frame, each point memorizing the two or more images from which it has been triangulated. An spatio-temporal 4D point cloud is obtained by adding time as the fourth dimension to all the 3D points; (ii) the Delaunay triangulation of the 4D point cloud is computed; (iii) the Delaunay pentatopes (full-dimensional simplices in R 4 ) are labeled inside or outside the spatio-temporal object minimizing some energy, a 4D oriented surface is then extracted as the set of 4D tetrahedra lying between inside and outside pentatopes; (iv) the 3D surface at a given time is obtained by intersecting this 4D hyper-surface with a temporal plane.

4D point cloud generation, 4D Delaunay triangulation

Given multiple video sequences of a dynamic scene, we first make a dense 3D point cloud for each time instant. Let normal to the optical ray of x t k,i , we evalute the normalized cross correlation (NCC) between a given neighborhood of x t k,i and its image by H around x t k ′ ,j . Then for each x t k,i , we keep the m best point x t k ′ ,j according to NCC, only if their NCC is higher than a given threshold, and add the corresponding X t kk ′ ,ij to the point cloud at time instant t (in practice m = 1 or 2). The final 3D point cloud is then obtained by merging close 3D points, so that a point of the cloud comes from possibly more than two images. Now we construct a "global" spatio-temporal point cloud by regarding time as a fourth dimension, and treating it similarly to the three spatial dimensions. At first sight, this is questionable since space is not homogenous to time regarding physical units. We obtain physical homogeneity of our 4D space by considering a scaling factor v between space and time dimensions. This scaling factor is homogeneous to a speed, and can be interpreted as a reference displacement per time unit. The global point cloud is obtained by taking a 4D point (X t i , vt) ∈ R 4 for the point X t i generated from the input images in time t. At the end, we compute the 4D Delaunay triangulation of the spatio-temporal cloud storing in each vertex the list of the views and keypoints from which it has been generated.

I t k , k ∈ {1, • • • , n}, t ∈ [0, T ]

4D hyper-surface extraction

In this step we compute a four-dimensional representation of the dynamic scene by extracting a 4D mesh from the Delaunay triangulation of the point cloud. This is done by labeling Delaunay pentatopes as inside or outside of the spatiotemporal scene. The final oriented hypersurface is then extracted as the set of facets between inside and outside pentatopes. The exact nature of these "facets" deserves clarification: they are tetrahedra with 4D coordinates, so they are indeed simplicial pieces of a hypersurface in R 4 . Now we make a graph of neighbor pentatopes which we will use to find the optimal label assignment. For that, we take Delaunay pentatopes as vertices and we add edges between every two pentatopes which are neighbor via a two or three dimensional face. In addition, we add a link between each vertex of the graph and the sink and the source. A globally optimal label assignment is then efficiently found by applying the graph cuts optimization method on this graph.

In the sequel, we note S the surface to be reconstructed. As discussed above, S is a union of 4D Delaunay facets. In order to find an optimal solution satisfying both spatial and temporal constraints, we minimize an energy composed of two terms, one dealing with visibility, and the other dealing with spatial and temporal smoothness,

E(S) = E vis (S) + E smooth (S) (1)
In the rest of this section, we give the exact definition of these two terms and we describe how they can be implemented in the graph cuts framework.

Visibility term

The visibility term that we propose for a spatio-temporal scene is a careful extension of the static case proposed in [START_REF] Labatut | Efficient multi-view reconstruction of largescale scenes using interest points, delaunay triangulation and graph cuts[END_REF]. The idea of their work is that if a point belongs to the final surface then it should be visible in the views from which it has been triangulated. This yields to the penalization of all the facets intersecting the ray between the point and the cameras from which it has been generated. In the dynamic case the same argument holds. A point which belongs to the final hypersurface should be visible in all its generating views. Consequently, all 4D pentatopes which intersect a 4D ray emanating from the point to the camera center of one of its generating views should be labeled as outside, and the pentatope behind the point should be labeled as inside. We remark that the spatio-temporal center of a camera at a given frame is its 3D center positioned in the temporal plane of that frame. Similarly to the static case, in order to make an energy which can be minimized via graph cuts, we take the number of intersections of a ray with the oriented hypersurface as the visibility term. At this point, there are several important remarks to be made. First, the visibility of a point at a given frame is defined only in the temporal plane corresponding to that frame. Therefore, the rays between the point and its generating views lie completely in the temporal plane which passes through that point.

Second, a 4D facet of the Delaunay triangulation passes generally through several consecutive frames. As a consequence, each intersection of a ray with a facet is considered as a penalizing "vote" for the facet. The final vote is then computed as the sum of all votes coming from different frames intersecting the facet. This is an important property since it makes a global visibility vote on every 4D facet taking in account the temporal coherence. Third, contrarily to the static case presented in [START_REF] Labatut | Efficient multi-view reconstruction of largescale scenes using interest points, delaunay triangulation and graph cuts[END_REF] where edges of the graph are only "full dimensional facets" (3D triangles) between Delaunay simplices, in the dynamic case, in addition to these facets we add an edge between every two pentatopes which are not neighbors via a full-dimensional facet (4D) but via a 3D facet. Figure 2 shows an example of this situation. For simplicity reasons we consider a lower dimensional scene: points are on 2D planes, time is the third dimension and the spatio-temporal object is extracted from the 3D Delaunay triangulation of the point cloud. Points a and b are on frame 1, and points c,d,e and f are on frame 2. Point o is the center of a camera from which f has been generated. The tetrahedra acde and bcdf have a one-dimensional intersection (the segment cd), but they should be connected in the graph since the ray f o intersects cd and therefore a penalization term should be added between them. It is important to note that despite the 3D intersection of the ray with the face acd, no penalization term should be added between the tetrahedra abcd and bcdf . That is because abcd does not appear in the static representation of the scene on frame 2.

The intersections of the ray with the triangulation can be computed in the four-dimensional space handling carefully the situation discussed above. However, as a ray always lies completely in a temporal plane, we propose to find these intersections more easily by intersecting the 3D ray with the 3D intersection of the triangulation with the temporal plane. Obviously only the pentatopes which make a full-dimensional (4D) temporal intersection should appear in the temporal slice. In this case, a 3D facet intersected by a ray will correspond to an edge of the graph, and the unnecessary intersections discussed in the example above will be omitted by definition.

camera center p 2 q 1 q 2 vertex surface ray p 0 λ ∞ p 0 p 1 s t p 1 λ out s t q 1 q 2
λ in s t Fig. 3. Top: A 3D slice of the 4D triangulation. A ray emanating from a vertex to a camera center intersects 3D cells. Bottom: The corresponding visibility-related energy term that penalizes the number of intersections with the ray and the edge weights of the crossed pentatope in the graph.

We should remark that the 3D intersection of a 4D Delaunay triangulation with a plane is a polyhedron which contains generally cells with more than four vertices. Figure 3 shows the 3D object and a ray intersecting the cells. The vertices p 0 , p 1 , q 1 , p 2 and q 2 shown in Figure 3(bottom) are the vertices of the graph which correspond to the pentatopes whose temporal intersections make the cells p 0 , p 1 , q 1 , p 2 and q 2 shown in Figure 3(top) respectively. Different visibility terms are added. The cell containing the camera should be labeled as outside: a term λ ∞ is added to the edge from source to p 0 . A 3D facet crossed by the ray from inside to outside should be penalized: a term λ out is added to the edge from p 1 to q 1 . The cell behind the origin of the ray should be labeled as inside: a term λ in is added to the edge from q 2 to the sink. The positive weights λ in , λ out and λ ∞ take into account the confidence in the reconstructed vertex yielding the ray. By summing up these visibility terms over all the frames, we make a complex distribution of "vote" for each pentatope taking in account the time coherence.

Spatio-temporal smoothness In order to take in account both spatial smoothness and temporal continuity, we propose to minimize the area of the 4D surface in R 4 . This yields to the sum of volumes over all the 4D tetrahedra between inside and outside pentatopes, E smooth (S) = A(S) = S dS = T ∈S A(T ) where S is the 4D surface to be reconstructed, T is a 4D tetrahedra, and A(T ) is the volume of the tetrahedra T in R 4 . Minimizing this term encourages smoothness in time and in space. As in the static case, this is trivially minimized in the graph cuts framework: for each pair of pentatopes (sharing a tetrahedra T ) represented by vertices p and q in our graph, a term w = A(T ) is added the edge p → q and to its opposite edge q → p.

3D surface extraction

The output 4D mesh cannot be used directly for rendering. Fortunately, the 3D scene at a given time instant is easily obtained by intersecting this 4D mesh with a temporal plane. This task can be performed very efficiently, even in real-time on GPUs (Graphics Processor Units), since it reduces to a marching tetrahedra algorithm [START_REF] Guéziec | Exploiting triangulated surface extraction using tetrahedral decomposition[END_REF] on the tetrahedra of the 4D mesh, with the temporal coordinate of vertices used as the scalar field for isocontouring. It produces one triangle or one quad per boundary tetrahedron intersected by the selected temporal plane.

Experimental results

We have implemented our method using CGAL (Computational Geometry Algorithms Library, homepage: www.cgal.org) [START_REF] Boissonnat | Triangulations in CGAL[END_REF]. CGAL defines most data structure and algorithms needed in our method. For the computation of 4D Delaunay triangulation we have used the Quickhull algorithm library (QHull) [START_REF] Barber | The quickhull algorithm for convex hulls[END_REF] (homepage: www.qhull.org). Fig. 5. A comparison between our method and the method of [START_REF] Labatut | Efficient multi-view reconstruction of largescale scenes using interest points, delaunay triangulation and graph cuts[END_REF] applied independly in each frame. Top: 3D meshes obtained by the method of [START_REF] Labatut | Efficient multi-view reconstruction of largescale scenes using interest points, delaunay triangulation and graph cuts[END_REF]. Bottom: corresponding 3D slices of the 4D representation obtained by our method. [START_REF] Labatut | Efficient multi-view reconstruction of largescale scenes using interest points, delaunay triangulation and graph cuts[END_REF] (Bottom): our method. The frame-by-frame reconstruction of the method [START_REF] Labatut | Efficient multi-view reconstruction of largescale scenes using interest points, delaunay triangulation and graph cuts[END_REF] makes no temporal continuity. In contrast, our method reconstructs correctly the trouser, avoids flicking artifacts and provides a much more continuous motion.

In our first experiment, we have tested our method on the first 60 frames of the "Trousers" sequence which is courtesy of R. White, K. Crane and D.A. Forsyth [START_REF] White | Capturing and animating occluded cloth[END_REF]. The sequence is acquired by 8 cameras at a 640 × 480 resolution. Figure 4 shows some images of this dataset. Despite their highly textured images, it is in fact a quite challenging dataset because of the very fast and complex motion of cloth and folds and severe self occlusion in various parts of the video. Regarding to an approximate size of the object we have chosen a spatio-temporal scaling factor v = 30 (space unit/frame). Figure 5 shows a comparison between our method and the method of [START_REF] Labatut | Efficient multi-view reconstruction of largescale scenes using interest points, delaunay triangulation and graph cuts[END_REF] applied independly in each frame. We have compared the 3D meshes output by their method for some frames of the sequence with the corresponding 3D slices of our 4D spatio-temporal scene. A total number of 793769 points have been generated for the initial point cloud. To provide a fair comparison, we have used the same point clouds in both methods.

We observe that the method of [START_REF] Labatut | Efficient multi-view reconstruction of largescale scenes using interest points, delaunay triangulation and graph cuts[END_REF] fails to separate correctly the two trouser legs when there is not enough distance between them. In addition, as shown in figure 6, their frame-by-frame reconstruction makes no temporal continuity. In contrast, relying on a global optimization, our method reconstructs correctly the trouser, avoids flicking artifacts and provides a much more continuous motion. This perfectly illustrates the capability of our approach to take advantage of temporal coherence in order to obtain more detailed and more continuous result. Please see supplemental material for a video illustrating better the result of our method. The computational times of our method and the method of [START_REF] Labatut | Efficient multi-view reconstruction of largescale scenes using interest points, delaunay triangulation and graph cuts[END_REF] for this experiment are 210 and 112 minutes respectively on a standard workstation. However, the most expensive part of our method is the computation of the 4D Delaunay triangulation. Fortunately, this can be strongly reduced using an optimized 4D Delaunay code.

In a second experiment, we have tested our method on the first 40 frames of the "Dancer" dataset which was made available to us by the 4Dviews company (http://4dviews.com). It is acquired by 14 calibrated and synchronized video cameras. Figure 7(top) shows some images of this dataset. The result shows that despite the lowly textured parts of the images, our method makes a correct 4D representation of the dancer. Figure 7(bottom) shows some 3D slices extracted from the 4D object.

Finally, we should remark that in order to have better visualization and to make better comparisons we have smoothed the results of our experiments. However, the output of our method might be used as a 4D compact representation, as a list of consecutive 3D meshes or as an initialization for variational spatio-temporal stereovision methods.

Discussion and Conclusion

We have presented a new method for multi-view reconstruction from videos adapted to dynamic cluttered scenes under uncontrolled imaging conditions. The main idea of our method is to regard time as the fourth dimension, and to extract a hyper-surface from the 4D Delaunay triangulation of the input points as the Fig. 7. Top: Some images of the "Dancer" dataset. Bottom: Some 3D slices extracted from the 4D representation of the "Danser" dataset, obtained by our method.

spatio-temporal representation of the scene. This is done by labeling Delaunay pentatopes as empty or occupied. A globally optimal assignment is efficiently found using graph cuts. We have validated our method on real video sequences. Our results prove the potential of what is, to our knowledge, the only globally optimal spatio-temporal multiview reconstruction method.

Fig. 1 .

 1 Fig. 1. (a) Voronoi diagram of a set of points in the plane. (b) Its dual Delaunay triangulation.

Fig. 2 .

 2 Fig.2. The tetrahedra acde and bcdf have a one-dimensional intersection (the segment cd), but they should be connected in the graph (refer to text for details).

Fig. 4 .

 4 Fig. 4. Some images of the "Trousers" dataset.

Fig. 6 .

 6 Fig.6. The first three consecutive frames of the Trousers dataset reconstructed by (Top): the method of[START_REF] Labatut | Efficient multi-view reconstruction of largescale scenes using interest points, delaunay triangulation and graph cuts[END_REF] (Bottom): our method. The frame-by-frame reconstruction of the method[START_REF] Labatut | Efficient multi-view reconstruction of largescale scenes using interest points, delaunay triangulation and graph cuts[END_REF] makes no temporal continuity. In contrast, our method reconstructs correctly the trouser, avoids flicking artifacts and provides a much more continuous motion.

  ) of the same type verifying the epipolar constraint up to a certain error (due to point extraction but also to calibration), we triangulate the corresponding 3D point X t kk ′ ,ij . Let H be the homography from I t k to I t k ′ induced by the plane at X t kk ′ ,ij

	denote the
	input video sequence. For each image we extract interest points x t k,i of several
	types (in practice Harris points and DOGs) without any scale information and
	with thresholds such that their number is high enough. For each image pair in
	a time instant t, (I t k , I t k ′ ) whose visual fields intersect, for all (x t k,i , x t k ′ ,j


Non-rigid Shape matching using Geometry and Photometry

Nicolas Thorstensen Abstract. In this paper, we tackle the problem of finding correspondences between three-dimensional reconstructions of a deformable surface at different time steps. We suppose that (i) the mechanical underlying model imposes time-constant geodesic distances between points on the surface; and that (ii) images of the real surface are available. This is for instance the case in spatio-temporal shape from videos (e.g. multiview stereo, visual hulls, etc.) when the surface is supposed approximatively unstretchable. These assumptions allow to exploit both geometry and photometry. In particular we propose an energy based formulation of the problem, extending the work of Bronstein et al. [START_REF] Bronstein | Generalized multidimensional scaling: a framework for isometry-invariant partial surface matching[END_REF]. On the one hand, we show that photometry (i) improves accuracy in case of locally elastic deformations or noisy surfaces and (ii) allows to still find the right solution when [START_REF] Bronstein | Generalized multidimensional scaling: a framework for isometry-invariant partial surface matching[END_REF] fails because of ambiguities (e.g. symmetries). On the other hand, using geometry makes it possible to match shapes that have undergone large motion, which is not possible with usual photometric methods. Numerical experiments prove the efficiency of our method on synthetic and real data.

Introduction

Most of the objects observed in the real world are non-rigid. This makes them particularly important for computer vision. What makes non-rigid shapes challenging is that their associated deformations exhibit a potentially infinite number of degrees of freedom. As a consequence they are hard to analyze. One typical example is the three-dimensional (3D) reconstruction of a person in a multiple cameras environment. More generally, matching 3D-reconstructed shapes have numerous applications, among which are reconstruction-based animation, motion analysis, shape recognition, physical phenomena analysis, etc.

A subject observed at different times usually results in meshes in different poses. In such a situation, image-based matching algorithms [START_REF] Vedula | Three-dimensional scene flow[END_REF][START_REF] Pons | Multi-view stereo reconstruction and scene flow estimation with a global image-based matching score[END_REF] tend to fail, as the motion in between the two poses is too large. Assuming the object surface is approximately unstretchable, Bronstein et al. [START_REF] Bronstein | Generalized multidimensional scaling: a framework for isometry-invariant partial surface matching[END_REF] designed a matching framework where geodesic distances between automatically chosen key points are preserved. Yet, such a criterion fails in case of ambiguities like symmetries. It also yields inaccuracy when the unstretchability hypothesis is violated.

Our contributions

Therefore we propose to take the best of both worlds and design an extension of the work of [START_REF] Bronstein | Generalized multidimensional scaling: a framework for isometry-invariant partial surface matching[END_REF] by adding a photometric term to their energy. We obtain a robust multi-resolution 3D surface mapping procedure that combines photometric and geometric information(c.f. Figure 1). We experiment it for non-rigid surface correspondence between two surfaces observed a different time steps in a multiple cameras environment and demonstrate its superiority.

The paper is organized as follows. Section 2 reviews previous work for nonrigid shape matching. Section 3 states the problem formulation and Section 4 explores the algorithmic implementation. Finally, numerical experiments on real data are reported in Section 5 and Section 6 concludes.

Related Work

The correspondence problem is one of the fundamental challenges in computer vision. Might it be in the context of optical flow, calibration or surface registration. For rigid surface and point cloud registration, Iterative Closest Point (ICP) and its variants [START_REF] Gelfand | Robust global registration[END_REF] are the standard algorithms. Operating on a purely geometric level, they rely on approximated differential quantities, e.g. curvature, or more robust surface descriptors [START_REF] Belongie | Shape matching and object recognition using shape contexts[END_REF]. Recently, several algorithms also address the problem of non rigid surface registration. They can be mainly divided into two categories: geometric and photometric. Whereas geometric methods assume the geometry known, photometric methods estimate structure and motion.

Geometric methods

A non-rigid counterpart to ICP was introduced in [START_REF] Amberg | Optimal step nonrigid icp algorithms for surface registration[END_REF]. Further in [START_REF] Brown | Global non-rigid alignment of 3-D scans[END_REF], the authors propose a non rigid registration method by piecewise rigid registration and local thin plate spline alignment of range scan images. Elad et al. [START_REF] Elad | On bending invariant signatures for surfaces[END_REF] take an intrinsic point of view of the matching problem. By isometrically embedding two surfaces into a common Euclidean space the problem is reduced to the simpler of rigid matching between surfaces. Following this idea, Bronstein et al. [START_REF] Bronstein | Efficient computation of isometryinvariant distances between surfaces[END_REF] take it one step further and solve the correspondence problem by embedding the surface directly into another. The major drawback of this approach is the difficulty to handle ambiguities like shape symmetries. As the method solely relies on the distance function of the surface, the algorithm tends to wrongly match pairs of symmetric points (see Figure 2).

Photometric methods

On the photometric side, non rigid registration is commonly known under the name of scene flow. Vedula et al. [START_REF] Vedula | Three-dimensional scene flow[END_REF] were the first to introduce the concept of scene flow which is the 3D extension of 2D optical flow. It is the simultaneous reconstruction of structure and motion between time frames. Further work on scene flow was undertaken by other teams [START_REF] Furukawa | Dense 3d motion capture from synchronized video streams[END_REF][START_REF] Huguet | A variational method for scene flow estimation from stereo sequences[END_REF][START_REF] Pons | Multi-view stereo reconstruction and scene flow estimation with a global image-based matching score[END_REF]. As pointed out by Starck et al. [START_REF] Starck | Correspondence labelling for wide-timeframe free-form surface matching[END_REF], such approaches are limited to small displacements.

Where geometry meets photometry

Lastly, one can take advantages of both approaches. Several recent papers address this problem. In [START_REF] Starck | Spherical matching for temporal correspondence of non-rigid surfaces[END_REF], Starck et al. provide a method to register shapes from silhouette. The method embeds the two surfaces into a common spherical domain. Using the analytic expression of geodesic distances on the sphere allows to minimize the geometric distance between correspondences and distance between the associated color feature. This work is the most related to ours in the sense that they perform an isometric mapping into a common metric space and use geodesic distances as regularization for the appearance minimization.

In order to be robust with respect to topological changes, hey then developed a matching algorithm based on a Markov random field optimization [START_REF] Starck | Correspondence labelling for wide-timeframe free-form surface matching[END_REF]. A much more heuristic approach is used by Ahmed et al. [START_REF] Ahmed | Dense correspondence finding for parametrization-free animation reconstruction from video[END_REF]. A set of SIFT features is computed which are then used for initialization of a refinement model. Lastly, the work of Varanasi et al. [START_REF] Varanasi | Temporal surface tracking using mesh evolution[END_REF] makes also use of SIFT. The 3D position of the interest points are used to compute a sparse estimate of the motion field. By applying a diffusion operator on the sparse motion field, they recover a dense motion field.

Most of the algorithms presented so far rely on heuristics to put two temporal meshes into correspondence. In this work we want to consider a more rigorous model and propose a variational refinement method using ideas from variational stereo [START_REF] Keriven | Complete dense stereovision using level set methods[END_REF] and geometric embedding methods [START_REF] Elad | On bending invariant signatures for surfaces[END_REF][START_REF] Bronstein | Efficient computation of isometryinvariant distances between surfaces[END_REF]. In this way, we recover temporal correspondences between meshes at different resolution and independently of the combinatorial configuration. Further, our method can handle local elastic deformations such that the correspondence is consistent with the observed images.

Variational matching

Let S 1 and S 2 be the two surfaces to match. Each of them is observed by a certain number of cameras. Although not required, we simplify notations and suppose that the number and positions of the cameras are constant, so that both surfaces are observed by n cameras defined by their respective projections Π 1 , . . . Π n . We denote by I k i the i th image of surface k. Following [START_REF] Bronstein | Efficient computation of isometryinvariant distances between surfaces[END_REF], we use a Lagrangian point of view where a set of m correspondences are constrained to move on the surface such that they minimize a given energy. We denote (P 1 i , P 2 i ) such a correspondence where P k i ∈ S k and Θ = {(P 1 i , P 2 i ), 1 ≤ i ≤ m} the set of all correspondences. m is a fixed integer that can be estimated during initialization (see Section 4.2). Our energy is a function of parameter Θ that writes:

E tot (Θ) = αE geom (Θ) + βE photo (Θ) (1)
The first term E geom is the geometric part, taken from the work of Bronstein et al. [START_REF] Bronstein | Generalized multidimensional scaling: a framework for isometry-invariant partial surface matching[END_REF] while E photo is our photometric attachment. As usual, α and β are positive constants that control the relative weights among these terms. Minimizing energy (1) with respect to Θ will position the correspondences on the mesh such that their projections in the images minimize a photometric dissimilarity measure while geodesic distances on the surfaces are respected.

Geometry

Bronstein et al. [START_REF] Bronstein | Efficient computation of isometryinvariant distances between surfaces[END_REF] propose to embed near-isometric surfaces one into another by minimizing the following energy:

E geom (Θ) = i>j (d S 1 (P 1 i , P 1 j ) -d S 2 (P 2 i , P 2 j )) 2 . ( 2
)
where d S k is the geodesic distance on surface S k . Again, this energy suffers mainly from two weaknesses: (i) symmetries yield ambiguities and (ii) if the object undergoes locally elastic deformations between shape S 1 and shape S 2 , geodesic distances are not exactly preserved. Note that this is also the case when the surface are 3D reconstructions, since they are unavoidably noisy.

Photometry

For our image matching term E photo , we chose the normalized cross-correlation to measure similarity between corresponding points. Its simplicity, robustness in the presence of varying lighting conditions and differentiability make it a common choice in variational methods. Each surface point is generally seen from several cameras and one might be tempted to correlate multiple pairs of images. However, in our experiments, the number of cameras is relatively small. Thus, using information from only one pair of camera for each surface point reveals to be enough. As a first step, for each surface S k , we associate to each point M ∈ S k an optimal image I k l k (M ) . Choice of labels l k might be guided by different criteria. Here, we compute partitions of the surfaces following [START_REF] Allène | Seamless image-based texture atlases using multi-band blending[END_REF]. This method assigns smoothly each point to a label corresponding to the camera from which it is best viewed. Using graph-cut optimization, the labeling is obtained by minimizing a weighted sum of two terms which realizes a good trade off between resolution and color continuity, while respecting occlusions (c.f. Figure 2b).

Let l k (i) be a short notation for l k (P k i ), our photometric energy then writes:

E photo (Θ) = m i=1 g[N CC(I 1 l 1 (i) • Π l 1 (i) , I 2 l 2 (i) • Π l 2 (i) )(P 1 i , P 2 i )] (3)
where g is a positive decreasing function and N CC(f 1 , f 2 )(M 1 , M 2 ) denote the normalized cross-correlation of functions f 1 and f 2 between two related neighborhoods of points M 1 and M 2 . Following the stereovision work of Keriven et al. [START_REF] Keriven | Complete dense stereovision using level set methods[END_REF], we approximate locally the surfaces by their tangent planes at points M k . In their case, only one surface is considered and M 1 and M 2 are the same point, with the same tangent plane, thus the same neighborhood. Their correlation boils down to correlating image regions related by a homography. In our case, we suppose that the tangent plane to S 1 at point M 1 and the tangent plane to S 2 at M 2 are related by a given two-dimensional isometry I M 1 ,M 2 sending M 1 to M 2 . Under this assumption, neighborhoods on the respective tangent planes are related and the correlation N CC(f 1 , f 2 )(M 1 , M 2 ) is correctly defined. Moreover, it (and its derivatives) remains easy to compute since corresponding image regions are still related by a homography. Introducing the isometry I M 1 ,M 2 are each point pair (M 1 , M 2 ) a) b)

Fig. 3. a) On the left, the front image of the first mesh, red dots being the projections of some key points. On the right, the front image of the second mesh, showing the reprojections of the corresponding points, in red obtained with the method in [START_REF] Bronstein | Efficient computation of isometryinvariant distances between surfaces[END_REF] (initialized with our method), in green obtained with our method. b) Shows the same correspondences directly on the meshes with the same color code. The inaccuracies of [START_REF] Bronstein | Efficient computation of isometryinvariant distances between surfaces[END_REF], here due to local elastic deformations, are corrected by our photometric + geometric criterion.

might be thought as problematic since one would require to match the surfaces to know it, yielding a chicken and egg problem. Practically, this is not the case.

As usual, we will minimize the energy by mean of a gradient descent starting from a coarse initialization (see Section 4). This approximate solution reveals to be sufficient to obtain a robust I M 1 ,M 2 . We proceed the following way: (i) each correspondence point P k i defines a geodesic distance map d S k (P k i , .) on S k ; (ii) the gradients of these distance maps at a given point M k defines local directions that should correspond from one surface to the other if M 1 corresponds to M 2 ; (iii) as a consequence, the best 1 isometry from the tangent plane at M 1 to the one at M 2 that sends both M 1 to M 2 , and the distance gradients directions at M 1 to the ones at M 2 , is a good estimate of I M 1 ,M 2 . Please note that computing distance maps is no extra cost since it is part of E geom . Note also that the isometries I are actually needed only for the pairs (P 1 i , P 2 i ) and that they will be refined during the gradient descent iterations as the pairs moves.

Optimization

Discretization

We suppose that the surfaces are both discretized as collections of triangles. Following [START_REF] Bronstein | Efficient computation of isometryinvariant distances between surfaces[END_REF], points P k i are taken as barycenters of triangle vertices. Θ consists in a choice of triangles and corresponding barycentric coordinates. The geodesic distances between all vertices of the mesh are computed using the Fast Marching algorithm for triangular surfaces [START_REF] Sethian | Fast marching methods[END_REF][START_REF] Kimmel | Computing geodesic paths on manifolds[END_REF]. Geodesic distance is then interpolated like in [START_REF] Bronstein | Efficient computation of isometryinvariant distances between surfaces[END_REF] (note that the problem is not only to interpolate the distance at some barycenter, but also the distance to some barycenter).

For the photometric part of the energy, discretization is not a particular issue: the labeling method [START_REF] Allène | Seamless image-based texture atlases using multi-band blending[END_REF] is designed for triangle meshes, and we use a standard normal interpolation method to estimate the tangent planes.

Minimizing the energy with respect to both the P 1 i 's and the P 2 i 's is obviously not well posed. Although different cases might happen, in our experiments we have no further constraint on the choice of the points to be matched. Thus, we fix points P 1 i to their initial position (see Section 4.2) and minimize the energy with respect to the positions of points P 2 i . As written above, we use a classical gradient descent. Properly minimizing it is not trivial because the problem is non-convex. In order to cope with local minima, we apply a multi resolution strategy, considering the problem at several scales. Once a solution is found at a coarse scale, it is used to initialize the problem at a finer scale. Our problem has two scalable dimensions. The first one is the number of correspondences and the second is the scale of the images. This leads to a two step multi resolution scheme. Starting with a small number of correspondences, we iteratively increase the number of points by interpolating the solution from the coarser level to the next finer level. This scheme is adapted from [START_REF] Bronstein | Efficient computation of isometryinvariant distances between surfaces[END_REF]. Then, at each level, we perform a gradient descent in a multi scale manner using a Gaussian pyramid of the images.

Initialization

We first have to initialize the correspondences. Copying [START_REF] Bronstein | Efficient computation of isometryinvariant distances between surfaces[END_REF], we take advantage of the geodesic distance maps and use the farthest point sampling (FPS) strategy [START_REF] Moenning | Fast marching farthest point sampling[END_REF] to get geometry-based feature points on the surfaces. For near isometric surfaces we can expect the sampling to be almost identically distributed on both surfaces S 1 and S 2 [START_REF] Bronstein | Efficient computation of isometryinvariant distances between surfaces[END_REF]. Taking photometry into account to avoid geometric ambiguities, we then reject points that have an autocorrelation score below a given threshold, thus corresponding to non textured regions. As in [START_REF] Bronstein | Efficient computation of isometryinvariant distances between surfaces[END_REF], points are then associated using branch and bound optimization [START_REF] Dasgupta | [END_REF], yielding m initial pairs. Here, to the initial geodesic distance based criterion, we add a photometric one in order to get rid of geometric ambiguities. Because no correlation is possible (at this stage, tangent planes cannot yet be related by isometries), we use SIFT descriptor based similarity, being thus invariant to scale and orientation. The results of the initialization of the correspondences can be viewed in Figure 2. Note how geometric ambiguities are solved.

Gradient descent

Optimization is performed at all scales until convergence is reached, i.e. norm of the gradient is below a given threshold. The expressions of the gradients of the geometric and the photometric parts of our energy can be found in [START_REF] Bronstein | Efficient computation of isometryinvariant distances between surfaces[END_REF] and [START_REF] Keriven | Complete dense stereovision using level set methods[END_REF] respectively. Remember that Θ consists in these coordinate but also in the choice of the triangles to which the barycenters are related. As in [START_REF] Bronstein | Efficient computation of isometryinvariant distances between surfaces[END_REF], the gradients are computed for a fixed such choice. However, when a point P k i gets out of its related triangle, we force it to stop at the reached edge and assign it to the triangle at the "other side" of this edge. Doing this way, points travel gently from one triangle to another if needed.

Results

In order to validate the proposed method, we run several experiments on real and synthetic data. First, we test it on a synthetic dataset. In a second experiment, we validate our algorithm on real images.

Validation

Our first experiment focuses on the validation of our energy by testing the algorithm on a synthetic dataset (8 cameras) and comparing to the result of [START_REF] Bronstein | Efficient computation of isometryinvariant distances between surfaces[END_REF]. This experiment aims at justifying the photometric part of the energy. The parameter α is set to 1.0 whereas β equals 8.0. We take 12 pairs of correspondences and a 3 level Gaussian pyramid. The advantages of our initialization having already be demonstrated on Figure 2, we rather launch the original method proposed by Bronstein et al. [START_REF] Bronstein | Efficient computation of isometryinvariant distances between surfaces[END_REF] with our initialization. The red dots in the left image of Figure 3a are some of the P 1 i projected on the front image of S 1 . The red dots in the right image of Figure 3a are the projections of the corresponding P 2 i obtained after running the optimization of [START_REF] Bronstein | Efficient computation of isometryinvariant distances between surfaces[END_REF]. The green dots correspond to the result obtained with our combined photometric-geometric optimization. One can clearly see, the green dots are consistent with the initial sampling in the left image although the zone around the knee and shoulder exhibit elastic deformation. Whereas the red dots in the right image ignore the image signal and are pushed away by the local elastic deformations. Figure 3b shows the same points on the meshes.

Real data

In order to see how our method performs, we run several experiments on image data courtesy of J. Starck2 [START_REF] Starck | Surface capture for performance based animation[END_REF], again with 8 cameras. In this experiment the number of correspondences is 150 and the number of image levels is set to 3. α and β are set to 0.9 and 1.5 respectively. The results are depicted in Figure 4. Notice how the method of [START_REF] Bronstein | Efficient computation of isometryinvariant distances between surfaces[END_REF] fails to solve the matching problem. Local elastic deformations are observed in both cases (Figure 4a and4c) and wrong matches occur because of symmetry. Nevertheless, our method can handle the symmetries and local elastic deformations as can be noticed the around the hair and the back in Figure 4b and in the zones located on the skirt and the hair in Figure 4d.

Conclusion

We have described a variational framework for the correspondence problem of non-rigid surfaces coming from multi-view reconstruction or from any setup that provides images of the scene. By using photometric and geometric information, our method improves the one by Bronstein et al. [START_REF] Bronstein | Efficient computation of isometryinvariant distances between surfaces[END_REF] and allows to cope with elastic stretches and symmetries of the shape. Like in the initial work, and contrary to usual photometric methods, we are not limited to small deformations. Numerical experiments show the success of our method. Our future work includes extending our framework to surfaces with varying topology. A first step in this direction could be the use of more robust embedding such as the one proposed in [START_REF] Bronstein | A gromovhausdorff framework with diffusion geometry for topologically-robust non-rigid shape matching[END_REF]. Another improvement is the use of second order information in the gradient descent in order to perform Newton steps. Finally, we notice that regularization is important in low-textured image regions. Therefore we currently investigate the effect of automatically setting the control weights as proposed in [START_REF] Vu | Towards high-resolution largescale multi-view stereo[END_REF] in another context.

1 and Renaud Keriven 1 Fig. 1 .

 11 Fig. 1. a) Multiple views of a person in two different poses; b) Mesh based representation of the reconstructed object and heat map of the geodesic distance to the lower left foot; c) Correspondence result after optimization

Fig. 2 .

 2 Fig. 2. a) Initialization using solely geometric information is sensible to symmetry in shapes. b) Color coded labeling of the meshes (see text) and our photometric plus geometric initialization, where ambiguities have been solved.

Fig. 4 .

 4 Fig.4. Results on a real-worl dataset[START_REF] Starck | Surface capture for performance based animation[END_REF]. a) and c) show the results obtained using the method in[START_REF] Bronstein | Efficient computation of isometryinvariant distances between surfaces[END_REF]. b) and d) are obtained using our method.

http://personal.ee.surrey.ac.uk/Personal/J.Starck/

Acknowledgement We would like to thank A.M. and M.M. Bronstein for their kind help and for providing their code.


• All available data on storms and storm impact gathered, backward modeling of wind/waves • 121 events listed between 1226 and 1868.

• 53 events analyzed between 1940 and 2008

Main conclusions

Objectives

Obtain an homogenous data at the regional scale, and share the various initiatives to monitor storms and storm impacts ;

Capitalize and disseminate the gathered information. Sur la base des connaissances historiques, une tempête de ce type peut occasionner des impacts sur le littoral et un suivi de la tempête peut être mis en oeuvre selon le protocole défini dans le cadre du réseau tempête.

Impacts

- -----------------------------------------------------------------------------------------------------------------------Pour rappel, les suivis suivants sont préconisés (pour plus de détails se référer au protocole) :

-3m < Hs < 4m : suivi léger avec relevé des réseaux de mesures (station météo ; bouées, marégraphes ; stations vidéo) sans déplacement sur le terrain ; -4m < Hs < 5 m : suivi amélioré avec en sus des observations de terrain simples (formulaire, photos, points GPS a main) ; -Hs > 5m : suivi complet avec en sus des levés complémentaires sur le terrain (photographies aériennes et levés topo-bathymétriques) et des sorties de modèle.

----------------------------------------------------------------------------------------------------------------------------s il ne vous est pas possible de mobiliser un suivi, merci de prendre contact avec l'équipe d'animation du réseau tempête (y.balouin@brgm.fr et y.delatorre@brgm.fr ).

  Gather an exhaustive information on storm impacts to: Improve our understanding of storm processes and validate the existing local and regional cartographies of coastal hazards and coastal risks ; Improve the management strategy of coastal hazard. Double observation : data and observation exist but hardly available ; very heterogeneous.

  information from already existing approaches no specific budget; no moral and financial obligation.

•

  No transfer of intellectual properties 7th EUREGEO, June 12-15th, Bologna, Italy > 7 Monitoring of specific places defined by the partners within their own data collection program Automatic message sent 3 days before an event Storm thresholds based on historical analysis 3m Ce message s'applique au littoral du Languedoc-Roussillon. Il est a vocation informative uniquement et s'adresse aux membres du réseau tempête du Languedoc-Roussillon.

  selected by field surveyers to characterise impacts, water levels, erosion, marine floodingAll partners have received the message indicating a possible event 27 defined stations were visited: more than 100 observations Data and observations: state of the beach, water levels, wave conditions, impacts, shoreline retreat Formulaires + photos + GPS points Evaluation of the storm event Wave under the 1y RT, and moderate storm surge Impacts increased by the storm group Moderate impact: erosion of beach and dune fronts, damage on coastal Sharing human and technical effort to obtain an exhaustive and quantitative information on storms and storm impacts Operational procedure successful Datasets, photos, formularies shared in an online web server (restricted to the partnership) Perspectives Implementation of an on-line tool for data transfer WebGis and automatic edition of storm reports Implementation of the historical database Use of datasets for calibration of models http://littoral.languedocroussillon.fr Conclusion