In Perfect Shape: Certifiably Optimal 3D Shape Reconstruction from 2D Landmarks
Heng Yang, Luca Carlone
Introduction
3D object detection and pose estimation from a single image is a fundamental problem in computer vision. Despite the progress in semantic segmentation , depth estimation , and pose estimation , reconstructing the 3D shape and pose of an object from a single image remains a challenging task .
A typical approach for 3D shape reconstruction is to first detect 2D landmarks in a single image, and then solve a model-based optimization to lift the 2D landmarks to form a 3D model . For the optimization to be well-posed, the unknown shape is assumed to be a 3D deformable model, composed by a linear combination of basis shapes, handcrafted or learned from a large corpus of training data . The optimization then seeks to jointly optimize the coefficients of the linear combination (shape parameters) and the camera pose to minimize the reprojection errors between the 3D model and the 2D landmarks. This model-based paradigm has been successful in several applications such as face recognition , car model fitting , and human pose estimation .
Despite its long history and broad range of applications, there is still no globally optimal solver for the non-convex optimization problem arising in 3D shape reconstruction. Therefore, most existing solutions adopt a local optimization strategy, which alternates between solving for the camera pose and the shape parameters. These techniques, as shown in prior works , require an initial guess for the solution and often get stuck in local minima. In addition, 2D landmark detectors are prone to produce outliers, causing existing methods to be brittle . Therefore, the motivation for this paper is two-fold: (i) to develop a certifiably optimal shape reconstruction solver, and (ii) to develop a robust reconstruction algorithm that is insensitive to a large amount of outlier 2D measurements (e.g., ).
Contributions. Our first contribution is to formulate the shape reconstruction problem as a polynomial optimization problem and apply Lasserre’s hierarchy of Sums-of-Squares (SOS) relaxations to relax the non-convex polynomial optimization into a convex semidefinite program (SDP). We show the SOS relaxation of minimum order 2 empirically solves the non-convex shape reconstruction problem exactly and provides a global optimality certificate. The second contribution is to apply basis reduction, a technique that exploits the sparse structure of the polynomial in the objective function, to reduce the size of the resulting SDP. We show that basis reduction significantly improves the efficiency of the SOS relaxation without compromising global optimality. To the best of our knowledge, this is the first certifiably optimal solver for shape reconstruction, and we name it Shape⋆. Our third contribution is to robustify Shape⋆ by adopting a truncated least squares (TLS) robust cost function and solving the resulting robust estimation problem using graduated non-convexity . The resulting algorithm, named Shape, is robust against outliers and does not require an initial guess.
The rest of this paper is organized as follows. Section 2 reviews related work. Section 3 introduces notation and preliminaries on SOS relaxations. Section 4 introduces the shape reconstruction problem. Section 5 develops our SOS solver (Shape⋆). Section 6 presents an algorithm (Shape) to robustify the SOS relaxation against outliers. Section 7 provides experimental results in both simulations and real datasets, while Section 8 concludes the paper.
Related Work
We limit our review to optimization-based approaches for 3D shape reconstruction from 2D landmarks. The interested reader can find a review of end-to-end shape and pose reconstruction using deep learning in .
Notation and Preliminaries
We now give a brief summary of SOS relaxations for polynomial optimization. Our review is based on . We first introduce the notion of SOS polynomial.
and is called the Gram matrix of .
Now consider the following polynomial optimization:
the ideal and the -th truncated ideal of , where is the degree of a polynomial. The ideal is simply a summation of polynomials with polynomial coefficients, a construct that will simplify the notation later on. We call
the quadratic module and the -th truncated quadratic module generated from . Note that the quadratic module is similar to the ideal, except now we require the polynomial coefficients to be SOS. Apparently, if , then is nonnegative on If , then , with and . For any , since , so ; since and , so . Therefore, . Putinar’s Positivstellensatz describes when the reverse is also true.
Based on Putinar’s Positivstellensatz, Lasserre derived a sequence of SOS relaxations that approximates the global minimum of problem (3) with increasing accuracy. The key insight behind Lasserre’s hierarchy is twofold. The first insight is that problem (3), which we can write succinctly as , can be equivalently written as (intuition: the latter pushes the lower bound to reach the global minimum of ). The second intuition is that we can rewrite the condition , using Putinar’s Positivstellensatz (Theorem 2), leading to the following hierarchy of Sums-of-Squares relaxations.
Lasserre’s hierarchy of order is the following SOS program:
which can be written as a standard SDP. Moreover, let be the global minimum of (3) and be the optimal value of (8), then monotonically increases and when . More recently, Nie proved that under Archimedeanness, Lasserre’s hierarchy has finite convergence generically (i.e., for some finite ).
In computer vision, Lasserre’s hierarchy was first used by Kahl and Henrion to minimize rational functions arising in geometric reconstruction problems, and more recently by Probst et al. as a framework to solve a set of 3D vision problems. In this paper we will show that the SOS relaxation as written in eq. (8) allows using basis reduction to exploit the sparsity pattern of polynomials and leads to significantly smaller semidefinite programs.
Problem Statement: Shape Reconstruction
Certifiably Optimal Shape Reconstruction
This section shows how to develop a certifiably optimal solver for problem (13). Our first step is to algebraically eliminate the translation and obtain a translation-free shape reconstruction problem, as shown below.
The shape reconstruction problem (13) is equivalent to the following translation-free optimization:
Further, let and be the global minimizer of the above translation-free optimization (14), then the optimal translation can be recovered as:
A formal proof of Theorem 4 can be found in the Supplementary Material. The intuition behind Theorem 4 is that if we express the landmark coordinates and 3D basis shapes with respect to their (weighted) centroids and , we can remove the dependence on the translation . This strategy is inspired by Horn’s method for point cloud registration , and generalizes to the weighted and non-centered case.
This section applies Lasserre’s hierarchy as described in Theorem 3 to solve the translation-free shape reconstruction problem (14). We do this in two steps: we first show problem (14) can be formulated as a polynomial optimization in the form (3); and then we add valid constraints to make the feasible set Archimedean.
then it becomes clear that is a polynomial function of with degree 4. Because the Lasso regularization is linear in , the objective function is a degree-4 polynomial.
In summary, the translation-free problem (14) is equivalent to a polynomial optimization with a degree-4 objective , constrained by 15 quadratic equalities (eq. (19)) and linear inequalities .
Now we can certify the Archimedeanness of :
with and (cf. Theorem 2).
Apply Lasserre’s Hierarchy. With Archimedeanness, we can now apply Lasserre’s hierarchy of SOS relaxations.
The SOS relaxation of order ()The minimum relaxation order is 2 because has degree 4. for the translation-free shape reconstruction problem (14) is the following convex semidefinite program:
where is the objective function defined in (14), are the inequality constraints , are the equality constraints defined in (19), and , , are the sizes of matrices and vectors.
While a formal proof of Proposition 6 is given in the Supplementary Material, we observe that (22) immediately results from the application of Lasserre’s hierarchy to (8), by parametrizing with monomial bases , and PSD matrices , (one for each ), and by parametrizing with monomial basis and coefficient vectors (one for each ). Problem (22) can be written as an SDP and solved globally using standard convex solvers (e.g. YALMIP ). We call the SDP written in (22) the primal SDP. The dual SDP of (22) can be derived using moment relaxation , which is readily available in GloptiPoly 3 .
Extract Solutions from SDP. After solving the SDP (22), we can extract solutions to the original non-convex problem (14), a procedure we call rounding.
Let and be the optimal solutions to the SDP (22) at order ; compute as the eigenvector corresponding to the minimum eigenvalue of , and then normalize such that the first entry is equal to 1. Then an approximate solution to problem (14) can be obtained as:
where is the true (unknown) global minimum of problem (14). We define the relative duality gap as:
which quantifies the quality of the SOS relaxation.
Certifiable Global Optimality. Besides extracting solutions to the original problem, we can also verify when the SOS relaxation solves the original problem exactly.
Let and be the optimal solutions to the SDP (22) at order . If (the corank is the dimension of the null space of ), then is the global minimum of problem (14), and the relaxation is said to be tight at order . Moreover, the relative duality gap and the solution extracted using Proposition 7 is the unique global minimizer of problem (14).
The proof of Theorem 8 is given in the Supplementary Material. Empirically (Section 7), we observed that the relaxation is always tight at the minimum relaxation order . Note that even when the relaxation is not tight, one can still obtain an approximate solution using Proposition 7 and quantify how suboptimal the approximate solution is using the relative duality gap .
2 Basis Reduction
Despite the theoretical soundness and finite convergence at order , the size of the SDP (22) is , which for becomes , implying that the size of the SDP grows quadratically in the number of bases . Although there have been promising advances in improving the scalability of SDP solvers (see for a thorough review), such as exploiting sparsity and low-rankness , in this section we demonstrate a simple yet effective approach, called basis reduction, that exploits the structure of the objective function to significantly reduce the size of the SDP in (22).
In a nutshell, basis reduction methods seek to find a smaller, but still expressive enough, subset of the full vector of monomials on the right-hand side (RHS) of eq. (23), to explain the objective function on the left-hand side (LHS). There exist standard approximation algorithms for basis reduction, discussed in and implemented in YALMIP . However, in practice we found the basis selection method in YALMIP failed to find any reduction for the SDP (22). Therefore, here we propose a problem-specific reduction, which follows from the examination of which monomials appear on the LHS of (23).
The SOS relaxation of order with basis reduction for the translation-free shape reconstruction problem (14) is the following convex semidefinite program:
Comparing the SDP (27) and (22), the most significant change is replacing the full monomial basis in (22) with a much smaller monomial basis that excludes degree-2 monomials purely supported in and . This reduction is motivated by analyzing the monomial terms in . Although a formal proof of the equivalence between (22) and (27) remains open, we provide an intuitive explanation in the Supplementary Material. After basis reduction, the size of the SDP (27) is , which is linear in and much smaller than the size of the original SDP (22) For , , while .. Section 7 numerically shows that the SDP after basis reduction gives the same (tight) solution as the original SDP.
3 Shape⋆: Algorithm Summary
To summarize the derivation in this section, our solver for the shape reconstruction problem (13), named Shape⋆, works as follows. It first solves the SDP (27) and applies the rounding described in Proposition 7 to compute an estimate of the shape parameters and rotation and possibly certify its optimality. Then, Shape⋆ uses the closed-form expression (17) to retrieve the translation estimate .
Robust Outlier Rejection
Section 5 proposed a certifiably optimal solver for problem (13). However, the least squares formulation (13) tends to be sensitive to outliers: the pixel measurements in eq. (9) are typically produced by learning-based or handcrafted detectors , which might produce largely incorrect measurements (e.g. due to wrong data association ), which in turn leads to poor shape reconstruction results. This section shows how to regain robustness by iteratively solving the weighted least squares problem (13) and adjusting the weights to reject outliers.
The key insight is to substitute the least square penalty in (13) with a robust cost function, namely the truncated least squares (TLS) cost . Hence, we propose the following TLS shape reconstruction formulation:
where (introduced for notational convenience), and implements a truncated least squares cost, which is quadratic for small residuals and saturates to a constant value for residuals larger than a maximum error .
Our second insight is that can be written as , by introducing extra slack binary variables . Therefore, we can write problem (29) equivalently as:
The final insight is that now we can minimize (30) by iteratively minimizing (i) over (with fixed weights ), and (ii) over the weights (with fixed ). The rationale for this approach is that step (i) can be implemented using Shape⋆ (since in this case the weights are fixed), and step (ii) can be implemented in closed-form. To improve convergence of this iterative algorithm, we adopt graduated non-convexity , which starts with a convex approximation of problem (30) and uses a control parameter to gradually increase the amount of non-convexity, till (for large ) one solves (30). The resulting algorithm named Shape is given in Algorithm 1. We refer the reader to the Supplementary Material and for a complete derivation of Algorithm 1 and for the closed-form expression of the weight update in line 1 of the algorithm.
Shape is deterministic and does not require an initial guess. We remark that the graduated non-convexity scheme in Shape (contrarily to Shape⋆) is not guaranteed to converge to an optimal solution of (30), but we show in the next section that it is empirically robust to outliers.
Experiments
Implementation details. Both Shape⋆ and Shape are implemented in Matlab, with both SOS relaxations (22) and (27) implemented using YALMIP and the resulting SDPs solved using MOSEK .
2 Shape⋆ for Outlier-Free Reconstruction
3 Shape##{\#} for Robust Reconstruction
Conclusions
We presented Shape⋆, the first certifiably optimal solver for 3D shape reconstruction from 2D landmarks in a single image. Shape⋆ is developed by applying Lasserre’s hierarchy of SOS relaxations combined with basis reduction to improve efficiency. Experimental results show that the SOS relaxation of order 2 always achieves global optimality. To handle outlying measurements, we also proposed Shape, which solves a truncated least squares robust estimation problem by iteratively running Shape⋆ without the need for an initial guess. We show that Shape achieves robustness against outliers on the FG3DCar dataset and outperforms state-of-the-art solvers.
Proof of Theorem 4
Here we prove Theorem 4 in the main document. Recall the weighted least squares optimization for shape reconstruction in eq. (13) and denote its objective function as , with :
In order to marginalize out the translation , we compute the derivative of w.r.t. :
and set it to , which allows us to write in closed form using and :
with and being the weighted centers of the 2D landmarks and the 3D basis shapes :
Then we can substitute the expression of in (A3) back into the objective function in (A1) and obtain an objective function without translation:
we can see the equivalence between the objective function in eq. (A5) and the objective function in eq. (14) of Theorem 4. The constraints remain unchanged because we only marginalize out the unconstrained variable . Therefore, the shape reconstruction problem (13) is equivalent to the translation-free problem (14), and the optimal translation can be recovered using eq. (A3). ∎
Proof of Proposition 6
Here we prove the SOS relaxation of order () for the translation-free shape reconstruction problem (14) is the semidefinite program in (22). First, let us rewrite the general form of Lasserre’s hierarchy of order in eq. (8) in Theorem 3 as the following:
In words, the constraints of (A8) ask the polynomial to be written as a sum of two polynomials and , with in the -th truncated ideal of , and in the -th truncated quadratic module of .
Next, we use the definition of the -th truncated ideal and the -th truncated quadratic module to explicitly represent and . First recall the definition of the -th truncated ideal in eq. (5), which states that must be written as a sum of polynomial products between the equality constraints ’s and the polynomial multipliers ’s:
with being the vector of unknown coefficients associated with the monomial basis . The size of is equal to the length of , which can be computed by , with being the number of variables in , and being the maximum degree of the monomial basis. Similarly for , we recall the definition of the -th truncated quadratic module in eq. (7), which states that must be written as a sum of polynomial products between the inequality constraints ’s and the SOS polynomial multipliers ’s:
and the degree of each polynomial product must be no greater than , i.e., . For our specific shape reconstruction problem, we have , , and . Since has degree 0, can have degree up to . All , have degree 1, so , can have degree up to . However, because SOS polynomials can only have even degree, can only have degree up to . For , they have degree 2, so their corresponding SOS polynomial multipliers can have degree up to . Now for each SOS polynomial , from the Gram matrix representation in eq. (2), we can associate a PSD matrix with it using corresponding monomial bases:
Finally, by inserting the expressions of in (A12) back to the expression of in (A11), and inserting the expression of in (A10) back to the expression of in (A9), we can convert the SOS relaxation of general form (A8) to the semidefinite program (22). ∎
Proof of Theorem 8
According to , the dual SDP of (22) is the following SDP:
and is in the null-space of . Therefore, the solution extracted using Proposition 7 is also the unique global minimizer of problem (14). ∎
Derivation of Proposition 9
Here we show the intuition for using the basis reduction in Proposition 9. In the original SOS relaxation (22), the parametrization of the SOS polynomial multipliers , and the polynomial multipliers , uses the vector of all monomials up to their corresponding degrees (cf. (A10) and (A12)), which leads to an SDP of size that grows quadratically with the number of basis shapes . In basis reduction, we do not limit ourselves to the vector of full monomials, but rather parametrize , and with unknown monomials bases , and , which allows us to rewrite (23) as:
with the hope that , and have much smaller sizes (we limit ourselves to the case of , at which level the relaxation is empirically tight).
As described, one can see that the problem of finding smaller , and , while keeping the relaxation empirically tight, is highly combinatorial in general. Therefore, our strategy is to only consider the following case:
Expressive: choose such that contains all the monomials in ,
Balanced: choose and such that the sum can only have monomials from .
In words, condition (i) ensures that the right-hand side (RHS) of (A20) contains all the monomials of the left-hand side (LHS). Condition (ii) asks the three terms of the RHS, i.e., , and , to be self-balanced in the types of monomials. For example, if contains extra monomials that are not in the LHS, then those extra monomials better appear also in and/or so that they could be canceled by summation. Under these two conditions, it is possible to have equation (A20) holdWhether or not these are sufficient or necessary conditions remains open. However, leveraging Theorem 8 we can still check optimality a posteriori..
The choices in both conditions depend on analyzing the monomials in . Recall the expression of in (14) and the expression in (18) for each term inside the summation, it can be seen that only contains the following types of monomials:
and the key observation is that does not contain degree-4 monomials purely in or , i.e., and , or any degree-3 monomials in and . Therefore, when choosing , we can exclude degree-2 monomials purely in and from , and set A more rigorous analysis should follow the rules of Newton Polytope , but the intuition is the same as what we describe here. as stated in Proposition 9. This will satisfy the expressive condition (i), because can have the following monomials:
and those in (A25) cover the monomials in . Replacing with is the key step in reducing the size of the SDP, because it reduces the size of the SDP from to , i.e., from quadratic to linear in .
In order to satisfy condition (ii), when choosing and , the goal is to have the product between , and , result in monomials that appear in , and ensure that monomials that do not appear in the latter can simplify our in the summation. For example, as stated in Proposition 9, we choose and will contain monomials 1, and . Because ’s have monomials , and , we can see that will contain the following monomials:
This still satisfies the balanced condition, because monomials of in (A27) balance with monomials of in (A25), and monomials of in (A27) balance with monomials of in (A26). Similarly, choosing makes have monomials , and , and because ’s have monomials , and , we see that contains the following monomials:
which balance with monomials in from (A25) and (A26).
We remark that we cannot guarantee that the SOS relaxation resulting from basis reduction can achieve the same performance as the original SOS relaxation and we cannot guarantee our choice of basis is “optimal” in any sense. Therefore, in practice, one needs to check the solution and compute and to check the optimality of the solution produced by (27). Moreover, it remains an open problem to find a better set of monomials bases to achieve better reduction (e.g., knowing more about the algebraic geometry of and could possibly enable using the standard monomials as a set of bases ).
Derivation of Algorithm 1
For a complete discussion of graduated non-convexity and its applications for robust spatial perception, please see .
In the main document, for robust shape reconstruction, we adopt the TLS shape reconstruction formulation:
where is called the residual, and implements a truncated least squares cost. Recalling that , we can rewrite the TLS shape reconstruction as a joint optimization of and the binary variables ’s, as in eq. (30) in the main document. However, as hinted in the main document, due to the non-convexity of the TLS cost, directly solving the joint problem or alternating between solving for and binary variables ’s would require an initial guess and is prone to bad local optima.
The idea of graduated non-convexity (GNC) is to introduce a surrogate function , governed by a control parameter , such that changing allows to start from a convex proxy of , and gradually increase the amount of non-convexity till the original TLS function is recovered. The surrogate function for TLS is stated below.
The truncated least squares function is defined as:
where is a given truncation threshold. The GNC surrogate function with control parameter is:
By inspection, one can verify is convex for approaching zero (() and retrieves in (A32) for . An illustration of is given in Fig. A1.
The nice property of the GNC surrogate function is that when is close to zero, is convex, which means the only non-convexity of problem (A31) comes from the constraints and can be relaxed using the SOS relaxations.
For the GNC surrogate function , the simple trick of introducing binary variables () would not work. However, Black and Rangarajan showed that this idea of introducing an outlier variable can be thought of an outlier variable: when , the measurement is an inlier, when , the measurement is an outlier. can be generalized to many robust cost functions. In particular, for the GNC surrogate function, we have the following.
The GNC surrogate TLS shape reconstruction:
with defined in (A33), is equivalent to the following optimization with outlier variables ’s:
where is the following outlier process:
The derivation of in (A36) follows the Black-Rangarajan procedure in Fig. 10 of . ∎
In words, the Black-Rangarajan duality allows us to rewrite the non-convex shape reconstruction problem as a joint optimization in and outlier variables ’s. The interested readers can find closed-form outlier processes for many other robust cost functions in the original paper .
Leveraging the Black-Rangarajan duality, for any given choice of the control parameter , we can solve problem (A35) in two steps: first we solve using Shape⋆ with fixed weights ’s, and then we update the weights with fixed . In particular, at each iteration (corresponding to a given control parameter ), we perform the following:
Variable update: minimize (A35) with respect to , with fixed weights :
where we have dropped the term because it is independent from . This problem is exactly the weighted least squares problem (13) and can be solved using Shape⋆ (cf. line 1 in Algorithm 1). Using the solutions , we can compute the residuals (cf. line 1 in Algorithm 1).
Weight update: minimize (A35) with respect to , with fixed residuals :
where we have dropped because it is a constant for the optimization. This optimization, fortunately, can be solved in closed-form. We take the gradient of the objective function with respect to :
and observe that when , and when . Therefore, the global minimizer is:
and this is the weight update rule in line 1 of Algorithm 1.
After both the variables and weights are updated using the shaperobust approach described above, we increase the control parameter to increase the non-convexity of the surrogate function (cf. line 1 of Algorithm 1). At the next iteration , the updated weights are used to perform the variable update. The iterations terminate when the change in the objective function becomes negligible (cf. line 1 of Algorithm 1) or after a maximum number of iterations (cf. line 1 of Algorithm 1). Note that all weights are initialized to (cf. line 1 in Algorithm 1), which means that initially all measurements are tentatively accepted as inliers, therefore no prior information about inlier/outlier is required.
FG3DCar Qualitative Results
Fig. 14 shows 9 full qualitative results comparing the performances of Altern+Robust , Convex+Robust and Shape on the FG3DCar dataset under to outlier rates. One can further see that the performance of Shape is sensitive to outliers, while the performances of Altern+Robust and Convex+Robust gradually degrade and fail at to outliers.
& Convex+Robust Shape Honda CRV Altern+Robust Convex+Robust Shape Mercedes-Benz GL450
& Convex+Robust Shape Volvo V70 Altern+Robust Convex+Robust Shape Saab 93