3-D Rigid Models from Partial Views—
—Global Factorization

Since the strongest cue to infer 3-D shape from video is the 2-D motion of the image brightness pattern, our problem has been often referred to as structure from motion (SFM). Since using a large number of views, rather than simply two consecutive frames, leads to more constrained problems, thus to more accurate 3-D models, current research has been focused on multi-frame SFM. The factorization method introduced in [38] overcomes the difficulties of multi-frame SFM—nonlinearity and large number of unknowns—by using matrix subspace projections. In [38], the trajectories of feature points are collected into an observation matrix that, due to the rigidity of the 3-D object, is highly rank deficient in a noiseless situation. The 3-D shape of the object and the 3-D motion of the camera are recovered from the rank deficient matrix that best matches the observation matrix. The work of [38] was extended in several ways, e.g., geometric projection models [37, 30], 3-D shape primitives [34, 31, 4], recursive formulation [27], uncertainty [28, 2, 6, 5], multibody scenario [15]. Besides recovering 3-D rigid structure from video sequences, recent approaches to several other image processing/computer vision problems require determining linear or affine low dimensional subspaces from noisy observations. These problems include object recognition [9, 39], applications in photometry [35, 29, 10], and image alignment [42, 43].

In general, such low dimensional subspaces are found by estimating rank deficient matrices from noisy observations of their entries. When the observation matrix is completely known, the solution to this problem is easily obtained from its Singular Value Decomposition (SVD) [16]. However, in practice, the observation matrix may be incomplete, i.e., some of its entries may be unknown (unobserved). When recovering 3-D structure from video, the observation matrix collects 2-D trajectories of projections of feature points [38, 37, 30, 27, 6, 5] or other primitives [34, 31, 28, 4]. In real life video clips, these projections are not visible along the entire image sequence due to the occlusion and the limited field of view. Thus, the observation matrix is in general incomplete. For the specific cases we focus in this paper—videos that show views all around (non-transparent) 3-D objects—this incomplete characteristic of the observation matrix is particularly noticeable.

The problem of estimating a rank deficient matrix from noisy observations of a subset of its entries deserved then the attention of several researchers. References [38] and [24] propose sub-optimal solutions. In [38], the missing values of the observation matrix are “filled in”, in a sequential way, by using SVDs of observed submatrices. In [24], the author proposes a method to combine the constrains that arise from the observed submatrices of the original matrix. More recently, several bidirectional optimization schemes were proposed, e.g., [25, 19, 12, 1, 23, 40, 14], and the problem was also addressed by using nonlinear optimization [13].

However, none of the methods above deal with “self-inclusion”, i.e., with the fact that a region that disappears due to self-occlusion may appear again later. When this happens, the re-appearing region is usually treated as a new region, i.e., as a region that was never seen before. This procedure has two drawbacks: i) the problem becomes less constrained than it should, thus leading to less accurate estimates of 3-D structure; and ii) further 3-D processing is needed to fuse the recovered multiple versions of the same real-world regions.

We propose a global approach to recover complete 3-D models from video. Global in the sense that our method computes the simplest 3-D rigid object that best matches the entire set of 2-D image projections. This way we avoid having to post-process several partial 3-D models, each obtained from a smaller set of frames, or an inaccurate 3-D model obtained from the entire set of frames without detecting re-appearing regions. We develop a global cost function that balances two terms—model fidelity and complexity penalization. The model fidelity term measures the error between the model—a 3-D shape and a set of re-appearing regions—and the observations, as in Maximum Likelihood (ML) estimation. This error is simply given by the distance of a re-arranged observation matrix to the appropriate space of rank deficient matrices. The penalty term measures the complexity of the 3-D model, which is easily coded by the number of feature points used to describe the observations. By minimizing this global cost, we get what statisticians usually call a Penalized Likelihood (PL) [18] estimate of the 3-D structure. Through PL estimation, re-appearing regions are then detected when the increase of the complexity of the 3-D model does not compensate a slightly better fit to the observations, meaning that a more complex 3-D model would fit the observation noise rather than the 3-D real-world object.

The minimization of the PL cost function requires the computation of rank deficient matrices from partial observations. In the paper, we study two iterative algorithms developed for this problem. The first algorithm is based in a well known method to deal with missing data—the Expectation-Maximization (EM) [26]. Although the authors don’t refer it, the bidirectional scheme in [25] is also EM-based. However, as detailed below, the EM algorithm we propose is more general and computationally simpler than the one in [25]. Our second two-step iterative scheme is similar to Wiberg’s algorithm [41] and related to the one used in [36] to model polyhedral objects. It computes, alternately, in closed form, the row space matrix and the column space matrix whose product is the solution matrix. We call this the Row-Column (RC) algorithm.

In the paper, we illustrate the behavior of both iterative algorithms with simple cases and evaluate their performance with more extensive experiments. In particular, we study the impact of the initialization on the algorithm’s behavior. From these experiments, we conclude that the RC algorithm is more robust than EM in what respects to the sensibility to the initialization. Furthermore, the number of iterations needed for good convergence and the computational cost of each iteration are both smaller for RC than for EM. Obviously, the performances of both EM and RC improve when the initial estimate provided to the algorithms is more accurate. Any sub-optimal method, e.g. the ones in [38, 24], can be used to compute such an initialization. Our experience shows that, with a simple initialization procedure, even for high levels of noise and large amount of missing data, both iterative algorithms converge: i) to the global optimum; and ii) in a very small number of iterations. Naturally, under our global factorization approach, the rank deficient matrices can also be computed by using nonlinear methods such as [13].

We use EM and RC in the context of recovering complete 3-D models from video, through the minimization of our global PL cost. Our experiments show that fitting the rank deficient matrix to the entire (re-arranged) observation matrix (which, naturally, misses several entries) leads to better 3-D reconstructions than those obtained by combining partial 3-D models estimated by fitting submatrices to smaller subsets of data (each corresponding to a subset of features that were visible in a subset of consecutive frames).

In section II we overview the matrix factorization approach to the recovery of 3-D structure. Section III details the problem addressed in this paper—the direct recovery of complete 3-D models from incomplete 2-D trajectories of features points. In section IV we describe our solution as the minimization of a global cost function. Section V describes the algorithms used to estimate rank deficient matrices from incomplete observations. In sections VI and VII, we present experiments and an application to video compression. Section VIII concludes the paper.

Preliminary versions of parts of this work appeared in [19, 20, 21, 17].

In our approach, re-appearing features are detected through the minimization of a global cost-function that takes into account both the accuracy of the model and its complexity. Thus, the re-arranged matrix ${\mathbf{W}}_{R}$ is given by

where $E_{{\mathcal{S}}_{4}}\left({\mathbf{W}}_{r}\right)$ measures the error of the model ${\mathbf{W}}_{r}$ as its distance to the space of rank 4 matrices, i.e., it is a likelihood term, and $P_{r}$ is the number of feature points of the model, i.e., it codes the model complexity (matrix ${\mathbf{W}}_{r}$ is $2F\times P_{r}$). The parameter $\alpha$ balances the two terms. Naturally, the minimization in (2) is constrained to the set of matrices that can be found by re-arranging the original observations matrix ${\mathbf{W}}$. We denote this set by ${\mathcal{R}}({\mathbf{W}})$.

The likelihood term $E_{{\mathcal{S}}_{4}}\left({\mathbf{W}}_{r}\right)$ in (2) models the object rigidity. As outlined in section II, in a noiseless situation, the matrix collecting the trajectories of points of a rigid object is rank 4. This term is then naturally given by the distance of the matrix ${\mathbf{W}}_{r}$ to the space of the $2F\times P_{r}$ rank 4 matrices. Since not all the entries of ${\mathbf{W}}_{r}$ are known, we measure this distance by masking the unknown entries:

where $\left\|.\right\|_{F}$ represents the Frobenius norm, which arises from assuming that the observations are corrupted by white Gaussian noise, ${\mathcal{S}}_{4}$ denotes the space of the $2F\times P_{r}$ rank $4$ matrices, $\odot$ represents the elementwise product, also known as the Hadamard product, and the matrix ${\mathbf{M}}_{r}$ is a binary mask that accounts for the known entries of the matrix ${\mathbf{W}}_{r}$, i.e., the entry $(i,j)$ of ${\mathbf{M}}_{r}$ is $m_{ij}=1$ if the entry $w_{ij}$ of ${\mathbf{W}}_{r}$ is known and $m_{ij}=0$ otherwise.

In Bayesian inference approaches, the term $\alpha P_{r}$ in (2) is considered to be a prior and the minimization (2) leads to the Maximum a Posteriori (MAP) estimate [11] (the cost function in (2) is the negative log posterior probability). Penalization terms, such as $\alpha P_{r}$ in (2), have also been introduced in the model selection literature through information-theoretic criteria like Akaike’s AIC [32] or the Minimum Description Length (MDL) principle [8]. Naturally, different basic principles lead to different choices for the balancing parameter $\alpha$ but the structure of the cost function is most often as in (2). This generic form (2) is usually called by the statisticians a penalized likelihood (PL) cost function [18]. In our case, to find a valid range for the weight parameter $\alpha$, we performed several experiences. By testing with pertinent values for the number of features, number of images, noise level, $\%$ of missing data, and number of re-appearing features, we concluded that there is a wide range of values for $\alpha$ that leads to a good balance between maximizing the number of correct detections of re-appearing features (probability of detection) and minimizing the number of incorrect detections (probability of false alarm).

To find candidates ${\mathbf{W}}_{r}\in{\mathcal{R}}({\mathbf{W}})$, our algorithm starts by selecting from the original observation matrix ${\mathbf{W}}$, pairs of columns that can be merged with each other. Naturally, these pairs correspond to disjoint tracking periods. For example, for the matrix ${\mathbf{W}}$ in the right side of Fig. 1, each one of the six last columns could be merged with the first column, therefore, in what respects to the feature corresponding to the first column, there are seven possible situations: it could either have re-appeared, generating one of the trajectories of the six last columns, or remain occluded for the remaining of the video.

To obtain a computationally feasible algorithm, we prune the search—our algorithm decides by comparing the costs of merging each selected pair of columns with the one of considering that there are not re-appearing features, i.e., of the model ${\mathbf{W}}_{r}={\mathbf{W}}$. The process is then repeated until the cost (2) does not decrease by merging any selected pair of columns. The search could be further pruned, e.g., using the local approach outlined in the previous section to guide the process, thus testing only pairs of columns corresponding to feature points that are close in 3-D.

What remains to be addressed is the computation of the likelihood term $E_{{\mathcal{S}}_{4}}({\mathbf{W}}_{r})$, given by expression (3). The minimization in (3) requires finding the rank 4 matrix that best matches the incompletely known matrix ${\mathbf{W}}_{r}$. There is not known closed-form solution for this problem. In the following section, we develop iterative algorithms that compute such rank deficient approximations in very few iterations. The likelihood term $E_{{\mathcal{S}}_{4}}({\mathbf{W}}_{r})$ is then computed as

where $\widehat{{\mathbf{W}}_{r}}$ is the rank 4 matrix that best matches the known entries of the matrix ${\mathbf{W}}_{r}$, i.e.,

The iterative algorithms we describe in the sequel converge in very few iterations to $\widehat{{\mathbf{W}}_{r}}$, when adequately initialized. In the appendix, we describe a strategy to compute such initialization. We reduce the computational burden of the global approach by performing this initialization step only once, when testing the model that corresponds to the original observation matrix, ${\mathbf{W}}_{r}={\mathbf{W}}$, and using the same initial guess when testing the other models.

The Expectation-Maximization (EM) approach to estimation problems with missing data works by enlarging the set of parameters to estimate—the data that is missing is jointly estimated with the other parameters. The joint estimation is performed, iteratively, in two alternate steps: i) the E-step estimates the missing data given the previous estimate of the other parameters; ii) the M-step estimates the other parameters given the previous estimate of the missing data, see [26] for a review on the EM algorithm.

In our case, given the initial estimate $\widehat{{\mathbf{W}}_{r}}^{(0)}$, the EM algorithm estimates in alternate steps: i) the missing entries of the input matrix ${\mathbf{W}}_{r}$; ii) the rank $4$ matrix $\widehat{{\mathbf{W}}_{r}}$ that best matches the “complete” data. The algorithm performs these two steps until convergence, i.e., until the error measured by the Frobenius norm in (5) stabilizes.

We now describe the Row-Column (RC) algorithm—another iterative approach, similar to Wiberg’s algorithm [41], to the estimation of a rank deficient matrix that best matches an incomplete observation. From our experience, summarized in section VI, the RC algorithm is not only computationally cheaper than EM, avoiding SVD computations and exhibiting faster convergence, but also more robust than EM to initializations far from the solution.

For the RC algorithm, we parameterize the $2F\times P_{r}$ rank $4$ matrix $\widetilde{{\mathbf{W}}}$ in (5) as the product

where $\widetilde{{\mathbf{A}}}$ determines the column space of $\widetilde{{\mathbf{W}}}$ and $\widetilde{{\mathbf{B}}}$ its row space. The rank 4 matrix $\widehat{{\mathbf{W}}_{r}}$ that best matches ${\mathbf{W}}_{r}$ is obtained by minimizing the cost function in (5) with respect to (wrt) the column space and row space matrices, i.e.,

By using the re-parameterization (11), we map the constrained minimization (5) wrt $\widetilde{{\mathbf{W}}}\!\in\!{\mathcal{S}}_{4}$ into the unconstrained minimization (12) wrt $\widetilde{{\mathbf{A}}}$ and $\widetilde{{\mathbf{B}}}$.

We minimize (12) in two alternate steps: i) the R-step assumes the column space matrix ${\mathbf{A}}$ is known and estimates the row space matrix ${\mathbf{B}}$; ii) the C-step estimates ${\mathbf{B}}$ for known ${\mathbf{A}}$. The algorithm is initialized by computing ${\mathbf{A}}^{(0)}$ from the initial estimate $\widehat{{\mathbf{W}}_{r}}^{(0)}$ and it runs until the value of the norm in (12) stabilizes.

When there is no missing data, i.e., when ${\mathbf{M}}_{r}\!=\!{\mathbf{1}}$ in (12), the R- and C- steps above are linear least squares (LS) problems whose solutions are obtained by using the Moore-Penrose pseudoinverse [16],

If we write steps R and C together as a recursion on one of the matrices ${\mathbf{A}}$ or ${\mathbf{B}}$, say, on the column space matrix ${\mathbf{A}}$, we get

which shows that, in the no missing data case, our RC algorithm is in fact implementing the application of the power method [16] to the matrix ${\mathbf{W}}_{r}{\mathbf{W}}_{r}^{T}$ (the factor $({\mathbf{A}}^{T}{\mathbf{W}}_{r}{\mathbf{W}}_{r}^{T}{\mathbf{A}})^{-1}{\mathbf{A}}^{T}{\mathbf{A}}$ is the normalization). The power method has been widely used to avoid the computation of the entire SVD when fitting rank deficient matrices to complete observations. We will see that, even when there is missing data, both R- and C- steps admit closed-form solution and the overall algorithm generalizes the power method in a very simple way.

We start by a simple case that allows an illustrative graphical representation—approximating $2\times 2$ rank 1 matrices from observations ${\mathbf{W}}_{r}$ that miss one of the entries. In this case, the error measured by the Frobenius norm in expressions (4) and (12), can be expressed in terms of a single parameter $\theta$. In fact, let the $2\times 2$ rank 1 matrix $\widetilde{{\mathbf{W}}}$ be written in terms of its column and row spaces as in (11),

Without loss of generality, impose that the row vector ${\mathbf{b}}$ has unit norm and write it in terms of a row angle $\theta$,

Now denote the minimum of (12) wrt the column space ${\mathbf{a}}$ for fixed row space ${\mathbf{b}}$, i.e., for fixed $\theta$, by ${{\mathbf{a}}}({\mathbf{W}}_{r},\theta)$, given by (19). The error in (4) and (12) is then rewritten as a function of $\theta$ as

Note that for any set of three entries of a $2\times 2$ matrix, there is always a value for the forth entry that makes the rank of the matrix equal to one, i.e., there is always a $2\times 2$ rank 1 matrix $\widetilde{{\mathbf{W}}}$ that fits exactly the observed entries of ${\mathbf{W}}_{r}$. Thus, we have $\min\mbox{error}(\theta)=0$.

The examples in the sequel illustrate the impact of the initialization on the behavior of the EM and RC algorithms with experiments that use the following ground truth rank 1 matrix $\widetilde{{\mathbf{W}}}$, observation mask ${\mathbf{M}}_{r}$, and observation ${\mathbf{W}}_{r}$:

where “$?$” represents the unobserved entry $w_{22}$ of the observation matrix ${\mathbf{W}}_{r}$.

Using the initial guess

we describe the evolution of the estimates $\widehat{{\mathbf{W}}_{r}}^{(k)}$ of the rank 1 matrix by plotting two equivalent representations of $\widehat{{\mathbf{W}}_{r}}^{(k)}$: i) its row space ${\mathbf{b}}=[b_{1},b_{2}]$; and ii) the corresponding angle $\theta$ as defined in (21). The left plot of Fig. 3 shows the level curves of the error function as function of the row vector ${\mathbf{b}}=[b_{1},b_{2}]$, superimposed with the evolution of the estimates of ${\mathbf{b}}$ for EM (dashed line) and RC (solid line). In this plot, the dotted line (optimum) are the row vectors that lead to zero estimation error. The right plot of Fig. 3 represents the same error function, now as a function of $\theta$, as defined in (22), (dotted line) superimposed with the locations of the $\theta$ estimates for EM (dashed line) and RC (solid line).

From the left plot of Fig. 3, we see that both EM and RC trajectories start at the same point (due to the equal initialization) and converge to points in the optimal line. As expected, the EM estimates of the row space vector have constant unit norm (due to the normalization in the SVD) while the RC estimates don’t. The good behavior of the algorithms is confirmed by the right plot of Fig. 3 that shows that both algorithms converge to a value of $\theta$ that makes $\mbox{error}(\theta)=0$, i.e., that minimizes $\mbox{error}(\theta)$.

To evaluate the convergence speed, we show in Fig. 4 the evolution of the estimation error along the iterative process for both algorithms (in the left plot, linear scale, in the right one, logarithmic scale). From Fig. 4, we see that RC converges in a smaller number of iterations than EM. Our experience with larger matrices in practical applications have confirmed the faster convergence of RC.

We now illustrate a drawback of the EM algorithm. Using the same data and the initial estimate

we get the first $1000$ iterations of Fig. 5, which is as described above for Fig. 3.

From the left plot of Fig. 5, we see that, while RC converges to the optimal line, the estimates given by the EM almost don’t change along the iterative process. This can also be seen in the right plot of Fig. 5, which shows that RC converges to the $\theta$ such that $\mbox{error}(\theta)=0$, while EM, after $1000$ iterations is still far from $\arg\min\mbox{error}(\theta)$. The left plot of Fig. 6 shows the evolution of the estimation errors. See that while the error of RC converges to zero in a few iterations, the error of EM almost doesn’t decrease during the first $100$ iterations.

The bad behavior of EM is due to the large initial guess for $w_{22}$ (note that $\widehat{w}^{(0)}_{22}$ is large when compared with the known entries of ${\mathbf{W}}_{r}$). Remember that EM starts by estimating the rank 1 matrix that best matches the initial guess $\widehat{{\mathbf{W}}_{r}}^{(0)}$. This rank 1 matrix is the matrix associated with the largest singular value of $\widehat{{\mathbf{W}}_{r}}^{(0)}$, which is highly constrained by the large spurious entry $\widehat{w}^{(0)}_{22}=22$ (note that while the singular values of the ground truth matrix $\widetilde{{\mathbf{W}}}$ are $\sigma_{1}(\widetilde{{\mathbf{W}}})\!\simeq\!4.9$ and $\sigma_{2}(\widetilde{{\mathbf{W}}})\!=\!0$, the singular values of the initial guess $\widehat{{\mathbf{W}}_{r}}^{(0)}$ are $\sigma_{1}(\widehat{{\mathbf{W}}_{r}}^{(0)})\!\simeq\!22.2$ and $\sigma_{2}(\widehat{{\mathbf{W}}_{r}}^{(0)})\!\simeq\!0.8$ due to the large entry $\widehat{w}^{(0)}_{22}=22$). Then, EM replaces the known entries of ${\mathbf{W}}_{r}$ in the new estimate (obtaining thus an estimate that is very close to the initial guess) and repeats the process. To better illustrate this very slow convergence of EM, we represent, in the right plot of Fig. 6, the evolution of the largest singular value of the estimate $\widehat{{\mathbf{W}}_{r}}^{(k)}$ for both algorithms. We see that while, as expected, the largest singular value of the RC estimates converges to the largest singular value of the solution $\widetilde{{\mathbf{W}}}$, $\sigma_{1}(\widetilde{{\mathbf{W}}})\!\simeq\!4.9$, the largest singular value of the first $100$ iterations of the EM estimates changes very slowly from its initial value $\sigma_{1}(\widehat{{\mathbf{W}}_{r}}^{(0)})\!\simeq\!22.2$.

The behavior just described also happens in situations other than the initial guesses of the unknown entries being too large when compared to the other entries. In fact, we observed the same behavior in situations where the observation matrix had large entries in rows or columns that contained missing entries. In these situations, due to those large entries, even small values for initial guess of the unknown entries had large impact on the row and column singular vectors associated with the large singular values that determined the best rank deficient approximations to the complete data matrices involved in EM. This lead to the same kind of very slow convergence illustrated in Fig. 5 and Fig. 6.

We now present results of approximating matrices likely to appear when processing real video sequences. We tested the EM and RC algorithms with noisy partial observations of rank 4 matrices. We used matrices with dimensions ranging from $2\times 2$ to $200\times 200$. The percentage of missing entries ranged from $10\%$ to $80\%$. We studied the behavior of the EM and RC algorithms in terms of the observation noise power and the amount of missing data, and the impact of the initialization.

We initialized both EM and RC with random values for the unknown entries. To illustrate the influence of the initialization on the convergence of the algorithms, Fig. 7 shows the percentage of experiments that converged (to an estimate close enough to the ground truth) in less than 100 iterations, as a function of the mean value of the random guesses for the unknown entries. Although representative for the entire range of experiments done, the results in the plots of Fig. 7 were obtained with noisy observations of $24\times 24$ rank 4 matrices that missed $70\%$ of their entries. The left plot of Fig. 7 shows three lines, each obtained by using EM with data generated with a ground truth matrix whose elements had mean $0.001$, $1$, and $1000$. The percentages of convergence for the RC algorithm are in the right plot.

The left plot of Fig. 7 shows that the EM algorithm converges almost always if the mean value of the initial guesses for the missing entries is smaller than the mean value of the observed entries. When we increase the values of the initial estimates of the missing entries, the percentage of convergence decreases abruptly, becoming close to zero when those values become much larger than the ones of the observed entries. This is in agreement with the behavior illustrated in the example of Fig. 5. In opposition, the right plot of Fig. 7 shows that the behavior of the RC algorithm is somewhat independent of the order of magnitude of the initialization. The experiments that lead to a non-convergent behavior of RC were such that the matrices whose inverse is computed in (19) and (18) were close to singular. We thus conclude that it is very important in practical applications to provide good initial estimates for both EM and RC algorithms.

Finally, we note that the relevance of a good initialization goes behind avoiding non-convergent behavior. In fact, we observed that both the amount of missing data and the noise level have strong impact on the algorithm’s convergence speed. Thus, when dealing with large percentages of missing entries and high levels of noise, as it may arise in practice, a better initialization not only improves the chance of a convergent behavior but also leads to a faster convergence.

We now use an initial guess $\widehat{{\mathbf{W}}_{r}}^{(0)}$ obtained by combining the column and row spaces given by the SVDs of the known submatrices of ${\mathbf{W}}_{r}$ as described in the appendix. In our tests, with this simple initialization procedure, $100\%$ of the runs of EM and RC converged to the ground truth matrix (mean square estimation error smaller than the observation noise variance) in a very small number of iterations, typically less than 10, even for high levels of noise.

To illustrate the impact of a better initialization, we use a noisy observation ${\mathbf{W}}_{r}$ of a $40\times 40$ rank deficient matrix $\widetilde{{\mathbf{W}}}$ with $900$ missing entries (observation noise variance $\sigma^{2}=1$). Figure 8 shows plots of the evolution of the estimation error, measured by the Frobenius norm, for both the EM and RC algorithms. While for the left hand side plot we used a random initialization, for the one in the right we used the complete procedure, i.e., we initialized the process by using the method described in the appendix. Since the error for the ground truth matrix $\widetilde{{\mathbf{W}}}$ is given by $\left\|(\widetilde{{\mathbf{W}}}-{\mathbf{W}}_{r})\odot{\mathbf{M}}_{r}\right\|_{F}\simeq\sigma\sqrt{40^{2}-900}\simeq 26.5$, we see from the plots of figure 8 that both the EM and RC algorithms provide good estimates. From the right hand side plot, we conclude that the initial guess provided by the heuristic method in the appendix enables a much faster convergence (in 2 or 3 iterations) to a solution with lower error.

The plot of Fig. 9 represents the average estimation error per entry after 20 iterations of EM and RC algorithms, for noisy observations of a $24\times 24$ rank 4 matrix $\widetilde{{\mathbf{W}}}$, with $70\%$ missing data, as a function of the observation noise standard deviation. We see that the average entry estimation errors after $20$ iterations are below $10^{-8}$ for noise standard deviation ranging from $10^{-2.5}$ to $10^{2.5}$ (the mean value of the entries of the ground truth matrix $\widetilde{{\mathbf{W}}}$ is $1$). This shows that, with a simple initialization procedure, both EM and RC algorithms converge to the optimal solution, even for very noisy observations. Furthermore, we conclude that the main impact of the observation noise is on the EM and RC convergence speeds—the slightly higher average error values on the right region of the plot of Fig. 9 indicates that the estimates were still converging to the optimal solution after $20$ iterations.

A relevant issue is the robustness of the rank reduction algorithms to the missing data. Our experience showed that the structure of the binary mask matrix ${\mathbf{M}}_{r}$ representing the known data is more important than the overall percentage of missing entries of ${\mathbf{W}}_{r}$. When recovering 3D structure from video, feature points enter and leave the scene in a continuous way, thus the typical structure of ${\mathbf{M}}_{r}$ is as represented in the left plot of Fig. 1, or, when considering re-appearing features, as represented in the right plot of the same Figure. We tested the EM and RC algorithms with several mask matrices ${\mathbf{M}}_{r}$ with this typical structure, with the heuristic initial guess. In these experiments, the algorithms converged always to the optimal solution in a very small number of iterations, independently of the percentage of missing data.

Similarly to the impact of the noise discussed above, the percentage of missing entries has impact on the convergence speed. In fact, we concluded that the larger is the portion of the matrix that is observed, i.e., the smaller is the percentage of missing data, the lower is the estimation error after a fixed number of iterations, i.e., the faster is the convergence of the iterative algorithms.

As referred above, the EM and RC algorithms converge in a very small number of iterations when initialized by the heuristic procedure described in the appendix. We now report an experimental evaluation of the computational costs of each iteration of EM and RC as functions of the observation matrix dimension.

We used $N\times 24$ observation matrices with missing data corresponding to a $(N-4)\times 20$ submatrix. The plots in Fig. 10 represent the number of MatLab^© floating point operations (FLOPS) and the computation time per iteration, as functions of $N$. From the left plot, we see that the number of FLOPS per iteration of the EM algorithm is larger than one of the RC algorithm. Furthermore, the FLOPS count for EM increases exponentially with $N$ (due to the SVD computation) while for RC it increases linearly with $N$. Thus, although the computation times in the right plot of Fig. 10 are smaller for EM than for RC (the reason being the very efficient MatLab^© implementation of the SVD), we conclude that RC is computationally simpler than EM. RC is even as simple as the methods to deal with complete matrices, since the most efficient way to compute the SVD is to use the power method [16] of which RC is a simple generalization.

We used a real-life video clip available at [22]. This clip shows a rotating ping-pong ball with painted dots. The left image of Fig. 11 shows the first of the $52$ video frames of the ball sequence. Due to the rotating motion, the ball occludes itself. However, it does not complete a turn, thus there are not re-appearing regions.

We used simple correlation techniques to track a set of $64$ feature points (these techniques are adequate for the short baseline video sequences we used but more sophisticated methods to cope with wide baselines are available in the literature [33]). Due to the camera-ball 3D rotation, the region of the ball that is visible changes along time, leading to an observation matrix with $\simeq 41\%$ missing entries. We used the RC algorithm to estimate the rank 4 matrix that best matched the observation and computed the 3-D structure by proceeding as in the factorization method of [38].

In the right image of Fig. 11, we represent the recovered 3-D shape, which shows that our method succeeded in recovering the spherical surface of the ball.

We now illustrate our method to recover complete 3-D models from video by minimizing a global cost function as described in section IV. We processed the observation matrix that lead to the left plot of Fig. 2. The 3-D shape recovered by our global method is shown in the right plot, where each of the pairs of re-appearing features have been correctly detected as representing the same 3-D point.

Fig. 12 describes an experiment with an object described by a larger number of 3-D feature points. We synthesized noisy versions of the 2-D trajectories of $372$ feature points located on the 3-D surface of a rotating cylinder. Then, we simulated occlusion and inclusion by removing significant segments of those trajectories. The left plot of Fig. 12 shows one of the $50$ synthesized frames. The small circles denote the noiseless projections and the points denote their noisy version, i.e., the data that is observed. Note that only an incomplete view portion of the cylinder is observed in each frame.

The data from the cylinder sequence was then collected on a $100\times 372$ observation matrix ${\mathbf{W}}$ with 9537 unknown entries ($\simeq 26\%$ of the total number). We used the global factorization algorithm to find the re-arranged matrix ${\mathbf{W}}_{R}$, i.e., to recover simplest 3-D rigid structure that could have generated the data. The right plot of Fig. 12 plots the final estimate of the 3D shape. We see that the complete cylinder is recovered. Due to the incorporation of the rigidity constraint, the 3D positions of the features points are accurately estimated, even in the presence of very noisy observations (compare the plots in Fig. 12).

We quantified the gain of using our method to recover SFM by measuring the 3-D reconstruction error. We synthesized noisy trajectories of 40 features randomly located in 3-D, 15 of them being “interrupted” to simulate object self-occlusion, leading to a $60\!\times\!55$ observation matrix ${\mathbf{W}}$ with 53.9% known entries ($x-$ and $y-$ coordinates in $[0,120]\!\times\![0,160]$, noise standard deviation $\sigma\!=\!3$). Our method generated a $60\times 40$ re-arranged matrix ${\mathbf{W}}_{R}$ with 74.2% known entries. By using ${\mathbf{W}}_{R}$ rather then ${\mathbf{W}}$ to recover SFM, we reduced the 3-D shape estimation error by approximately 50% and the 3-D motion error by approximately 70%.

In order to evaluate the sensitivity of our method to the observation noise, we plotted in Fig. 13 the probabilities of correctly detecting re-appearing features and the probability of incorrect detections (false alarms) as functions of the noise level, for the experimental setup described above. These probabilities were estimated from 100 runs for each level of noise. Although re-appearing features become harder to detect as the noise level increases, our method correctly detects more than 90% of them when the noise level is below $\sigma=5$. The plot of Fig. 13 also shows that our method has approximately zero false alarms for noise levels below $\sigma=7$. Note that the $x-$ and $y-$ coordinates are in $[0,120]\!\times\![0,160]$, so this level of performance seems to indicate that our method is able to cope with the noise in trajectories likely to happen when processing the majority of real videos.

This video—see two representative frames on the left and middle images of Fig. 14—shows a Rubik’s cube. The video sequence was obtained by rotating a hand-held camera around a table with the cube. In the leftmost image of Fig. 14, we superimposed with the video frame the visible features and the initial parts of their trajectories. Due to the cube self-occlusion, feature points enter and leave the scene. In particular, along the video sequence, all features disappear and several of them re-appear because the camera performs a complete turn around the cube.

To emphasize the advantages of using the algorithms described in this paper, we first applied to a segment of the Rubik’s cube video clip the original factorization method of [38] for complete matrices, obtaining the 3D shape represented on the left side of Fig. 15. This model was obtained with $28$ features and $18$ frames.

We then collected the entire set of the visible parts of the trajectories of $64$ features across $85$ frames in a $170\!\times\!64$ incomplete observation matrix ${\mathbf{W}}$. The structure of the missing part of ${\mathbf{W}}$ is coded by the $170\times 64$ binary mask ${\mathbf{M}}$ represented on the right side of Figure 14. The number of missing entries in ${\mathbf{W}}$ was about $62\%$. We computed the re-arranged matrix ${\mathbf{W}}_{R}$ by using our global factorization. The right image of Fig. 15 shows the texture-mapped 3-D shape recovered from ${\mathbf{W}}_{R}$.

This example illustrates how our global approach trivializes the usually hard task of merging partial estimates of a 3D model like the one in the left image of Fig. 15. Note that the top face is missing also in the right image of Fig. 15 because the position the cube model is shown in that image was not seen in the video sequence. The advantage of using our approach to recover 3-D rigid structure is two-fold. First, while recovering a complete 3-D model by fusing partial models as the one on the left side of Fig. 15 is a complex task, our method recovers directly the complete model shown on the right side of Fig. 15. Second, rather than processing subsets of the set of features and frames at disjoint steps, our method uses all the information available in a global way, leading to more accurate 3-D shapes as illustrated by the 3-D models in Fig. 15.

The 3D models recovered by our method can be used to represent in an efficient way the original video sequence as proposed in [3]—the video sequence is represented by the 3-D shape, texture, and 3-D motion of the objects. This leads to significant bandwidth saving since once the 3-D shape and texture of the objects have been transmitted, their 3-D motion is transmitted with a few bytes per frame.

We used this methodology to compress the entire sequence of $2161$ frames of the Rubik’s cube video clip. The compression ratio (relative to the original J-PEG compressed frames) was approximately $10^{3}$. Fig. 16 shows sample original frames (top row) and the corresponding compressed frames (bottom row). The differences of lighting between the top and bottom images are due to fact that the constancy of the texture of the 3-D model does not model the different light reflections that happen in the real scenario.

Before addressing the general case, we consider the simpler case where a number of columns of ${\mathbf{W}}_{r}$ are entirely known, i.e., a number of feature points are present in all frames, and a number of rows of ${\mathbf{W}}_{r}$ are entirely known, i.e., and a number of frames contain all feature projections. We collect those known columns in a submatrix ${\mathbf{W}}_{A}$ and those known rows in a submatrix ${\mathbf{W}}_{B}$. From the data in ${\mathbf{W}}_{A}$, we estimate the column space matrix ${\mathbf{A}}$ as

From ${\mathbf{W}}_{B}$ and the column space ${\mathbf{A}}$, the LS estimate of the row space matrix ${\mathbf{B}}$ is, see [16],

where ${\mathbf{A}}_{B}$ collects the rows of ${\mathbf{A}}$ that correspond to the rows of ${\mathbf{W}}_{B}$. We see that the matrix ${\mathbf{A}}_{B}^{T}{\mathbf{A}}_{B}$ must be nonsingular so the matrix ${\mathbf{W}}_{A}$ must have at least 4 linearly independent columns and the matrix ${\mathbf{W}}_{B}$ must have at least 4 linearly independent rows.

In general, it may be impossible to find 4 entire columns and 4 entire rows without missing elements in the matrix ${\mathbf{W}}_{r}$. We must then estimate the column and row space matrices ${\mathbf{A}}$ and ${\mathbf{B}}$ by combining the spaces that correspond to smaller submatrices of ${\mathbf{W}}_{r}$. We describe the algorithm for combining two column/row space matrices. The process is then repeated until the entire matrix ${\mathbf{W}}_{r}$ has been processed.

Select from ${\mathbf{W}}$ two submatrices ${\mathbf{W}}_{1}$ and ${\mathbf{W}}_{2}$ that have at least 4 columns and 4 rows without missing elements. We factorize ${\mathbf{W}}_{1}$ and ${\mathbf{W}}_{2}$ using (27) to obtain the corresponding column space matrices ${\mathbf{A}}_{1}$ and ${\mathbf{A}}_{2}$. If the observation matrix ${\mathbf{W}}$ is in fact well modelled by a rank 4 matrix, the submatrices of ${\mathbf{A}}_{1}$ and ${\mathbf{A}}_{2}$ that correspond to common rows in ${\mathbf{W}}$, denoted by ${\mathbf{A}}_{12}$ and ${\mathbf{A}}_{21}$, are related by a linear transformation,

We compute ${\mathbf{N}}$ as the LS solution of (29) and assemble a combined column space matrix ${\mathbf{A}}$ for the rows corresponding to ${\mathbf{W}}_{1}$ and ${\mathbf{W}}_{2}$:

where ${\mathbf{A}}_{2\setminus 1}$ denotes the submatrix of ${\mathbf{A}}_{2}$ that collects the rows that do not correspond to rows of ${\mathbf{W}}_{1}$. We compute the combined row space matrix ${\mathbf{B}}$ by using an analogous procedure.