Robust Non-Rigid Registration with Reweighted Position and Transformation Sparsity

TL;DR

This paper introduces a robust non-rigid registration method that employs reweighted sparsity on position and transformation, effectively handling noise and outliers in 3D shape alignment.

Contribution

It proposes a novel double sparsity model with reweighting scheme for improved robustness and guarantees convergence through an optimized subproblem approach.

Findings

Outperforms state-of-the-art methods in accuracy.

More robust to noise and outliers.

Effective on both public and real datasets.

Abstract

Non-rigid registration is challenging because it is ill-posed with high degrees of freedom and is thus sensitive to noise and outliers. We propose a robust non-rigid registration method using reweighted sparsities on position and transformation to estimate the deformations between 3-D shapes. We formulate the energy function with position and transformation sparsity on both the data term and the smoothness term, and define the smoothness constraint using local rigidity. The double sparsity based non-rigid registration model is enhanced with a reweighting scheme, and solved by transferring the model into four alternately-optimized subproblems which have exact solutions and guaranteed convergence. Experimental results on both public datasets and real scanned datasets show that our method outperforms the state-of-the-art methods and is more robust to noise and outliers than conventional non-rigid registration methods.

Tables

Algorithm 1. Algorithm of reweighting non-rigid registration

1.  Input: template $𝒱$, target $𝒰$.

2.   While not converged do

3.    Find correspondence mapping $f^{\left(l\right)}:𝒱\mapsto 𝒰$;

4.    Update ${𝐖}_{\text{D}}^{\left(l\right)}$ and ${𝐖}_{\text{S}}^{\left(l\right)}$ acco. to (10) and (11), resp.

5.    Solve transformations ${𝐗}^{\left(l\right)}$ via Algorithm (2);

6.   End while

7.  Output: $𝐗$

Algorithm 2. ADMM algorithm to solve (9)

1.  Input: ${\widetilde{𝐔}}_{f^{\left(l\right)}}\in {𝐑}^{N\times 3}$, $𝐕\in {𝐑}^{N\times 4N}$, $𝐁\in {𝐑}^{\left|ℰ\right|\times 4\left|𝒱\right|}$;

2.   Initialize: ${𝐗}^{(l,0)}=X^{\left(l-1\right)}$, ${𝐘}_1^{\left(0\right)},{𝐘}_2^{\left(0\right)}=0$;

${\mu }_1,{\mu }_2>0$, ${\rho }_1,{\rho }_2>1$;

3.   While not converged do

4.    Solve ${𝐂}^{(l,k+1)}$ by (15);

5.    Solve ${𝐀}^{(l,k+1)}$ by (17);

6.    Solve ${𝐑}_i^{(l,k+1)}$ by (18);

7.    Solve ${𝐗}^{(l,k+1)}$ by (LABEL:equ:LDL_decomposition_WS_WD)$\sim$(21);

8.    Update ${\mu }_1^{\left(k+1\right)}$, and ${\mu }_2^{\left(k+1\right)}$ according (14);

9.    Update $Y_1^{\left(k+1\right)}$, and $Y_2^{\left(k+1\right)}$ according (14);

10.   End while

11.  Output: ${𝐗}^{\left(l\right)}$.

Table 3. Table I: Comparison on running times

Num. vertexes	1000	2000	5000	10000
${\ell }_2$-norm	1.23s	3.51s	12.88s	29.78s
SNR	8.05s	17.36s	52.48s	119.06s
Group sparsity	7.39s	24.83s	59.96s	126.58s
Ours	7.17s	22.13s	55.68s	122.85s

Equations

$$E\left(\mathbf{X};f\right)=E_{data}\left(\mathbf{X};f\right)+\alpha E_{smooth}\left(\mathbf{X}\right)+\beta E_{orth}\left(\mathbf{X}\right),$$

$$\begin{array}{rl}{E}_{data}\left(\mathbf{X};f\right)& =\sum _{{\mathbf{v}}_i\in \mathcal{V}}{w_i\Vert \begin{array}{c}{\mathbf{X}}_i{\mathbf{v}}_i-{\tilde{\mathbf{u}}}_{f(i)}\end{array}\Vert }_1,\end{array}$$

$$\begin{array}{rl}\mathbf{W}& =\text{diag}\left(\sqrt{w_1},\cdots ,\sqrt{w_N}\right),\\ \mathbf{V}& =\text{diag}\left({\mathbf{v}}_1^{\top },\cdots ,{\mathbf{v}}_N^{\top }\right),\\ {\tilde{\mathbf{U}}}_f& ={\left[\begin{array}{ccc}{\tilde{\mathbf{u}}}_{f(1)}& \cdots & {\tilde{\mathbf{u}}}_{f(N)}\end{array}\right]}^{\top },\end{array}$$

$$E_{data}\left(\mathbf{X};f\right)={\Vert \begin{array}{c}\mathbf{W}\left(\mathbf{VX}-{\tilde{\mathbf{U}}}_f\right)\end{array}\Vert }_1.$$

$$\begin{array}{rl}{E}_{smooth}\left(\mathbf{X}\right)& =\sum _{{\mathbf{v}}_i\in \mathcal{V}}\sum _{{\mathbf{v}}_j\in {\mathcal{N}}_i}{\Vert \begin{array}{c}{\mathbf{X}}_i{\mathbf{v}}_i-{\mathbf{X}}_j{\mathbf{v}}_i\end{array}\Vert }_1.\end{array}$$

$$E_{smooth}\left(\mathbf{X}\right)={\Vert \begin{array}{c}\mathbf{BX}\end{array}\Vert }_1.$$

$$\begin{array}{rl}& E_{orth}\left(\mathbf{X}\right)=\sum _{i=1}^N{\Vert \begin{array}{c}{\mathbf{S}}_i{\mathbf{X}}_i-{\mathbf{R}}_i\end{array}\Vert }_F^2,\\ & \text{leavevmode}\text{leavevmode}\text{s.t.}{\mathbf{R}}_i^T{\mathbf{R}}_i=\mathbf{I},\text{leavevmode}det({\mathbf{R}}_i)>0,\end{array}$$

$$\begin{array}{rl}& \underset{\mathbf{X}}{\min }{\Vert \begin{array}{c}\mathbf{W}\left(\mathbf{VX}-{\tilde{\mathbf{U}}}_f\right)\end{array}\Vert }_1+\alpha {\Vert \begin{array}{c}\mathbf{BX}\end{array}\Vert }_1+\beta \sum _{i=1}^N{\Vert \begin{array}{c}{\mathbf{S}}_i{\mathbf{X}}_i-{\mathbf{R}}_i\end{array}\Vert }_F^2,\\ & \text{s.t.}{\mathbf{R}}_i^T{\mathbf{R}}_i=\mathbf{I},\text{leavevmode}det({\mathbf{R}}_i)>0.\end{array}$$

$$\begin{array}{rl}& \underset{\mathbf{X}}{\min }{\Vert \begin{array}{c}{\mathbf{W}}_{\text{D}}\left(\mathbf{VX}-{\tilde{\mathbf{U}}}_f\right)\end{array}\Vert }_1+\alpha {\Vert \begin{array}{c}{\mathbf{W}}_{\text{S}}\mathbf{BX}\end{array}\Vert }_1+\beta \sum _{i=1}^N{\Vert \begin{array}{c}{\mathbf{S}}_i{\mathbf{X}}_i-{\mathbf{R}}_i\end{array}\Vert }_F^2,\\ & \text{s.t.}{\mathbf{R}}_i^T{\mathbf{R}}_i=\mathbf{I},\text{leavevmode}det({\mathbf{R}}_i)>0.\end{array}$$

$${\mathbf{W}}_{\text{D}}^{(l)}(i,i)=\{\begin{array}{ll}\frac{1}{{\Vert {\mathbf{X}}_i^{(l-1)}{\mathbf{v}}_i-{\tilde{\mathbf{u}}}_{f(i)}^{(l)}\Vert }_1+{ϵ}_{\text{D}}},& f(i)\ne 0,\\ 0,& f(i)=0,\end{array}$$

$${\mathbf{W}}_{\text{S}}^{(l)}(r,r)=\frac{1}{{\Vert {\mathbf{X}}_i^{(l-1)}{\mathbf{v}}_i-{\mathbf{X}}_j^{(l-1)}{\mathbf{v}}_i\Vert }_1+{ϵ}_{\text{S}}},$$

$$\begin{array}{rl}\underset{\mathbf{X},\mathbf{C},\mathbf{A}}{\min }& {\Vert \begin{array}{c}\mathbf{C}\end{array}\Vert }_1+\alpha {\Vert \begin{array}{c}\mathbf{A}\end{array}\Vert }_1+\beta \sum _{i=1}^N{\Vert \begin{array}{c}{\mathbf{S}}_i{\mathbf{X}}_i-{\mathbf{R}}_i\end{array}\Vert }_F^2,\\ & \text{s.t.}\mathbf{C}={\mathbf{W}}_{\text{D}}\left(\mathbf{VX}-{\tilde{\mathbf{U}}}_f\right),\\ & \mathbf{A}={\mathbf{W}}_{\text{S}}\mathbf{BX},\\ & {\mathbf{R}}_i^T{\mathbf{R}}_i=\mathbf{I},\text{leavevmode}det({\mathbf{R}}_i)>0.\end{array}$$

$$\begin{array}{rl}& L(\mathbf{X},\mathbf{C},\mathbf{A},\{{\mathbf{R}}_i\},{\mathbf{Y}}_1,{\mathbf{Y}}_2,{\mu }_1,{\mu }_2)={\Vert \begin{array}{c}\mathbf{C}\end{array}\Vert }_1+\alpha {\Vert \begin{array}{c}\mathbf{A}\end{array}\Vert }_1\\ & +⟨{\mathbf{Y}}_1,\mathbf{C}-{\mathbf{W}}_{\text{D}}\left(\mathbf{VX}-{\tilde{\mathbf{U}}}_f\right)⟩\\ & +\frac{{\mu }_1}{2}{\Vert \begin{array}{c}\mathbf{C}-{\mathbf{W}}_{\text{D}}\left(\mathbf{VX}-{\tilde{\mathbf{U}}}_f\right)\end{array}\Vert }_F^2\\ & +⟨{\mathbf{Y}}_2,\mathbf{A}-{\mathbf{W}}_{\text{S}}\mathbf{BX}⟩+\frac{{\mu }_2}{2}{\Vert \begin{array}{c}\mathbf{A}-{\mathbf{W}}_{\text{S}}\mathbf{BX}\end{array}\Vert }_F^2\\ & +\beta \sum _{i=1}^N{\Vert \begin{array}{c}{\mathbf{S}}_i{\mathbf{X}}_i-{\mathbf{R}}_i\end{array}\Vert }_F^2,\\ & \text{s.t.}{\mathbf{R}}_i^T{\mathbf{R}}_i=\mathbf{I},\text{leavevmode}det({\mathbf{R}}_i)>0,\end{array}$$

$$\{\begin{array}{ll}{\mathbf{C}}^{(k+1)}=& \arg {\min }_{\mathbf{C}}\text{leavevmode}\Vert \mathbf{C}{\Vert }_1\\ & +⟨{\mathbf{Y}}_1^{(k)},\mathbf{C}-{\mathbf{W}}_{\text{D}}\left(\mathbf{V}{\mathbf{X}}^{(k)}-{\tilde{\mathbf{U}}}_f\right)⟩\\ & +\frac{{\mu }_1^{(k)}}{2}{\Vert \begin{array}{c}\mathbf{C}-{\mathbf{W}}_{\text{D}}\left(\mathbf{V}{\mathbf{X}}^{(k)}-{\tilde{\mathbf{U}}}_f\right)\end{array}\Vert }_F^2,\\ {\mathbf{A}}^{(k+1)}=& \arg {\min }_{\mathbf{A}}\text{leavevmode}\alpha \Vert \mathbf{A}{\Vert }_1+⟨{\mathbf{Y}}_2^{(k)},\mathbf{A}-{\mathbf{W}}_{\text{S}}\mathbf{B}{\mathbf{X}}^k⟩\\ & +\frac{{\mu }_2^{(k)}}{2}{\Vert \begin{array}{c}\mathbf{A}-{\mathbf{W}}_{\text{S}}\mathbf{B}{\mathbf{X}}^{(k)}\end{array}\Vert }_F^2,\\ {\mathbf{R}}_i^{(k+1)}=& \arg {\min }_{{\mathbf{R}}_i}\text{leavevmode}\beta {\sum }_{i=1}^N{\Vert \begin{array}{c}{\mathbf{S}}_i{\mathbf{X}}_i^{(k)}-{\mathbf{R}}_i\end{array}\Vert }_F^2\\ & \text{s.t.}{\mathbf{R}}_i^T{\mathbf{R}}_i=\mathbf{I},\text{leavevmode}det({\mathbf{R}}_i)>0\\ {\mathbf{X}}^{(k+1)}=& \arg {\min }_{\mathbf{X}}\text{leavevmode}⟨{\mathbf{Y}}_1^{(k)},{\mathbf{C}}^{(k+1)}-{\mathbf{W}}_{\text{D}}\left(\mathbf{VX}-{\tilde{\mathbf{U}}}_f\right)⟩\\ & +\frac{{\mu }_1^{(k)}}{2}{\Vert \begin{array}{c}{\mathbf{C}}^{(k+1)}-{\mathbf{W}}_{\text{D}}\left(\mathbf{VX}-{\tilde{\mathbf{U}}}_f\right)\end{array}\Vert }_F^2\\ & +⟨{\mathbf{Y}}_2^{(k)},{\mathbf{A}}^{(k+1)}-{\mathbf{W}}_{\text{S}}\mathbf{BX}⟩\\ & +\frac{{\mu }_2^{(k)}}{2}{\Vert \begin{array}{c}{\mathbf{A}}^{(k+1)}-{\mathbf{W}}_{\text{S}}\mathbf{BX}\end{array}\Vert }_F^2\\ & +\beta {\sum }_{i=1}^N{\Vert \begin{array}{c}{\mathbf{S}}_i{\mathbf{X}}_i-{\mathbf{R}}_i^{(k+1)}\end{array}\Vert }_F^2,\\ {\mathbf{Y}}_1^{(k+1)}=& {\mathbf{Y}}_1^{(k)}+{\mu }_1^{(k)}\left({\mathbf{C}}^{(k+1)}-{\mathbf{W}}_{\text{D}}\left(\mathbf{V}{\mathbf{X}}^{(k+1)}-{\tilde{\mathbf{U}}}_f\right)\right),\\ {\mathbf{Y}}_2^{(k+1)}=& {\mathbf{Y}}_2^{(k)}+{\mu }_2^{(k)}\left({\mathbf{A}}^{(k+1)}-{\mathbf{W}}_{\text{S}}\mathbf{B}{\mathbf{X}}^{(k+1)}\right),\\ \text{leavevmode}{\mu }_1^{(k+1)}=& {\rho }_1{\mu }_1^{(k)},\text{leavevmode}\text{leavevmode}{\rho }_1>1,\\ \text{leavevmode}{\mu }_2^{(k+1)}=& {\rho }_2{\mu }_2^{(k)},\text{leavevmode}\text{leavevmode}{\rho }_2>1.\end{array}$$

$$\begin{array}{rl}& {\mathbf{C}}^{(k+1)}=\\ & \text{shrink}\left({\mathbf{W}}_{\text{D}}\left(\mathbf{V}{\mathbf{X}}^{(k)}-{\tilde{\mathbf{U}}}_f\right)-\frac{1}{{\mu }_1^{(k)}}{\mathbf{Y}}_1^{(k)},\frac{1}{{\mu }_1^{(k)}}\right),\end{array}$$

$$\text{shrink}\left(x,\tau \right)=\text{sign}(x)\max (|x|-\tau ,0).$$

$${\mathbf{A}}^{(k+1)}=\text{shrink}\left({\mathbf{W}}_{\text{S}}\mathbf{B}{\mathbf{X}}^{(k)}-\frac{1}{{\mu }_2^{(k)}}{\mathbf{Y}}_2^{(k)},\frac{\alpha }{{\mu }_2^{(k)}}\right).$$

$$\begin{array}{rl}\mathbf{UD}{\mathbf{V}}^T& =\text{svd}({\mathbf{S}}_i{\mathbf{X}}_i^k),\\ {\mathbf{R}}_i^{k+1}& =\mathbf{U}\mathbf{missing}{V}^T.\end{array}$$

$$\begin{array}{rl}& \left({\mu }_1^{(k)}{\mathbf{V}}^{\top }{\mathbf{W}}_{\text{D}}^{\top }{\mathbf{W}}_{\text{D}}\mathbf{V}+{\mu }_2^{(k)}{\mathbf{B}}^{\top }{\mathbf{W}}_{\text{S}}^{\top }{\mathbf{W}}_{\text{S}}\mathbf{B}+\beta \sum _{i=1}^N{\mathbf{S}}_i^T{\mathbf{S}}_i\right)\mathbf{X}\\ & ={\mathbf{B}}^{\top }{\mathbf{W}}_{\text{S}}^{\top }\left({\mathbf{Y}}_2^{(k)}+{\mu }_2^{(k)}{\mathbf{A}}^{(k+1)}\right)\\ & +{\mathbf{V}}^{\top }{\mathbf{W}}_{\text{D}}^{\top }\left({\mathbf{Y}}_1^{(k)}+{\mu }_1^{(k)}\left({\mathbf{C}}^{(k+1)}+{\mathbf{W}}_{\text{D}}{\tilde{\mathbf{U}}}_f\right)\right)\\ & +\beta \sum _{i=1}^N{\mathbf{S}}_i^T{\mathbf{R}}_i^{(k+1)}.\end{array}$$

$$\begin{array}{rl}& \left(\mathbf{L},\mathbf{D}\right)=\\ & \text{ldl}\left({\mu }_1^{(k)}{\mathbf{V}}^{\top }{\mathbf{W}}_{\text{D}}^{\top }{\mathbf{W}}_{\text{D}}\mathbf{V}+{\mu }_2^{(k)}{\mathbf{B}}^{\top }{\mathbf{W}}_{\text{S}}^{\top }{\mathbf{W}}_{\text{S}}\mathbf{B}+\beta \sum _{i=1}^N{\mathbf{S}}_i^T{\mathbf{S}}_i\right),\end{array}$$

$$\begin{array}{rl}\mathbf{LQ}& ={\mathbf{V}}^{\top }{\mathbf{W}}_{\text{D}}^{\top }\left({\mathbf{Y}}_1^{(k)}+{\mu }_1^{(k)}\left({\mathbf{C}}^{(k+1)}+{\mathbf{W}}_{\text{D}}{\tilde{\mathbf{U}}}_f\right)\right)\\ & +{\mathbf{B}}^{\top }{\mathbf{W}}_{\text{S}}^{\top }\left({\mathbf{Y}}^{(k)}+{\mu }_2^{(k)}{\mathbf{A}}^{(k+1)}\right)+\beta \sum _{i=1}^N{\mathbf{S}}_i^T{\mathbf{R}}_i^{(k+1)}\\ \mathbf{DZ}& =\mathbf{Q},\\ {\mathbf{L}}^{\top }\mathbf{X}& =\mathbf{Z}.\end{array}$$

$$\underset{\mathbf{X}}{\min }\text{leavevmode}\text{leavevmode}{\Vert \begin{array}{c}\mathbf{W}(\mathbf{VX}-{\tilde{\mathbf{U}}}_f)\end{array}\Vert }_F^2+\alpha {\Vert \begin{array}{c}\mathbf{BX}\end{array}\Vert }_F^2.$$

Full text

1 Introduction

Non-rigid registration is a hot research topic in computer graphics and computer vision [32, 27, 14, 19], and is a key technique for dynamic 3-D reconstruction using a depth camera. Commodity depth sensors, e.g., Microsoft Kinect, become cheaper and more widely used, but depth images and reconstructed point clouds captured by such devices contain much noise. Hence, non-rigid registration methods robust to noise and outliers are highly desirable to scan dynamic scenes with deformable objects.

Given two input 3-D shapes, one as the template shape and the other as the target shape, non-rigid registration aims to find a suitable transformation that when applied deforms the template shape to be aligned with the target shape. Non-rigid registration is often formulated as an optimization problem. Most methods formulate some energy functional with both position and transformation constraints. The position constraint measures the closeness of the transformed template shape and the target shape, and the transformation constraint measures the fitness to model, which always includes the smoothness, namely the total energy of transformation differences of all the local neighbors. Most work uses the classic squared $\ell_{2}$-norm in the position constraint and the transformation constraint [16],[3],[28]. However, the quadratic energy functional is more easily affected by noise and outliers. To address this problem, Yang $et\leavevmode\nobreak\ al.$ [35] propose a sparse non-rigid registration (SNR) method with an $\ell_{1}$-norm regularized model for the transformation constraint. However, their position constraint is still based on the $\ell_{2}$-norm. In practice, e.g. for near piece-wise rigid deformation, which is common for real-world deformable objects, the positional error tends to concentrate on small regions. This cannot be modeled well using the $\ell_{2}$-norm.

In this paper, we propose a non-rigid registration method with sparsity-regularized position and transformation constraints. The distribution of positional errors and transformation differences for typical non-rigid deformation can be well modeled using the Laplacian distribution, or equivalently, the $\ell_{1}$-norm should be used to measure both the positional errors and transformation differences, which is therefore called dual sparsities. To promote the sparsity, we adopt a reweighted sparse model, which is solved by the alternating direction method of multipliers (ADMM). The proposed method is evaluated on public datasets [9, 33] and real datasets captured by a RGB-D depth sensor. The results demonstrate that the proposed method obtains better results than the state-of-the-art non-rigid registration methods.

The main contributions of this work are summarized as:

• •

We propose a dual-sparsity based non-rigid registration method on both position and transformation constraints. The proposed model is robust against outliers as the sparsity terms allow a small fraction of regions with larger deviations.

• •

We incorporate orthogonality constraints in the dual-sparsity non-rigid registration framework to promote locally rigid transformations.

• •

We equip the dual-sparsity based non-rigid registration model with a reweighted scheme to iteratively enhance sparsity in the series of alternating optimization subproblems.

2 Related Work

3-D shape registration consists of rigid registration and non-rigid registration. Rigid registration aims to find a global rigid-body transformation, while non-rigid registration needs to find a set of local transformations that align two shapes.

In rigid registration, the 3-D shapes are assumed to be aligned by a Euclidean transformation, including rotation and translation. Iterative Closest Point (ICP) and its variants [5] are the dominant algorithms for rigid registration. This kind of methods alternates between two steps: 1) finding closest points and 2) solving the optimal transformation. As an improved method of ICP, Chen et al. [10] minimize the shortest distance between a point in the template and the tangent plane of the closest point on the target. Pottmann et al. [23] propose a registration method with quadratic convergence, which gives faster and more stable convergence than the standard ICP  [22]. Bouaziz $et\leavevmode\nobreak\ al.$ [7] propose a new variant of the ICP algorithm, which uses sparsity-inducing norms to represent the positional constraint, and achieve better results for the situation with noise and outliers. Their work focuses on rigid registration with low degrees of freedom, and hence regularization is not necessary.

When shapes have large deformations from template to target, automatic non-rigid registration is necessary. It is more challenging due to its high degrees of freedom, and an appropriate deformation model is the key for an efficient and robust algorithm.

Some methods compute global rigid transformations for bones and local non-rigid transformations near joints, which is essentially a piecewise rigid transformation model. Allen et al. [1] place markers on the object to help reconstruct the pose of scan and use it as a basis for modeling deformation. Pekelny et al. [21] use predefined bone information to find bone transformations.

Some models take more generic deformations into consideration. Chui et al. [12] use the thin-plate spline (TPS) as the non-rigid transformation model. Papazov et al. [20] allow points to move freely and use an additional uniform distribution to limit noise and outliers, and propose an ordinary differential equation (ODE) model. Local affine transformations [2] are also frequently used in non-rigid registration. Liao et al. [17] use differential coordinates as local affine transformations with smoothness constraints. Amberg et al. [3] use a stiffness term to ensure similarity of adjacent transformations. Rouhani et al. [24] model non-rigid deformation as an integration of locally rigid transformations. In our work, we use local affine transformation with an orthogonality constraint as it allows more flexibility to capture fine surface details while keeping local shapes.

Non-rigid registration is often formulated as an energy functional with data and regularization terms. Most of the non-rigid registration work models the data term in $\ell_{2}$-norm in a least-squares sense [29],[3].

Regularization terms help to preserve smoothness, making the optimization more robust to noise and outliers, and $\ell_{2}$-norm is also widely used in regularization terms. Süßmuth et al. [30] use a generalized as-rigid-as-possible energy [28] to promote smoothness. Liao et al. [17] define a transformation model using the TPS [12], and use graduated assignment for non-rigid registration and optimization. Wand et al. [34] take a set of time-varying point data as input, and reconstruct a single shape and a deformation field that fit the data. To improve robustness, Li et al. [16] solve correspondences, confidence weights, and a deformation field within a single optimization framework using $\ell_{2}$-norm. Hontani et al. [15] propose a statistical shape model (SSM) which is incorporated into the nonrigid ICP (NICP), and outliers can be detected based on their sparseness. Yang $et\leavevmode\nobreak\ al.$ [35] propose a sparse non-rigid registration (SNR) method with an $\ell_{1}$-norm regularized model for the smoothness. However, their $\ell_{2}$-norm position constraint cannot model the concentricity of positional errors well.

In this paper, based on the observation that the deformations of 3-D surfaces vary smoothly and the positional distances and transformation differences are sparse, we propose a non-rigid registration method with sparse position and transformation constraints. The model is efficiently solved by the alternating direction method under the augmented Lagrangian multiplier framework.

3 Motivation

The traditional quadratic data term assumes the Gaussian distribution of positional errors. However, transformations are piecewise smooth signals residing on 3D surfaces, resulting in larger positional errors for geometric details and joints. This suggests that the positional errors are sparse, and should be modeled by a heavy-tailed distribution, rather than being dense and modeled by a rapidly vanishing Gaussian distribution. This is verified in Fig. 1(b). We uniformly pick up $10\%$ ground truth matchings (vertices) as correspondences, and solve the transformations using the SNR method [35] which measures the positional errors in the standard quadratic term to avoid bias towards the $\ell_{1}$-norm. The Laplacian distribution fits the histogram significantly better than the Gaussian distribution, suggesting the use of sparsity-promoting $\ell_{1}$-norm in the data term.

Our $\ell_{1}$-norm sparsity measures equally coordinate differences of each dimension. Another possibility is to use the sum of Euclidean distances (group sparsity) between corresponding points, which also well fits the distribution of positional errors as shown in Fig. 1(c). The group sparsity advocates sparsity for each Euclidean distance as a whole, while the $\ell_{1}$-norm allows a large distance along a particular dimension although the Euclidean distance is not significant. In this sense, $\ell_{1}$-norm is more flexible to preserve large non-rigid deformation along some dimensions. Such an advantage is also observed in the anisotropic total variation (TV) [13] that applies the $\ell_{1}$-norm on the image gradient over the isotropic TV [25] that measures TV as the sum of $\ell_{2}$-norm (not squared). Birkholz [6] showed that anisotropic TV achieves better denoising performance in preserving the geometries of corners in images. We choose the $\ell_{1}$-norm to measure the positional errors for its potential flexibility, and also for its easier and faster implementation with an element-wise shrinkage (cf. Table 1 for statistics of running times).

4 The Proposed Method

4.1 Iterative Framework

We iteratively compute the deformation between the template shape and the target shape. Each iteration consists of two steps. In the first step, the correspondences between template and target are estimated using the registration result obtained from the last iteration. At the beginning of the iteration, we use a technique based on local geometric similarity and diffusion pruning of inconsistent correspondence [31] as it often provides reliable correspondences. Alternative correspondence techniques or manual specification of a few correspondences may instead be used (an example is shown in Fig. 5). These computed correspondences are referred to as the correspondence mapping $f$. Then, we use the closest points between template and target shapes to find additional correspondences similar to ICP. In the second step (Sec. 4.2), we propose an energy-minimization approach based on dual-sparsity representation to estimate the non-rigid transformations using the correspondences obtained from the first step.

4.2 Deformation Estimation

Let $\mathbf{v}_{i}\triangleq[x_{i},y_{i},z_{i},1]^{\top}$ be a 3D point in the homogenous coordinate. Denote by $\mathcal{V}\triangleq\{{\mathbf{v}_{1}},\cdots,\mathbf{v}_{N}\}$ a template set of 3D points and by $\mathcal{U}\triangleq\{\mathbf{u}_{1},\cdots,\mathbf{u}_{M}\}$ a target set of 3D points, where $N$ and $M$ are the numbers of points. Denote by $\mathbf{u}_{f(i)}\in\mathcal{U}$ the correspondence of $\mathbf{v}_{i}\in\mathcal{V}$. Define $f:\{1,\cdots,N\}\mapsto\{0,1,\cdots,M\}$ as the index mapping from the template points to the target points, where $f(i)=0$ means the corresponding vertex cannot be found for the $i$-th vertex. Denote by $\mathbf{X}_{i}$ the $3\times 4$ transformation matrix for point $\mathbf{v}_{i}$. Define $\mathcal{X}\triangleq\left\{\mathbf{X}_{1},\cdots,\mathbf{X}_{N}\right\}$ as the set of non-rigid transformations. For compact notation, we define $\mathbf{X}\triangleq\left[\mathbf{X}_{1},\cdots,\mathbf{X}_{N}\right]^{\top}$ as a matrix containing the $N$ transformation matrices to be solved. The proposed method is to find non-rigid transformations $\mathbf{X}$ that transforms the template $\mathcal{V}$ into the target $\mathcal{U}$ as accurately as possible, given a correspondence mapping $f$.

The non-rigid registration is formulated as the minimization of the following energy function:

[TABLE]

where $E_{data}\left(\mathbf{X}\right)$, $E_{smooth}\left(\mathbf{X}\right)$ and $E_{orth}\left(\mathbf{X}\right)$ are data term, smoothness term, and orthogonality constraint, respectively. $\alpha$ and $\beta$ adjust the importance of different terms. The data term measures the position accuracy, the smoothness term imposes a smoothness constraint so that the original ill-posed problem (defined by only the data term) is now well-posed, and the orthogonality constraint promotes locally rigid transformations, which is particularly needed for underconstrained scenarios such as partial meshes.

Data term: We measure the accuracy of deformation as the closeness of the transformed points to their corresponding target points. We assign a weight, denoted by $w_{i}$, for each point. The weight $w_{i}$ is one if there is a corresponding point on the target shape for $\mathbf{v}_{i}$, and zero otherwise. Hence, we propose the following data term

[TABLE]

where $\tilde{\mathbf{u}}_{f(i)}$ is the Cartesian coordinate of $\mathbf{u}_{f(i)}$.

For the compact representation in algorithm derivation, we define the following matrix/vector form of the variables to reformulate data term (2):

[TABLE]

where $\textrm{diag}(\cdot)$ is a diagonal matrix containing the input elements as diagonal entities. Then, the data term can be rewritten as

[TABLE]

Smoothness term: In the smoothness term, local rigidity is assumed: for vertex $\mathbf{v}_{i}$, the transformations of neighboring vertices $\mathbf{v}_{j}\in\mathcal{N}_{i}$ should have very close transformed positions when applied to $\mathbf{v}_{i}$. Therefore, we define the following smoothness term:

[TABLE]

Define a graph $\mathcal{G}\triangleq\left(\mathcal{V},\mathcal{E}\right)$, where the vertices of the graph are the 3D points in $\mathcal{V}$, and the edges of the graph are denoted by $\mathcal{E}$. For a 3D mesh, edges of the graph are simply defined by the edges of the mesh; for 3D point clouds, edges can be defined by connecting each vertex with its $K$-nearest neighbors ($K$ is typically set to 6). Denote the neighborhood of vertex $\mathbf{v}_{i}$ by $\mathcal{N}_{i}$, and an edge $e_{ij}$ is defined between each neighboring vertex $\mathbf{v}_{j}$ and $\mathbf{v}_{i}$. So, we have $\mathcal{E}=\left\{e_{ij}\mid\mathbf{v}_{j}\in\mathcal{N}_{i},\mathbf{v}_{i}\in\mathcal{V}\right\}$. Similar to the data term, we define a differential matrix $\mathbf{K}\in\{-1,1\}^{|\mathcal{E}|\times|\mathcal{V}|}$ on the graph $\mathcal{G}$ for concise presentation. Concretely, each row of $\mathbf{K}$ corresponds to an edge in $\mathcal{E}$ and each column corresponds to a vertex in $\mathcal{V}$. Each row in $\mathbf{K}$ has only two nonzero entries. For example, assuming the $r^{\textrm{th}}$ row in $\mathbf{K}$ associates with edge $e_{ij}$, then the entry related to the reference vertex $\mathbf{v}_{i}$ is set at 1, while the one related to the neighboring vertex $v_{j}$ is set at -1, i.e. $k_{ri}=1$ and $k_{rj}=-1$. Let $\mathbf{k}_{i:}$ denote the $i^{\textrm{th}}$ row of $\mathbf{K}$. We introduce a matrix $\mathbf{B}\in R^{|\mathcal{E}|\times 4|\mathcal{V}|}$, where the $i^{\textrm{th}}$ row of $\mathbf{B}$ is defined as $\mathbf{b}_{i:}:=\mathbf{k}_{i:}\otimes\mathbf{v}_{i}^{\top}$. Therefore, the cost of transformation smoothness is rewritten as

[TABLE]

Orthogonality constraint: Especially for partial meshes with large motions, the problem may be underconstrained leading to large distortions. In this case, orthogonality constraint is effective in better preserving local shapes and making the solution more reasonable.

[TABLE]

where $\mathbf{R}_{i}$ is a $3\times 3$ rotation matrix, and $\mathbf{S}_{i}$ is a constant $3\times 4$ matrix that extracts the rotation component of $\mathbf{X}_{i}$. $det(\mathbf{R}_{i})>0$ ensures that $\mathbf{R}_{i}$ is a rotation matrix, not a mirrored matrix.

The final energy function has the following compact form with matrix-vector notations:

[TABLE]

Reweighting: To further promote sparsity, both the data term and the smoothness term are weighted, and the weighting matrices are updated at each iteration of non-rigid registration. The weighted version of the dual-sparsity model (8) is defined as follows:

[TABLE]

where $\mathbf{W}_{\textrm{D}}$ and $\mathbf{W}_{\textrm{S}}$ are diagonal weighting matrices for the data term and smoothness term, respectively. The weighting matrices are updated according to the $\ell_{1}$-norm of the corresponding entries. For the data term, the weights are updated as

[TABLE]

where $l$ represents the index of iteration, $\epsilon_{\textrm{D}}$ is a constant to avoid the division-by-zero issue, and is set as 0.01 in the experiments. Similarly, the weights for the smoothness term are updated as

[TABLE]

where $\epsilon_{\textrm{S}}$ is a constant which is set as 0.01 in the experiments, and the $r^{\textrm{th}}$ row of matrix $\mathbf{B}\mathbf{X}$ is associated with edge $e_{ij}$ between $\mathbf{v}_{i}$ and $\mathbf{v}_{j}$.

To solve the problem, we first transform the minimization (9) into the following form with auxiliary variables $\mathbf{A}$ and $\mathbf{C}$:

[TABLE]

Then, we solve the constrained minimization (12) using the augmented Lagrangian method (ALM) [4]. The ALM method converts the original problem (12) to iterative minimization of its augmented Lagrangian function:

[TABLE]

where ($\mu_{1}$, $\mu_{2}$) are positive constants, ($\mathbf{Y}_{1}$,$\mathbf{Y}_{2}$) are Lagrangian multipliers, and $\langle\cdot,\cdot\rangle$ denotes the inner product of two matrices considered as long vectors. Under the standard ALM framework, ($\mathbf{Y}_{1}$, $\mathbf{Y}_{2}$ ) and ($\mu_{1}$, $\mu_{2}$) can be efficiently updated. However, each iteration has to solve $\mathbf{A}$, $\mathbf{C}$, $\{\mathbf{R}_{i}\}$ and $\mathbf{X}$ simultaneously, which is difficult yet computationally demanding. Hence, we resort to the alternate direction method (ADM) [8] to optimize $\mathbf{A}$, $\mathbf{C}$, $\{\mathbf{R}_{i}\}$ and $\mathbf{X}$ separately at each iteration:

[TABLE]

The $\mathbf{C}$-subproblem has the following closed solution:

[TABLE]

where shrink($\cdot$,$\cdot$) is the shrinkage function applied on the matrix element-wise:

[TABLE]

The $\mathbf{A}$-subproblem is solved in a similar way:

[TABLE]

The $\mathbf{R}_{i}$-subproblem can be explicitly solved using Procrustes projection:

[TABLE]

If the obtained matrix has a negative determinant, take $\mathbf{R}_{i}$ with the opposite sign to turn the matrix into a rotation matrix.

Being quadratic, the $\mathbf{X}$-subproblem can be readily solved by using the first-order optimality condition:

[TABLE]

However, the straightforward matrix inversion in solving (19) is inefficient or even practically impossible for large-scale problems, e.g., registration of tens of thousands of points. This can be relieved by using the LDL decomposition:

[TABLE]

where $\mathbf{L}$ and $\mathbf{D}$ are the lower triangular matrix and the diagonal matrix of the LDL decomposition. Then, the linear equations in (19) is solved by solving the following much easier linear systems:

[TABLE]

The iterative non-rigid registration with reweighting is summarized in Algorithm 1, and the algorithm for minimization (9) is summarized in Algorithm 2.

5 Experimental Results

In this section, we evaluate the performances of the proposed method on clean datasets (Section 5.1), noisy datasets (Section 5.2), and real scans (Section 5.3). Running times of our method are reported in Section 5.4.

5.1 Results on Clean Datasets

We evaluate the proposed method on two datasets: TOSCA high-resolution dataset [9] and a human motion dataset [33]. Fig. 2 and Fig. 3 give the registration results on cat and jumping datasets, compared with the classic $\ell_{2}$-norm regularized non-rigid ICP method and the SNR method [35]. The results are shown as the overlap of the deformed template shape (blue) and the target shape (gray) and the fitting errors are color-coded on the reconstructed mesh for visual inspection. Denote $\mathbf{g}_{i}$ as the ground-truth correspondence of $\mathbf{v}_{i}$. For a vertex $\mathbf{v}_{i}$, the registration error is defined as $\|\mathbf{X}_{i}\mathbf{v}_{i}-\mathbf{g}_{i}\|^{2}_{2}$. The compared classic $\ell_{2}$-norm based non-rigid ICP method is formulated as optimizing:

[TABLE]

The smoothness constraint of this kind of methods is imposed on the transformation differences. To ensure fair comparison, we adjust the weight $\alpha$ until we get the most accurate registration without loss of smoothness for each method. The result shows that our method achieves the best results with less fitting errors in the areas with intensive deformations than the SNR method [35] and the classic $\ell_{2}$-norm regularized non-rigid ICP method, such as the tail of the cat and the wrinkles around the waist of the person highlighted in rectangles.

We also compare our method with state-of-the-art non-rigid registration [16] in Fig. 4. Obvious registration errors can be seen in the result of the method in [16], especially in the right foot (top) and head (bottom), while the methods with sparse representation (SNR [35] and our method) achieve better registration results. The method in [16] works effectively when the template and target shapes are close so that good initial correspondences can be obtained, but the pose changes substantially in this example. Moreover, our result is more accurate and better-distributed for the whole body than the SNR method [35], due to the sparse constraint on the position.

To evaluate the robustness of the proposed method, we manually assign 35 correspondences on Jumping dataset, and compare the result of our method with the SNR method [35] and the $\ell_{2}$-regularized method. As shown in Fig. 5, our method achieves the best result, especially around the places with substantial deformation, e.g., the right knee.

To evaluate the effectiveness of the proposed reweighting scheme, we compare the registration results with and without reweighting on Bouncing dataset in Fig. 6. The parameters $\epsilon_{\textrm{D}}$ and $\epsilon_{\textrm{S}}$ are set as 0.01. As shown in the figure, the reweighting scheme significantly improves the registration results.

5.2 Results on Noisy Datasets

1) Correspondences with partially incorrect matchings:

It is common to include incorrect correspondences using established methods. We simulate this in two cases. In the first case, we obtain two thirds of correspondences using diffusion pruning [31] and the remaining one third using local geometric feature matching based on SHOT signatures [26]. The majority of correspondences from the former are correct while many correspondences from the latter are incorrect due to the ambiguity of local features. In the second case, we generate all the correspondences using SHOT signatures. Fig. 7 gives the results for the two cases in a difficult situation which involves very complex transformations from template to target. As shown in the figure, our method is significantly more robust than the SNR method [35] with respect to incorrect correspondences.

2) Target shapes with noise or outliers:

In the first case, 3-D shapes of targets are polluted with dense noise along the norm directions of the associated vertices. All the target vertices are perturbed with Gaussian noise. The standard deviation of the noise $\sigma$ is normalized by $\bar{l}$, where $\bar{l}$ is the average length of triangle edges on the associated target mesh, and chosen in the range of $[0.1,1]$. Fig. 8 gives the registration results compared with the SNR method [35] and the $\ell_{2}$-norm regularization method. The results show that our method is more robust to noise, performing significantly better for models with high noise levels.

In the second case, 3-D shapes of targets are polluted with sparse outliers along the normal directions of the associated vertices. Fig. 9 gives the results for the situations when $1\%$, $2\%$, $5\%$ of target vertices are perturbed with Gaussian noise. The results show that our method is more robust than the other two methods, particularly for cases with larger proportion of outliers.

To evaluate the effectiveness of the proposed reweighting scheme, we also compare the registration results with and without reweighting for noise and outlier cases on Bouncing dataset in Fig. 10 and Fig. 11. The parameters $\epsilon_{\textrm{D}}$ and $\epsilon_{\textrm{S}}$ are set as 0.01. The standard deviation of the noise $\sigma$ is set as 1, and the percentage of outliers is set as $50\%$. It can be seen that the reweighting scheme contributes significantly to improve the registration results for the dataset with noise and outliers.

We compare the registration results with different parameter settings for the reweighting scheme on Bouncing dataset with $50\%$ outliers in Fig. 12 to evaluate the influence of the paremeters $\epsilon_{\textrm{D}}$ and $\epsilon_{\textrm{S}}$. To make experiments more tractable, we adjust both parameters consistently (i.e. $\epsilon_{\textrm{D}}=\epsilon_{\textrm{S}}=\epsilon$). It can be seen that the best setting is 0.006 for this case, which has the smallest fitting errors. However, the performance is quite close, and 0.01 is a generally good choice (found in experiments).

5.3 Results on Real Scans

Fig. 13 presents the results on real scans generated by Kinect Fusion [18] using Kinect V2.0. The real scans are very challenging, because they have much noise and a large number of outliers. Moreover, each mesh is incomplete and the topology between the template and the target is inconsistent. Hence, it is difficult to obtain sufficient and reliable correspondences. The overlap of the deformed template and the target show that the $\ell_{2}$-norm regularization method and the SNR method present misalignments around the hands, arms and some other joints which have large deformations, while the result of our method is well-distributed and better registered.

Fig. 14 gives an example of generating a complete color mesh for a human head. A base mesh is scanned by Kinect Fusion using Kinect V2.0, and four partial color meshes are registered to the base mesh using our method. The textures are blended by solving the Poisson equation over the surface of mesh [11]. As shown in the figure, our method correctly registers the input view surfaces with better registration than alternative methods, and successfully generates a watertight color mesh.

5.4 Running times

We compare the running times of the proposed method with the $\ell_{2}$-norm regularized method, SNR method, and group sparsity method on Crane dataset. We downsample the meshes into smaller meshes with 1K to 10K vertices. The number of nICP registration iterations for each method is set as 20, and $\ell_{1}$-norm has extra 20 inner iterations for each outer iteration. All the experiments are performed on a desktop computer with Intel i5 3.2GHz CPU and 8GB RAM. The comparison results are shown in Table 1. Our method has similar time complexity as SNR.

6 Conclusions

This paper proposes a non-rigid registration method with reweighted sparse position and transformation constraints. We formulate the energy function with dual sparsity on both the data term and the smoothness term, and define the smoothness constraint using local rigidity. The dual-sparsity based non-rigid registration model is equipped with a reweighting scheme, and solved by the alternating direction method under the augmented Lagrangian multiplier (ADM-ALM) framework which have exact solutions and guaranteed convergence. Experimental results on both public datasets and real scans show that our method provides significantly improved results over alternative methods, especially for more challenging cases, and is more robust to noise and outliers.

Acknowledgments

The authors would like to thank Ke Li for her help with some experiments, and thank Shuai Lin for help with comparative experiments with [16].