1 Introduction
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 [20]. In this paper, we propose an efficient, fully probabilistic method for surface registration based on MetropolisHastings (MH) sampling. Efficiency is gained by the use of a specialised proposal based on a nonrigid Iterative Closest Point (ICP) algorithm. With this, we can benefit from geometryaware 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 pointtopoint 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 nonrigid surface registration [5]. 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 [19]. The Coherent Point drift method [17] is a probabilistic alternative to nonrigid ICP. Yet, it does not provide an uncertainty estimate for the registration result. In [13] 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 nonrigid registration.
We formulate the surface registration problem as approximating the posterior distribution over all possible instances of pointtopoint 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 nonrigid image registration. MH has also been used in [16] to fit an Active Shape Model to images and in [23] 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 randomwalk approach, which suffers from very long convergence times when working in highdimensional parameter spaces. To overcome this problem, informed proposal distributions can be designed to perform directed sample updates. In [11]
, 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
[9] they propose to combine the local randomwalk 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 nonrigid ICP algorithm can be incorporated into the MH algorithm to take geometric information into account. We introduce a novel ICPbased 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 [6], 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 nonrigid 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 nonrigid 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 ICPbased proposal leads to better and more robust registration results than the standard nonrigid 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.
2 Background
In this section, we formally introduce the GPMM as presented in [14]. We show how the analytic posterior, which we use in our ICPproposal, is computed. Furthermore, we formulate the classical nonrigid ICP algorithm within the GPMM setting.
2.1 Gaussian Process Morphable Model (GPMM)
GPMMs are a generalisation of the classical point distribution models (PDMs) [14]. 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
(1) 
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:
(2) 
for a given . Classical PDMs are obtained by defining the mean function as:
(3) 
and the covariance function as: k_PDM(x,x’) = 1n1 ∑_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, bsplines, or thinplatesplines. To model smooth deformations, we choose a zeromean Gaussian process with the following covariance function:
(4) 
where is a Gaussian kernel
(5) 
The model, as stated above, is a possibly infinitedimensional nonparametric model. Lüthi et al.
[14], propose the use of a truncated KarhunenLoève expansion, in order to obtain a lowrank approximation of the Gaussian process. In this representation, the Gaussian process is approximated as:(6) 
where is the number of basis functions used in the approximation and
are the ith eigenvalue and eigenfunction of the covariance operator associated with the covariance function
. Consequently, any deformationis uniquely determine by a coefficient vector
. Hence, we can easily express any surface as:(7) 
with associated probability
(8) 
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
(9) 
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:
(10) 
Here we defined , the kernel matrix and a vector of the target deformations as .
2.3 Nonrigid Iterative Closest Point (ICP)
ICP can be used to perform surface registration by iteratively applying a nonrigid transformation to a template surface, as introduced in [5]. Here, we describe our nonrigid 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 [14]. 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 [14]. 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).
3 Method
In our framework we formulate the registration problem as Bayesian inference, where the posterior distribution of the parameters
given the target surface is obtained as:(11) 
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
(12) 
The distance between the th point and its closest point on the surface
is rated using a zeromean normal distribution with the expected standard deviation for a good registration. The variance
is 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 MetropolisHastings algorithm [22] to generate samples from the posterior distribution in the form of a Markovchain. The MH algorithm is summarised in algorithm 1.
A general way to explore the parameter space is to use a randomwalk proposal, i.e. a Gaussian update distribution in the parameter space
(13) 
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.
3.2 ICPProposal
A randomwalk in the parameter space of the model is timeconsuming as it usually is highdimensional. 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 ICPproposal, 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 .

We generate
(14) with being a steplength. If a small steplength 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:
(15) 
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 ICPproposal is exactly the same as in the deterministic ICP section 2.3. In the ICPproposal, 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 ICPproposal 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 ICPproposal. 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.
4 Experiments
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 ICPproposal in comparison to randomwalk is shown and we demonstrate how the likelihood function can be changed to allow for different registration properties. We compare the standard nonrigid ICP with our ICPproposal 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 Scalismo^{1}^{1}1https://github.com/unibasgravis/scalismo.
For the femur experiments, we use 50 healthy femur meshes extracted from computer tomography (CT) images^{2}^{2}2Available via the SICAS Medical Image Repository [10].. 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) [6], 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 ICPproposal and the standard randomwalk. The ICPproposal configuration is mentioned together with its definition in section 3.2. For the randomwalk, 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 randomwalk and the ICPproposal is shown. The experiment is performed with a GPMM with a lowrank 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 runtimes from a similar femur registration experimental setup, now with a GPMM with 100 components. The method with our ICPproposal 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 ICPproposal converges faster while the randomwalk 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 ICPproposal. 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 100200 samples as shown in fig. 1, but we continue sampling for the exploration of the posterior.
We observe that the ICPproposal 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 [2] by changing the likelihood to the function:
(16) 
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 [23] 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
(17) 
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.
5 Discussion
As it is common in nonrigid 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 [16], or texture, illumination and camera position as in [23]. In the burnin phase, the rigid parameter integration would allow us to adjust for a bad initial alignment and after the burnin 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 ICPproposal or the randomwalk 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 realtime applications, offline analysis can still profit form our approach.
Our integration of the ICP algorithm and the ICPproposal has been made publicly available online^{3}^{3}3https://github.com/unibasgravis/icpproposal.
6 Conclusion
In this paper, we presented our robust, probabilistic registration framework. Our main contribution is the informed proposal for the MetropolisHastings 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 highdimensional model space that would be difficult to explore with pure randomwalk. 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 nonrigid 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.
References

[1]
B. Amberg, S. Romdhani, and T. Vetter.
Optimal step nonrigid icp algorithms for surface registration.
In
2007 IEEE Conference on Computer Vision and Pattern Recognition
, pages 1–8. IEEE, 2007.  [2] N. Aspert, D. SantaCruz, 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.
 [3] P. J. Besl and N. D. McKay. Method for registration of 3D shapes. In Sensor Fusion IV: Control Paradigms and Data Structures, volume 1611, pages 586–607. International Society for Optics and Photonics, 1992.
 [4] Y. Chen and G. Medioni. Object modelling by registration of multiple range images. Image and vision computing, 10(3):145–155, 1992.
 [5] J. Feldmar and N. Ayache. Rigid, affine and locally affine registration of freeform surfaces. International journal of computer vision, 18(2):99–119, 1996.
 [6] T. Gerig, A. MorelForster, C. Blumer, B. Egger, M. Lüthi, S. Schönborn, and T. Vetter. Morphable face modelsan open framework. arXiv preprint arXiv:1709.08398, 2017.

[7]
H. Hufnagel, X. Pennec, J. Ehrhardt, N. Ayache, and H. Handels.
Generation of a statistical shape model with probabilistic point correspondences and the expectation maximizationiterative closest point algorithm.
International journal of computer assisted radiology and surgery, 2(5):265–273, 2008.  [8] H. Hufnagel, X. Pennec, J. Ehrhardt, H. Handels, and N. Ayache. Pointbased statistical shape models with probabilistic correspondences and affine emicp. In Bildverarbeitung für die Medizin 2007, pages 434–438. Springer, 2007.
 [9] 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.
 [10] 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.
 [11] A. Kortylewski, M. Wieser, A. MorelForster, 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.
 [12] 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.
 [13] H. Liu, T. Liu, Y. Li, M. Xi, T. Li, and Y. Wang. Point cloud registration based on mcmcsa icp algorithm. IEEE Access, 2019.
 [14] 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.
 [15] 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.
 [16] A. MorelForster, 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.
 [17] A. Myronenko and X. Song. Point set registration: Coherent point drift. IEEE transactions on pattern analysis and machine intelligence, 32(12):2262–2275, 2010.
 [18] Y. Pan, B. Yang, F. Liang, and Z. Dong. Iterative global similarity points: A robust coarsetofine integration solution for pairwise 3d point cloud registration. In 2018 International Conference on 3D Vision (3DV), pages 180–189. IEEE, 2018.
 [19] 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.
 [20] P. Risholm, A. Fedorov, J. Pursley, K. Tuncali, R. Cormack, and W. M. Wells. Probabilistic nonrigid 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.
 [21] P. Risholm, S. Pieper, E. Samset, and W. M. Wells. Summarizing and visualizing uncertainty in nonrigid registration. In International Conference on Medical Image Computing and ComputerAssisted Intervention, pages 554–561. Springer, 2010.
 [22] C. Robert and G. Casella. Monte Carlo statistical methods. Springer Science & Business Media, 2013.
 [23] S. Schönborn, B. Egger, A. MorelForster, and T. Vetter. Markov chain monte carlo for automated face image analysis. International Journal of Computer Vision, 123(2):160–183, 2017.
Comments
There are no comments yet.