Fabio Cuzzolin
email: fabio.cuzzolin@brookes.ac.uk

Diana Mateus
email: mateus@in.tum.de

Radu Horaud
email: radu.horaud@inria.fr

Robust Temporally Coherent Laplacian Protrusion Segmentation of 3D Articulated Bodies

Keywords: Unsupervised segmentation, 3D shape analysis, 3D shape tracking, Human motion analysis, Spectral embedding

In motion analysis and understanding it is important to be able to fit a suitable model or structure to the temporal series of observed data, in order to describe motion patterns in a compact way, and to discriminate between them. In an unsupervised context, i.e., no prior model of the moving object(s) is available, such a structure has to be learned from the data in a bottom-up fashion. In recent times, volumetric approaches in which the motion is captured from a number of cameras and a voxel-set representation of the body is built from the camera views, have gained ground due to attractive features such as inherent view-invariance and robustness to occlusions. Automatic, unsupervised segmentation of moving bodies along entire sequences, in a temporally-coherent and robust way, has the potential to provide a means of constructing a bottom-up model of the moving body, and track motion cues that may be later exploited for motion classification. Spectral methods such as locally linear embedding (LLE) can be useful in this context, as they preserve "protrusions", i.e., high-curvature regions of the 3D volume, of articulated shapes, while improving their separation in a lower dimensional space, making them in this way easier to cluster. In this paper we therefore propose a spectral approach to unsupervised and temporally-coherent body-protrusion segmentation along time sequences. Volumetric shapes are clustered in an embedding space, clusters are propagated in time to ensure coherence, and merged or split to accommodate changes in the body's topology. Experiments on both synthetic and real sequences of dense voxel-set data are shown. This supports the ability of the proposed method to cluster body-parts consistently over time in a totally unsupervised fashion, its robustness to sampling density and shape quality, and its potential for bottom-up model construction.

Introduction, Background and Contributions

In motion analysis and understanding it is crucial to be able to fit a suitable model to a temporal series of data observed in association with an object's evolution, in order to describe motions in a compact way, discriminate between them, or infer the pose, however described, of the moving object. In an unsupervised context in which no prior model of the moving object(s) is available, such a structure has to be learned from the data in a bottom-up fashion. In recent times, multi-view stereo methods, e.g., [START_REF] Furukawa | Carved visual hulls for high accuracy image-based modeling[END_REF][START_REF] Hernandez | Silhouette and stereo fusion for 3d object modeling[END_REF], using a set of synchronized and calibrated cameras and yielding volumetric or surface representations, have gained ground, due to attractive features such as inherent view-invariance and relative robustness to occlusion. Moreover, it is also possible to extract 3D motion flows from this kind of data, e.g., [START_REF] Pons | Multi-view stereo reconstruction and scene flow estimation with a global image-based matching score[END_REF].

Meshes, in particular, are widely used in computer vision, computer graphics and computeraided design to describe objects for a variety of tasks such as object recognition, deformation, animation, visualization, etc. Often, these tasks require the objects to be decomposed into sub-parts: mesh segmentation, therefore, plays a crucial role in many applications, and is a well-studied problem in computer graphics [START_REF] Huang | Shape decomposition using modal analysis[END_REF][START_REF] Katz | Mesh segmentation using feature point and core extraction[END_REF][START_REF] Lien | Approximate convex decomposition of polyhedra[END_REF][START_REF] Shamir | A survey on mesh segmentation techniques[END_REF][START_REF] Shapira | Consistent mesh partitioning and skeletonisation using the shape diameter function[END_REF], as it links skeletal representations useful for animation with mesh representations necessary for rendering. Most approaches, however, deal with smooth meshes that are acquired from artists or from laserscans of real actors. In opposition, multiple-stereo methods for 3D reconstruction, while they are true to reality at every instant and allow the texture information to be mapped in a straightforward way to the reconstructed geometry, may generate topological and geometric artifacts when assuming no prior knowledge of the scene being observed [START_REF] Zaharescu | Topology-adaptive mesh deformation for surface evolution, morphing, and multi-view reconstruction[END_REF]. For example, the two legs of a human actor may get clubbed together, or holes can be introduced in the resulting 3D model because of silhouette extraction problems [START_REF] Franco | Efficient Polyhedral Modeling from Silhouettes[END_REF]. In addition, the obtained meshes are not usually smooth and suffer from local surface distortions due to image and silhouette noise. As they often make use of surface metrics, such as curvature, which are extremely sensitive to noise, or they assume genus zero shapes [START_REF] Katz | Mesh segmentation using feature point and core extraction[END_REF], the most common segmentation approaches in computer graphics are limited in their applicability towards visually acquired meshes. The same remarks hold for volumetric representations of the object(s) of interest as sets of voxels.

More recently, temporal segmentation of visually reconstructed meshes has been addressed as well. [START_REF] Varanasi | Temporally coherent segmentation of 3D reconstructions[END_REF] proposes a method that performs both convex segmentation of a static mesh and temporally-coherent segmentation of a mesh sequence. [START_REF] Franco | Learning Temporally Consistent Rigidities[END_REF] proposes a generative probabilistic method that alternates between segmenting a mesh into rigid parts and estimating the rigid motion of each part along time. These methods perform temporal segmentation in the mesh/voxel domain and do not take advantage of the properties of spectral embeddings of articulated shapes.

Mathematically, both mesh and voxel representations can be thought of as graphs describing either the surface or the volume of a shape, in which each vertex of the graph represents a 3D point, while each graph edge connects vertices representing adjacent 3D points. Therefore, mesh/voxel segmentation can be addressed in the framework of graph partitioning. Spectral graph theory [START_REF] Cvetkovic | Spectra of Graphs: Theory and Applications[END_REF][START_REF] Chung | Spectral Graph Theory[END_REF] provides an extremely powerful framework allowing to cast the graph partitioning problem into a spectral clustering problem [START_REF] Belkin | Laplacian eigenmaps for dimensionality reduction and data representation[END_REF][START_REF] Luxburg | A tutorial on spectral clustering[END_REF]. In particular, it has been shown that the Laplacian embedding of a graph is well suited for representing the graph's vertices into an isometric space spanned by a few eigenvectors of the Laplacian matrix. The problem has been thoroughly investigated. Spectral methods allowing geometric representations of graphs are well established, and have been applied to automatic circuit partitioning for VLSI design [START_REF] Alpert | Spectral partitioning with multiple eigenvectors[END_REF], image segmentation [START_REF] Shi | Normalized cuts and image segmentation[END_REF], 2D point data [START_REF] Ng | On spectral clustering: analysis and an algorithm[END_REF], docu-ment mining [START_REF] Lafon | Diffusion maps and coarse-graining: A unified framework for dimensionality reduction, graph partitioning, and data set parameterization[END_REF], web data [START_REF] Fouss | Random-walk computation of similarities between nodes of a graph, with application to collaborative recommendation[END_REF], and so forth. Standard methods assume that the Laplacian matrix characterizing the data is a perturbed version of an ideal case in which groups of data-point are infinitely separated from each other. Links between Laplacian eigenmaps and linear locally embedding (LLE) [START_REF] Belkin | Laplacian eigenmaps for dimensionality reduction and data representation[END_REF], kernel PCA [START_REF] Bengio | Learning eigenfunctions links spectral embedding and kernel PCA[END_REF], kernel K-means [START_REF] Dhillon | Kernel kmeans: spectral clustering and normalized cuts[END_REF], random walks on graphs [START_REF] Meila | A random walks view of spectral segmentation[END_REF], PCA [START_REF] Saerens | The principal components analysis of a graph, and its relationships to spectral clustering[END_REF][START_REF] Fouss | Random-walk computation of similarities between nodes of a graph, with application to collaborative recommendation[END_REF] and diffusion maps [START_REF] Lafon | Diffusion maps and coarse-graining: A unified framework for dimensionality reduction, graph partitioning, and data set parameterization[END_REF] have been since established.

To the best of our knowledge [START_REF] Liu | Segmentation of 3D meshes through spectral clustering[END_REF] is the first attempt to apply spectral clustering to meshes using [START_REF] Ng | On spectral clustering: analysis and an algorithm[END_REF]. A specific feature of both meshes and voxel sets is, however, that the associated graphs are very sparse and that vertex connectivity is relatively uniform over the graph. Under these circumstances, no obvious partitioning of this kind of graph into strongly connected components that are only weakly interconnected exists: this makes mesh segmentation unfeasible via standard graph-partitioning or spectral-clustering algorithms [START_REF] Chung | Spectral Graph Theory[END_REF][START_REF] Ng | On spectral clustering: analysis and an algorithm[END_REF], in which the multiplicity of the zero eigenvalue of the Laplacian corresponds to the number of connected components. Firstly, in the case of meshes there is no "eigengap" [START_REF] Luxburg | A tutorial on spectral clustering[END_REF]: this makes difficult to estimate the dimension of the spectral embedding, and hence the number of clusters, in a completely unsupervised way. Secondly, the eigenvalues of any large semi-definite positive symmetric matrix can be estimated only approximately; this means that it is not easy to study the eigenvalue multiplicities which play a crucial role in the analysis of the Laplacian eigenvectors [START_REF] Biyikoglu | Laplacian Eigenvectors of Graphs[END_REF]. Finally, ordering the eigenvectors based on these estimated eigenvalues can be less than reliable [START_REF] Mateus | Articulated shape matching using laplacian eigenfunctions and unsupervised point registration[END_REF].

In [START_REF] Sharma | Mesh segmentation using Laplacian eigenvectors and Gaussian mixtures[END_REF], these issues were addressed based on an analysis of the nodal domains of a graph [START_REF] Biyikoglu | Laplacian Eigenvectors of Graphs[END_REF], and the interpretation of the eigenvectors as the principal components of a mesh [START_REF] Saerens | The principal components analysis of a graph, and its relationships to spectral clustering[END_REF]. The theory of nodal domains can be viewed as an extension of the analysis of the Fiedler vector [START_REF] Chung | Spectral Graph Theory[END_REF] to the other eigenvectors.

In order to be applicable to the scenario addressed in this paper, in which time series of measurements are exploited to estimate the pose of a moving object or discriminate between different types of motion, the obtained mesh/voxel segmentation has to be consistent throughout the temporal sequence. In particular, we are interested in surfaces or volumes representing articulated objects, i.e., sets of rigid bodies linked by articulations. Several approaches have been proposed which rely on the availability of sets of correspondences between pairs of volumes/surfaces. For instance, geodesic distances on both meshes and volumes are invariant to pose and articulation, and can be used for matching surfaces which differ by a non-rigid deformation [START_REF] Bronstein | Generalized multidimensional scaling: A framework for isometry-invariant partial surface matching[END_REF]. However, they are very sensitive to topological changes (which, as we recalled, often occur in visually reconstructed meshes). [START_REF] Bronstein | Topology-invariant similarity of nonrigid shapes[END_REF] uses a metric based on diffusion over the surface through a heat kernel, which is less sensitive in this regard. Functions defined on the mesh-volume are more resistant to topological changes. [START_REF] Shapira | Consistent mesh partitioning and skeletonisation using the shape diameter function[END_REF] uses the local "thickness" of the mesh-volume to consistently partition different surfaces over articulations and pose changes. These and other works in the geometry processing community use surface curvature to derive the final segmentation, which is fine for smooth 3D models but is not reliable for visually reconstructed meshes. An example of these matching based approaches is [START_REF] Golovinskiy | Consistent segmentation of 3d models[END_REF] which consistently segments sets of objects such as chairs and airplanes, starting from point-wise correspondences via global ICP registration, without exploiting the temporal information.

Other matching based approaches rely on the use of a graphical model to explicitly match two surfaces differing from a non-rigid deformation [START_REF] Anguelov | The correlated correspondence algorithm for unsupervised registration of non-rigid surfaces[END_REF][START_REF] Starck | Correspondence labeling for widetimeframe free-form surface matching[END_REF]. These point correspondences can be later used to recover the articulated structure of the shape, as proposed in [START_REF] Anguelov | Recovering articulated object models from 3D range data[END_REF]. [START_REF] Chang | Automatic registration for articulated shapes[END_REF] proposes an interesting alternative to the graphical model approach: it computes a putative set of rigid transformations between surface points on two articulated shapes, and treats these transformations (instead of the points themselves) as the labels for the graphical model.

A general limitation of model-based approaches to mesh/voxel segmentation is when dealing with unknown scenes, as in the case of multiple interacting actors. These limitations can be addressed by exploiting the temporal consistency between individual reconstructions. Preliminarily computing sets of correspondences between surface or volume points can greatly help to enforce such consistency. Unfortunately, surface/volume matching is a very hard problem, its results suffering from the presence of surface occlusions, sensitivity to topological noise and disjoint objects.

In this paper we focus on volumetric representations, and propose an approach to automatic, unsupervised spectral segmentation of moving bodies (represented as sets of stereoreconstructed voxels) along entire sequences, in a temporally coherent, robust way, as a means of both inferring a bottom-up model of an unknown articulated body and providing motion cues that can be later used for motion/action recognition, e.g., Figure 1 (Throughout the paper, the voxel representations are visualized as point clouds.) In this endeavor we rely neither on any prior model of the moving body nor on sets of previously computed point correspondences, while accounting for topological artifacts.

Recent attempts to solve the problem at hand proposed manifold learning techniques to process spatiotemporal data, e.g., [START_REF] Jenkins | A spatio-temporal extension to isomap nonlinear dimension reduction[END_REF][START_REF] Lin | Learning nonlinear manifolds from time series[END_REF]; these methods rely on enforcing temporal relationships when embedding time sequences, a hard task when handling dense volume representations. We propose instead a mechanism to enforce temporal consistency of segments obtained by clustering shapes in a spectral-embedded space, in which clusters are both collinear and remarkably stable under articulated motions, and propagate them over time. Within the class of spectral embedding techniques that are available today we favor local linear embedding (LLE) [START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF] for a number of desirable features LLE exhibits in the specific scenario of unsupervised segmentation: while it conserves shape protrusions in virtue of its local isometry, their separation is increased and their dimensionality reduced as an effect of the covariance constraint. In pose estimation, learning-based or example-based approaches [START_REF] Brand | Shadow puppetry[END_REF][START_REF] Grauman | Inferring 3D structure with a statistical image-based shape model[END_REF][START_REF] Elgammal | Inferring 3D body pose from silhouettes using activity manifold learning[END_REF] directly relate visual information (features) to learned body configurations without the need for accurate a priori models of the studied shape. These methodologies are not affected by the initialization issue, but are limited by the use of training sets of examples. In contrast, techniques have been proposed that directly infer body poses from multiple image cues or volume sequences: when the articulated structure of the moving body is not known at all, the automated recovery of at least a skeleton-like model can be useful prior to any pose estima-tion step. Skeletonization methods recover the intrinsic articulated structure of 3D shapes, either directly in 3D [START_REF] Brostow | Novel skeletal representation for articulated creatures[END_REF], or in an embedded space [START_REF] Chu | Markerless kinematic model and motion capture from volume sequences[END_REF][START_REF] Sundaresan | Segmentation and probalistic registration of articulated body models[END_REF]. As it is able to automatically recover the basic structure of an unknown body in terms of protrusions, and to do so consistently along a sequence, the spectral segmentation framework that we propose can be seen as a preliminary step in which both a time-invariant stick-like structure and its configuration in the different time instants are recovered.

In action recognition, on the other hand, cues need to be extracted from a video sequence in order to recognize the class of motion being performed. While many successful approaches do not rely on dynamics but favor extracting clues from spatiotemporal volumes [START_REF] Blank | Actions as space-time shapes[END_REF][START_REF] Kim | Canonical correlation analysis of video volume tensors for action categorization and detection[END_REF], others are based on extracting features, instant-by-instant, and explicitly encode motion dynamics via a graphical model [START_REF] Bissacco | Classification and recognition of dynamical models: The role of phase, independent components, kernels and optimal transport[END_REF][START_REF] Chaudhry | Histograms of oriented optical flow and Binet-Cauchy kernels on nonlinear dynamical systems for the recognition of human actions[END_REF][START_REF] Gupta | Understanding videos, constructing plots learning a visually grounded storyline model from annotated videos[END_REF] such as, for instance, a hidden Markov model [START_REF] Shi | Discriminative human action segmentation and recognition using semi-Markov model[END_REF][START_REF] Carter | Supplementing Markov chains with additional features for behavioural analysis[END_REF]. In this second case, one has to decide whether to extract features in 2D from individual images, or use multiple-stereo techniques to produce volumetric/surface reconstructions of the moving body. The availability of systems of multiple, synchronized cameras has made this approach both feasible and promising [START_REF] Cuzzolin | Action modeling with volumetric data[END_REF][START_REF] Cuzzolin | Invariant action classification with volumetric data[END_REF]. As the recovered body protrusions are consistently tracked over time along a sequence, the framework proposed here can be easily exploited as a preprocessing stage in which, for example, the 3D locations of the centroids of the recovered segments are used as observations to be fed to a parameter estimation algorithm (such as Baum-Welch for HMMs) in order to identify the performed motion class.

This paper is an extended version of [START_REF] Cuzzolin | Robust spectral 3D-bodypart segmentation along time in Workshop on Human Motion -Understanding[END_REF] and [START_REF] Cuzzolin | Coherent Laplacian 3-D protrusion segmentation In Computer Vision and Pattern Recognition[END_REF]. Here we outline the LLE method, which is essential for understanding our temporal-coherent segmentation technique, and we provide an interpretation of LLE in terms of geometric features. A physical interpretation of LLE can be found in [START_REF] Cuzzolin | Coherent Laplacian 3-D protrusion segmentation In Computer Vision and Pattern Recognition[END_REF]. The segmentation algorithm, initially proposed in [START_REF] Cuzzolin | Robust spectral 3D-bodypart segmentation along time in Workshop on Human Motion -Understanding[END_REF] and sketched in [START_REF] Cuzzolin | Coherent Laplacian 3-D protrusion segmentation In Computer Vision and Pattern Recognition[END_REF], is described below in more detail. The experimental section contains extensive quantitative evaluations based on synthetic data, which were segmented automatically, as well as more results with real sequences and more comparisons with competing methods. Moreover, with respect to [START_REF] Cuzzolin | Robust spectral 3D-bodypart segmentation along time in Workshop on Human Motion -Understanding[END_REF][START_REF] Cuzzolin | Coherent Laplacian 3-D protrusion segmentation In Computer Vision and Pattern Recognition[END_REF] we provide a sensitivity analysis with respect to the LLE free parameters and, in particular, we discuss the algorithm robustness with respect to these parameters, to topological changes in the moving shape, and to the resolution of the voxel grid. We propose an interpretation of protrusion segmentation in terms of the nodal domains of the Laplacian eigenvectors. This in turn provides insights on how to choose the dimension of the embedded space. An application to bottom-up model recovery is proposed and illustrated as well.

Paper Organization

The remainder of this paper is organized as follows. Section 2 starts by outlining the proposed approach, its objectives and the problem constraints. We motivate our approach by the interesting geometric features of LLE. Section 3 illustrates the proposed algorithm in detail, step by step: the use of k-wise clustering to segment the embedded cloud, starting from detected branch terminations (section 3.2), followed by how to propagate clusters over time to ensure consistency (section 3.3) and merge/split them to fit the topology of the body (3.4). The whole procedure is summarized in section 3.5. Section 4 describes in detail an extensive set of experimental validations that provides evidence, on both synthetic and real data, on the performance of the algorithm and in comparison with competing methods. The results are benchmarked with ground-truth data: therefore we thoroughly test the way our method and other methods can cope with topology transitions. Moreover we study the robustness of the algorithm. Section 5 outlines potential applications to bottom-up model recovery. Section 6 concludes the paper.

Problem Formulation and Proposed Solution

The main objective of this paper is a method for performing unsupervised and temporally coherent protrusion clustering. The proposed methodology should be able to produce an automatic segmentation of an unknown moving body. Such a body can be represented indifferently as a set of surface points or a set of voxels obtained by volumetric reconstruction from a set of synchronized cameras: here, nevertheless, we will pay special attention to volumetric representations. Without loss of generality, the body is supposed to be composed by a collection of linked rigid parts (a condition that could be possibly relaxed). The produced segments are "meaningful", in the sense that they are linked to actual segments of the moving body. The segmentation process takes place in a completely unsupervised way, i.e., neither with previously annotated data nor with human intervention. The segmentation algorithm processes an entire sequence of surface/volumetric shapes in a consistent way, i.e., the shape is segmented into the "same" elements throughout the sequence. The methodology is able to adapt to topological changes due to imperfections or self-contacts. Throughout the paper, the following constraints are assumed:

no a priori model is available for the moving shape, but instead some model of the latter can be recovered a posteriori, and no correspondences between sets of surface/volume points are necessary to ensure the consistency of the segmentation (freeing the approach from the limitations of matching algorithms).

The remainder of this section explains how spectral segmentation, once LLE is applied, yields geometric features that allow us to design such a system, to be later used for model/pose estimation or action recognition (Section 5), and outline an approach to the above defined problem based on it. Then we will detail the various steps of the process, and validate it in a variety of scenarios.

Locally Linear Embedding

LLE [START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF] is a graph Laplacian [START_REF] Zhou | Large mesh deformation using the volumetric graph Laplacian[END_REF] algorithm which computes the set of d-dimensional embeddings Y = {Y i , i = 1, .., N} of a set of input points X = {X i , i = 1, .., N}, while preserving their local structure. For each data point X i , the algorithm computes the weights {W i j , j = 1, ..., N(i)} that best linearly reconstruct X i from its k nearest neighbors {X j , j = 1, .., k} by solving the following constrained least-square problem: arg min

∑ N i=1 |X i -∑ k j=1 W i j X j | 2 .
{Y i } N ∑ i=1 Y i - k ∑ j=1 W i j Y j 2 .
(

) 1
The embedded cloud Y is constrained to be centered at the origin

0 of R d (∑ i Y i = 0) and to have unit covariance: 1 N N ∑ i=1 Y i Y T i = I. ( 2
)
The objective function to minimize ( 1) is in fact a quadratic form ∑ N i, j=1 M i j Y i ,Y j (where . denotes the usual scalar product of vectors) involving the "affinity" matrix:

M . = (I - W ) T (I -W ), W = [w i j ].
The optimal embedding (up to a global rotation) is found by computing the bottom d + 1 eigenvectors of M, and discarding the bottom (unit) one.

Geometric Features of LLE

Despite having been originally introduced as an unsupervised dimensionality reduction technique, LLE displays a number of attractive geometric features for unsupervised segmentation as well, as it is illustrated by Figure 3. Given a cloud of 3D points representing, for instance, the voxels occupying the volume of an articulated shape, for a large interval of its parameters k (scale of the neighborhood to preserve) and d (dimensionality of the embedding space), the corresponding embedded cloud after LLE exhibits two interesting properties:

1. it is (roughly) lower dimensional, i.e, the embedded points live (approximately) on a manifold of dimension lower than the original one (a 1-D tree in Figure 3-a); 2. the protrusions (chains of smaller rigid links) in the original shape are mapped to roughly equally spaced, widely separated branches of the embedded cloud.

Intuitively, this is due to the fact that the covariance constraint (2) pulls outwards and stretches the chain of links formed by all local neighborhoods, redistributing them on a (roughly) lower dimensional manifold. As the force pulls those chains in a radial direction, their separation increases in the embedding space. Clustering/segmenting such protrusions becomes much easier there: compare in Figure 3-a)-right the output of the ISOMAP embedding on the same shape. However, our purpose is to cluster moving objects in a temporally coherent way: we need the segments obtained at different time instants to be consistent, i.e., to cover the same body-parts. Some embedding schemes, such as ISOMAP, are inherently pose-invariant under articulated motion (as, while the articulated body evolves, the geodesic distances between all pairs of points do not vary). This is not true, in a strict sense, for LLE and other graph Laplacian embeddings [START_REF] Belkin | Laplacian eigenmaps for dimensionality reduction and data representation[END_REF]. However, for articulated bodies formed by a number of rigid parts linked by rotational joints, all local neighborhoods incident on a rigid part are preserved along the motion, while only the few neighborhoods interested by the evolving joint(s) are affected (Figure 3

The Proposed Formulation

Shape protrusions are then conserved in the LLE spectral domain, while their separation is increased and their intrinsic dimensionality reduced. The effect is dramatic in situations in which body-parts are physically close to each other in the original space, and makes clustering significantly easier after embedding. Moreover, the embedded shape of an articulated body is to a large extent pose-invariant, making the propagation in time of the obtained clusters much easier.

Here we are going to exploit these geometric features of LLE to devise an unsupervised, time-consistent method for protrusion segmentation of 3D articulated shapes. The proposed segmentation technique does not assume any prior model (not even a weak topological or graphical one). The number of clusters themselves is unknown while it is inherently determined by the lower dimensional structure of the embedded cloud. The proposed method follows the following steps:

1. At each time-instant the current 2D voxel-cloud is mapped via LLE to a suitable embedding space; 2. The most suitable number of clusters is estimated by counting the number of branches in the embedded cloud; 3. The k-wise clustering technique devised in [START_REF] Agarwal | Beyond pairwise clustering[END_REF] is used to segment such branches, due to their intrinsic lower dimensionality; 4. Such a embedded-space clustering induces protrusion segmentation in the original 3D space, and finally 5. The resulting segmentation is propagated to the next time-instant to ensure temporal coherence of the segmented body.

In the next section we are going to provide a detailed description of each step.

Method

As articulated shapes are mapped by LLE onto an embedded space, it would make no sense to employ generic k-means [START_REF] Macqueen | Some methods for classification and analysis of multivariate observations[END_REF] to segment embedded-space protrusions (as it happens in classical spectral clustering [START_REF] Ng | On spectral clustering: analysis and an algorithm[END_REF]). We therefore adopt the k-wise clustering introduced in [START_REF] Agarwal | Beyond pairwise clustering[END_REF] to segment protrusions in the embedding space. Given a sequence of 3D shapes, the initial segmentation is consistently propagated along time, and the number of clusters estimated in an automatic way to fit the changing topology of the moving articulated body.

K-wise Clustering of the Embedded Shape

Consider Figure 3. For a dimension d = 3 of the LLE embedding space, the embedded cloud is (for a wide interval of values of k) a tree-like one-dimensional curve. When trying to segment its branches, then, it is natural to look for clusters formed by sets of roughly collinear points. Traditional clustering algorithms, such as k-means [START_REF] Macqueen | Some methods for classification and analysis of multivariate observations[END_REF], are based on measuring pairwise distances between data points. In opposition, it is possible to define measures of similarity between triplets of points, which signal how close these triplets are to being collinear (Figure 5-left) [START_REF] Heiser | Triadic distance models: Axiomatization and least squares representation[END_REF][START_REF] Hayashi | Two dimensional quantification based on the measure of dissimilarity among three elements[END_REF]). More generally, the problem of clustering points based on similarity between k-tuples of points is called k-wise clustering. An interesting approach to k-wise clustering has been proposed in [START_REF] Agarwal | Beyond pairwise clustering[END_REF], based on the notion of "hypergraph". A weighted undirected hyper-graph H is a pair (V, h), where V is the set of vertices of H. Subsets of V of size k are called hyper-edges (just as subsets of two elements in a normal graph are potential edges). The function h associates nonnegative weights h(z) with each hyper-edge (k-tuple) z, and measures the affinity of each hyper-edge (Figure 5-right), in analogy with the weight of edges in conventional graphs.

Consider then a set of data-points Y = {Y i , i = 1, ..., N}.

1. First, an affinity hyper-graph H is built by measuring the affinity of all the k-tuples of points in Y ; 2. Next, a weighted graph G with the same vertices Y approximating the hyper-graph H is constructed by constrained least square optimization; 3. Finally, the approximating graph G is partitioned into n parts via a spectral clustering algorithm [START_REF] Shi | Normalized cuts and image segmentation[END_REF][START_REF] Ng | On spectral clustering: analysis and an algorithm[END_REF].

In the problem of interest to us, the hyper-graph to approximate has as vertices the embedded points {Y i }. If the embedding is d-dimensional, d-tuples of embedded points are considered as hyper-edges: in particular, hyper-edges are triads of points when d = 3. A natural choice for the affinity of these triads of embedded points is the area of the triangle they form, e.g., Figure 5-left, or the volume of the (d -1)-dimensional hyper-edge in the general case. The outcome of the above clustering algorithm is a segmentation of the embedded cloud Y , which can be trivially transferred back to the original 3D space using the known correspondences between the 3D points and the embedded-space points. A mechanism for automatically estimating the number n of clusters is the next step in the algorithm.

Branch Detection and Number of Clusters

The fact that embedded clouds typically appear as one-dimensional strings formed by a number of branches (corresponding to the extremities of the moving body) provides us with a simple method to estimate, at each time instant, the most suitable number of clusters (Figure 6).

Each point of the embedded cloud is tested to decide whether or not it is a branch termination. This test is performed by looking for its nearest neighbors, for a certain threshold distance empirically learned from the data. The best interpolating line (in blue) for all the neighbors (in green in Figure 6) of the considered point is found, and all such neighbors are projected onto it. A point of the embedded cloud (red star) is considered a branch termination if the projections of all its neighbors on the interpolating line lay on one side of its own projection (red square), e.g., Figure 6-left. It is not considered a termination when its projection has neighbors on both sides (Figure 6-right). This algorithm proves to work extremely well on embedded clouds generated through LLE. It becomes then possible to detect transitions in the topology of the moving body whenever they happen, and modify the number of clusters accordingly. We will return on this in Section 3.4.

Temporal Consistency and Seed Propagation

When considering entire sequences of 3D clouds we need to ensure the temporal consistency of the segmentation: in normal situations (no topology changes due to contact of different body-parts) the cloud has to be decomposed into the "same" groups in all instants of the sequence. We propose a propagation scheme in which centroid clusters at time t are used to generate initial seeds for clustering at time t + 1 (Figure 7). Let n be the number of clusters.

The Seed Propagation Algorithm

1. The embedded cloud {Y i (t), i = 1, ..., N(t)} at time t is clustered using d-wise clustering (Section 3.1, Figure 7-bottom-left); 2. for each cluster centroid c j (t) ( j = 1, ..., n) in the embedding space (a square in Figure 7-bottom-left), the original datapoint, denoted by X i j (t), whose embedding Y i j (t) is the closest neighbor of c j (t) (and which we call "3D centroid") is found (a square of the same color in Figure 7-top-left):

i j (t) = arg min i=1,...,N Y i (t) -c j (t) 2 ;
3. at time t + 1 the dataset of 3D input points X(t + 1) = {X i (t + 1), i = 1, ..., N(t + 1)} is augmented with the positions of the old "3D centroids" at time t (colored circles), yielding the set (Figure 7-top-right): X ′ (t + 1) = X(t + 1) ∪ {X i j (t), j = 1, ..., n}; 4. LLE is applied to the extended dataset X ′ (t + 1), obtaining (Figure 7-bottom-right) 1 Y (t + 1) ∪ {c ′ j (t + 1), j = 1, ..., n}, where c ′ j (t + 1) is the embedded version of the old 3D centroid X i j (t) in the new embedded cloud at time t + 1 (a colored circle in Figure 7-bottom-right); 5. such embedded points c ′ j (t + 1) are then used as seeds from which to cluster via k-wise the new embedded cloud Y (t + 1).

Topology Changes and Dynamic Clustering

The question of how to initialize the seeds at t = 0 naturally arises. Besides, even though working in an embedding space helps to dramatically contain the issue with segmenting body-parts which get dangerously close to each other, instants in which different parts of the articulated body come into contact still have important effects on the shape of the embedded cloud. In fact, in an unsupervised context in which we do not possess any prior knowledge about the number of rigid parts which form the body or the way they are arranged, but only 1 The embeddings c ′ j (t + 1) of the previous 3D centroids X i j (t) in the new embedded cloud can also be computed by out of sample extension [START_REF] Bengio | Out-of-sample extensions for LLE, Isomap, MDS, eigenmaps, and spectral clustering[END_REF].

the fact that they belong to an articulated shape, there is no reason to distinguish touching body-parts. It appears more sensible to adapt the number of clusters to the number of actually distinguishable protrusions.

The branch detection algorithm of Section 3.2 provides, on one side, a suitable tool for initializing the embedded clustering machinery, while allowing at the same time the framework to adapt the number and location of the clusters when a topology change occurs.

The Cluster Merging-Splitting Algorithm

1. At each time instant t all the branch terminations of the embedded cloud Y (t) are detected; if t = 0 they are used as seeds for d-wise clustering; 2. otherwise (t > 0), first standard k-means with k = d is performed on Y (t) using branch terminations as seeds, yielding a rough partition of the embedded cloud into distinct branches; 3. if two or more propagated seeds c ′ j (t) (see Section 3.3) fall inside the same element of this preliminary partition obtained via k-means, they are replaced by the branch termination of the element of the preliminary partition; 4. for each element of the preliminary partition of Y (t) which does not contain any propagated seed, a new seed is defined as the related branch termination; 5. finally, d-wise clustering is applied to the resulting set of seeds.

Step 3 takes place when previously separated protrusions get too close (in the original 3D space) to be distinguished: it makes then sense to merge the corresponding clusters. Step 4 represents the opposite event in which a body-part which was previously impossible to distinguish becomes well separated, requiring the introduction of a new cluster. As a result, clusters are allowed to merge and/or split according to topological changes in the moving articulated body. We will validate this approach to topology change management in Section 4.7.

Summary of the Method

Let us at this point summarize our approach for unsupervised, robust segmentation of parts of moving articulated bodies in a consistent way along a time sequence, by integrating the separate algorithms of Sections 3.1, 3.2, 3.3, 3.4 into a coherent whole. Given the values d, k of the parameters of LLE, for each time instant t: if t = 0, we use as seeds the branch terminations; if t > 0, seeds are obtained from the old propagated centroids {c ′ j (t), j = 1, ..., n(t -1)} and the new branch terminations via the splitting/merging procedure described in Section 3.4; 4. this yields a new set of centroids {c j (t), j = 1, .., n(t)} in the new embedding space Y ′ (t); 5. the labeling of the embedded points induces a segmentation in the original 3D shape (Figure 8-d); 6. the new cluster centroids c j (t) are re-mapped to 3D (Section 3.3), the corresponding 3D centroids X i j (t) are added to the data-set X(t + 1) at time t + 1 (Figure 7); 7. we go back to step 1, until the last time instant of the sequence is attained.

1. the current data-set X ′ (t) = X(t) ∪ {X i j (t -1)} for t > 0, X ′ (0) = X(0) = {X i (0), i = 1

Experiments

Experimental Setup

We tested the algorithm of Section 3.5 on both synthetic and real data, in order to have both qualitative and quantitative assessments on its performances.

In the first experiment (Section 4.3) we tested our algorithms on a set of synthetic sequences depicting a moving person, for which a "ground truth" segmentation of the person's body into rigid links was available, e.g., Figure 9. The sequences were generated by simulating the evolution of a human body model formed by a kinematic model [START_REF] Knossow | Human motion tracking with a kinematic parameterization of extremal contours[END_REF], plus the volume elements representing its rigid links. As it can be seen in Figure 9, the grid sampling was such that the height of a standing person was 200 voxels, corresponding to a real-world resolution of around 1cm.

In a second experiment (Section 4.4) we used a multiple-camera setup composed by eight synchronized and cross-calibrated video cameras to capture images sequences in a constrained environment, such that the silhouette of a moving person could be easily extracted from each image. Using the resulting silhouettes, a standard space-carving algorithm was employed to generate at each time instant a voxel-based representation of the person. The frame rate of the cameras was 30 frames per second. The original voxel size was 1cm. Figure 10 shows as an example a few frames from the "dancer", one of the real-world voxelset sequences, and a number of images for the same sequence 1 . As it can be appreciated, the original resolution of the voxel-set grid is quite high. In the case of sequences for which motion capture estimates were not available, and could not be used to provide ground truth on the "true" segmentation into rigid links, we could nevertheless provide visual results which provide qualitative evidence on the effectiveness of the approach.

In a third experiment (Section 4.5) we were able to calculate performance scores for real-world sequences for which motion capture data was available, and a ground truth segmentation of the moving person could be built. In all these experiments we compared our results with those of similar schemes in which seeds are also passed to the next frame to ensure time consistency, but clustering is either performed in 3D on the original data-set using a Gaussian mixture and the EM algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the em algorithm[END_REF] (time-consistent EM clustering) or in the ISOMAP space [START_REF] Tenenbaum | A global geometric framework for nonlinear dimensionality reduction[END_REF] using k-means (time-consistent ISOMAP clustering).

In Section 4.6 we further analyzed the sensitivity of the algorithm to changes in the values of its basic parameters: the set of eigenvectors selected to compute the desired LLE embedding, the size k of the neighborhood, and the resolution of the grid of voxels. Finally, in Section 4.7 the robustness with respect to topological changes in the moving shape was assessed.

Fig. 11 Left: Ground-truth labels for voxels are automatically generated from an articulated model of a virtual human. Each voxel is assigned the label of the closest generating link. Right: different performance scores can be calculated by comparing the unsupervised segmentation produced by an algorithm with the three a-priori segmentations shown here (black lines in the plots).

Ground Truth and Performance Scores

In order to quantitatively assess the performance of our segmentation algorithm we produced ground-truth segmentation for the analyzed sequences. In the case of synthetic sequences we used as ground-truth the labels automatically generated by associating each (synthetic) voxel with the closest link of the model which had generated the voxel set in the first place, e.g., Figure 11-left.

We worked out three distinct performance measures. Firstly, for a body composed by a number of rigid parts, a segmentation is valid if it does not cut in half a rigid link, i.e., the segmentation is a coarsening of the partition of all rigid links. We defined as coarsening score the average percentage of points across all rigid links which belong to the majority cluster (the unsupervised cluster containing the largest number of points, whichever it is) for the link.

Secondly, the obtained segmentation can be compared with three different "natural" subdivisions of the human body, e.g., Figure 11-right, by counting for each a priori segment s the histogram distribution hist s (l) of the different unsupervised labels l, and retain the percentage of "majority" voxels, i.e., those associated with the unsupervised cluster with the highest frequency: score(s) = max l hist s (l). A segmentation score can then be obtained as the simple average score for all segments: segm = 1 #s ∑ s score(s). When the obtained unsupervised clusters correspond to the a-priori segments the largest cluster contains 100% of the points for each segment, and the score is equal to 1.

Finally, the consistency along time of the obtained unsupervised segmentation is measured by storing at t = 0 (as a reference) for each link of the body the percentage of voxels which belong to each unsupervised cluster, and measure for each t > 0 and for each link l the similarity of the current label distribution and the initial one, as 1 minus the L 1 distance between the two histograms. This consistency score tends to one when the proportion of unsupervised cluster labels is constant in time for all the ground-truth segments of the articulated body.

Experiment 1: Quantitative Performance on Synthetic Data

Figure 12 shows the segmentation scores for six different synthetic sequences, of length comprised between 25 and 70 frames. We can appreciate how, typically, the obtained unsupervised segmentation turns out to be very consistent along time, as witnessed by a consistency score (in red) between 95% and 100% at all times, even for fairly long sequences Fig. 12 Segmentation scores obtained by comparing the labeling generated by our segmentation algorithm with ground truth labels automatically provided for a number of synthetic sequences of articulated motion. Scores obtained over six different sequences of different length from 25 to 70 frames are shown. From top left to bottom right: "walk", "walk on mars", "waking up", "laughter", "surf1", and "surf2". Red: consistency score. Blue: coarsening score. Solid, dash-dot and dashed black lines plot the segmentation scores respectively associated with the three a-priori segmentations of Figure 11, taken in the same order. picturing quite different types of motion. In all cases the boundaries between unsupervised clusters normally lie in correspondence of actual articulations, as witnessed by the values of the coarsening score (in blue). As for the a-priori segmentation scores (the three different black curves) we can notice that, at least for the given collection of sequences, the designed unsupervised spectral clustering algorithm seems to favor the second a-priori partition of Figure 11-right, i.e., it tends to highlight the outermost rigid links rather than whole protrusions (such as legs and arms). This unexpected result can be explained by pointing out that lower density regions (such as articulated joints) cause the embedded cloud to bend.

Figure 13 shows the corresponding results produced by time-consistent EM clustering in the original 3D space. Not only the absolute segmentation performance (measured by the black plots) is consistently, dramatically worse than that of the proposed algorithm, but the obtained segmentation is not at all consistent along time (red curves). Incidentally, EM segmentation appears to relatively favor the first a-priori segmentation (whole legs and arms), as the dominance of the solid black curve attests. Clusters drift inside the shape during the motion, and usually span different distinct body-parts. This phenomenon is made clear in Figure 14, where the irregular trajectories of the obtained 3D clusters for two sub-sequences of "walk" and "surf" are plot, while the obtained segmentation is visually rendered for two key-frames of the two sequences, showing its lack of relation with articulated segments.

Experiment 2: Qualitative Assessment on Real Data

In a second series of tests we applied the algorithm of Section 3.5 to several different realworld, high-resolution sequences of voxelsets, generated by our multi-camera acquisition system. Among those, a 200-frame-long sequence capturing a dancer who moves and swirls all around the scene, and voxelset resolution of 1cm.

Figure 15 illustrates some typical results of the dynamic segmentation algorithm described in Section 3.5, for a number of real-world sequences. It can be appreciated how the resulting unsupervised segmentation turns out to be pretty consistent in time, yielding very smooth cluster centroid trajectories, not only in situations where body-parts are well separated ( a), b)) but also during walking gaits (c) or even extremely complicated motions like the dance performed in d), e). Figure 15-f) shows how the algorithm copes with transitions in the body's topology, in particular the arm of a walking person touching their torso. Cluster evolution displays remarkable stability before and after the transition.

Experiment 3: Quantitative Assessment on Real Sequences

When motion capture data are available it is possible to measure quantitatively the performance of the three competing methods (the proposed time-consistent d-wise LLE segmentation, time consistent EM clustering, and time-consistent ISOMAP embedding with K-means clustering) on real-world sequences as well, by assigning as ground-truth to each voxel the label of the closest estimated link.

Voxelset sequences actually captured through a multi-camera system suffer from a number of unpleasant phenomena, like presence of large gaps or holes in the grid of voxels (e.g. the gap in the wrists in Figure 16-middle), noise, disconnected components, not to We can still notice a couple of drops in the scores (top-left), mirrored by a brief sudden glitch in clusters trajectories, (top-right). They correspond to frames (such as t = 408, whose voxelset segmentation is shown in the middle) in which gaps are so wide that they affect the quality of the segmentation, even though the shape of the embedding cloud remains stable.

Figure 17 shows the scores obtained by all three competing methods over two other challenging sequences of real voxelsets. It is apparent how our method exhibits strong resilience to data of very poor quality, easily outperforming both EM clustering in 3D or k-means clustering in the geodesic-based ISOMAP space. At times glitches due to extremely corrupted data (t = 8,t = 17, right) appear, but the topology adaptation algorithm of Section 3.4 brings swiftly the segmentation back on track.

Sensitivity analysis

Estimating the Optimal Number of Neighbors

The proposed segmentation methodology critically relies on the "desirable" geometric properties of LLE discussed in Section 2.1, which in turn depend on its two basic parameters: the size k of the neighborhoods and the dimension d of the embedding space. The former, in particular, affects both the stability of the embedded shape along time and its lowerdimensionality, from which the estimation of the number of clusters itself depends. It can be noticed empirically that, while the embedded shape shows a remarkable stability for some values of k, this is not in general the case for arbitrary such values. Figure 18-a) shows how the embedded shape varies for two frames (top, bottom) in the "dancer" sequence as k ranges from 10 to 30. Fig. 18 a) This shows how the lower-dimensionality, separability and stability of the embedded cloud along a sequence depend on the size of the neighborhood k in the LLE algorithm. b) Unsuitable values of k are characterized by "anomalous" neighborhoods which span distinct body-parts (middle), in opposition to the case of admissible values (left). Right: plots of the distance between the farthest point of the neighborhood and all the others in the two cases (top: "regular" neighborhood; bottom: "anomalous" neighborhood).

For higher values of k some neighborhoods of points in a given body-part comprise regions of a different body-part, e.g., Figure 18-b). However, in those "anomalous" neighborhoods the farthest element (as it belongs to another, distinct link) is relatively distant from all others. If we plot the distance between this point and all its fellows we can notice a large jump, e.g.,Figure 18-b) (bottom right). This is not the case for "regular" neighborhoods spanning a single rigid part (top right). We can then set as acceptable value for k any of those which yield only "regular" neighborhoods.

Fig. 19 As stationary functions on the shape approximated by a cloud of points, Laplacian eigenfunctions are associated with symmetries and protrusions of the underlying domain. Here the level sets (in blue) of the first five eigenfunctions of the depicted "monster" are shown to partition the toy along directions that follow its natural axes of symmetry: horizontal (leftmost image), left/right (second picture), or radiating from the barycenter (last image). .

Laplacian Methods, Selection of Eigenfunctions and Body-part Resolution

Spectral methods have long been used in the computer graphics community as a tool to analyze and process meshes representing surfaces of objects inside virtual scenes. Generally speaking, they share a common framework in which an affinity matrix M which reflects the structure of the input set of points is defined and an eigen decomposition of this matrix which yields its eigenvalues and eigenvectors is performed. They can be roughly classified into two categories: methods in which M encodes adjacency information (such as LLE or ISOMAP), and methods based on the graph Laplacian, an operator L : f → L f mapping functions f : X → R, X i → f (X i ) = f i defined on sets of points (vertices) X = {X i , i = 1, ..., N} forming a graph, of the form:

(L f ) i ∝ ∑ j∈N(i) w i j • ( f i -f j ) (3)
where N(i) is the set of neighbors of the point X i (the vertices connected to it by an edge) in X, and w i j is the weight of the edge joining X i with j-th neighbor. Functions on a set X of N points can be trivially represented as the vector of their values on the N points: f = [ f (X i ) = f i ] T . Therefore, the N by N affinity matrix M computed by LLE can also be seen as an operator on the space of functions defined on the cloud of N points X, as the matrix multiplication f ′ = M × f yields another vector f ′ of N components, representing another function on X. As it can be proven that the affinity matrix operator can be approximated by the square Laplacian, M f = (I -W ) T (I -W ) f ≈ 1 2 L 2 f , Locally Linear Embedding can be considered a form of Laplacian embedding. Now, the level sets of Laplacian eigenfunctions are related to the geometry of the underlying shape X. In particular nodal sets [START_REF] Jakobson | Geometric properties of eigenfunctions[END_REF], i.e., zero-level sets of eigenfunctions, partition X into a set of "nodal domains" strongly related to protrusions and symmetries of the underlying grid of points. Figure 19 shows a clear example of how the level sets of different eigenfunctions of the Laplacian operator defined on (the discretized version of) any given shape follow the symmetries of the shape itself [START_REF] Levi | Laplace-Beltrami eigenfunctions: Towards and algorithm that understands geometry[END_REF]. As the n-th eigenfunction has at most n nodal domains, the first eigenfunctions are associated with relatively "coarse" partitions of the shape of interest.

As the LLE affinity matrix is linked to the graph Laplacian operator, its eigenvectors are also associated with specific symmetries of the underlying cloud of points. In other words, the choice of which eigenvectors of M we select after SVD determines the structure of the embedded cloud. While in standard LLE always the bottom d are chosen, d being the dimension of the resulting embedding, LLE is used here just as a tool to make the structure of the 3D shape in terms of protrusions clearly emerge. It makes then more sense to look for an appropriate, arbitrary selection of eigenfunctions of the affinity matrix able to facilitate the clustering of the shape at hand.

To visually render these concepts, we can visualize functions f (and in particular eigenvectors of M) on the original cloud X of 3D points as "colored" versions of the same cloud, in which the color of each point X i is determined by the value f (X i ) of the function in that particular point.

We used this technique to depict in Figure 20-top the top four eigenvectors of the LLE affinity matrix for frame 209 of the "Mars" sequence. It is striking (but to be expected, following the brief discussion on graph Laplacian methods) to notice how the positive (red) and negative (dark blue) peaks of the different eigenfunctions of M are located exactly on the protrusions of the underlying shape, i.e., its high-curvature regions. Selecting a specific set of eigenfunctions is then equivalent to highlighting the associated peaks. The resulting embedded clouds may even possess a different number of distinguishable branches. In Figure 20-bottom-left, for instance, eigenfunctions 1,2,3 determine an embedding not able to resolve the head as a separate protrusion. The head appears instead as an additional cluster when selecting eigenfunctions 2,3,4 (right).

Influence of Voxel-grid Resolution

As we just argued, the quality of the segmentation depends on the "good" features of the embedding, namely its lower dimensionality and improved branch separation. In turn, those depend (amongst other factors) on the assumption that the graph Laplacian of the actual cloud of points representing the shape is a good approximation of the (unknown) Laplace-Beltrami operator.

Figure 21 illustrates how performances degrade as the number of voxels (seen as samples of an underlying continuous shape) decreases. To allow a fair comparison we kept roughly constant at all runs the product s • k between number of neighbors k and sampling factor s (as the original voxel sets contain as many as 25000 voxels, it is convenient to sub-sample them to limit computation time). In Figure 21 the consistency score and the segmentation score for the "natural" segmentation of the body into lower limbs, head, and bulk are shown for two subsequences of "space surf" and "reveil" (top and bottom, respectively). For more static sequences (bottom) the segmentation score (black) is more stable as a function of the product k • s, even though subject to sudden jumps. The consistency score (red) is more sensitive, while performance degrades consistently when grid resolution is reduced for sequences in which the body moves significantly (top). This is due to two different effects. On one side, for sparser clouds of points the graph Laplacian is a worse approximation of the Laplace-Beltrami operator. On the other, as the list of points is randomly sampled to produce the reduced dataset, the fewer the samples the less they happen to be distributed along the nodes of a regular 3D grid, as in the original voxel set. As this deforms the local structure of the neighborhoods, it causes irregular instability in the embedded cloud along the sequence which influences in particular the consistency score (red).

Robustness with Respect to the LLE Parameters

It is important to assess the sensitivity of the algorithm to the main parameters k and d (or, better, the indices of the selected eigenvectors). Figure 22 illustrates how the segmentation scores vary when different values of the parameters k, d are used to compute the embedding (assuming for sake of simplicity that we want to select the first d eigenfunctions). The consistency score, the segmentation score (w.r.t. the second a-priori segmentation of Figure 11-right) and their average are plotted as functions of k and d, for k = {14, 16, 18, 20, 22} and d = {3, 4, 5, 6, 7}, for two synthetic subsequences of "walk" and "surf". The stability of both consistency in time and quality of the segmentation in a fairly large region of the parameter space is apparent, and speaks for the robustness of the approach.

Robustness to Topological Changes

No matter how robust to pose variation LLE embedding may be, instants in which different parts of the articulated body come into contact still have important effects on the shape of the embedded cloud. These events have truly dramatic consequences on embeddings based on measuring geodesic distances along the body, since new paths appear, affecting in general the distance between all pairs of points in the cloud.

Figure 23 illustrates how the splitting-merging algorithm of Section 3.4 copes with such changes. Full consistency of the segmentation cannot be preserved, as the number of protrusions in the embedding shape changes, and we do not make use of point-to-point correspondences (which would themselves have trouble coping with such situations). Nevertheless, most protrusion end up been preserved through the event changing the body's topology, while the others can be recovered as soon as they can be distinguished once again.

Figure 24, instead, compares the segmentation scores of methods based on the local (LLE) and global (ISOMAP) structure of the shape in similar situations, in which the topology of the body changes, for two significant synthetic sequences ("wake up", left, and "clap", right). As the plots on top show, propagating clusters in the LLE space exhibits superior performance and robustness, as the algorithm smoothly adapts to topology changes in virtue of the geometric features of LLE embedding.

The bottom plots in Figure 24 show, as an example, the segmentations produced by kmeans in the ISOMAP space for the final frame of the two sequences. As ISOMAP embed- dings do not exhibit any detectable lower-dimensionality or separation between body-parts / protrusions, k-wise clustering is just not feasible there. Clustering is then performed in the ISOMAP space by k-means. Indeed, Figure 25 illustrates how ISOMAP (as a representative of all geodesic-based spectral methods) copes with the same transitions of Figure 23. The comparison clearly shows the advantage of doing clustering after Laplacian embedding.

An Application to Bottom-up Model/Pose Recovery

In our view, the unsupervised segmentation algorithm we propose can be seen as a building block of a wider motion analysis framework. We mentioned in the introduction two possible applications of time-consistent unsupervised segmentation: action recognition and bottom- up body-model recovery. As for the first one, a look at Figures 14,15, and 16 clearly shows the smoothness and coherence of the tracks associated with the centroids of the unsupervised clusters generated by our framework. Just as in the case of image features, we can assume that these tracks are generated by a graphical model, such as for instance an HMM [START_REF] Elliot | Hidden Markov models: estimation and control[END_REF] (see Figure 26), or possibly a sophisticated hierarchical, non-linear [START_REF] Ralaivola | Dynamical modeling with kernels for nonlinear time series prediction[END_REF], or even chaotic [START_REF] Ali | Chaotic invariants for human action recognition[END_REF] model. Classical algorithms can then be employed to estimate the parameters of this model in a generative approach to recognition: each test sequence is then classified according to the label of the training model which is the closest to the one that more likely generated it.

Another natural application of our unsupervised, time-consistent segmentation algorithm is recovering and fitting simple stick models to the segmented clusters along the sequence. As it provides a coherent protrusion segmentation along a sequence, and protrusions typically correspond to chains of rigid links, its output is suitable to reconstruct rough models of the moving body as a first step, for instance, of a model-free motion capture algorithm. However, here we do not aspire at providing a full solution to this problem. We can illustrate with a simple example of how this can be done. Ellipsoids can be easily fitted to the segmented protrusions, for instance by aligning the moments or principal axes, or position sticks along the main axis of the segmented body-part.

Figure 27 shows the resulting model fitted to a sequence of voxel-sets representing a hand and its fingers. In an augmented reality environment the 3D reconstruction of a hand could be used to interact with virtual objects in a virtual environment. From this ellipsoidbased representation one can easily infer an implicit-surface representation which can be used for real-time interactive applications. Such a parametric description of the object is complementary to data, e.g., voxel, representations. It can be noted by looking at Figure 27 that our methodology does not (by construction) necessarily identify all articulations. It does, however, identify remarkable protrusions in the moving body. As the hand evolves and more fingers become clearly separated from the palm, they are isolated as separated clusters and the rough model of the hand can be updated. Even though limb articulations do not show up in the embedding space, these can be later identified by refining the obtained segments in a more comprehensive model recovery framework.

Conclusions

In this paper we presented a novel dynamic, unsupervised spectral segmentation scheme in which moving articulated bodies are clustered in an embedding space, and clusters are propagated in time to ensure temporal consistency. By exploiting some desirable geometric characteristics of LLE, in particular, we can estimate the optimal number of clusters in order to merge/split them when topology transitions occur. We compared the performance of our algorithm versus (time-consistent) EM clustering in 3D, k-means clustering in ISOMAP space, and ground truth labeling provided through motion capture or synthetic generation.

An extension of the proposed approach to widely separated, non contiguous poses can be devised in a quite straightforward manner by resorting for cluster propagation to methods that match different poses of the same articulated object by aligning their embedded images [START_REF] Jain | Robust 3d shape correspondence in the spectral domain[END_REF][START_REF] Bronstein | Generalized multidimensional scaling: A framework for isometry-invariant partial surface matching[END_REF][START_REF] Mateus | Articulated shape matching using laplacian eigenfunctions and unsupervised point registration[END_REF]. Though the methodology has been proposed for articulated objects, we can also hope to extend it to certain classes of deformable objects. This could be done by exploiting the property of graph Laplacian methods which states that the eigenvalues of the graph Laplacian remain stable under deformations of the object to segment which preserve its volume [START_REF] Zhou | Large mesh deformation using the volumetric graph Laplacian[END_REF]. Much further work needs to be done in this sense, and will be the aim of our research in the near future.

Fig. 1

 1 Fig.1Sample results of the proposed approach to unsupervised spectral segmentation of moving bodies (represented as sets of stereo-reconstructed voxels) along entire sequences. Image sequences from multiple cameras (a) are used to generate a volumetric reconstruction of the moving body at each time instant (b); the final result is a consistent unsupervised segmentation of the moving body into segments closely related to the actual body's links, and their centroids' trajectories.

  Then the low-dimensional embedded vectors {Y i } associated with the original cloud of points {X i } are computed by minimizing the sum of the squared reconstruction errors (Figure 2): arg min

Fig. 2

 2 Fig. 2 Graphical description of the LLE algorithm.

Fig. 3 a

 3 Fig. 3 a) How LLE (middle) and ISOMAP [68] (right) map the same 3D cloud (left) for the same number of neighbors k = 13. LLE increases the separation of protrusions like legs and arms, unlike methods based on geodesic distances. b) The number of neighborhoods affected by articulated motion is relatively small..

  -b).

  Figure 4 illustrates two different poses of a dancer performing a ballet, and how the related embedded clouds obtained through LLE for d = 3 and k = 10 are stable (up to an orthogonal transformation).

Fig. 4

 4 Fig. 4 Stability of LLE under articulated motion. Two different poses (left and middle) of the same person are mapped onto the same embedded shape (right), for a large interval of parameter values and up to an orthogonal transformation.

Fig. 5

 5 Fig. 5 Left: rationale of 3-wise (k-lines) clustering. The areas of the triangles defined by triads of points measure their collinearity. As the points tend towards collinearity, the area of the triangle tends to 0. Right: affinity of k-tuples in a hypergraph H.

Fig. 6

 6 Fig. 6 Termination (left) and internal (right) points of the embedded cloud are characterized by their projection (red square) on the line (in blue) interpolating their neighborhoods (in green) being an extremum of the interval of all projections.

Fig. 7

 7 Fig. 7 Seed propagation for consistent clustering along time in the embedding space. The anti-images of the embedded centroids at time t are added to the 3D cloud X(t + 1) at time t + 1. Their own embeddings c ′ j (t + 1) are the seeds from which to start clustering the new embedded cloud Y (t + 1).

  , ..., N(0)} for t = 0 (Figure 8-a), is mapped to an embedding space of dimension d, yielding an embedded cloud Y (t) = LLE[k, d](X ′ (t)); 2. all the branch terminations of the embedded cloud Y ′ (t) are detected (Section 3.2): the number of clusters n(t) at time t is then set to the number of branches (plus one for the torso), Figure 8-b; 3. the embedded cloud Y ′ (t) is clustered into n(t) groups by d-wise clustering (Section 3.1) starting from n(t) seeds (Figure 8-c):

Fig. 8

 8 Fig. 8 Graphical illustration of the proposed coherent spectral segmentation algorithm.

Fig. 9 A

 9 Fig.9A few sample frames from the synthetic sequence "reveil" used in the first experiment. The automatically generated ground-truth segmentation into rigid links is shown in colors.

Fig. 10

 10 Fig. 10 Top: sample frames of the "dancer" sequence. Bottom: frames 20-26 of the real voxel-set sequence "dancer" used in the second experiment.

Fig. 13

 13 Fig.13Segmentation scores obtained by time-consistent EM clustering on the same synthetic sequences of Figure12. For sequence 5, the scores produced by two runs of EM are shown (as there is a random component to it): its performance is consistently not good.

Fig. 14

 14 Fig. 14 Segmentation results (left) and centroid trajectories (right) obtained by our algorithm (top) and time-EM clustering (bottom) for the sequences "walk" and "surf1". The voxelset segmentation for the last frame of the sequence is shown as an example.

Fig. 15

 15 Fig. 15 Examples of results produced by the dynamic algorithm of Section 3.5 on a number of real-world sequences: a) "fly" sequence of 11 frames; b) "arm-waving" sequence of length 50; c) "walking" motion, length 10; d) a subsequence of "danceuse", 16 frames; e) another "danceuse", length 10. Both the evolution of whole unsupervised clusters in 3D and their centroids' trajectories are shown. f) Example of topology transition management: after 11 frames in which arms are spread out (left) the left arm gets in contact with the torso: the algorithm adapts the number of clusters accordingly and proceeds to segment in a smooth way for other 7 frames (right, from a different viewpoint).

Fig. 16

 16 Fig. 16 Segmentation scores (left) and centroid trajectories (right) obtained by our algorithm for the realworld sequence "mars", from t = 389 to t = 415. Top: our algorithm. Bottom: time-consistent EM clustering. The actual voxelset segmentation corresponding to the critical frame t = 408 of the sequence is also shown (middle).

Fig. 17

 17 Fig. 17 Performance scores for two more real-world sequences of length 29 and 40 respectively, and for all the competing methods. Only the coarsening (blue) and the (second) segmentation (in black) scores are plotted. Solid lines: our method, time-consistent d-wise clustering in LLE space; dashed lines: time-consistent EM clustering in 3D. Dash and dotted lines: time-consistent k-means clustering in ISOMAP space. The performance scores of the proposed approach clearly dominate all the others'.

  Figure 16-top illustrates how, for adequate values of the parameters (in the plots d = 4, k = 25), the segmentation scores achieved by our algorithm are once again very high, with the method showing remarkable resilience to unreliable data capture (red line). This compares very favorably with the behavior of time-consistent EM clustering (bottom).

Fig. 20

 20 Fig.20Different eigenfunctions capture different aspects of the shape's geometry. Top: the first four eigenfunctions for the frame t = 209 of the "mars" sequence. Each eigenfunction corresponds to a 1 by N eigenvector of the affinity matrix: its value of the i-th point of the grid is the i-th entry of the eigenvector and is rendered here in color. Bottom-left: segmented embedding cloud and corresponding 3D segmentation for the same frame, obtained by selecting eigenvectors 1,2, and 3. Bottom-right: segmentation associated with eigenvectors 2,3, and 4.

Fig. 21

 21 Fig.21How the discrete sampling of a 3D shape influences the performance of the algorithm. Top: sequence "space surf" from t = 300 to t = 319. Bottom: "reveil" from t = 150 to t = 200. Left: consistency (red) and segmentation (black) scores throughout each sequence, for different values of the product k by s: 3 × 14 (solid), 4 × 10 (dashed), 5 × 8 (dash-dotted), 6 × 7 (dotted) for "space surf"; 5 × 18, 6 × 15 and 7 × 14 for "reveil". Center: example segmentation (t = 319 for s = 3, k = 14 for the top sequence; t = 218 for the bottom one). Right: corresponding centroid trajectories.

Fig. 22

 22 Fig. 22 Consistency (left), segmentation (middle), and average (right) scores obtained over two sample synthetic sequence (top and bottom) for different values of the parameters k = {14, 16, 18, 20, 22} (on the abscissa) and d = {3, 4, 5, 6, 7} (on the ordinate) of the algorithm. Top performances are achieved for a wide range of the parameters. Top: "walk", from t = 71 to t = 80. Bottom: "surf", from t = 1 to t = 25.

Fig. 23

 23 Fig.23How the splitting/merging algorithm deals with topology changes in the embedded space. At t = 16 the left arm of the dancer touches her scarf (left), and a single cluster covers body, arm, and scarf. Then (t = 18) the left arm becomes visible and a new cluster is assigned to it, while the dancer's feet get too close to each other to be distinguished (middle). Finally (t = 20) her legs pull again apart, generating a separate cluster for each of them (right).

Fig. 24

 24 Fig. 24 Measuring the relative performance of "local" methods based on graph Laplacians (represented in this paper by LLE), EM clustering, and "global" embeddings based on geodesic distances (represented by ISOMAP) for sequences affected by topology changes. Left column: "wake up" sequence, from t = 95 to t = 114. Right column: "clap" sequence. Top: consistency (blue) and segmentation (black) scores throughout the sequence: solid -our method; dashed -time-consistent EM clustering; dash-dotted -time-consistent kmeans clustering in ISOMAP space. Middle: some examples of how our algorithm copes with topology transitions are given in terms of 3D segmentations before, during and after the transition. Bottom: ISOMAP clusterings for the final frames of the two sequences.

Fig. 25

 25 Fig. 25 Behavior of ISOMAP along the sequence of topology transitions of Figure 23. The shape of the embedded cloud (we chose d = 2 for sake of readability) changes dramatically, gravely affecting the segmentation in the original 3D space.

Fig. 26

 26 Fig.26The segments' centroid trajectories produced by our unsupervised segmentation algorithm are consistent enough to be used as feature vectors for a parameter estimation algorithm generating a dynamical model (for instance, a hidden Markov model) which represent the input video sequence. Such a model can later be used for action classification.

Fig. 27

 27 Fig.27Fitting a rough model formed by a number of ellipsoids to the obtained clusters is straightforward and yields remarkable results. Top: a sequence of voxel-sets capturing a hand's motion in an augmented reality environment where the user interacts with virtual objects. Bottom: corresponding rough articulated model fitting.

These data are available on-line at http://4drepository.inrialpes.fr/public/datasets.