Simultaneous completion and spatiotemporal grouping of corrupted motion tracks

The Visual Computer https://doi.org/10.1007/s00371-021-02238-8 ORIGINAL ARTICLE Simultaneous completion and spatiotemporal grouping of corrupted motion tracks Antonio Agudo1·Vincent Lepetit2·Francesc Moreno-Noguer1 Accepted: 15 June 2021 © The Author(s) 2021 Abstract Given an unordered list of 2D or 3D point trajectories corrupted by noise and partial observations, in this paper we introduce a framework to simultaneously recover the incomplete motion tracks and group the points into spatially and temporally coherent clusters. This advances existing work, which only addresses partial problems and without considering a unified and unsupervised solution. We cast this problem as a matrix completion one, in which point tracks are arranged into a matrix with the missing entries set as zeros. In order to perform the double clustering, the measurement matrix is assumed to be drawn from a dual union of spatiotemporal subspaces. The bases and the dimensionality for these subspaces, the affinity matrices used to encode the temporal and spatial clusters to which each point belongs, and the non-visible tracks, are then jointly estimated via augmented Lagrange multipliers in polynomial time. A thorough evaluation on incomplete motion tracks for multiple-object typologies shows that the accuracy of the matrix we recover compares favorably to that obtained with existing low-rank matrix completion methods, specially under noisy measurements. In addition, besides recovering the incomplete tracks, the point trajectories are directly grouped into different object instances, and a number of semantically meaningful temporal primitive actions are automatically discovered. Keywords Point track completion ·Spatiotemporal clustering ·Augmented Lagrangian multiplier 1 Introduction Motion visual tracking is an important and essential component of perception that has been an active research area in computer vision for past two decades. The developments of 2D and 3D visual tracking algorithms have shown rapid progress thanks to the explosive growth of video data which in turn creates high demand for accuracy and speed of tracking methods. Current approaches are motivated to design faster and better methods in spite of the challenges that exist in this topic, especially robustness to large occlusions, drasBAntonio Agudo [email protected] Vincent Lepetit [email protected] Francesc Moreno-Noguer [email protected] 1Institut de Robòtica i Informàtica Industrial, CSIC-UPC, 08028 Barcelona, Spain 2Université de Bordeaux, cours de la Libération, 33405 Talance, France tic scale change, accurate localization, multi-object tracking, and recovery from failure [25,27]. Despite the success in addressing numerous challenges under a wide range of scenarios, a number of core problems still remain unsolved. A major challenge in real scenarios is handling missing entries of the data, due to ad hoc data collection, presence of outliers, sensor failure, or partial knowledge of relationships in a dataset. For instance, to recover object motions and deformations from video, the tracking algorithm may lose the track of features in some image frames due to lack of visibility or mismatches. In a similar manner, for 3D tracking, multi-camera systems (such as motion capture systems) [26,46] are applied to obtain the time-varying evolution of a scenario. While these systems are now capable of recovering most of the observations, they can fail on real-world scenarios, such as those formed by multiple objects while performing different activities, deforming, moving, and even interacting between them. In these cases, missing tracks continually appear, as either self-occlusions or occlusions between objects. It is worth mentioning that this is especially relevant in outdoors scenarios, where current algorithms to estimate motion tracks 123 A. Agudo et al. often produce partial solutions with a wide amount of missing entries. In many fields, an underlying tenet is that the data may contain certain type of structure that enables intelligent processing and representation, and they can be characterized by using parametric models. Assuming that, the visual tracking completion problem can be addressed as a matrix completion one. To this end, one can use the well-known linear subspaces since they are easy to estimate, and often effective in many real-world applications. For instance, these models have been successfully used to characterize several types of visual data, such as motion [39,49], shape [32], and texture [35]. Maybe, the most common choice it is to use the principal component analysis method that is based on the hypothesis that the data are approximately drawn from a low-rank subspace. Unfortunately, real data from complex scenarios can rarely be well described by a single low-rank subspace. For these cases, a more reasonable model is to assume data are lying near several subspaces, i.e., the data are considered as samples approximately drawn from a union of several low-rank subspaces. The generality and importance of subspaces naturally lead to a challenging problem of subspace clustering, whose goal is to group data into clusters with every cluster corresponding to a different subspace. Solving clustering and finding low-dimensional representations of data are important unsupervised learning problems in machine learning with numerous applications, including image segmentation, system identification, data visualization, and collaborative filtering to name just a few. The problem becomes even more complex if the data are partially observed, due to either sensor failure or visual occlusions [1,2,4,24]. In this paper, we propose, to the best of our knowledge, the first attempt to approximate high-dimensional data using a dual union of low-dimensional subspaces, accounting for two distinct criteria. Additionally, input data are assumed to be corrupted by partial observations and noise. We apply our approach to the specific case in which the input data encode 2D or 3D point trajectories of multiple dynamic objects with large percentages of missing entries, and we aim at hallucinating these missing tracks while simultaneously approximating the data using spatial and temporal subspaces, as well as filtering the noisy measurements. We will formulate the problem as a matrix completion one. Input data will be arranged into a matrix with the missing entries set to zero. To encode data similarities, we introduce two affinity matrices to be learned. We will then devise an optimization scheme based on augmented Lagrangian multipliers (ALM) to simultaneously and efficiently estimate the missing entries, and the bases and dimensionality for of each low-rank subspace. The proposed algorithm is unsupervised, does not depend on the initialization nor relies on training data at all, and can be solved in polynomial time. An important corollary of our approach is that applying off-the-shelf state-of-the-art spectral clustering on the estimated affinity matrices results in consistent temporal and spatial segmentations of the input data. We evaluate the proposed algorithm on 2D and 3D incomplete motion capture and real sequences of several objects performing complex actions and interacting with each other. We will show that the accuracy of the completed tracks we obtain improves that of state-of-the-art methods by a considerable margin, while we additionally provide a spatiotemporal clustering of the data, which in most cases has a direct physical interpretation (either the object identity or the type of motion it is performing). 2 Related work The most standard approach to perform matrix completion is to assume the underlying data lies in a single lowdimensional subspace. Early works [30,44] enforced this constraint based on expectation maximization strategies to optimize non-convex functions of the model parameters and the missing entries. Other attempts constrain the solution space using trajectory [24], or spatiotemporal models [5]. Nevertheless, all these methods require a good initialization, and most importantly, they need to set the rank of the subspace a priori, performing poorly when the dimension of the subspace increases. Additionally, trajectory-based methods normally use a predefined basis, making them very problem specific. To address these limitations, another family of low-rank matrix completion techniques has been recently proposed [7,9–11,14,15]. These methods estimate missing entries by optimizing the convex surrogate of the rank, i.e., by they enforce the nuclear norm of the complete matrix. These ideas were also applied in problems where the matrix directly includes visual tracking information, imposingsmooth[42,43,54] and sparse [47] representations. When the underlying subspace is not consistent with standard basis components and missing track locations are spread uniformly at random, these approaches are guaranteed to recover missing entries. Unfortunately, matrix completion techniques based on a single low-rank subspace cannot handle the challenging and more general scenario in which input data lie in a union of low-rank subspaces (e.g., when dealing with simultaneous and incomplete tracks of multiple objects performing complex motions). Data segmentation from full annotations was proposed by assuming a union of subspaces by means of a subspace clustering based on sparse representation [18] or seeking the lowest rank one [34]. Going back to the completion problem, the objective would extend to recovering the missing entries together with the clustering of the data according to the subspaces. Mixture of factor analyzers [22], mixture of probabilistic principal component 123 Simultaneous completion and spatiotemporal grouping of corrupted motion tracks Table 1 Qualitative comparison of our approach against competing techniques Approaches Missing Automatic Temporal Spatial Unified entries rank clustering clustering framework Category #1  Category #2  Category #3   Category #4    Ours     We split the methods into four incremental categories, depending on a number of desirable properties, namely robustness to missing entries, automatic estimation of rank, computation of temporal/spatial clusters, and unified formulation. Note that the method proposed in this paper is the only that can simultaneously offers all these properties. Some examples for every category are: #1 methods to solve completion using a single low-rank formulation [7,9–11,14,15,43], #2 methods to solve completion by means of multiple subspaces where the rank is known a priori [6,45], #3 methods able to automatically retrieve the matrix rank, such as [36,40,50], and finally #4, methods that can solve the problem in a unified manner [17,20]. It is worth noting that [3] used a spatiotemporal constraint for shape reconstruction, but not for shape completion as we do in this paper. Self-expressive models [18,34] solved for one type of clustering but considering full data analysis [45], and incremental matrix completion algorithms with K-subspaces [6] are some early methods used to address grouping and completion of multi-subspace data. Again, the performance of these methods highly depends on the initialization and degrades for large subspace ranks. A polynomial number of data points in the ambient space dimension are required in [19] which often cannot be met in highdimensional datasets. Ma et al. [36] proposed an algebraic approach to model data drawn from a union of subspaces based on generalized principal component analysis. Yet, due to the difficulty of estimating the polynomials from data, the method is sensitive to noise and is computationally very demanding. This strategy was extended in [40], yielding more robust solutions but onlyfor low-dimensional input data and a reduced number of subspaces. A Lipschitz monotonic function was assumed to model the low-rank matrix in [21], even though this cannot cover the case of multiple subspaces. Another family of solutions proposed solving completion and clustering as a two-stage problem [50], by first obtaining a similarity graph for clustering and then applying low-rank matrix completion to each cluster. While this is an interesting direction, the solution proposed in [50]isproneto fail when subspaces intersect or when the initial grouping is incorrect. To solve this limitation, Elhamifar [17] has proposed self-expressive models for simultaneous clustering and completion of incomplete data. Along the same line, Fan and Chow [20] have recently presented a sparse representation to solve the problem. However, these approaches can only cluster the data based on one single criterion. In parallel, some works have relied on neural networks to learn temporal clustering [52] and infer missing entries [37,55], but solving just a single problem. In all cases, these approaches propose to exploit a loss function as we do in this paper, but they require a large amount of training data to learn the model and demand a specific hardware to complete the training step. Unfortunately, this cannot be assumed for generic scenarios, where an unknown number of unknown object typologies can deform, move, and even interact between them, doing the process of simultaneously obtaining training data for track completion, spatial groups, and temporal ones very hard and expensive in practice. In contrast, our formulation can solve the problem in just few seconds in a commodity computer, without requiring sophisticated hardware, nor prior knowledge about the scenario to be solved. Moreover, none of them simultaneously solve multiple clustering and completion as we propose in this paper. 2.1 Our contributions We go beyond previous works by proposing an efficient and robust method that does not require initialization, and it can jointly perform two types of clustering (spatial and temporal), while recovering missing entries and filtering the rest. To the best of our knowledge, no previous approach has jointly addressed the three problems in a unified and unsupervised framework. To this end, we assume the input data to lie in a dual union of low-rank subspaces, where no a priori knowledge about the dimensionality of the subspaces or which data points belong to which subspace is required. It is worth noting that our approach does not require any training data at all. Additionally, the proposed solution can handle situations with complex motion patterns, affected by large degrees of overlapping and percentage of missing entries, in a completely unsupervised manner. Table 1summarizes a qualitative comparison of our approach and the aforementioned techniques to jointly solve completion and clustering. 123 A. Agudo et al. 3 Preliminaries and problem statement 3.1 Notation Matrices are represented with boldface uppercase letters, e.g., X. In particular, IAis used to denote the identity matrix of size A×A, and 1Aa column vector of ones of size A×1. The entries of matrices are denoted by means of subscripts [·]. For instance, X[: j]corresponds to the jth column of the matrix X,X[i:] is the ith row of the matrix X, and X[ij]indicates its (i,j)th entry. We also define two types of products: X⊗Zto denote the Kronecker product, and XZto denote the Hadamard (or element-wise) one. The negative of a binary matrix Xis denoted as ¯ X. We also define several norms on matrices: The l∞-norm is defined as X∞=max(i,j)|X[ij]|, and the l2,1-norm as X2,1=ΣjX[: j]2, where the l2-norm of a vector is denoted by X[: j]2. The Frobenius and nuclear norms are represented as XFand X∗=Σiσi(X), respectively, with σi(X)being the ith singular value of the matrix X. Finally, the Euclidean inner product between two matrices is denoted as X,Z=tr(XZ), where tr(·)represents the trace of a matrix. 3.2 Problem Formulation Let us consider Ftemporal subspaces {Sf}F f=1of dimension {df>0}F f=1in a C-dimensional space, and Gspatial subspaces {Sg}G g=1of dimension {dg>0}G g=1in a Hdimensional space. Let Y∈RC×TbeamatrixofTdata points lying on the union of the temporal subspaces, and ˆ Y∈RH×NamatrixofNdata points lying on the union of spatial subspaces. If we assume that both dimensions can be factorized by a factor D(i.e., C=DN and H=DT), the two matrices Yand ˆ Ycan then contain exactly the same number of values but in a different arrangement. Additionally, we will assume that only some entries of these matrices are observed, i.e., some locations can include null values. To denote this, we include the matrix ˜ Y, a sparse version of Y, in which non-observed entries are set to zero. Our problem consists in, given an incomplete and noisy matrix ˜ Yof data points from motion tracking, retrieving the full matrix, Yor ˆ Y, and clustering the data into the underlying temporal and spatial subspaces. To this end, we will encode the spatial and temporal subspaces using affinity matrices. It is worth noting that both the bases ({Sf}F f=1and {Sg}G g=1) and dimensions of each subspace (dfand dg, respectively) are not known a priori, nor to which cluster each data point belongs to. The incomplete and noisy input matrix can be provided by any tracking algorithm, by considering, for instance, optimization [25,27] or deep learning approaches (Joo et al. [28]). We next describe our unsupervised and unified approach that can solve the problem without requiring any training data at all. 4 Spatiotemporal subspace clustering Drawing inspiration on the ideas of [3] for reconstructing non-rigid shapes, we next generalize a spatiotemporal constraint for joint motion track matrix completion and clustering. Note that this constraint was not used previously in the literature for completing missing entries as we present here. We first introduce the two types of interpretations of the tracking matrices we shall use. After that, and considering the previous interpretations, we will introduce the temporal and spatial constraints, extending our formulation to handle missing tracks. 4.1 Motion tracking matrix interpretations Let us consider a dynamic set of ND-dimensional points tracked along Ttime instances. For the particular case of D=3, i.e., a tridimensional space, we shall denote by xt i= [xt i,yt i,zt i]the spatial coordinates of the ith point at time instant t. All acquired point coordinates can be collected into the matrix Y∈RDN×Tin an unordered manner in terms of any type of grouping that stores the x,y, and zcoordinates in a block matrix form as: Y=⎡ ⎢ ⎣ x1 1... x1 Ny1 1... y1 Nz1 1... z1 N . . ..... . .. . ..... . .. . ..... . . xT 1... xT NyT 1... yT NzT 1... zT N ⎤ ⎥ ⎦  . We could assume the previous motion tracking matrix admits a low-rank decomposition of rank K(K=1for rigid objects), where Krepresents the number of bases in a single subspace. We know from the structure from motion theory this matrix is of low rank (Dai et al. [16], Xiao et al. [48]), but since no information about the motion is assumed, only a low-rank constraint can be considered. However, as discussed above, the single low-rank assumption may not have sufficient expressiveness power to model complex motion patterns of multi-object tracks. It is worth mentioning that if we know some kind of clustering or grouping of the Tdata points, we might handle this situation by enforcing the low-rank assumption to every particular cluster. In this work, however, the number and type of clusters is not known a priori, making the problem more challenging and generic. Consequently, we need to jointly solve for completion and clustering, without assuming any information about the dimensionality of the subspaces. Since each column of the matrix Yencodes all points at a time instant, this matrix cannot be directly used to retrieve 123 Simultaneous completion and spatiotemporal grouping of corrupted motion tracks spatial similarities. To address this limitation, we consider a new DT ×Nmatrix ˆ Y, for which each column stores the point tracks. Following the previous case of D=3, this matrix can be written as1: ˆ Y= ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ x1 1x1 2... x1 N y1 1y1 2... y1 N z1 1z1 2... z1 N . . .. . ..... . . xT 1xT 2... xT N yT 1yT 2... yT N zT 1zT 2... zT N ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ that is also low rank as Ybut differing in value. Both matrices use two different matrix arrangements of the data points, but they include exactly the same information. We can map from Ymatrix to ˆ Yusing the relation: ˆ Y=(ID⊗Y)A,(1) where Ais a (DDN)×Nbinary matrix. The inverse mapping can be written as: Y=(ˆ Y⊗ID)B,(2) where Bis also a (DDT)×Tbinary matrix. Both Aand B matrices are known a priori, and they can be easily obtained by considering the data structure in data matrices ˆ Yand Y. 4.2 Dual union of spatiotemporal subspaces The arrangement of the point tracks through the matrices Y or ˆ Ygives two different interpretations, and each of it can be associated with a distinct subspace clustering process. For instance, when analyzing the temporal domain using Y,we can define an affinity matrix to capture the temporal similarities between instances at different time steps. This relation can be written as: Y=YT +Et,(3) where Tencodes a temporal affinity T×Tmatrix and Etis a DN ×Tresidual noise. In this context, the temporal affinity matrix Tmeasures the similarities between D-dimensional poses along time. Using this relation, we enforce Yto lie in a union of Sftemporal subspaces, each of them with rank df. We could say that the matrix Tis the lowest rank representation of the data Ywith respect to itself. It is worth 1Note that other dimensions could be used, such as 2D motion tracking where D=2. In this case, we only need to eliminate the rows of the matrices Yand ˆ Ycorresponding to the zcomponent. Theoretically, other dimensions can be also used, such as 1D image sensors. noting that Twill be block-diagonal when the data samples have been grouped together in Yaccording to the subspace memberships. This block pattern is lost for random entries, obtaining null entries when no affinities are provided. This type of self-expressive model was previously used by [18,34] in the context of subspace clustering. Similarly, we can analyze the spatial domain through the matrix ˆ Y, by introducing an affinity matrix associated with a union of spatial subspaces in the presence of noise. In this case, we can write: ˆ Y=ˆ YS +Es,(4) where Sencodes a spatial affinity N×Nmatrix and Esis a DT ×Nresidual noise. In this case, we are enforcing ˆ Yto lie in a union of Sgspatial subspaces of rank dg, respectively, measuring the similarities between D-dimensional points in a same time instant. Basically, Sand Tare made of low-rank coefficients that define the union of subspaces in every domain, respectively. Once these affinity matrices are learned from data, off-the-shelf spectral clustering algorithms like [13] can be applied on each of them to discover the grouping in every domain. The temporal clustering splits the data into motion primitives and the spatial one into different object instances. Nevertheless, the previous formulation requires full measurements on the tracking matrices Yor ˆ Y, which is not often the case in real applications. Previous subspace clustering algorithms assume the observation matrices to be complete [3,18,29,34,56]. As mentioned above, other approaches [17,20] proposed an algorithm to jointly estimate missing entries and build a similarity graph for clustering, when considering a single union of temporal subspaces. The algorithm we present in the following section goes beyond these approaches and allows solving the matrix completion problem when considering the data to be spanned by two different union of subspaces. Our approach can handle high levels of missing entries and noisy measurements, and solve the problem by means of a one-stage optimization algorithm. This means our approach can produce more accurate solutions than competing techniques, while being more general. 5 Motion tracking completion and spatiotemporal subspace clustering We next present our algorithm to simultaneously recover missing entries and estimate two similarity matrices for computing the spatial and temporal grouping. Note that neither prior information nor training data is used at all. The input to our algorithm are incomplete motion tracks of NDdimensional points observed along Ttime instances that are arranged into the matrix ˜ Y. In addition, we also introduce 123 A. Agudo et al. an observation matrix O∈RN×Twith binary entries that indicate whether the coordinates of a point at a specific time instant are observed or not. 5.1 Proposed formulation Let us denote by Θ≡{Y,ˆ Y,T,S,Et,Es}the set of model parameters we have to learn from the input data Γ≡{˜ Y,O}. We introduce an optimization framework ruled by a cost function that accounts for the spatiotemporal clustering constraints of Eqs. (3) and (4), and enforces the similarity matrices Tand Sto be spanned by low-rank subspaces. Consequently, the combination of both constraints enforces the data in order to lie in a dual union of subspaces. Indeed, the single union of subspaces model can be seen as a degenerate case of our model (see Remark 1). Since rank minimization is a non-convex NP-hard problem [41], the nuclear norm is approximated by its convex relaxation [12,14]. Additionally, in order to be able to deal with data corrupted by noise and outliers, we use l2,1-norm regularization, as the convex relaxation of the l2,0-norm [33]. The objective function can therefore be written as: arg min Θ (1D⊗O)˜ Y−Y2 F+φ(T∗+S∗) +γY∗+λtEt2,1+λsEs2,1(5) subject to Y=YT +Et ˆ Y=ˆ YS +Es (ID⊗Y)A=ˆ Y where {φ,γ,λt,λ s}are predefined penalty term parameters. Remark 1 When the data points are not connected in the spatial domain, it means that the affinity matrix Sbecomes the identity IN(we assume the data points are clean in this domain, i.e., Es=0), and hence, our formulation degenerates to a union of temporal subspaces. On the other hand, when this occurs in the temporal domain (T=ITand Et=0), our formulation degenerates to a union of spatial subspaces. 5.2 Efficient augmented Lagrangian multiplier optimization The optimization problem in Eq. (5) can be efficiently solved in a unified manner via an ALM method [8,31]. Without loss of generality, we set λ≡λt≡λs. In order to reduce the number of parameters and the complexity of the problem while improving convergence, we choose to bring the clustering constraints into the energy function using several Lagrange multipliers with a unique penalty weight β>0. In addition, we introduce three support matrices Y≡M,T≡J, and S≡K, to obtain the corresponding augmented Lagrangian function that can be written as: arg min ΘL (1D⊗O)˜ Y−Y2 F+φ(J∗+K∗) +γM∗+λEt2,1+Es2,1 +L1,Y−YT−Et+β 2Y−YT−Et2 F +L2,ˆ Y−ˆ YS−Es+β 2ˆ Y−ˆ YS−Es2 F +L3,(ID⊗Y)A−ˆ Y+β 2(ID⊗Y)A−ˆ Y2 F +L4,Y−M+β 2Y−M2 F +L5,T−J+β 2T−J2 F +L6,S−K+β 2S−K2 F(6) where ΘL≡{M,Y,J,T,K,S,ˆ Y,Es,Et}includes the tracking completion, spatiotemporal similarity parameters, and residual noises. The Lagrange multipliers are defined as {L1,L4}∈RDN×T,{L2,L3}∈RDT×N,L5∈RT×T, and L6∈RN×N. Recall that we do not need to know the dimensions nor the bases of the temporal and spatial subspaces a priori, since Eq. (5) automatically selects the appropriate number of data points from every spatiotemporal subspace. We propose to solve the problem in Eq. (6) by minimizing each variable individually and in closed form, while keeping fixed the rest of model parameters. Algorithm 1explains the details. The expressions for estimating Y,T,S, and ˆ Y(steps 4, 6, 8, and 9) are obtained by computing the derivatives of Eq. (6)inY,T,S, and ˆ Y, respectively, and equating to zero. The subproblems to recover M,J,K,Et, and Esare convex and have closed-form solutions. Particularly, for steps 2, 5, and 7, we apply a singular value thresholding minimization [10] with a ‘shrinkage operator’ S∗ β(x)=max(0,x−∗ β) where ∗={φ,γ}. In order to optimize the noise terms Et and Es(steps 10 and 11, respectively), we apply the Lemma 4.1 in [51]. After each iteration, the Lagrange multipliers are updated according to standard rules as shown in lines 1213. Additionally, we also update the penalty weight β(step 14) to guarantee the convergence of our algorithm, following the upper bounded requirement of the alternating direction methods. Particularly, we apply a factor of 1.1 to increase β every iteration. The theoretical convergence of our algorithms is not easy to proof, as the method is based on nine different blocks. However, we have empirically observed that for all experiments reported in the following section, the algorithm always converged in about 190 −220 iterations. Additionally, we observe the optimality gap obtained in every iteration to monotonically decrease. An example of this analysis is dis123 Simultaneous completion and spatiotemporal grouping of corrupted motion tracks Algorithm 1: Algorithm for optimizing Eq. (6). Input : Incomplete trajectories ˜ Yon a D-dimensional space. Parameters {φ,γ,λ},andβ Output: Matrix completion Y(or ˆ Y), temporal Tand spatial Saffinity matrices for clustering 1while not converged do /* Update Model Parameters */ 2M=min γ βM∗+1 2M−(Y+L4 β)2 F 3D=M−L4 β+(ˆ Y−L3 β)⊗IDB+(Et− L1 β)(IT−T)(IT−T)(IT−T)+2IT−1 4Y=(1D⊗O)1 2+β(2˜ Y+βD)+1D⊗¯ OD 5J=min φ βJ∗+1 2J−(T+L5 β)2 F 6T=(YY+IT)−1(Y(Y−Et)+J+YL1 −L5 β) 7K=min φ βK∗+1 2K−(S+L6 β)2 F 8S=(ˆ Yˆ Y+IN)−1(ˆ Y(ˆ Y−Es)+K+ˆ YL2 −L6 β) 9ˆ Y=(Es−L2 β)(IN−S)+L3 β+((ID⊗ Y)A(IN−S)(IN−S)+IN−1 10 Et=min λ βEt2,1+1 2Et−(Y−YT +L1 β)2 F 11 Es=min λ βEs2,1+1 2Es−(ˆ Y−ˆ YS +L2 β)2 F /* Update Lagrange Multipliers */ 12 L1=L1+β(Y−YT −Et); 13 L2=L2+β(ˆ Y−ˆ YS −Es) 14 L3=L3+β(ID⊗Y)A−ˆ Y 15 L4=L4+β(Y−M) 16 L5=L5+β(T−J) 17 L6=L6+β(S−K) /* Update Penalty Weight */ 18 β=min(1.1·β, βmax) /* Check Convergence */ 19 Y−YT −Et∞< 20 ˆ Y−ˆ YS −Es∞< 21 (ID⊗Y)A−ˆ Y∞< 22 Y−M∞< 23 T−J∞< 24 S−K∞< 25 Setting: In our experiments, we use =10−8,and βmax =1012 since the values in Yare normalized within the range [-1,1]. played in Fig. 1, where both constraints and full errors in Eq. (6) are represented for a specific case. As it can be seen, after around 50 iterations all constraints are almost perfectly satisfied and the overall energy converges. 5.3 Complexity analysis The most computationally demanding parts of Algorithm 1 are the steps 2, 5 and 7, which require computing several SVD operations over matrices of size DN ×T,T×Tand N×N, respectively. Hence, our problem can be solved in polynomial time with a computational complexity of at Fig. 1 Convergence analysis: energy reduction as a function of the number of iterations. Evolution of the error for the six constraints (denoted as Cc, with c={1,...,6}) and the full energy in Eq. (6) as a function of the number of iterations until convergence (corresponding to the Jump scenario described in “Experimental results” section). Note that two different scales are used to represent the errors of the constraints (left axis) and the full error (right axis). For visualization purposes, we plot the full energy scaled by a factor of 0.1 Fig. 2 Computation time as a function of the number of frames, points and iterations. Computation time versus number of iterations until convergence on the mocap sequences described in “Experimental results” section, for two (red dots) and four (blue dots) people. Next to each dot are indicated the number of images of the sequence. In all cases, the number of iterations until convergence always remains within reasonable bounds. The corresponding computation time depends on the number of frames and points most of O(N2T+T3+N3)[23]. Note that this complexity could be easily reduced by orthogonalizing the columns of the matrices Yand ˆ Y. The computation times (in unoptimized MATLAB code) on a commodity laptop with an Intel Core i7 processor at 2.4GHz for motion capture sequences for two and four people are displayed in Fig. 2. On average, the median computation time in experiments with sequences between 277 and 652 frames, and two people (N=82 123 A. Agudo et al. points) was of 51 seconds. Processing between 214 and 432 frames, and four people (N=164 points) required a median time of 44 seconds. In any case, to handle larger datasets we could use current results Yao et al. [53] on the use of SVD operations on large datasets to address large-scale low-rank problems. This could really help to reduce the reported complexity. Moreover, our formulation could be extended to be employed in a sequential manner, being this a part of our future work. 6 Experimental results In this section we report the performance of our algorithm to solve motion tracking completion, as well as temporal and spatial clustering on several challenging datasets. For all cases, we denote by ρthe fraction of missing entries in the input data. In all experiments, we set φ=1.0, γ=2.0 and λ=0.03. It is worth pointing out that we do not need fine-tuning these parameters, as the results were stable for wide range of values for φ∈[0.1,10] and γ∈[0.2,20]. Regarding the competing approaches, we will compare our algorithm, denoted as spatiotemporal track completion (ST2C), with the low-rank matrix completion (LRMC) [10], and the bilinear factorization matrix completion (BFMC) [9], two approaches where the rank is automatically estimated. We do not include [6,22]asthese methods require knowing the rank of every subspace a priori. Unfortunately, neither can we report the results of [17,20]as its source code is not publicly available. Recall, however, that both approaches did only consider a single union of subspaces. To establish a quantitative evaluation, we will compute three types of errors: the temporal eTC and spatial eSC clustering error as well as the motion tracking completion eMTC (this error is equivalent to a matrix completion evaluation) that are defined as: eTC =#Misclassified frames #All frames ,(7) eSC =#Misclassified points #All points ,(8) eMTC =Y−YGT 2 F YGT2 F .(9) where YGT and Ydenote the ground truth and the recovered matrices, respectively. For the temporal clustering error, we have obtained the ground truth segmentation over noisefree and complete measurement matrices by applying [34] to compute the similarity matrices and [13] to obtain the clusters. Spatial ground truth were annotated by hand. This means the evaluation we propose for temporal clustering is Fig. 3 Patterns of missing entries. Vpatterns used to simulate missing entries in the Jump sequence. White and black cells denote non-visible and visible points, respectively. Top: ρ=0.4 of random missing entries. Bottom: ρ=0.4 of structured missing entries actually an implicit comparison with respect to the competing approach [34] by assuming clear measurements. 6.1 Real experiments on motion capture data We evaluate the proposed approach on the CMU MoCap dataset. We consider several scenarios with either two or four people interacting and performing complex motions in 3D. On average, the sequences we consider are 433 frames long, and the number of points per frame is either 82 (two people) or 164 (four people). Specifically, we select eight sequences with two people: 23_16 (Jump): subjects alternating synchronized jumping jacks; 19_05 (Pull): a subject pulls the other by the elbow; 22_20 (Violence): a subject picks up high stool and threatens to strike the other; 20_06 (Soldiers): subjects follow a soldiers march; 23_19 (Stares Down):a subject stares down the other and leans with hands on high stool; 22_12 (Stumbles): a person stumbles into the other; 20_09 (Nursery): people follow a nursery rhyme; and 22_10 (Shelters): a person shelters the other from harm. A total of four sequences with four people are considered, synthetically generated by combining pairs of sequences with two people. All sequences are corrupted in three different ways: 1) randomly removing a fraction ρ={0.1,...,0.8}of entries of the measurement matrix Y; 2) removing a structured fraction ρ={0.1,...,0.4}of entries of the measurement matrix Ywhere we emulate temporal self-occlusions or lack of visibility, by including patterns with 50% of structured missing entries per frame; and 3) adding noise to the observed points, according to a Gaussian distribution with standard deviation σnoise =τ 100 ψ, where τcontrols the amount of noise, and ψ represents the maximum distance of a point to the centroid of all the points. An example of these artifacts is shown in Fig. 3, for both random and structured missing entries. Figures 4and 5summarize the results for two and four people, respectively. Each graph depicts the results of all 3 methods for one specific sequence, at increasing levels or 123 Simultaneous completion and spatiotemporal grouping of corrupted motion tracks Fig. 4 Motion completion errors of different algorithms as a function of the missing entries rate ρon motion capture sequences with two subjects. Each algorithm is evaluated under noise-free (τ=0) and noisy (τ={1,2}) data. For visualization purposes, the error of LRMC has been divided by a factor 3.5 in all graphs. Top: Random missing entries. Bottom: Structured missing entries missing data for the two types of cases we propose. Solid and dashed lines represent results for noise-free and noisy measurements, respectively. Our approach and BFMC [9] show similar error patterns, even though ours being always consistently better. A breaking point is achieved earlier by BFMC [9], showing our superiority in terms of robustness against this type of artifacts. As it can be seen, our solution by assuming noise can even provide better solutions than the competing approaches for clean annotations. The performance of LRMC [10] is far below the other two algorithms. We hypothesize this is due to the pseudo-block structure of the missing data, as each missing point does indeed represent three—recall that for this experiment, D=3—adjacent null elements in Y. This is especially relevant when the missing entries are structured, as it can be seen in the bottom part of Figs. 4and 5. Note also that BFMC [9] and LRMC [10] are specifically designed for matrix completion. These algorithms do not provide any kind of affinity measure, that allows subsequent clustering. Some instances for several scenarios when the missing entries are random are displayed in Fig. 6. Moreover, our algorithm is faster than the competing approaches, producing an speed up of 2.7×when BFMC [9] is considered. As we have commented above, our approach also estimates spatial and temporal clustering. Tables 2and 3 summarize the mean error for each sequence and all levels of missing data for the random and structured cases, respectively, for noiseless (τ=0) and noisy (τ={1,2}) measurements. As it can be seen, our approach produces very good results for most of the sequences, especially in terms of spatial clustering where we obtain an almost negligible clustering error. In fact, our algorithm produces better spa123 A. Agudo et al. 56. Zhu, Y., Huang, D., de la Torre, F., Lucey, S.: Complex non-rigid motion 3D reconstruction by union of subspaces. In: Computer Vision and Pattern Recognition, pp. 1542–1549 (2014) Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Antonio Agudo received the M.Sc. degree in industrial engineering and electronics in 2010, M.Sc. degree in computer science in 2011, and the Ph.D. degree in computer vision and robotics in 2015, from University of Zaragoza. He was a visiting student at vision group of Queen Mary University of London in 2013 and with the vision and imaging science group of University College London in 2014. He was also a visiting fellow at Harvard University in 2015, and at Universitè de Bordeaux in 2018 and 2019. After two years as a postdoctoral fellow at Institut de Robòtica i Informàtica industrial, CSIC-UPC, in Barcelona, he joined as an associate researcher of the Consejo Superior de Investigaciones Científicas (CSIC) in 2017. His research interests include non-rigid structure from motion, machine learning, and deformation analysis to medical and robotics applications. He is a recipient of the 2018 Best Paper Award Honorable Mention from the European Conference on Computer Vision. Vincent Lepetit is a Full Professor at the LaBRI, University of Bordeaux, and a senior member of the Institut Universitaire de France (2018-2025). He also supervises a research group in Computer Vision for Augmented Reality at the Institute for Computer Graphics and Vision, TU Graz. His research interests include vision-based Augmented Reality, 3D camera tracking, machine learning, object recognition, and 3D reconstruction. He is an editor for the IEEE Transactions on Pattern Analysis and Machine Intelligence, the International Journal of Computer Vision, and the Computer Vision and Image Understanding journal. Francesc Moreno-Noguer received the MSc degrees in industrial engineering and electronics from the Technical University of Catalonia (UPC) and the Universitat de Barcelona in 2001 and 2002, respectively, and the PhD degree from UPC in 2005. From 2006 to 2008, he was a postdoctoral fellow at the computer vision departments of Columbia University and the Ecole Polytechnique Fédérale de Lausanne. In 2009, he joined the Institut de Robòtica i Informàtica industrial in Barcelona as an associate researcher of the Spanish Scientific Research Council. His research interests include retrieving rigid and non-rigid shape, motion, and camera pose from single images and video sequences. He received UPC’s Doctoral Dissertation Extraordinary Award for his work. 123