Reconstructing a dynamic world: What is next?

Title: Reconstructing a dynamic world: What is next? Author: Fernando Gastón Codony Advisor: Antonio Agudo Martínez Department: Institut de Robòtica i Informàtica industrial Academic year: 2022/2023 Master of Science in Advanced Mathematics and Mathematical Engineering Abstract In this thesis, we study the problem denoted as Non-Rigid Structure from Motion and tackle the two distinct parts of the problem: motion and shape estimation together with the corresponding temporal segmentation into actions of the body. For the motion estimation, we implement a Single Rotation algorithm (SRA) based on the Weiszfeld algorithm for the median of points in Rn. We use Brand’s method to obtain a full corrective matrix and several estimates for the rotation matrices but find that the first column triplet obtains better results than any of the triplets found by Brand’s method, making SRA pointless. For the shape factor, we make an assumption that shapes lie in a temporal union of subspaces. We perform Sparse Subspace Clustering to jointly reconstruct the shape matrix while computing an affinity matrix that can be used to cluster frames of the input data according to which subspace they belong to. Keywords Non-Rigid Structure from Motion, Geometric Computer Vision, Single Rotation Averaging, Sparse Subspace Clustering. 1 Contents 1 Introduction 3 1.1 Rigid Structure from Motion ............................. 3 1.2 Non-Rigid Structure from Motion .......................... 6 1.3 About the under-constraint nature of NRSfM .................... 7 2 Single Rotation Averaging 9 2.1 Representations of rotation ............................. 9 2.2 SOp3qand sop3q................................... 10 2.3 Distances in SOp3q................................. 12 2.4 l2averaging ..................................... 12 2.5 l1averaging ..................................... 13 2.6 SRA in NRSfM .................................... 15 3 Sparse subspace clustering 19 3.1 SSC with ADMM .................................. 20 3.2 Joint shape reconstruction and SSC ......................... 21 4 Conclusion 27 2 1. Introduction In this thesis, we study a classical problem in geometric computer vision: Non-Rigid Structure from Motion (NRSfM). The task consists in reconstructing 3D shapes of deforming 3D objects from a sequence of 2D measurements obtained using regular cameras. NRSfM has many applications in various fields such as virtual reality (VR), medical surgery or Computer-Generated Imagery (CGI) so there is been a lot of research devoted to this area of computer vision. This problem has been well studied since the seminal work by Bregler et al. [5] was published in the year 2000, generalizing the framework of Rigid Structure from Motion to the non-rigid case. Following this publication, many papers were published trying to improve the original method and studying theoretical aspects of the problem. In what follows, we will introduce the mathematical formulation of the problem, which can be split into two phases: motion and shape estimation. Motion estimation refers to the retrieval of the camera position in each frame and shape estimation refers to the retrieval of the actual 3D structure in each frame. In this thesis, we will tackle both problems separately. On the one hand, we will work on trying to improve the motion estimation part of NRSfM by using Single Rotation Averaging. On the other, we will tackle the shape estimation part of the problem by making an assumption that shapes lie in a temporal union of subspaces and we will use to Sparse Subspace Clustering to jointly learn the shape and cluster frames. There are four papers which shape the direction in which the study of NRSfM has evolved over time and we think are important to understand: 1. Tomasi and Kanade (1992) propose the factorization approach for Structure from Motion (SfM). 2. Bregler et al. (2000) introduce the shape basis constraint for Non-Rigid Structure from Motion (NRSfM). 3. Xiao et al. (2004) proved that the problem is under-constrained and using the orthonormality constraint alone does not recover the shape coefficients uniquely due to an ambiguity in the shape basis. 4. Akhter et al. (2009) prove that even though the solution is not unique any of the solutions will be a valid solution to the structure from motion problem. Next, we will discuss these four papers in detail and explain why they are important in this study. 1.1 Rigid Structure from Motion Before NRSfM, we will first look at Rigid Structure from Motion (SfM). This is a simpler problem since the 3D shape we want to reconstruct is rigid, i.e., it does not change over time. This makes the number of unknowns be greatly reduced, the 3D shape we want to retrieve has 3Punknowns, in the non-rigid case this number increases to 3PF, since we want to know the 3D structure at 3 Figure 1: Non-Rigid Structure from Motion setting. A shape deforms between two consecutive frames (Sfto Sf`1) and the camera moves making the orientation vectors of the camera shift. each frame. The input data is a set of Pfeature points tracked through Fframes, tpuf p,vf pq|f“1 ... F,p“ 1 ... Pu.1The goal is to recover the camera pose and the 3D points corresponding to the tracked features. The camera pose at each frame is given by the coordinates of its origin tfand its orientation, that is a pair of unit vectors pif,jfq. Under the assumption of an orthogonal camera the projection rays are parallel to kf“ifˆjf. Each of the Pfeatures that are tracked will correspond to a certain position in the 3D space, sp“ pxp,yp,zpqJfor p“1 ... P. We set the origin of the world reference system at the centroid of these points 1 PřP p“1sp. Given these definitions, then under orthogonal camera assumption the projection of a particular 3D point can be written as2: uf p“ifpsp´tfqvf p“jfpsp´tfq, (1) Given the tracked features we can define the following matrices: U“¨ ˚ ˝ u1 1... u1 P . . ..... . . uF 1... uF P ˛ ‹ ‚V“¨ ˚ ˝ v1 1... v1 P . . ..... . . vF 1... vF P ˛ ‹ ‚W“ˆU V˙, (2) 1This can be obtained using, for example, the KLT (Kanade-Lucas-Tomasi) feature tracker or modern deeplearning features over a sequence of images of the object we want to reconstruct. 2We will always assume the cameras to behave 4 WPR2FˆPis the measurement or observation matrix. We take the mean of the xand y coordinates at every frame as: af“1 P P ÿ p“1 uf p,bf“1 P P ÿ p“1 vf p, (3) and center the measurement matrix, constructing it instead using ˜uf p“uf p´afand ˜vf p“vf p´bf. The corrected/registered measurement matrix is ˜ W. This is helpful because we can rewrite Eq. (1) as: ˜uf p“uf p´af“ifpsp´tfq´ 1 P P ÿ q“1 uf q. (4) Now, because we set the origin of the world reference system to be at the centroid of the points we track: P ÿ q“1 uf q“ P ÿ q“1 ifpsq´tfq “ if P ÿ q“1 sq´ P ÿ q“1 iftf“ ´Piftf, (5) and ˜uf p“ifpsp´tfq`iftf“ifspand similarly ˜vf p“jfsp. Then, ˜ Wcan be expressed as: ˜ W“R¨S“ ¨ ˚ ˚ ˚ ˚ ˚ ˚ ˚ ˝ i1 . . . iF j1 . . . jF ˛ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‚`s1... sP˘“ˆ˜ U ˜ V˙, (6) where Ris the matrix with the orientation vectors of the camera at every frame (R2Fˆ3) and Sis the a matrix with the 3D positions of the tracked features (R3ˆP). rankpABq ď minprankpAq, rankpBqq. Since the ranks of Sand Rare at most 3, the rank of ˜ Wis at most 3 and we can use a truncated SVD of ˜ Wto decompose it: ˜ W“UΣVJ“ pUΣ1{2qpΣ1{2VJq. (7) Using only the first 3 columns of U, rows of VJand singular values in Σ: ˆ W“ pˆ Uˆ bSigma1{2qp ˆ bSigma1{2ˆ VJq “ RS “ˆ RGG´1ˆ S. (8) However, the SVD is not unique and any 3ˆ3 invertible matrix Gwill give an additional valid decomposition of ˜ Was shown in Eq. (8). This matrix Gis called corrective matrix in the literature and to fins it the orthonormality conditions of the orientation vectors is imposed, such that: iJ fif“ˆ iJ fGGJˆ if“1jJ fjf“ˆ jJ fGGJˆ jf“1iJ fjf“ˆ iJ fGGJˆ jf“0. (9) A solution for such a matrix Gcan be found for example through a least squares approach. 5 1.2 Non-Rigid Structure from Motion However, the problem we study in this work is a little bit different in nature. In NRSfM the object of study are deformable objects that change shape through time [2,3,1], so instead of having a fixed shape matrix Sfor all frames, we instead have a shape Sffor every frame. Similarly to what Tomasi and Kanade did for Rigid reconstruction, for a certain frame fP t1 ... Fu, we have Ptracked points with pu,vqcoordinates in the image reference: ˜upfq 1upfq 2... upfq P vpfq 1vpfq 2... vpfq P¸“Rpfq¨Spfq, (10) with RpfqPR2ˆ3a rotation matrix3and SpfqPR3ˆPa shape matrix with the 3D position of the tracked points in the world reference in that frame. However, NRSfM is ill-posed in the sense that if the shape could arbitrarily deform the problem would simply not be solvable, there would simply be too many parameters to infer. That’s why Bregler et al. [5] assume that the shape Scan be expressed as a linear combination of a set of K basis shapes tS1... SKu, where SiPR3ˆP, so we can write the expression above as follows: R¨S“R¨ K ÿ i“1 ciSi“ K ÿ i“1 ciRSi“`c1Rc2R... cKR˘¨ ˚ ˚ ˚ ˝ S1 S2 . . . SK ˛ ‹ ‹ ‹ ‚ , (11) for some configuration weights ciPR.We are restricting the non-rigid shape to those that can be generated by linear combinations of the basis shapes. Taking the full observation matrix WPR2FˆP: W“¨ ˚ ˚ ˚ ˚ ˚ ˝ up1q 1... up1q P vp1q 1... vp1q P ... upFq 1... upFq P vpFq 1... vpFq P ˛ ‹ ‹ ‹ ‹ ‹ ‚ “¨ ˝ cp1q 1Rp1q... cp1q KRp1q ... cpFq 1RpFq... cpFq KRpFq˛ ‚¨ ˚ ˝ S1 . . . SK ˛ ‹ ‚“M¨B, (12) with MPR2Fˆ3Kand BPR3KˆP. Eq. (12) shows that the tracking matrix has rank 3K. We assume Kăă F,P. We can decompose Winto ˆ Mand ˆ M, by using only the first 3Ksingular values and vectores of its (truncated) SVD, W“UDVJ, writing: W“MB “ˆ MGG´1ˆ B, (13) with ˆ MPR2Nˆ3Kand ˆ BPR3KˆP. Like before the decomposition of Wis not unique and any invertible matrix GPR3Kˆ3Kcan be inserted into the equation as shown above to obtain a new 3ˆr1r2r3 r4r5r6˙ 6 possible decomposition. Again we find this matrix Gby enforcing the orthonormality constraint of the rotation matrices Rpfqfor fP t1 ... Fu. However, reconstructing the entire corrective matrix Gis not necessary. In fact only one of the Kcolumn triplet of the matrix Gis needed. Let ˆ M2i´1:2ibe the ith pair of rows of matrix ˆ M, and let Gkbe the kth column triplet of G4, then from Eq. (12): ˆ M2i´1:2iGk“cik Rii“1, ... , F k “1, ... , K. (14) So, for any kP t1, ... , Kuwe can obtain the rotation matrices for any frame. However, something we will not be discussing throughout this work is how to choose this Kgiven some data. To obtain the camera rotations and the shape coefficients from M, we do another SVD. First take the rows of M: qpfq“ pcpfq 1Rpfq... cpfq KRpfqq “ ˆc1r1c1r2c1r3... cKr1cKr2cKr3 c1r4c1r5c1r6... cKr4cKr5cKr6˙, (15) we can reorder the terms in the matrix as: qpfq“¨ ˝ c1r1c1r2c1r3c1r4c1r5c1r6 ... cKr1cKr2cKr3cKr4cKr5cKr6˛ ‚“¨ ˚ ˚ ˚ ˝ c1 c2 . . . cK ˛ ‹ ‹ ‹ ‚`r1r2... r6˘, (16) which can be obtained via an SVD by enforcing a rank equal to 1. However, obtaining Ris just the first step. What we have done now is just retrieving the orientation of the cameras, a problem which we call motion estimation. However, we still have to retrieve the shape matrix S, this step is called shape estimation. One way of doing so could be to use the pseudoinverse of Rand obtain Sby: S“R:W. (17) In Section 3we will discuss a more sophisticated way of tackling this second part of the problem. 1.3 About the under-constraint nature of NRSfM An important theoretical result in the study of NRSfM is the uniqueness of the solution to the problem and the validity of the solutions. In [5] only the orthogonality constraint of the rotation matrix is enforced. However, Xiao et al. [19] argue these results in ”ambiguous and invalid solutions” because of the non-uniqueness of the shape bases (”any linear transformation of a set of shape bases yields a new set of eligible bases”). In fact, the orthonormality constraint (Eq. (14)) can be rewritten as follows: ˆ M2iGkGJ kˆ MJ 2i´ˆ M2i´1GkGJ kˆ MJ 2i´1“0 , i“1...F(18) ˆ M2i1GkGJ kˆ MJ 2i“0 , i“1...F(19) 4That is GkPR3Kˆ3with the first column being column number 3kof G, the second column being column 3k`1 and the last one column 3k`2. 7 Method 0% outliers 10% outliers 50% outliers Hartley’s l1(with outlier rejection) 0.1975 0.2136 0.3526 Hartley’s l1(without outlier rejection) 0.1962 0.2359 1.1322 Robust SRA 0.1968 0.2149 0.3547 Our implementation (Hartley’s l1) 0.1893 0.2303 1.1307 Table 1: Average angular distance between ground truth and l1-norm mean of SRA methods for N“100 and different ratios of outliers. Figure 2: Steps in Weiszfeld algorithm for SRA In order to check the correctness of our implementation, we will compare its results against the implementation of SRA algorithms published by the authors of ”Robust Single Rotation Averaging” [15]. The authors of the original paper shared their implementation of a new method to perform SRA as well as their implementation of the original algorithm proposed by Hartley [10] but with two key differences. Firstly, they change the initialization of the algorithm. Secondly, they add an outlier rejection schema so at each iteration the rotations that are very far from the current estimate of the average are not taken into account. The code used in the paper is publicly available so we compared our implementation with the ones provided by the authors. To do so, we genereate a random rotation which will be the ground truth that we want to retrieve. Then, we generate Nperturbations of up to 5 degrees of rotation which we will average. We then compute a distance metric (in this case the angle metric) between the original rotation and our average to see how far we are from it. Furthermore, we also check the robustness to outliers by letting a fraction of those Nrotations be completely random rotations, see Fig. 3. The results we obtain are shown in Table 1. These results agree with the findings of the original paper [15], we can clearly see that both their Robust SRA and their implementation of Hartley’s original method with outlier rejection perform better when the fraction of outliers increases. As expected, all methods get worse results as the number of outliers increases. 14 Figure 3: Example of the data used to evaluate SRA. Rotations are shown as the unit vector ˆvfrom their angle-axis representation. The ground truth rotation is shown in green, the perturbations are shown in black and the outliers are shown in red. However, we can see that when we remove the outlier rejection schema in their implementation of SRA it behaves very similarly to my implementation both in the presence and absence of outliers. 2.6 SRA in NRSfM The idea to use rotation averaging for motion estimation is not new. For example, in 2004, already Govindu uses the l2average for this task [9]. Furthermore, in another paper by the same author titled ”Robustness in motion estimation” , they are concerned with the impact that outliers have on the performance on the l2average algorithm they used for motion estimation. To minimize their impact in this paper the author uses a RANSAC approach to identify and discard outliers prior to performing the average. Therefore, when Hartley [11] introduced the l1average, it was also used in order to perform NRSfM. The goal of this section it to replicate the results of ”Organic Priors in Non-Rigid Structure from Motion” by Kumar et al. [14]. The idea is the one we stated at the start of this section: compute Kestimates of the rotation matrices at each frame, Rfand obtain an average of those Kestimates using the algorithms we presented in this section. As we saw in Eq. (14) one does not need the entire corrective matrix to obtain the rotations and the shape for each frame, instead only a column triplet is required. So, the idea is to obtain the full corrective matrix and from each of the Kcolumn triplets obtain Krotation matrices for each frame and then perform SRA. See Fig. 4for a schematic view of the procedure. The paper by Kumar states that they use Dai’s method [6] which is based on a semidefinite programming (SDP) approach to obtain the column triplets. To understand it, let us define some concepts. Let Qk“GkGJ kPR3Kˆ3Kfor some column triplet Gkand qk“vecpQkqwhere vecp¨q 15 Figure 4: Diagram showing how to obtain a single rotation from the full corrective matrix. Picture taken from [14]. is simply the vectorization operation. This matrix Qkis the Gram matrix of a column-triplet from which we can retrieve Gkusing a truncated SVD of rank 3. Then we can define a matrix A: Ai“ˆpˆ Mibˆ Miqp1, :q´pˆ Mibˆ Miqp4, :q pˆ Mibˆ Miqp2, :q˙, (30) A“rAJ 1,AJ 2, ..., AJ FsJ, (31) where ˆ Midenotes the double columns of ˆ M, i.e., ˆ M2i:2i´1. Now, we can rewrite Eq. (14) as: Aqk“0. (32) Then, the algorithm proposed by Dai et al. [6] follows from the fact that Qkhas low rank (rank 3) and is positive semidefinite since it is defined as a Gram matrix. The idea is to formulate an optimization problem whose solution is the precisely a matrix satisfying the above conditions to obtain as a solution Qk, from which the Gkcan be obtained. Instead of imposing Qkto have rank 3, we relax this condition and define a rank minimization problem. However, rank minimization is an NP-hard problem, so instead a nuclear norm minimization problem is solved. The nuclear norm is defined as the sum of the the eigenvalues which for the case of symmetric and positive semidefinite matrices is simply the trace. All of this conditions motivate the authors to formulate the following optimization problem to retrieve Qk: min tracepQkq(33) s.t.Qkě0 AvecpQkq “ 0 16 Once we solve this optimization using an SDP solver and obtain Qk,Gkis retrieved using SVD. Additionally after finding Gkthe authors also use a non-linear refinement procedure, via an unconstrained minimization: min Gk F ÿ i“1» –˜1´ˆ M2iGkGJ kˆ MJ 2i ˆ M2i´1GkGJ kˆ MJ 2i´1¸2 `˜2ˆ M2i´1GkGJ kˆ MJ 2i ˆ M2i´1GkGJ kˆ MT 2i´1¸2fi fl. (34) Therefore, this method returns just one column triplet and Kumar [14] does not specify how to obtain the full corrective matrix. However, the most well-known method to find the full corrective matrix from a single column-triplet is Brand’s method. Brand’s method Now, having found a column triplet Gkby Dai’s method [6], we have to find the full corrective matrix. Brand’s method is used to obtain the full corrective matrix from a single one of its columntriplets. Assume one of the column triplets, G1for example, has already been estimated and used to find the truncated rotation matrix for every frame, Rf“ pif,jfqJ, has been obtained using Eq. (14). Then, take the cross product of the orientation vectors to obtain the third row of the rotation matrix zfand notice that for every fP t1...Fu: Rfzf“0 c1fRfzf“0ÝÑ ˆ MfG1zf“0 c2fRfzf“0ÝÑ ˆ MfG2zf“0 ... cKf Rfzf“0ÝÑ ˆ MfGKzf“0 Using the vectorization operation, we obtain that: pzJ fbˆ MfqvecpGkq “ 0@f,k, (35) which means that all column triplets vecpGkqbelong to the null space of the matrix: N“¨ ˝ zJ 1bˆ M1 ... zJ Fbˆ MF ˛ ‚PR9KˆK. (36) Therefore, Brand’s method consists in using a previously obtained column triplet to then compute Nand a basis for its null space, which will give the full corrective matrix G. However, after coding the procedure in Matlab, we do not get very good results. We tested the code on the ”yoga” sequence kith K“10 and found that the column triplets that were obtained through Brand’s method perform worse than the original triplet found via Dai’s method [6] (see 17 Figure 5: Rotation error obtained from each column triplet obtained by Brand’s method in yoga sequence. The original triplet that was used to obtain the full corrective matrix obtained a lower error, as shown by horizontal line (Rotation error “0.088) Fig. 5). Hence, applying SRA using the rotations found with the full corrective matrix did not give a better estimate of the ground truth rotation than the original triplet found with Dai’s method [6]. Another option that we also explored was to use the unconstrained optimization problem used for refining the estimate of Gkin Eq. (34). The idea was to use random initialization of Gkand obtain several estimates of Gkwith the idea of doing the same thing we planned to do with the full corrective matrix: obtain several estimates of the rotations and perform SRA. However, we found that the optimization was very sensitive to the initialization and the obtained Gkalso gave pretty bad results. 18 3. Sparse subspace clustering In this section, we explore sparse subspace clustering in the shape estimation part of the NRSfM problem. Remember that the basic assumption we have been working under throughout this thesis is that the shape at every frame is a linear combination of some basis shapes. That is, we assume the shapes at each frame all belong to the same subspace. Now, we alter that assumption to be somewhat more general and consider the shapes to belong to union of subspaces. The basic idea is to jointly learn the shape but also cluster frames in the input data that belong to the same ”group” or ”subspace”. For a video of a human performing several tasks one could understand a subspace as the action that is being performed. If the video contains a section where a person is walking up to a table and then picking up a glass of water and drinking it, the video could be grouped into two distinct sections. Let XPRnˆNbe a collection of Ndata points in Rn, that is, every column xiis a data point. Furthermore, assume each column is samples from a union of Kindependent linear subspaces tSkuK k“1each of dimension dkăă n. The subspace clustering problem refers to the problem of finding the number of subspaces, their dimensions, a basis for each subspace, and the segmentation of the data from X. Now, assume we want to represent each xias a linear combination of the other columns. Naturally, we would expect that the contribution of the columns in the same subspace be high. Conversely, the contribution of the columns in other (independent) subspaces should be zero. In fact, if for every xiPSkthe rest of the columns of Xthat belong to Skspan the entire subspace we could find a matrix CPRNˆNsuch that: X“XC (37) s.t. Cii “0 and such that if xiand xjbelong to different subspaces then Cij ‰0. However, due to noisy measurement this will not be the case in general. In this case, what we want to do instead is finding the solution to the following optimization problem: min CPRNˆNp||C||l`||X´XC||l1q(38) s.t. diagpCq “ 0, where Cis the matrix with the coefficient of the linear combinations of the columns and X´XC is the ”error” due to outlying entries [12]. Different choices of norms can be used for both || ¨ ||land || ¨ ||l1. A natural choice of the first one is, for example, the l1-norm since it is a well-known fact that it encourages sparsity. For the second norm, one could use the Frobenius norm penalizing big errors [8]. We will be using these two throughout the thesis The algorithms to perform sparse subspace clustering generally follow the following steps: 19 - Finding coefficient matrix Cby solving the optimization problem in Equation 38 - Construct a weighted graph with affinity matrix W“ |C|`|CJ|where each node corresponds to one of the columns of X(one node for each data point). - Finally, we use spectral clustering on Wto cluster nodes (data points) and obtain the grouping we desire. Spectral clustering is a well-known algorithm so the only part of this entire process that requires some work is finding the coefficient matrix C. Let us see now look at a classic algorithm to solve the SSC clustering problem. 3.1 SSC with ADMM To solve this problem Elhamifar et al. [8] propose the following procedure. The idea is to introduce an auxiliary matrix S: min C,A||C||1`µ 2||X´XA||2 F(39) s.t. A“C´diagpCq AJ1“1 Then, you can use the Alternating Direction Multiplier Method (ADMM) by minimizing the Augmented Lagrangian with respect to Cand Aindependently and iteratively at each step: LpC,Aq “||C||1`µ 2||X´XA||2 F(40) `ρ 2p||A´pC´diagpCqq||2 F`||AJ1´1||2 2q ` ă λ1,pAJ1´1q ą ` ă λ2,A´pC´diagpCqq ą, where λ1PRNˆNand λ2PRNare the Lagrange multipliers and the inner product between matrices is defined as ăA,Bą“ trpABJq. To obtain the update rules at each iteration we derive the Lagrangian with respect to Cand Aand obtain the following update rules: Api`1q“ pµXJX`ρId`ρ11Jq´1pµXJX`ρp11J`Cpiqq´1λT 1´λ2q, (41) Ci`1“J´diagpJqwhere J“T1{ρpAi`1`λpiq 2{ρq, (42) with T1{ρpXqapplying signpxij qmaxp0, xij ´1{ρqto all the elements of the matrix X. 20 Figure 6: Left: 3 subspaces in R3: a plane (2D) and two lines (1D). The circles indicate some samples in those subspaces. Right: Example of the resulting affinity matrix when using SSC on a synthetic dataset. We take 3 randomly generated subspaces of R10 of dimensions 3, 5 and 7 and we sample 50 points from each of the subspaces and apply SSC. The horizontal bar at the top shows the results of the clustering algorithm. 3.2 Joint shape reconstruction and SSC The idea now is to use SSC to reconstruct the shape matrix S. To do so we have to modify the optimization problem slightly to add the constraints on Simposed by the NRSfM problem, that is W“RS: min S,C,A||C||1`µ 2||S#´S#A||2 F(43) s.t. A“C´diagpCq AJ1“1 W“RS where S#PR3PˆFis simply a reordering of the elements of SPR3FˆPso that each columns f contains the x,yand zcoordinates of all points in frame f. To handle this relationship during the optimization, we add another hard constraint to encode the relationship between Sand S#. We additionally add a term to penalize the rank of the shape matrix. This is inspired by a work from Kumar et al. [13] where a similar problem is tackled with SSC, but instead of clustering frames their goal was to cluster points together (a cluster of points represents an individual body). We do this by adding the nuclear norm of S#to the objective function since it is a known fact that it is a convex relaxation of the rank function. 21 min S,S#,C,A||C||1`µ 2||S#´S#A||2 F`γ||S#||˚(44) s.t. A“C´diagpCq AJ1“1 W“RS fpSq “ S# where ||S#||˚is the nuclear norm of S#and fpSqand its inverse are defined as: S#“fpSq“pSJbI3qP S“f´1pS#q“pI3bpS#qJqT for binary matrices Pand T. Now, the corresponding Augmented Lagrangian is: LpS,S#C,Aq “||C||1`µ 2||S#´S#A||2 F`(45) ρ 2p||A´pC´diagpCqq||2 F`||AJ1´1||2 2`||W´RS||2 Fq`||fpSq´S#||2 Fq ăλ1,A´pC´diagpCqq ą ` ă λ2,pAJ1´1qą` ăλ3,W´RS ą`ăλ4,fpSq´S#ą where λ3PR2FˆP. After adding these new constraints the update rules for Cand Aremain unchanged, however we need to derive the update rules for Sand S#. First, to obtain the update rule for Swe have to minimize the Lagrangian with respect to Sat each iteration. However, to make it easier to derive the expression above with respect to S, we can use f´1to rewrite the Lagrangian and not have to deal with fpSq: min S ρ 2p||W´RS||2 F`||fpSq´S#||2 Fq` ă λ3,W´RS ą`ăλ4,fpSq´S#ą“ “min S ρ 2p||W´RS||2 F`||S´f´1pS#q||2 Fq` ă λ3,W´RS ą`ăf´1pλ4q,S´f´1pS#q ą Now, we can simply derive with respect to S: BL BS“B BS”ρ 2p||W´RS||2 F`||S´f´1pS#q||2 Fq` ă λ3,W´RS ą`ăf´1pλ4q,S´f´1pS#q ąı “´ρRJpW´RSq`ρpS´f´1pS#qq´RJλ3`f´1pλ4q Once again, setting the expression above to 0 and isolating Swe obtain the update rule: 22 S“ pId`RJRq´1„f´1pS#q` RJλ3 ρ´f´1pλ4q ρ`RJWȷ Now, in order to update S#, we have to solve the following minimization problem: min S#γ||S#||˚`µ 2||S#´S#A||2 F s.t.fpSq “ S# In order to solve this, we compute the derivative of the Lagrangian with respect to S#, disregarding the nuclear norm term, which can easily be done considering that ||X||2 F“trpXXJqand the following two identities from matrix calculus B BXtrpAXJq “ Aand B BXtrpXAXJq “ XAJ`XA [17]: BL BS#“µ 2B||S#´S#A||2 F BS#`ρ 2B||fpSq´S#||2 F BS#`B ă λ4,fpSq´S#ą BS# “µS#pI´AqpI´AqJ´ρpfpSq´S#q´λ4 Setting the expression above equal to 0 gives the following update rule for S#: S#“ pλ4`ρfpSqqpµpI´AqpI´AqJ`ρIq´1 Finally, we just apply the Singular Value Thresholding algorithm to S#. The final algorithm to perform the reconstruction of the shape matrix using SSC is presented in Algorithm 3. 23