The ability to quantify the uncertainty of a surface registration is important in many areas of shape analysis. It is especially useful for the reconstruction of partial data, or when analysing data where the exact correspondence is unclear, such as smooth surfaces. Uncertainty estimates can be used to build robust Statistical Shape Models (SSMs) by assigning an uncertainty measure to each landmark [15, 7]. Within the medical area, knowing the uncertainty of surface reconstructions is important to make surgical decisions based on registered data . In this paper, we propose an efficient, fully probabilistic method for surface registration based on Metropolis-Hastings (MH) sampling. Efficiency is gained by the use of a specialised proposal based on a non-rigid Iterative Closest Point (ICP) algorithm. With this, we can benefit from geometry-aware properties of ICP, while at the same time obtaining an uncertainty estimate for the registration result.
The ICP algorithm is arguably one of the most popular surface registration methods. It was originally developed for rigid alignment of point sets [4, 3]. ICP iteratively estimates the point-to-point correspondences between two surfaces and then computes a transform based on the established correspondence. Finally, this transform is applied to the reference surface. The algorithm has later been modified for non-rigid surface registration . In general, standard ICP is very efficient and produces good results. However, a bad initialisation may lead the algorithm to local optima from which it is unable to recover. Moreover, it is impossible to estimate the correspondence uncertainty of the registration result. Multiple extensions have been proposed to make the algorithm more robust [1, 8, 18]. An extensive review of different ICP methods can be found in . The Coherent Point drift method  is a probabilistic alternative to non-rigid ICP. Yet, it does not provide an uncertainty estimate for the registration result. In  they improve the ICP method based on Simulated Annealing and MCMC. Their method is a robust version of the ICP algorithm, which is able to find the global optimal rigid registration of point clouds. They do, however, not measure the registration uncertainty, nor are they able to perform non-rigid registration.
We formulate the surface registration problem as approximating the posterior distribution over all possible instances of point-to-point correspondence given a target surface. With this approach, the registration uncertainty can directly be obtained from the posterior. The MH algorithm is used to sample surface registrations from the posterior distribution. Our method can avoid local optima and aims to capture the full posterior distribution of registrations. A similar approach has previously been applied in [21, 12] to estimate the uncertainty in non-rigid image registration. MH has also been used in  to fit an Active Shape Model to images and in  to fit a Morphable Face Model to an image, both of which only make use of the framework to avoid local optima and to easily integrate different likelihood terms. The main problem with the MH approach is the commonly used random-walk approach, which suffers from very long convergence times when working in high-dimensional parameter spaces. To overcome this problem, informed proposal distributions can be designed to perform directed sample updates. In 
, a Bayesian Neural Network is estimated to produce informed samples for the MH framework in the case of 3D face reconstruction from a 2D image. This, however, requires training a neural network for each class of shapes to be registered. In they propose to combine the local random-walk with an image dependent global proposal distribution. This image dependent distribution is, however, not directly transferable to the problem of surface registration.
In this paper, we show how the standard non-rigid ICP algorithm can be incorporated into the MH algorithm to take geometric information into account. We introduce a novel ICP-based proposal that can be used within the MH algorithm. Our proposal is able to make informed updates while still maintaining the beneficial convergence properties of the MH algorithm. To propose probabilistic surface deformations, we make use of a Gaussian Process Morphable Model (GPMM) as our prior model. In , GPMMs are applied to face surface registration. However, the authors formulate the registration problem as a parametric registration problem which does not provide an uncertainty measure for the established correspondence.
We validate our method on several surface registration tasks. In our first experiment, the registration targets are reconstructions of femur bones. The biggest challenge is establishing correspondence along the long smooth surface of the femur shaft. We show that in this case, the standard non-rigid ICP algorithm consistently ends up in local optima, while our method is robust against a bad initialisation and in addition can quantify the uncertainty in this region. In the second set of experiments, we use our method to reconstruct missing data. To that end, we compute the posterior distribution of registrations for faces where the nose has been removed. Unlike ICP, we are able to give an uncertainty estimate for each individual point in the surface reconstruction.
The main contributions of this paper are:
We introduce an informed proposal for the MH framework, based on the non-rigid ICP algorithm, section 3.2.
We show that different likelihood functions (independent L2, Hausdorff, exponential, etc.) can be used with our approach, section 4.2.1.
We verify that the MH algorithm with our ICP-based proposal leads to better and more robust registration results than the standard non-rigid ICP algorithm does, section 4.2.
We show how the posterior distribution for the reconstruction of missing data can be obtained with our method, section 4.3.
In this section, we formally introduce the GPMM as presented in . We show how the analytic posterior, which we use in our ICP-proposal, is computed. Furthermore, we formulate the classical non-rigid ICP algorithm within the GPMM setting.
2.1 Gaussian Process Morphable Model (GPMM)
GPMMs are a generalisation of the classical point distribution models (PDMs) . The idea is to model the deformations, which relate a given reference shape to the other shapes within a given shape family, using a Gaussian process. More formally, let be a reference surface. We obtain a probabilistic model of possible target surfaces by setting
where the deformation field is distributed according to a Gaussian process with mean function and covariance function . Let be a set of surfaces for which correspondence to the reference surface is known. This means that we can express any surface as:
for a given . Classical PDMs are obtained by defining the mean function as:
and the covariance function as: k_PDM(x,x’) = 1n-1 ∑_i=1^n (u_i(x)-μ_PDM(x)) (u_i(x’)-μ
_PDM(x’))^T. However, GPMMs also allow us to define the mean and covariance function analytically. Different choices of covariance functions lead to different well known deformation models, such as radial basis functions, b-splines, or thin-plate-splines. To model smooth deformations, we choose a zero-mean Gaussian process with the following covariance function:
where is a Gaussian kernel
The model, as stated above, is a possibly infinite-dimensional non-parametric model. Lüthi et al., propose the use of a truncated Karhunen-Loève expansion, in order to obtain a low-rank approximation of the Gaussian process. In this representation, the Gaussian process is approximated as:
where is the number of basis functions used in the approximation and. Consequently, any deformation
is uniquely determine by a coefficient vector. Hence, we can easily express any surface as:
with associated probability
2.2 Posterior Model
GPMMs make it simple and efficient to constrain a model to match known correspondences, such as user annotations or the estimated correspondence from an ICP step. Indeed, the corresponding posterior model is again a Gaussian process, whose parameters are known in closed form.
Let be a GPMM and be the certainty of each known landmark. Every landmark on the reference surface can then be matched with its corresponding landmark on the target. The set consists of the reference landmarks and its expected deformation to match the target
with being subject to Gaussian noise . Using Gaussian process regression, we obtain the posterior model , which models the possible surface deformation that are consistent with the given landmarks. Its mean and covariance function are given by:
Here we defined , the kernel matrix and a vector of the target deformations as .
2.3 Non-rigid Iterative Closest Point (ICP)
ICP can be used to perform surface registration by iteratively applying a non-rigid transformation to a template surface, as introduced in . Here, we describe our non-rigid ICP implementation which we then compare our proposed method to.
To register a reference surface to a target surface using ICP, two steps are iterated a predefined number of times or until convergence. To start with, the current surface is initialised as the reference surface :
Establish correspondence: for each landmark in the current surface , find the closest point on the target surface , .
We estimate a new surface as the mean of the posterior model described in the former section following . The noise term is the uncertainty of the correspondence established in step 1. The covariance for the noise term needs to be set such that is close to , as the ICP algorithm has no notion of correspondence uncertainty.
We can use a statistical model learned from data or an analytically defined model with a manually specified smoothness level as described in . In practice, the ICP algorithm often finds wrong correspondences if the size of and greatly differs. This is especially problematic in the case of elongated thin structures. It is, therefore, useful to also establish the correspondence from to from time to time. In our implementation, we randomly switch the direction (target to reference or reference to target).
In our framework we formulate the registration problem as Bayesian inference, where the posterior distribution of the parametersgiven the target surface is obtained as:
The prior defined in eq. 8 pushes the solution towards a more likely shape given the GPMM space, by penalising unlikely shape deformations. The likelihood term can easily be customised with different distance measures and probability functions, depending on the application goal at hand. We are usually interested in modelling the distance between two surfaces, for which we can use the independent point evaluator likelihood
The distance between the -th point and its closest point on the surfaceis the observation noise of the points of our target surface.
3.1 Approximating the Posterior Distribution
The posterior distribution defined in eq. 11 can unfortunately not be obtained analytically. Yet, we can compute the unnormalised density value for any shape described by . This allows us to use the Metropolis-Hastings algorithm  to generate samples from the posterior distribution in the form of a Markov-chain. The MH algorithm is summarised in algorithm 1.
A general way to explore the parameter space is to use a random-walk proposal, i.e. a Gaussian update distribution in the parameter space
We usually combine differently scaled distributions, each with a specified , to allow for both local and global exploration of the posterior distribution. For each proposal, one distribution is chosen at random.
A random-walk in the parameter space of the model is time-consuming as it usually is high-dimensional. Instead, we propose to accelerate convergence by using informed proposals. In order for the proposal to reach a unique stationary distribution, we have to be able to compute the transition probability, which requires the proposal to be stochastic. We propose the ICP-proposal, which takes the geometry of the model and the target into account to guide the proposed change. Internally, a posterior model based on estimated correspondences is used to propose randomised informed samples.
From a current state , we propose an update by executing the following steps:
Sample points on the current model instance .
For every point , find the closest point on the target surface .
Construct the list of corresponding landmark pairs according to eq. 9 and define the noise term for each correspondence pair.
Compute the closed form analytic posterior (eq. 10) based on the corresponding landmarks.
Take a random sample from the posterior model .
with being a step-length. If a small step-length is chosen, the current proposal is only slightly adjusted in the direction of the surface, resulting in a locally restricted step. With a step size of 1.0, the proposed solution is an independent sample from the posterior in eq. 10.
The noise in step 3 is modelled as a multivariate normal distribution with low variance along the normal direction and high variance along the surface. The variance at each point in is computed by:
where is the surface normal at the position in the mesh and and are perpendicular vectors to the normal. The variance along the vectors is set to and . This noise term ensures that the posterior model from step 4 takes into account that the correspondence, especially in flat regions, is not well defined along the surface.
Step 1 and 2 of the ICP-proposal is exactly the same as in the deterministic ICP section 2.3. In the ICP-proposal, we also randomly chose the sampling direction to escape local optima faster.
3.2.1 Computing the transition probability
For each new proposal from the ICP-proposal distribution, we need to compute the transition probability as part of the acceptance threshold (see algorithm 1 step 4). The transition probability is equal to the probability of sampling from the analytically computed posterior model , computed in step 4 of the ICP-proposal. For , the transition probability is computed in the same way. We solve eq. 14 for after swapping and and evaluate its likelihood under the posterior distribution , based on . The likelihood computed in the acceptance threshold is evaluated as stated in eq. 12.
In the following, we perform a series of registration experiment on surfaces of femur bones as well as a reconstruction experiment of missing data on face scans. The convergence properties of the ICP-proposal in comparison to random-walk is shown and we demonstrate how the likelihood function can be changed to allow for different registration properties. We compare the standard non-rigid ICP with our ICP-proposal for the registration of femurs. Finally, we show how the posterior of a missing part can be computed without the need for defining point correspondence. We perform all of our experiments with the open source library Scalismo111https://github.com/unibas-gravis/scalismo.
For the femur experiments, we use 50 healthy femur meshes extracted from computer tomography (CT) images222Available via the SICAS Medical Image Repository .. Each surface is complete, i.e. no holes or artefacts. This setting is optimal for the standard ICP algorithm described in section 2.3 and therefore serves as a fair comparison to our probabilistic implementation. For the face experiment, we make use of the publicly available Statistical Face Model (BFM) , together with 10 available face scans that are excluded from the model.
4.1 Convergence Comparison
In this experiment, we compare the convergence of the sampling method with our ICP-proposal and the standard random-walk. The ICP-proposal configuration is mentioned together with its definition in section 3.2. For the random-walk, we use a mixture of the proposals defined in eq. 13, with being set to six different levels, from to , and all six proposal distributions equally sampled from. In fig. 1, the convergence time of the standard random-walk and the ICP-proposal is shown. The experiment is performed with a GPMM with a low-rank approximation of 50 components. We randomly initialise the model parameters and start 5 registrations in parallel. As expected, our proposal leads to much faster convergence.
In fig. 2, we see run-times from a similar femur registration experimental setup, now with a GPMM with 100 components. The method with our ICP-proposal converges in the same number of iterations as for the 50 component experiment but takes 2.5 times longer in time. Again, our method with the ICP-proposal converges faster while the random-walk does not reach the same registration quality even after 1M samples.
4.2 Registration Comparison - ICP vs Sampling
In this experiment, we compare the ICP method with our probabilistic method using the ICP-proposal. We define a femur GPMM from a single femur bone (not contained in the test set), with both local and global smooth deformation kernels, and approximate it with 50 basis functions, see section 2.1. From the model, we sample 100 random meshes as an initial starting point for the registration. From the 50 femurs used as targets, we visualise a representative part of our registration results in fig. 4. The naming scheme combines the number of the femur (0 to 49), the registration method (ICP or MH with ICP) as well as the likelihood function which was used (). The box plot shows the variation of the average surface distances from all 100 registrations of each target. For the ICP method, we perform 1000 iterations. For our sampling based method, we similarly draw 1000 samples from the posterior and choose the one with highest unnormalised posterior value as the registration. The chain usually converges within 100-200 samples as shown in fig. 1, but we continue sampling for the exploration of the posterior.
We observe that the ICP-proposal is consistently better than even the best ICP registration and at the same time provides much less fluctuation in the quality. For timing comparison, the sampling method takes approximately 2.5 times longer per iteration than the standard ICP. The time increase is due to the calculation of the likelihood and the transition probability for every proposed sample.
In fig. 3, we show the uncertainty of the established correspondence of a result from the comparison experiment. Depicted is the uncertainty of individual points from the posterior distribution from a single registration. Note the high values of the uncertainties along the shaft, which indicate that the established correspondence is less reliable in that region. On the other hand, the uncertainty is very low outside the shaft region. We also observe that there is no variance along the normals, which indicates that the uncertainty is purely in the shift of correspondence.
4.2.1 Alternative Hausdorff distance likelihood
Instead of computing the posterior distribution based on the distance metric, we can register for a better Hausdorff distance  by changing the likelihood to the function:
with being the Hausdorff distance between the two meshes and
being the exponential distribution with pdf. The rate parameter is in our experiment set to . In fig. 5, we compare the registration results based on their Hausdorff distance, and we compare results from sampling using the likelihood, eq. 12, and the Hausdorff likelihood, eq. 16. As expected we can successfully focus the registration to avoid large Hausdorff distances. The distance difference for the likelihood experiment is shown in fig. 6 and shows that while optimising for the Hausdorff distance the average surface distance is increasing only slightly.
4.2.2 Discussion of drawbacks from ICP
The main problem with ICP is that it cannot recover from local optima. If the algorithm finds the closest point on the wrong part of the target, we end up with a wrong registration. In fig. 7, we show an example of the target mesh, the reference mesh and the final registrations. The ICP method has no problem to get the overall length of the bone correct but fails to correctly establish correspondence along the shaft. The sampling method does a better job when registering the reference mesh which is much different from the target mesh.
4.3 Posterior Estimation for Missing Data
In this experiment, we use the BFM model to estimate the posterior distribution of noses as well as the correspondence uncertainty of the available parts. For each of the face scans, we clip away the nose, which is then the part to be reconstructed. In fig. 8, we see that there is no perfect overlap between the face model and the scan as in the femur experiment. Therefore, we need to change our likelihood function to be able to adjust for a possible changing number of correspondence points during the registration. We use the Collective average likelihood introduced in  and extend it with the Hausdorff likelihood. The idea behind this, is to be able to obtain a close fit on average, while also penalising any far away individual points P(Γ_T—→α) = N(d_CL(Γ_T, Γ[→α]); 0, σ^2) ⋅Exp( d_H(Γ_T, Γ[→α]), λ), where
Here has the closest point on the target and being the number of landmarks in the reference mesh . On a technical note, we need to filter away all the points which have a closest point lying on a surface boundary. The filtering of the correspondences is necessary as we register a partial surface. For example, the points from the model’s nose would be mapped to the closest points on the border of the region of missing data in the target, and hence disturb the result.
In fig. 9, we show the correspondence uncertainty from the posterior registration distribution. A larger uncertainty is noticed around the outer region of the face model as there is no anatomical distinctive point correspondence, which prevents the model to explore different correspondence possibilities. Likewise, we also see uncertainty in the missing part, the nose. If we only look at the uncertainty along the normal of the face surface, we see that the uncertainty is only high in the nose regions. The uncertainty maps can for instance be used to detect missing areas, or to build better SSMs by incorporating the uncertainty of the registration process.
As it is common in non-rigid registration, we do not infer the rigid alignment, and therefore assume that a good initial alignment can be obtained prior to the registration process. Fortunately, MH allows for easy integration of proposal distributions of parameters other than the shape. The proposal distribution eq. 13 can therefore easily be extended to include translation, rotation, scaling as shown in , or texture, illumination and camera position as in . In the burn-in phase, the rigid parameter integration would allow us to adjust for a bad initial alignment and after the burn-in phase, it would allow us to explore a broader posterior distribution, where the registration difference in the rigid difference is included.
In our experiments, we have either purely made use of the ICP-proposal or the random-walk proposal. A mixture of different levels of these proposals can easily be combined in the case a broader proposal distribution is wanted while still maintaining the awareness of the geometry.
The main drawback of the sampling approach is its speed. While it is not usable within real-time applications, offline analysis can still profit form our approach.
Our integration of the ICP algorithm and the ICP-proposal has been made publicly available online333https://github.com/unibas-gravis/icp-proposal.
In this paper, we presented our robust, probabilistic registration framework. Our main contribution is the informed proposal for the Metropolis-Hastings algorithm together with a thorough evaluation. Our informed proposal integrates ICP in the update step, which results in faster convergence. This makes it possible to work in a high-dimensional model space that would be difficult to explore with pure random-walk. Hence, we combine the convergence properties of standard ICP with the probabilistic features of MCMC sampling methods, allowing us to obtain an estimate for the correspondence uncertainty along with the registration result. Moreover, our method allows for different likelihood terms to be used, while the choice is restricted to the norm in standard ICP. We experimentally showed that our approach is robust against bad initialisation and leads to better registration results than non-rigid ICP. In the case of missing data, the posterior over possible reconstructions can be estimated with our method. Thus our method can provide robust registration results for critical tasks such as surface reconstruction in the medical domain where an estimate of correspondence uncertainty is required for surgical decision making.
-  B. Amberg, S. Romdhani, and T. Vetter. Optimal step nonrigid icp algorithms for surface registration. In , pages 1–8. IEEE, 2007.
-  N. Aspert, D. Santa-Cruz, and T. Ebrahimi. Mesh: Measuring errors between surfaces using the hausdorff distance. In Proceedings. IEEE International Conference on Multimedia and Expo, volume 1, pages 705–708. IEEE, 2002.
-  P. J. Besl and N. D. McKay. Method for registration of 3-D shapes. In Sensor Fusion IV: Control Paradigms and Data Structures, volume 1611, pages 586–607. International Society for Optics and Photonics, 1992.
-  Y. Chen and G. Medioni. Object modelling by registration of multiple range images. Image and vision computing, 10(3):145–155, 1992.
-  J. Feldmar and N. Ayache. Rigid, affine and locally affine registration of free-form surfaces. International journal of computer vision, 18(2):99–119, 1996.
-  T. Gerig, A. Morel-Forster, C. Blumer, B. Egger, M. Lüthi, S. Schönborn, and T. Vetter. Morphable face models-an open framework. arXiv preprint arXiv:1709.08398, 2017.
H. Hufnagel, X. Pennec, J. Ehrhardt, N. Ayache, and H. Handels.
Generation of a statistical shape model with probabilistic point correspondences and the expectation maximization-iterative closest point algorithm.International journal of computer assisted radiology and surgery, 2(5):265–273, 2008.
-  H. Hufnagel, X. Pennec, J. Ehrhardt, H. Handels, and N. Ayache. Point-based statistical shape models with probabilistic correspondences and affine em-icp. In Bildverarbeitung für die Medizin 2007, pages 434–438. Springer, 2007.
-  V. Jampani, S. Nowozin, M. Loper, and P. V. Gehler. The informed sampler: A discriminative approach to bayesian inference in generative computer vision models. Computer Vision and Image Understanding, 136:32–44, 2015.
-  M. Kistler, S. Bonaretti, M. Pfahrer, R. Niklaus, and P. Büchler. The virtual skeleton database: an open access repository for biomedical research and collaboration. Journal of medical Internet research, 15(11):e245, 2013.
-  A. Kortylewski, M. Wieser, A. Morel-Forster, A. Wieczorek, S. Parbhoo, V. Roth, and T. Vetter. Informed mcmc with bayesian neural networks for facial image analysis. arXiv preprint arXiv:1811.07969, 2018.
-  L. Le Folgoc, H. Delingette, A. Criminisi, and N. Ayache. Quantifying registration uncertainty with sparse bayesian modelling. IEEE transactions on medical imaging, 36(2):607–617, 2016.
-  H. Liu, T. Liu, Y. Li, M. Xi, T. Li, and Y. Wang. Point cloud registration based on mcmc-sa icp algorithm. IEEE Access, 2019.
-  M. Lüthi, T. Gerig, C. Jud, and T. Vetter. Gaussian Process Morphable Models. IEEE Transactions on Pattern Analysis and Machine Intelligence, pages 1–1, 2018.
-  J. Ma, F. Lin, J. Honsdorf, K. Lentzen, S. Wesarg, and M. Erdt. Weighted robust pca for statistical shape modeling. In International Conference on Medical Imaging and Augmented Reality, pages 343–353. Springer, 2016.
-  A. Morel-Forster, T. Gerig, M. Lüthi, and T. Vetter. Probabilistic fitting of active shape models. In International Workshop on Shape in Medical Imaging, pages 137–146. Springer, 2018.
-  A. Myronenko and X. Song. Point set registration: Coherent point drift. IEEE transactions on pattern analysis and machine intelligence, 32(12):2262–2275, 2010.
-  Y. Pan, B. Yang, F. Liang, and Z. Dong. Iterative global similarity points: A robust coarse-to-fine integration solution for pairwise 3d point cloud registration. In 2018 International Conference on 3D Vision (3DV), pages 180–189. IEEE, 2018.
-  F. Pomerleau, F. Colas, R. Siegwart, et al. A review of point cloud registration algorithms for mobile robotics. Foundations and Trends® in Robotics, 4(1):1–104, 2015.
-  P. Risholm, A. Fedorov, J. Pursley, K. Tuncali, R. Cormack, and W. M. Wells. Probabilistic non-rigid registration of prostate images: modeling and quantifying uncertainty. In 2011 IEEE International Symposium on Biomedical Imaging: From Nano to Macro, pages 553–556. IEEE, 2011.
-  P. Risholm, S. Pieper, E. Samset, and W. M. Wells. Summarizing and visualizing uncertainty in non-rigid registration. In International Conference on Medical Image Computing and Computer-Assisted Intervention, pages 554–561. Springer, 2010.
-  C. Robert and G. Casella. Monte Carlo statistical methods. Springer Science & Business Media, 2013.
-  S. Schönborn, B. Egger, A. Morel-Forster, and T. Vetter. Markov chain monte carlo for automated face image analysis. International Journal of Computer Vision, 123(2):160–183, 2017.