Rigid and Articulated Point Registration with Expectation Conditional Maximization

In plain words

Registration means moving one set of 3-D points so that it lines up with another. Here one set is a model, such as the points on a hand, and the other is a set of noisy measurements in which some points are junk. The classic ICP method pairs each measurement with its nearest model point and repeats. This paper instead treats each model point as the centre of a soft cloud of uncertainty, and each cloud may be stretched in its own direction. Every measurement is considered to belong to one of those clouds, or to a background of random junk.

The method alternates two steps. First it estimates how likely each measurement is to belong to each cloud. Then it updates the pose and the shape of the clouds, using those soft weights. A rigid object is aligned in one go; its rotation is found with a convex relaxation, because the exact problem is hard. A jointed object such as a hand is aligned one part at a time: the first part takes the measurements it explains, those are removed, and the next part is fitted to what remains.

The authors report lower errors than ICP-style baselines on their simulated trials and on one stereo example. Their hand tracking is shown qualitatively, with no accuracy numbers.

TL;DR

ECMPR is an ECM variant of EM for registering model points to noisy, outlier-laden observations. Each model point is a Gaussian centre with an anisotropic covariance, and a uniform class absorbs outliers. The rotation comes from an SDP relaxation, and articulated objects are fitted one rigid part at a time. The paper reports lower errors than TriICP and ICP on its own examples. No official code was found. The outlier sphere radius r and the regularization ε are not stated.

Problem

• Matching two 3-D point sets needs point-to-point correspondences, and these are unknown §1.
• ICP alternates hard nearest-point assignments with a transform estimate. It is fast but can be trapped in local minima and is sensitive to the initial pose and to the match threshold §1.
• Trimmed and RANSAC-style variants must examine many random subsets before a good one is found, and articulated objects need one trimmed subset per part §1.
• Gaussian mixtures with one isotropic covariance give a closed-form rotation, but they cannot represent anisotropic noise §1, §5.1. With full anisotropic covariances, the rotation subproblem has no closed form §1, §5.2.
• Articulated methods that fit parts separately, or that enforce joint constraints afterwards, only satisfy those constraints approximately §1. The paper's rigid links are exact by construction §1, §6.1.
• Assumed setting: model points lie on one rigid object, or on the parts of one open kinematic chain. Observations are the inliers, which are noisy copies of the model, plus outliers drawn uniformly. A partition of the observations is not needed in advance §2.1, §6.1.
• The number of outliers does not need to be known in advance §7.2.

Method

Model and notation §2.1, §3

1. Model points \mathcal{X} = \{X_i\}_{i=1}^{n} and observations \mathcal{Y} = \{Y_j\}_{j=1}^{m}, in 3-D. A rigid pose maps a model point with \mu(X_i;\Theta) = R X_i + t, where \Theta = \{R, t\} eq. 1, §2.1.
2. Each model point is the centre of one Gaussian component with covariance \Sigma_i. Observation Y_j is drawn from component i or from an outlier class n+1 with uniform density 1/V over the working volume V eq. 4, eq. 5, §3.
3. Each observation has a hidden label Z_j \in \{1,\dots,n+1\}. The priors are p_{in} = v/V for each inlier component and p_{out} = (V - nv)/V for the outlier class, where v = 4\pi r^3/3 is the volume of a small sphere of radius r around each model point. The paper assumes nv \ll V eq. 11, §4.

E-step §4, eq. 12–14

4. For each observation and each model point, compute the inlier posterior. The Gaussian normalization is folded into the outlier constant, which is the \emptyset_{3D} term below.

\alpha_{ji} = \frac{|\Sigma_i|^{-1/2}\exp\!\left(-\tfrac12\|Y_j-\mu(X_i;\Theta)\|^2_{\Sigma_i}\right)} {\sum_{k=1}^{n}|\Sigma_k|^{-1/2}\exp\!\left(-\tfrac12\|Y_j-\mu(X_k;\Theta)\|^2_{\Sigma_k}\right)+\emptyset_{3D}} eq. 12

\emptyset_{3D} = 1.5\sqrt{2\pi}\,r^{-3} eq. 13

5. The outlier posterior is one minus the sum of the inlier posteriors eq. 14, §4.

CM-step 1: pose §4, eq. 19–24

6. Choose the pose that minimizes the weighted residuals, with the covariances held fixed at their current values eq. 19, §4.
7. Equivalently, collapse the observations into one virtual observation per model point. Its weight \lambda_i is the sum of the posteriors for that point, and its position W_i is their weighted mean eq. 22, eq. 23, §4. The pose then minimizes \tfrac12 \sum_i \lambda_i \|W_i - \mu(X_i;\Theta)\|^2_{\Sigma_i} eq. 24, §4.

CM-step 2: covariances §4, eq. 20–21

8. With the new pose, update each covariance from the posterior-weighted outer products of the residuals eq. 20, §4.
9. Alternatively use one covariance shared by all components eq. 21, §4. The paper prefers the shared form when the number of data points is small §4. Its Figure 2 example uses the shared form fig. 2.
10. Before use, add \varepsilon I to each covariance. The paper argues that this keeps the shape of the ellipsoid, and that it prevents collapse onto a single point §4.

Rigid pose §5

11. Given the rotation, the translation has a closed form that weights each virtual observation by its inverse covariance eq. 27, §5.2.
12. Isotropic case: the translation is a weighted mean of the residuals eq. 29, §5.1. The rotation is then found in closed form by one of the standard absolute-orientation methods §5.1, eq. 30.
13. Anisotropic case: the rotation is written as the vector r = \mathrm{vec}(R) and the rank-one matrix \rho = r r^\top eq. 31, §5.2. The criterion becomes a quadratic in r with six orthonormality constraints eq. 32, §5.2.
14. The rank-one condition \rho = r r^\top is nonconvex. The paper replaces it by the convex condition \rho \succeq r r^\top eq. 33, eq. 34, §5.2. The solution of this relaxed problem is used directly as the new rotation §5.2, step 3a.
15. The matrices A (9×9) and b (9×1) are given in Appendix A, in terms of the virtual observations, weights and covariances app. A.

\min_{r}\ \tfrac12\left(r^\top A\, r + 2\, b^\top r\right) \quad\text{s.t.}\quad r^\top \Delta_{kl}\, r = \delta_{kl} eq. 32

ECMPR-rigid §5.2

# ECMPR-rigid  [§5.2]
# Input:  model points X (n, 3); observations Y (m, 3); covariances (n, 3, 3)
# Output: rotation R (3, 3), translation t (3,), labels (m,)

def ecmpr_rigid(X, Y, sigmas):
    # Start from the identity rotation and the zero translation  [§5.2]
    R = identity(3)
    t = zeros(3)
    while not converged:
        # E-step: inlier posteriors plus the outlier column  [eq. 12, 14]
        alpha = e_step(X, Y, R, t, sigmas)              # (m, n+1)
        # Virtual observations and their weights  [eq. 22, 23]
        W, lam = virtual_observations(alpha, Y)         # (n, 3), (n,)
        # CM-step 1: rotation by SDP relaxation  [eq. 32-34]
        R_new = sdp_rotation(X, W, lam, sigmas)
        # CM-step 2: translation in closed form  [eq. 27]
        t_new = translation(X, W, lam, sigmas, R_new)
        # CM-step 3: covariances, per point or shared  [eq. 20 or 21]
        sigmas = covariance_update(X, Y, alpha, R_new, t_new)
        # Convergence is judged on the rotation change only  [§5.2 step 3]
        converged = squared_norm(R_new - R) < epsilon  # epsilon not stated
        R, t = R_new, t_new
    # Maximum a posteriori labels; index n means outlier  [§5.2 step 4]
    labels = argmax(alpha, axis=1)
    return R, t, labels

Articulated motion §6

16. The articulated model is a chain of transforms T_p = Q_0 Q_1 \cdots Q_p eq. 36, §6.1. The root Q_0 has a free rotation and a free translation. Each later Q_p has only a rotation eq. 37, §6.1.
17. Parts are expanded incrementally: T_p = T_{p-1} Q_p eq. 38, §6.1. Written out, the rotation of part p is R_{0,p-1} R_p and its translation is t_{0,p-1} eq. 39, §6.1.

T_p = \begin{bmatrix} R_{0,p-1} R_p & t_{0,p-1} \\ 0^\top & 1 \end{bmatrix} eq. 39

18. The root is registered as a rigid part, with the whole observation set eq. 40, §6.2.
19. For part p \ge 1, the substitutions U_p = R_{0,p-1} R_p R_{0,p-1}^\top and V_i = R_{0,p-1} X_i turn the problem into a rigid one in U_p with the translation fixed to t_{0,p-1} eq. 42, eq. 43, eq. 44, §6.2. The fit is therefore rotation-only for every part after the root.
20. After solving for U_p, recover R_p = R_{0,p-1}^\top U_p R_{0,p-1} eq. 42, §6.2.
21. After each part is fitted, its inliers are removed from the observation set before the next part. The observations are not partitioned in advance §6.2.

ECMPR-articulated §6.2

# ECMPR-articulated  [§6.2]
# Input:  parts: list of (n_p, 3) model points, root first
#         observations: (m, 3) data points
# Output: rotations R_0..R_P (relative to the parent), root translation t_0

def ecmpr_articulated(parts, observations):
    # Root part: free rigid motion over all observations  [eq. 40]
    R_0, t_0 = ecmpr_rigid(parts[0], observations)
    # Remove the root's inliers before the next part  [§6.2 step 1]
    remaining = remove_inliers(observations, parts[0], R_0, t_0)
    rotations = [R_0]
    for p in range(1, len(parts)):
        # Parent rotation R_{0,p-1}, the product of earlier rotations
        parent = product(rotations)
        # Model points expressed in the parent frame  [eq. 43]
        in_parent = parts[p] @ parent.T
        # Rotation-only fit with the translation fixed at t_0  [eq. 44]
        U_p = ecmpr_rotation_only(in_parent, remaining, t_0)
        # Undo the parent rotation to get the relative rotation  [eq. 42]
        R_p = parent.T @ U_p @ parent
        rotations.append(R_p)
        # Pose of part p is R_{0,p} with translation t_0  [eq. 39]
        remaining = remove_inliers(remaining, parts[p], parent @ R_p, t_0)
    return rotations, t_0

Hyperparameters

Results

Author-reported. Setup: simulated rigid trials with 15 model points and 25 observations, plus one real stereo pair of a walking person (223 model points, 249 data points) §7.1, fig. 5. The baseline is TriICP, and ICP for the stereo pair.

Table 1 summarizes the simulated rigid experiments. Rows 2 and 3 both add anisotropic noise and differ only in the covariance model tab. 1:

Figure 4 plots the same measures against the true rotation angle, over 1,000 trials at each angle. Its curves show the mean and the mean plus or minus one standard deviation fig. 4. The text reports that the anisotropic model beats the isotropic one on both pose error and correct assignments §7.1.

Table 2 compares ECMPR and ICP on the stereo pair tab. 2, fig. 5:

Articulated hand results are qualitative. Simulated sequences of 120 poses are tracked with 240 inliers and 72 outliers. The paper reports that the tracker recovers the poses when it starts from large covariances, and that small starting covariances lose the finger motion fig. 6, fig. 7, fig. 8. On real stereo sequences of 100 frames each, the tracker keeps about 250 inliers per frame on average §7.2, fig. 9, fig. 10, fig. 11. The paper gives no accuracy numbers for the real hand sequences.

Paper vs code

No official code.

The paper and its arXiv metadata link no repository and no project page. A GitHub repository search for "ECMPR", "expectation conditional maximization point registration", "articulated point registration", "rigid articulated registration ECM" and "point registration org:INRIA" returned no matching repository.

Pitfalls

• The outlier radius is missing. The uniform-class term \emptyset_{3D} = 1.5\sqrt{2\pi} r^{-3} eq. 13 depends on a sphere radius r that the paper never gives, in the text or in the experiments §4. Choose r and state the choice. The working volume V cancels from the E-step and needs no value §4.
• Gaussian normalization is absorbed into \emptyset_{3D} (derived). Equation 12 drops the (2\pi)^{-3/2} factor from every Gaussian, and the \emptyset_{3D} term includes it eq. 4, eq. 12, eq. 13. If your Gaussians carry the normalization, the outlier term must become 1/v = 3/(4\pi r^3) instead. Use one convention only.
• The assumption nv \ll V is not checked §4. It fails with many model points or a large radius r, and then the outlier prior is no longer small.
• The εI covariance regularization is not analysed. The paper's claim that it leaves eccentricity and orientation unchanged §4 is not exact: adding εI shifts every eigenvalue by ε, which changes the ratio between the axes. The value of ε is not stated. The introduction says that the convergence proofs carry over §1, but the paper gives no argument for the εI step or for the SDP approximation.
• Convergence uses the rotation only §5.2 step 3. Translation, covariances and the objective do not enter the test, and no iteration cap is given.
• The SDP is a relaxation, and the paper does not say how to recover a rotation from it. The relaxed solution is used directly as R^{q+1} §5.2 step 3a. The solver is not named. The rank-one condition is dropped, so nothing forces \rho to be rank one or R to be orthonormal. The constraints list only R R^\top = I; the requirement |R| = +1 §5 is not in the relaxation. The abstract calls the method "optimal" [abstract], but the paper proves no optimality for the original nonconvex problem.
• Six constraints, not nine (derived). The paper writes the constraint for each pair (k, l) with k, l = 1, 2, 3 eq. 32, but says there are six \Delta_{kl} matrices §5.2. The matrix equation is symmetric, so there are six distinct constraints. Implement the pairs with k \le l.
• Dimension typo. The text calls r = \mathrm{vec}(R) a "9×2" vector §5.2. It is 9×1.
• Vectorization and Kronecker conventions. Appendix A defines \mathrm{vec} with the row-major order, and defines the Kronecker product in block form app. A. Check A and b numerically against the direct criterion eq. 28 before use. A column-major \mathrm{vec} gives a different A.
• Articulated links are rotations about one point (as written). Every Q_p with p \ge 1 has zero translation eq. 37. Under eq. 39, every part's translation equals the root translation t_0. The model then has all joints pivoting at the same point, unless the part points are shifted into place. The paper's source comments out the per-link translation term that would give general joints §6.1. The paper does not say how the finger model points are placed relative to joint centres §7.2.
• Parts after the root estimate rotation only. The translation of each part is fixed to t_0 eq. 44, so only the rotation is free in the rigid routine for those parts §6.2. The paper does not say how the translation is then checked.
• Sequential fitting makes errors sticky. The root's inliers are removed before the fingers are fitted §6.2. The paper does not state the order of the parts beyond the chain order, and it does not test the sensitivity to that order.
• Starting covariance matters. The paper says small starting covariances leave four of the five fingers unbent §7.2, fig. 6. It also reinitializes large covariances at each frame, but the size is not stated §7.2, fig. 11.
• Costs the abstract hides. Each E-step evaluates all n \times m pairs of model points and observations eq. 12. The paper lists KD-trees and classification EM as future speed-ups §8. The timings in Table 1 have no hardware named tab. 1. No timing is given for the articulated tracker.
• The "correct matches" metric is not comparable. The paper says that ECMPR's MAP labels and the nearest-point rule of ICP count matches differently §7.1. Compare the pose errors, as the paper does §7.1.
• ICP baselines are underspecified. The ICP thresholds are not listed tab. 2, fig. 5, and the "best" and "worst" rows are chosen after the fact. The TriICP sampling count is not stated §7.1.

Lineage

• Builds on: expectation maximization [Dempster, Laird & Rubin 1977] and its conditional variant ECM [Meng & Rubin 1993], both cited in §2.2, §2.3. The uniform outlier component follows [Banfield & Raftery 1993], as cited in §1.
• Related to the isotropic Gaussian-mixture registration family, cited in §1: Chui and Rangarajan, Granger and Pennec, and Myronenko et al. The paper's rigid rotation step for the isotropic case uses the closed-form methods of Arun et al., Horn, and Umeyama §5.1.
• Compared against: ICP [Besl & McKay 1992; Zhang 1994] §1; TriICP [Chetverikov et al. 2005] §7.1; and the articulated ICP of Munderman et al. §1.
• Refinement: the paper notes that its relaxed rotation is a good starting point for a nonlinear refiner such as Williams and Bennamoun's §5.2. The paper does not use that refiner in its algorithm.
• Superseded by: not stated.