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 -norm in the position constraint and the transformation constraint ,,. However, the quadratic energy functional is more easily affected by noise and outliers. To address this problem, Yang  propose a sparse non-rigid registration (SNR) method with an -norm regularized model for the transformation constraint. However, their position constraint is still based on the -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 -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 -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  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.  minimize the shortest distance between a point in the template and the tangent plane of the closest point on the target. Pottmann et al.  propose a registration method with quadratic convergence, which gives faster and more stable convergence than the standard ICP . Bouaziz  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.  place markers on the object to help reconstruct the pose of scan and use it as a basis for modeling deformation. Pekelny et al.  use predefined bone information to find bone transformations.
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 -norm in a least-squares sense ,.
Regularization terms help to preserve smoothness, making the optimization more robust to noise and outliers, and -norm is also widely used in regularization terms. Süßmuth et al.  use a generalized as-rigid-as-possible energy  to promote smoothness. Liao et al.  define a transformation model using the TPS , and use graduated assignment for non-rigid registration and optimization. Wand et al.  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.  solve correspondences, confidence weights, and a deformation field within a single optimization framework using -norm. Hontani et al.  propose a statistical shape model (SSM) which is incorporated into the nonrigid ICP (NICP), and outliers can be detected based on their sparseness. Yang  propose a sparse non-rigid registration (SNR) method with an -norm regularized model for the smoothness. However, their -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.
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 ground truth matchings (vertices) as correspondences, and solve the transformations using the SNR method  which measures the positional errors in the standard quadratic term to avoid bias towards the -norm. The Laplacian distribution fits the histogram significantly better than the Gaussian distribution, suggesting the use of sparsity-promoting -norm in the data term.
Our -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 -norm allows a large distance along a particular dimension although the Euclidean distance is not significant. In this sense, -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)  that applies the -norm on the image gradient over the isotropic TV  that measures TV as the sum of -norm (not squared). Birkholz  showed that anisotropic TV achieves better denoising performance in preserving the geometries of corners in images. We choose the -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  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 . 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 be a 3D point in the homogenous coordinate. Denote by a template set of 3D points and by a target set of 3D points, where and are the numbers of points. Denote by the correspondence of . Define as the index mapping from the template points to the target points, where means the corresponding vertex cannot be found for the -th vertex. Denote by the transformation matrix for point . Define as the set of non-rigid transformations. For compact notation, we define as a matrix containing the transformation matrices to be solved. The proposed method is to find non-rigid transformations that transforms the template into the target as accurately as possible, given a correspondence mapping .
The non-rigid registration is formulated as the minimization of the following energy function:
where , and are data term, smoothness term, and orthogonality constraint, respectively. and 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 , for each point. The weight is one if there is a corresponding point on the target shape for , and zero otherwise. Hence, we propose the following data term
where is the Cartesian coordinate of .
For the compact representation in algorithm derivation, we define the following matrix/vector form of the variables to reformulate data term (2):
where is a diagonal matrix containing the input elements as diagonal entities. Then, the data term can be rewritten as
Smoothness term: In the smoothness term, local rigidity is assumed: for vertex , the transformations of neighboring vertices should have very close transformed positions when applied to . Therefore, we define the following smoothness term:
Define a graph , where the vertices of the graph are the 3D points in , and the edges of the graph are denoted by . 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 -nearest neighbors ( is typically set to 6). Denote the neighborhood of vertex by , and an edge is defined between each neighboring vertex and . So, we have . Similar to the data term, we define a differential matrix on the graph for concise presentation. Concretely, each row of corresponds to an edge in and each column corresponds to a vertex in . Each row in has only two nonzero entries. For example, assuming the row in associates with edge , then the entry related to the reference vertex is set at 1, while the one related to the neighboring vertex is set at -1, i.e. and . Let denote the row of . We introduce a matrix , where the row of is defined as . Therefore, the cost of transformation smoothness is rewritten as
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.
where is a rotation matrix, and is a constant matrix that extracts the rotation component of . ensures that is a rotation matrix, not a mirrored matrix.
The final energy function has the following compact form with matrix-vector notations:
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:
where and are diagonal weighting matrices for the data term and smoothness term, respectively. The weighting matrices are updated according to the -norm of the corresponding entries. For the data term, the weights are updated as
where represents the index of iteration, 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
where is a constant which is set as 0.01 in the experiments, and the row of matrix is associated with edge between and .
To solve the problem, we first transform the minimization (9) into the following form with auxiliary variables and :
Then, we solve the constrained minimization (12) using the augmented Lagrangian method (ALM) . The ALM method converts the original problem (12) to iterative minimization of its augmented Lagrangian function:
where (, ) are positive constants, (,) are Lagrangian multipliers, and denotes the inner product of two matrices considered as long vectors. Under the standard ALM framework, (, ) and (, ) can be efficiently updated. However, each iteration has to solve , , and simultaneously, which is difficult yet computationally demanding. Hence, we resort to the alternate direction method (ADM)  to optimize , , and separately at each iteration:
The -subproblem has the following closed solution:
where shrink(,) is the shrinkage function applied on the matrix element-wise:
The -subproblem is solved in a similar way:
The -subproblem can be explicitly solved using Procrustes projection:
If the obtained matrix has a negative determinant, take with the opposite sign to turn the matrix into a rotation matrix.
Being quadratic, the -subproblem can be readily solved by using the first-order optimality condition:
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:
where and 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:
The iterative non-rigid registration with reweighting is summarized in Algorithm 1, and the algorithm for minimization (9) is summarized in Algorithm 2.
|Algorithm 1. Algorithm of reweighting non-rigid registration|
|1. Input: template , target .|
|2. While not converged do|
|3. Find correspondence mapping ;|
|4. Update and acco. to (10) and (11), resp.|
|5. Solve transformations via Algorithm (2);|
|6. End while|
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  and a human motion dataset . Fig. 2 and Fig. 3 give the registration results on cat and jumping datasets, compared with the classic -norm regularized non-rigid ICP method and the SNR method . 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 as the ground-truth correspondence of . For a vertex , the registration error is defined as . The compared classic -norm based non-rigid ICP method is formulated as optimizing:
The smoothness constraint of this kind of methods is imposed on the transformation differences. To ensure fair comparison, we adjust the weight 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  and the classic -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  in Fig. 4. Obvious registration errors can be seen in the result of the method in , especially in the right foot (top) and head (bottom), while the methods with sparse representation (SNR  and our method) achieve better registration results. The method in  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 , 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  and the -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 and 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  and the remaining one third using local geometric feature matching based on SHOT signatures . 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  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 noiseis normalized by , where is the average length of triangle edges on the associated target mesh, and chosen in the range of . Fig. 8 gives the registration results compared with the SNR method  and the -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 , , 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 and are set as 0.01. The standard deviation of the noise is set as 1, and the percentage of outliers is set as . 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 outliers in Fig. 12 to evaluate the influence of the paremeters and . To make experiments more tractable, we adjust both parameters consistently (i.e. ). 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  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 -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 . 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 -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 -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.
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.
The authors would like to thank Ke Li for her help with some experiments, and thank Shuai Lin for help with comparative experiments with .
-  B. Allen, B. Curless, and Z. Popović. Articulated body deformation from range scan data. ACM Trans. Graph., 21(3):612–619, 2002.
-  B. Allen, B. Curless, and Z. Popović. The space of human body shapes: reconstruction and parameterization from range scans. ACM Trans. Graph., 22(3):587–594, 2003.
B. Amberg, S. Romdhani, and T. Vetter.
Optimal step nonrigid icp algorithms for surface registration.
IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 1–8, 2007.
-  D. P. Bertsekas. Constrained optimization and lagrange multiplier methods. Computer Science and Applied Mathematics, Boston: Academic Press, 1, 1982.
-  P. J. Besl and N. D. McKay. Method for registration of 3-D shapes. In Robotics-DL tentative, pages 586–606, 1992.
-  H. Birkholz. A unifying approach to isotropic and anisotropic total variation denoising models. JCAM, 235(8).
-  S. Bouaziz, A. Tagliasacchi, and M. Pauly. Sparse iterative closest point. Computer Graphics Forum, 32(5):113–123, 2013.
S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein.
Distributed optimization and statistical learning via the alternating
direction method of multipliers.
Foundations and Trends® in Machine Learning, 3(1):1–122, 2011.
-  A. M. Bronstein, M. M. Bronstein, and R. Kimmel. Numerical geometry of non-rigid shapes. Springer Science & Business Media, 2008.
-  Y. Chen and G. Medioni. Object modeling by registration of multiple range images. In IEEE International Conference on Robotics and Automation, pages 2724–2729. IEEE, 1991.
-  M. Chuang, L. Luo, B. J. Brown, S. Rusinkiewicz, and M. Kazhdan. Estimating the Laplace-Beltrami operator by restricting 3D functions. Computer Graphics Forum, 28(5):1475–1484, 2009.
-  H. Chui and A. Rangarajan. A new point matching algorithm for non-rigid registration. Computer Vision and Image Understanding, 89(2):114–141, 2003.
-  S. Esedolu et al. Decomposition of images by the anisotropic Rudin-Osher-Fatemi model. CPAM, 57(12).
-  K. Guo, F. Xu, Y. Wang, Y. Liu, and Q. Dai. Robust non-rigid motion tracking and surface reconstruction using L0 regularization. IEEE International Conference on Computer Vision (ICCV), 2015.
-  H. Hontani, T. Matsuno, and Y. Sawada. Robust nonrigid icp using outlier-sparsity regularization. In IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 174–181, 2012.
-  H. Li, R. W. Sumner, and M. Pauly. Global correspondence optimization for non-rigid registration of depth scans. Computer graphics forum, 27(5):1421–1430, 2008.
-  M. Liao, Q. Zhang, H. Wang, R. Yang, and M. Gong. Modeling deformable objects from a single depth camera. In IEEE International Conference on Computer Vision (ICCV), pages 167–174, 2009.
-  R. A. Newcombe, A. J. Davison, S. Izadi, P. Kohli, O. Hilliges, J. Shotton, D. Molyneaux, S. Hodges, D. Kim, and A. Fitzgibbon. KinectFusion: Real-time dense surface mapping and tracking. In 10th IEEE international symposium on mixed and augmented reality (ISMAR), pages 127–136. IEEE, 2011.
-  R. A. Newcombe, D. Fox, and S. M. Seitz. Dynamicfusion: Reconstruction and tracking of non-rigid scenes in real-time. In IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 343–352. IEEE, 2015.
-  C. Papazov and D. Burschka. Deformable 3d shape registration based on local similarity transforms. In Computer Graphics Forum, volume 30, pages 1493–1502. Wiley Online Library, 2011.
-  Y. Pekelny and C. Gotsman. Articulated object reconstruction and markerless motion capture from depth video. Computer Graphics Forum, 27(2):399–408, 2008.
-  H. Pottmann, Q.-X. Huang, Y.-L. Yang, and S.-M. Hu. Geometry and convergence analysis of algorithms for registration of 3d shapes. International Journal of Computer Vision, 67(3):277–296, 2006.
-  H. Pottmann, S. Leopoldseder, and M. Hofer. Registration without ICP. Computer Vision and Image Understanding, 95(1):54–71, 2004.
-  M. Rouhani, E. Boyer, and A. D. Sappa. Non-rigid registration meets surface reconstruction. In Proc. 3D Vision, 2014.
-  L. Rudin et al. Nonlinear total variation based noise removal algorithms. PDNP, 60(1).
-  S. Salti, F. Tombari, and L. Di Stefano. SHOT: unique signatures of histograms for surface and texture description. Computer Vision and Image Understanding, 125:251–264, 2014.
-  K. Sidorov, S. Richmond, D. Marshall, et al. Efficient groupwise non-rigid registration of textured surfaces. In IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 2401–2408. IEEE, 2011.
-  O. Sorkine and M. Alexa. As-rigid-as-possible surface modeling. In Proc. Symposium on Geometry Processing, pages 109–116, 2007.
-  R. W. Sumner and J. Popović. Deformation transfer for triangle meshes. ACM Trans. Graph., 23(3):399–405, 2004.
-  J. Süßmuth, M. Winter, and G. Greiner. Reconstructing animated meshes from time-varying point clouds. Computer Graphics Forum, 27(5):1469–1476, 2008.
-  G. K. Tam, R. R. Martin, P. L. Rosin, and Y.-K. Lai. Diffusion pruning for rapidly and robustly selecting global correspondences using local isometry. ACM Trans. Graph., 33(1):4, 2014.
-  J. Tong, J. Zhou, L. Liu, Z. Pan, and H. Yan. Scanning 3d full human bodies using kinects. IEEE Transactions on Visualization and Computer Graphics, 18(4):643–50, 2012.
-  D. Vlasic, I. Baran, W. Matusik, and J. Popović. Articulated mesh animation from multi-view silhouettes. ACM Trans. Graph., 27(3):97, 2008.
-  M. Wand, B. Adams, M. Ovsjanikov, A. Berner, M. Bokeloh, P. Jenke, L. Guibas, H.-P. Seidel, and A. Schilling. Efficient reconstruction of nonrigid shape and motion from real-time 3D scanner data. ACM Trans. Graph., 28(2):15, 2009.
-  J. Yang, K. Li, K. Li, and Y.-K. Lai. Sparse non-rigid registration of 3D shapes. In Computer Graphics Forum, volume 34, pages 89–99, 2015.