FilterReg: Robust and Efficient Probabilistic Point-Set Registration using Gaussian Filter and Twist Parameterization

Introduction

Point-set registration is the task of aligning two point clouds by estimating their relative transformation. This problem is an essential component for many practical vision systems, such as SLAM , object pose estimation , dense 3d reconstruction , and interactive tracking of articulated and deformable objects.

The ICP algorithm is the most widely used method for this task. ICP alternatively establishes nearest-neighbor correspondences and minimizes the point-pair distances. With spatial indices such as the KD-tree, ICP provides relatively fast performance. The literature contains many variants of the ICP algorithm; and provide a thorough review and comparison.

Despite its popularity, the ICP algorithm is susceptible to noise, outliers and occlusions. These limitations have been widely documented in the literature . Thus, a great deal of research has been done on the use of probabilistic models for point-set registration , which can in principle provide better outlier-rejection. Additionally, if each point is given a Gaussian variance, the point cloud can be interpreted as a Gaussian Mixture Model (GMM). Most statistical registration methods are built on the GMM and empirically provide improved robustness . However, these methods tend to be much slower than the ICP and can hardly scale to large point clouds, which severely limits their practical usability.

In this paper, we present a novel probabilistic registration algorithm that achieves state-of-the-art robustness as well as substantially faster computational performance than modern ICP implementations. To achieve it, we propose a computationally-efficient probabilistic model and cast the registration as a maximum likelihood estimation, which can be solved using the EM algorithm. With a simple augmentation, we formulate the E step as a filtering problem and solve it using advances in efficient Gaussian filters . We also present a customized permutohedral filter with improved efficiency while retaining sufficient accuracy for our task. Empirically our method is as robust as state-of-the-art GMM-based methods, such as . In terms of the speed, our method with CPU is 3-7 times faster than modern ICP implementations and orders of magnitude faster than typical robust GMM-based methods. Furthermore, the proposed method can be GPU-parallelized and is 7 times faster than the CPU implementation.

Additionally, we propose a simple and efficient twist parameterization that extends our method to articulated and node-graph deformable objects. Our method is easy to implement and achieves substantial speedup over direct parameterization. For articulated objects, the complexity of our method is almost independent of the DOFs, which makes it highly efficient even for high-DOF systems. Combining these components, we present a robust, efficient and general registration method that outperforms many competitive baselines on a variety of registration tasks. The video demo, supplemental document and source code are available on our project page.

Related Work

The problem of point set registration is extensively pursued in computer vision and an exhaustive review is prohibitive. In the following text, we limit our discussion to GMM-based probabilistic registration and review them roughly according to their underlying probabilistic models.

The earliest statistical methods implicitly assumed the model points, which is controlled by the motion parameters (such as the rigid transformation or joint angles), induce a GMM distribution over the 3d space. The observation points are independently sampled from this distribution. Later, several contributions derived the EM procedure rigorously from the aforementioned probabilistic model. This formulation has also been applied to the registration of multi-rigid , articulated and deformable objects.

Another type of algorithms is known as the correlation-based methods . These algorithms treat both observation points and model points as probabilistic distributions. The point-cloud registration can be interpreted as minimizing some distance between distributions, for instance the KL-divergence. To improve the efficiency, techniques such as voxelization or Support Vector Machine are used to create compact GMM representations.

In this paper, we assume that the observation points induce a probabilistic distribution over the space. Intuitively, the registration is to move the model points to positions with large posterior probability, subject to kinematic constraints. This formulation is related to several existing works , and a more technical comparison is presented in Sec. 3.2. In addition to the formulation, the key contribution of our work includes the introduction of the filter-based correspondence and twist parameterization built on the probabilistic model, as mentioned in Sec. 1. Combining these components, the proposed method is general, robust and efficient that outperform various competitive baselines.

Probabilistic Model for Registration

In this subsection, we present our probabilistic model for point-set registration. We use $X,Y$ to denote the two point sets, $x_{1},x_{2},..x_{M}$ and $y_{1},y_{2},...,y_{N}$ are points in $X$ and $Y$. We define the model $X$ as the point set that is controlled by the motion parameter $\theta$. Another point set $Y$ is defined as the observation, which is fixed during the registration.

We are interested in the joint distribution $p(X,Y,\theta)$. We assume given model geometry $X$, the observation $Y$ is independent of $\theta$, and the joint distribution $p(X,Y,\theta)$ can be factored as

where $\phi_{\text{geometric}}(X,Y)$ is the potential function that encodes the geometric relationship, and the potential $\phi_{\text{kinematic}}(X,\theta)$ encodes the kinematic model. The $\phi_{\text{kinematic}}(X,\theta)$ can encode hard constraints such as $X=X(\theta)$ and/or soft motion regularizers, for instance the smooth terms in and the non-penetration term in .

We further assume the kinematic model $\phi_{\text{kinematic}}(X,\theta)$ has already captured the dependency within model points $X$. Thus, conditioned on the motion parameter $\theta$, the points in $X$ are independent of each other. The distribution can be further factored as

A factor graph representation of our model is shown in Fig. 1. With these factorization schemes, the conditional distribution can be written as

Following several existing work , we let the geometric distribution of each model point $\phi_{\text{geometric}}(x_{i}|Y)$ be a GMM,

where $p(x_{i}|y_{j})=\mathcal{N}(x_{i};y_{j},\Sigma_{xyz})$ is the Probability Density Function (PDF) of the Gaussian distribution, $y_{j}$ is the Gaussian centroid and $\Sigma_{xyz}=\text{diag}(\sigma_{x}^{2},\sigma_{y}^{2},\sigma_{z}^{2})$ is the diagonal covariance matrix. An additional uniform distribution $p(x_{i}|y_{N+1})=\frac{1}{M}$ is added to account for the noise and outliers. Similar to , we use equal membership probabilities $P(y_{j})=\frac{1}{N}$ for all GMM components, and introduce a parameter $0\leq w\leq 1$ to represent the ratio of outliers.

We estimate the motion parameter $\theta$ and model points $X$ by maximizing the following log-likelihood,

here we restrict ourselves to the kinematic model $X=X(\theta)$ and leave the general case to supplemental materials. We use the EM algorithm to solve this optimization. The EM procedure is

M step: minimize the following objective function

where $M^{0}_{x_{i}}$ and $M^{1}_{x_{i}}$ are computed in the E step (6), $c=\frac{w}{1-w}\frac{N}{M}$ is a constant, and $w$ is the parameter that represents the ratio of outliers.

The EM procedure is conceptually related to ICP. The weight-averaged target point $({M^{1}_{x_{i}}}/{M^{0}_{x_{i}}})$ replaces the nearest neighbour in ICP, and each model point is weighted by $\frac{M^{0}_{X_{i}}}{M^{0}_{X_{i}}+c}$. Intuitively, the averaged target provides robustness to noise in observation, while the weight for each model point should reject outliers in the model. Please refer to supplemental materials for the complete derivation.

2 Discussion and Comparison

At a high level, the proposed formulation can be viewed as an “inverse” of Coherent Point Drift (CPD) and many similar formulations , as shown in Fig. 1. CPD assumes the observation points are independently distributed according to a GMM introduced by model points, while the proposed formulation directly assumes the observation points induce a GMM over the space. Empirically, both methods are very robust to noise and outliers and significantly outperform ICP.

On the perspective of computation, the proposed method is much more simple and efficient than CPD and similar formulations . The proposed method only requires sum over $Y$ (6), while CPD requires sum over both $Y$ and $X$. Moreover, if a spatial index is used to perform this sum, CPD must rebuild the index every EM iteration as the model points $X$ are updated. In our formulation, we only need to build the index once if the variance is fixed during EM iterations, which is sufficient for many applications .

Several existing works also build a GMM representation of the observation points. Compared with our method, they do not explicitly account for the outlier distribution and miss the weight $\frac{M^{0}_{X_{i}}}{M^{0}_{X_{i}}+c}$. Furthermore, these methods assume each model point is only correlated with one or several “nearest” GMM centroids, while conceptually we assume each model point is correlated with all observation GMM centroids. Additionally, combined with the filter-based correspondence and twist parameterization in Sec. 4 and Sec. 5, our method tends to be much faster than these works, as demonstrated by our experiments.

3 Several Extensions

The presented probabilistic formulation can be extended to incorporate many well-established GMM-based registration techniques. Additionally, these extensions can be efficiently computed in a unified framework using the filter-based E step in Sec. 4 and the twist-based M step in Sec. 5. We select the optimized variance proposed in , feature correspondence in and point-to-plane residual in as three practically important examples, although many other methods can also be integrated in a very similar way.

Features: Currently in the E step (6), only the 3d position is used to measure the similarity between the model and observation points. Similar to , the E step can be extended to incorporate features such as normal, SHOT , learned features or their concatenation. The E step for arbitrary feature is

where $f_{x_{i}}$ and $f_{y_{k}}$ are the feature value for point $x_{i}$ and $y_{k}$, $\Sigma_{f}$ is the diagonal covariance for the feature.

Optimized Variance: In our previous formulation, the variance of Gaussian kernel $\Sigma_{xyz}$ is used as a fixed parameter. Similar to CPD , if $\Sigma_{xyz}=\text{diag}(\sigma^{2},\sigma^{2},\sigma^{2})$, the variance $\sigma$ can be interpreted as a decision variable and optimized analytically. Please refer to supplemental materials for the detailed formula and derivation.

Point-to-Plane Distance: The objective in our M step (7) is similar to the point-to-point distance in ICP, which doesn’t capture the planar structure. A simple solution is to compute a normal direction to characterize the local planar structure in the vicinity of the target $({M^{1}_{x_{i}}}/{M^{0}_{x_{i}}})$

where $N_{y_{k}}$ is the normal of the observation point $y_{k}$. The objective in the M step is then a point-to-plane error

E Step: Filter-based Correspondence

In this section, we discuss the method to compute the E step (6) and several extensions (8 and 9). These specific E steps can be written into the following generalized form

where $v_{y_{k}}$ generalizes the 3d position $y_{k}$ and the unit weight in (6, 8) and the normal $N_{y_{k}}$ in (9). The $G(f_{x_{i}})$ generalizes $M_{x_{i}}^{0}$ and $M_{x_{i}}^{1}$ in (6, 8) and the normal $N_{x_{i}}$ in (9). The features $f_{x_{i}}$ and $f_{y_{k}}$ generalize 3d positions $x_{i}$ and $y_{k}$ in the Gaussian PDF $\mathcal{N}(x_{i};y_{k},\Sigma_{xyz})$. The features $f_{x_{i}}$ and $f_{y_{k}}$ are normalized to have identity covariance. We also omit the normalization constant $\text{det}(2\pi\Sigma_{xyz})^{-\frac{1}{2}}$ of the Guassian PDF $\mathcal{N}(x_{i};y_{k},\Sigma_{xyz})$ for notational simplicity.

Equ. (11) is known as the general Gaussian Transform and the Improved Fast Gaussian Transform (IFGT) is proposed for it. However, IFGT uses a k-means tree internally and there would be too many k-means centroids for typical parameters in the registration. Practically, is not much faster than brute-force evaluation for our task.

We instead propose to compute (11) using Gaussian filtering algorithms , which demonstrate promising accuracy and efficiency on image processing. The filtering operation that these algorithms accelerated is

which is a subset of the general Gaussian transform: the feature $f_{y_{i}}$ used to retrieve the filtered value $G(f_{y_{i}})$ must be included in the input point set $Y$.

In our case, we would like to retrieve the value $G(f_{x_{i}})$ using feature $f_{x_{i}}$ from another point cloud $X$, which cannot be directly expressed in (12). To resolve it, we propose the following augmented input:

where $F_{X}=[f_{x_{1}},f_{x_{2}},...,f_{x_{M}}]$, $F_{Y}=[f_{y_{1}},f_{y_{2}},...,f_{y_{N}}]$ and $V_{Y}=[v_{y_{1}},v_{y_{2}},...,v_{y_{N}}]$. The new input feature $F_{\text{filter-input}}$ and value $V_{\text{filter-input}}$ are suitable for these filtering algorithms , and the filtered output is

With this augmentation, we can apply these filtering algorithms as a black box to our problem. However, by exploiting the structure of these methods, we can make them more efficient for our tasks. In the following text, the permutohedral lattice filter is discussed as an example, which is used in our experiments.

2 Permutohedral Lattice Filter

We briefly review the filtering process of , an illustration is shown in Fig. 2. The $d$-dimension feature $f$ is first embedded in $(d+1)$-dimensional space, where the permutohedral lattice lives. In the embedded space, each input value $v$ Splats onto the vertices of its enclosing simplex with barycentric weights. Next, lattice points Blur their values with nearby lattice points. Finally, the space is Sliced at each input position using the same barycentric weights to interpolate output values.

Although the permutohedral filter has demonstrated promising performance on a variety of tasks, it is still not optimal for our problem. In particular, the index building in can be inefficient when the variance $\Sigma_{xyz}$ is too small. Additionally, naively apply to the E step (6) requires rebuilding the index every EM iteration as the model point $X$ is updated. To resolve these problems, we propose a customization of the permutohedral filter that is more efficient while retaining sufficient accuracy for our task. The detailed method is presented in the supplemental material.

M Step: Efficient Twist Parameterization

In this section, we present methods to solve the optimizations (7, 10) with the twist parameterization. We first discuss the twist in the general kinematic model, then specialize it to articulated and node-graph deformable objects.

We focus on the following general kinematic model,

where $T_{i}(\theta)\in SE(3)$ is a rigid transformation, $x_{i\_\text{reference}}$ is the fixed reference point for the $x_{i}$. Note that $T_{i}(\theta)$ depends on $i$ and the kinematic model is not necessarily a global rigid transformation.

Twist is a 6-vector that represents the locally linearized “change” of $SE(3)$. Let the twist $\zeta_{i}=(w_{i},t_{i})=(\alpha_{i},\beta_{i},\gamma_{i},a_{i},b_{i},c_{i})$ be the local linearization of $T_{i}$, we have

Thus, the Jacobian $\frac{\partial x_{i}}{\partial\zeta_{i}}=[\text{skew}(x_{i}),I_{3\times 3}]$ is a $3\times 6$ matrix, where $I_{3\times 3}$ is identity matrix, and $\text{skew}(x_{i})$ is a $3\times 3$ matrix such that $\text{skew}(x_{i})b=\text{cross}(x_{i},b)$ for arbitrary $b\in{R}^{3}$.

The optimization (7, 10) are least squares problems, and we focus on the following generalized form of them

where $r_{x_{i}}$ is the concatenated least-squares residuals that depends on $x_{i}$. We use the Gauss-Newton (GN) algorithm to solve (17). In each GN iteration we need to compute the following $A$ and $b$ matrices by

and the update of the motion parameters is $\Delta\theta=-A^{-1}b$. Thus, the primary computational bottleneck is to assemble the matrices $A$ and $b$. In the following text, we only discuss the computation of the $A$ matrix, while the computation of the $b$ vector is similar and easier. The computation of the $A$ matrix can be written as

where $\frac{\partial\zeta_{i}}{\partial\theta}$ is the Jacobian that maps the change of motion parameter $\theta$ to the change of the rigid transformation $T_{i}$, while the change of $T_{i}$ is expressed as its twist. Note that the term $\frac{\partial{r_{x_{i}}}}{\partial\zeta_{i}}=\frac{\partial{r_{x_{i}}}}{\partial x_{i}}\frac{\partial x_{i}}{\partial\zeta_{i}}$ is very easy to compute, as both $\frac{\partial{r_{x_{i}}}}{\partial x_{i}}$ and $\frac{\partial x_{i}}{\partial\zeta_{i}}$ are only dependent on $x_{i}$.

If the kinematic model is a global rigid transformation, we have $\frac{\partial\zeta_{i}}{\partial\theta}=I_{6\times 6}$ and $A=\sum_{x_{i}}((\frac{\partial{r_{x_{i}}}}{\partial\zeta_{i}})^{T}\frac{\partial{r_{x_{i}}}}{\partial\zeta_{i}})$. In the following subsections, we proceed to the articulated and node-graph deformable kinematic models.

Articulated objects consist of rigid bodies connected through joints in a kinematic tree. A broad set of real-world objects, including human bodies, hands and robots are articulated objects. If the kinematic model (15) is an articulated model, the motion parameter $\theta\in R^{N_{\text{joint}}}$ would be the joint angles, where $N_{\text{joint}}$ is the number of joints. The $T_{i}(\theta)$ is the rigid transform of the rigid body that the point $x_{i}$ is attached to. The computation of the $A$ matrix can be factored as

where $\zeta_{j}$ is the twist of rigid body $j$, and we have exploited $\frac{\partial\zeta_{i}}{\partial\zeta_{j}}=I_{6\times 6}$ if point $i$ is on rigid body $j$. Importantly, $\frac{\partial\zeta_{j}}{\partial\theta}$ is known as the spatial velocity Jacobian and is provided by many off-the-shelf rigid body simulators . The algorithm that uses (20) is shown in Algorithm 1.

The lines 1-4 of Algorithm 1 dominates the overall performance and the complexity is O($6^{2}M$), where $M$ is the number of model points and usually $M\gg N_{\text{joint}}$. Thus, the complexity of this algorithm is almost independent of $N_{\text{joint}}$. As a comparison, previous articulated registration methods need O($N_{\text{joint}}^{2}M$) time to assemble the $A$ matrix, and $N_{\text{joint}}$ is usually much larger than 6. Furthermore, lines 1-4 of Algorithm 1 is very simple to implement and can be easily GPU parallelized. Combined with an off-the-shelf simulator, the overall pipeline can achieve promising efficiency. On the contrary, previous methods typically need a customized kinematic tree implementation for real-time performance, while requires substantial software engineering effort to realize.

2 Node-Graph Deformable Model

To capture the motion of objects such as rope or cloth, we need a kinematic model which allows large deformation while preventing unrealistic collapsing or distortion. In this paper, we follow to represent the general deformable kinematic model as a node graph. Intuitively, the node graph defines a motion field in the 3D space and the reference vertex in Equ. (15) is deformed according to the motion field. More specifically, the node graph is defined as a set $\{[p_{j}\in R^{3},T_{j}\in SE(3)]\}$, where $j$ is the node index, $p_{j}$ is the position of the $j$th node, and $T_{j}$ is the $SE(3)$ motion defined on the $j$th node. The kinematic equation (15) can be written as

where $N_{i}(x_{i\_\text{reference}})$ is the nearest neighbor nodes of model point $x_{i\_\text{reference}}$, and $w_{ki}$ is the fixed skinning weight. The interpolation of the rigid transformation $T_{k}$ is performed using the DualQuaternion representation of the $SE(3)$.

The $A$ matrices for this kinematic model can be constructed using an algorithm very similar to Algorithm 1. The detailed method is provided in supplemental materials.

Results

We conduct a variety of experiments to test the robustness, accuracy and efficiency of the proposed method. Our hardware platform is an Intel i7-3960X CPU except for Sec. 6.5, where the proposed method is implemented with CUDA on a Nivida Titan Xp GPU. The video demo and the source code are available on our project page.

We follow CPD to setup an experiment on synthetic data. We use a subsampled Stanford bunny with 3500 points. The initial rotation discrepancy is 50 degrees with a random axis. The proposed method is compared with two baselines: CPD , a representative GMM-based algorithm; TrICP , a widely used robust ICP variant. Parameters for all methods are well tuned and provided in supplemental materials. We use the following metric to measure the pose estimation error

where $T_{\text{gt}}$ is the known ground truth pose, $x_{i\_\text{reference}}$ defined in (15) is the reference position. We terminate the algorithm when the twist (change of transformation) is less than a threshold. In this way, the final alignment error (22) is about 1 [mm] for all methods. All of the statistical results are the averaged value of 30 independent runs.

Fig. 3 shows the robustness of different algorithms with respect to outliers in the point sets. We add different number of points randomly to both the model and observation clouds. An example of such point sets with initial alignment is shown in Fig. 3 (a), the converged alignment by the proposed method and TrICP are shown in Fig. 3 (b) and Fig. 3 (c), respectively. The proposed method and CPD significantly outperform the robust ICP.

Fig. 4 shows the robustness of different algorithms with respect to noise in the point sets. We corrupt each point in both model and observation clouds with a Gaussian noise. The unit of the noise is the diameter of the Bunny. An example of such point sets with initial alignment are shown in Fig. 4 (a). Fig. 4 (b) and (c) are the final alignment by the proposed method and TrICP initialized from (a). Note that we use clean point clouds for better visualization. Our method and CPD are more accurate than the robust ICP.

Table. 1 summarizes the computational performance of each method. The running time is measured on clean point cloud. Our method is about 7 times faster than TrICP and two orders of magnitude faster than CPD . The proposed method with fixed $\sigma$ is faster per iteration, but need more iterations to converge. Overall the proposed method is as robust as the state-of-the-art statistical registration algorithm CPD , and runs substantially faster than the modern ICP implementation.

2 Rigid Registration on Real-World Data

We follow to setup this experiment: the algorithm is used to compute the frame-to-frame rigid transformation on the Stanford Lounge dataset . We register every 5th frame for the first 400 frames, each downsampled to about 5000 points. The average Euler angle deviation from the ground truth is used as the estimation error.

Fig. 5 (a) shows an example registration by the proposed method. Fig. 5 (b) shows the accuracy and speed of various algorithms. The results of baseline methods are from Our CPU (i7-3960X) is slightly inferior to (i7-5920K), and we observe similar accuracy and slightly worse speed using CPD and TrICP . Thus, we think our speed result are comparable to despite hardware difference. except for CPD . For CPD we use $\sigma_{\text{init}}=20\text{ [cm]}$ instead of the data-based initialization of , with which we observed improved performance. As the point-to-point error doesn’t capture the planar structure, the point-to-point version of the proposed method as well as many other point-to-point algorithms are less accurate on this dataset. The proposed method with point-to-plane error achieves state-of-the-art accuracy. On the perspective of computation, the proposed method significantly outperforms all the baselines, including GMMTree and EMICP that rely on a high-end Titan X GPU.

3 Global Pose Estimation using Learned Features

We demonstrate global pose estimation using motion-invariant features. The task is to align a pre-built geometric model to observation clouds from RGBD images, where both the model and observation clouds are colored by the learned feature . We use the proposed method with feature correspondence in Sec. 3.3 and fixed $\sigma=0.05$ as the feature has unit norm. The proposed method is compared with a modified TrICP : the nearest neighbour is searched in feature space (instead of 3d-space). After feature-based registration, we apply 3d-space local refinement to get the final alignment.

Fig. 6 shows an example registration. Note that we treat the background as outliers. As shown in Fig. 6 (a), the observation (RGBD cloud) is under severe occlusion and contains very ambiguous outliers. Fig. 6 (b) and (c) show feature-based registration by our method and the “feature” TrICP. The proposed method is more robust to the outliers and occlusion. Fig. 6 (d) and (e) show the final alignment using local refinement initialized from (b) and (c). The proposed method converges to correct pose while the baseline is trapped to bad alignment. Table. 2 summaries the success rate of both methods on 30 RGBD images with different view points and lighting conditions. Our method has a higher success rate and is more efficient than the baseline.

4 Articulated Tracking

The proposed method with articulated kinematic model is used to track a robot manipulating a box. The robot and box model has 20 DOFs (12 for the floating bases of the box and the robot, 8 for robot joints). We use drake for the kinematic computation in (20). We use fixed $\sigma=1$ [cm] and set the maximum EM iterations to be 15. Our template has 7500 points and the depth cloud has about 10000 points.

Fig. 7 (a) shows the snapshots of the tracked manipulation scenario. Fig. 7 (b) shows the live geometric model and the observation clouds. Points from observation are in black, while the geometric model is colored according to rigid bodies. Note that points from the table are treated as outliers. Fig. 7 (c) summaries the averaged per-frame performance of various algorithms. The proposed twist parameterization is an order of magnitude faster than direct parameterization. Combining the filter-based correspondence and twist parameterization leads to a real-time tracking algorithm and substantial performance improvement over articulated ICP and .

5 Application to Dynamic Reconstruction

The proposed method with node-graph deformable kinematic is implemented on GPU and used as the internal non-rigid tracker of DynamicFusion (our implementation). The proposed method is compared with the projective ICP, the original non-rigid tracker of . We use fixed $\sigma=2$ [cm]. Fig. 8 shows both methods operate on a RGBD sequence with relative fast motions. The proposed method tracks it correctly, while the projective ICP fails to track the right hand of the actor. The proposed method is more robust to fast and tangential motion than the projective ICP.

To test the efficiency of the proposed twist parameterization on node-graph deformable objects, we compare it with Opt , a highly optimized GPU least squares solver using direct parameterization. The per-frame computational performance of various algorithms is summarized in Table. 3. The GPU parallelization of our filter-based E step achieves 8 times speedup over the CPU version, and the proposed twist parameterization is about 20 times faster than .

Conclusion

To conclude, we present a probabilistic registration method that achieves state-of-the-art robustness, accuracy and efficiency. We show that the correspondence search can be formulated as a filtering problem, and employ advances in efficient Gaussian filtering methods to solve it. In addition, we present a simple and efficient twist parameterization that generalizes our method to articulated and deformable objects. Extensive empirical evaluation demonstrates the effectiveness our method.

Acknowledgments This work was supported by NSF Award IIS-1427050 and Amazon Research Award. The views expressed in this paper are those of the authors themselves and are not endorsed by the funding agencies.

References