A Closest Point Proposal for MCMC-based Probabilistic Surface Registration

TL;DR

This paper introduces a probabilistic approach to non-rigid surface registration using a novel closest-point proposal for MCMC, enabling better shape modeling, uncertainty visualization, and improved convergence over traditional methods.

Contribution

It presents a new closest-point proposal for MCMC-based surface registration, enhancing convergence and allowing flexible likelihood modeling in a probabilistic framework.

Findings

Improved convergence speed over random-walk strategies

Ability to visualize uncertainty in surface correspondences

Effective estimation of missing data without fixed point correspondences

Abstract

We propose to view non-rigid surface registration as a probabilistic inference problem. Given a target surface, we estimate the posterior distribution of surface registrations. We demonstrate how the posterior distribution can be used to build shape models that generalize better and show how to visualize the uncertainty in the established correspondence. Furthermore, in a reconstruction task, we show how to estimate the posterior distribution of missing data without assuming a fixed point-to-point correspondence. We introduce the closest-point proposal for the Metropolis-Hastings algorithm. Our proposal overcomes the limitation of slow convergence compared to a random-walk strategy. As the algorithm decouples inference from modeling the posterior using a propose-and-verify scheme, we show how to choose different distance measures for the likelihood model. All presented results are fully reproducible using publicly available data and our open-source implementation of the registration framework.

Equations

$$\Gamma =\{x+u(x)|x\in {\Gamma }_R\}$$

$${\Gamma }_i=\{x+u_i(x)|x\in {\Gamma }_R\}.$$

$$\displaystyle\begin{array}[]{rccl}\mu_{\text{PDM}}(x)&=&\frac{1}{n}&\sum_{i=1}^{n}u_{i}(x)\\ \\ k_{\text{PDM}}(x,x^{\prime})&=&\frac{1}{n-1}&\sum_{i=1}^{n}(u_{i}(x)-\mu_{\text{PDM}}(x))(u_{i}(x^{\prime})-\mu_{\text{PDM}}(x^{\prime}))^{T}.\end{array}$$

$$k(x,x^{\prime })$$

$$g(x,x^{\prime })$$

$$u[\vec{\alpha }](x)=\mu (x)+\sum _{i=1}^r{\alpha }_i\sqrt{{\lambda }_i}{\varphi }_i(x),{\alpha }_i\sim \mathcal{N}(0,1)$$

$$\Gamma [\vec{\alpha }]=\{x+\mu (x)+\sum _{i=1}^r{\alpha }_i\sqrt{{\lambda }_i}{\varphi }_i(x)|x\in {\Gamma }_R\}$$

$$p(\Gamma [\vec{\alpha }])=p(\vec{\alpha })={(2\pi )}^{-\frac{r}{2}}{\exp (-\Vert \vec{\alpha }\Vert )}^2.$$

$$L=\{(l_R^1,l_T^1-l_R^1),\dots ,(l_R^n,l_T^n-l_R^n)\}=\{(l_R^1,{\hat{u}}^1),\dots ,(l_R^n,{\hat{u}}^n)\},$$

$$\begin{array}{rl}{\mu }_p(x)& =\mu (x)+K_X(x{)}^T(K_{XX}+ϵ{)}^{-1}\hat{U}\\ k_p(x,x^{\prime })& =k(x,x^{\prime })+K_X(x{)}^T(K_{XX}+ϵ{)}^{-1}{K}_X(x^{\prime }).\end{array}$$

$$P(\vec{\alpha }|{\Gamma }_T)=\frac{P({\Gamma }_T|\vec{\alpha })P(\vec{\alpha })}{\int P({\Gamma }_T|\vec{\alpha })P(\vec{\alpha })d\vec{\alpha }}.$$

$$P({\Gamma }_T|\vec{\alpha })=\prod _{i=1}^n\mathcal{N}(d_{l^2}({\Gamma }_T^i,\Gamma [\vec{\alpha }{]}^i);0,{\sigma }_{l^2}^2),$$

$$P({\Gamma }_T|\vec{\alpha })=Exp(d_H({\Gamma }_T,\Gamma [\vec{\alpha }]);{\lambda }_H)$$

$$Q({\vec{\alpha }}^{\prime }|\vec{\alpha })\sim \mathcal{N}(\vec{\alpha },{\sigma }_l).$$

$${\vec{\alpha }}^{\prime }=\vec{\alpha }+d({\vec{\alpha }}_o-\vec{\alpha })$$

$${\Sigma }_{s_i}=[\vec{n},{\vec{v}}_1,{\vec{v}}_2]\left[\begin{array}{ccc}{\sigma }_n^2& 0& 0\\ 0& {\sigma }_v^2& 0\\ 0& 0& {\sigma }_v^2\end{array}\right][\vec{n},{\vec{v}}_1,{\vec{v}}_2{]}^T,$$

$$P({\Gamma }_T|\vec{\alpha })\propto \mathcal{N}(d_{CL}({\Gamma }_T,\Gamma [\vec{\alpha }]);0,{\sigma }_{CL}^2)\cdot Exp(d_H({\Gamma }_T,\Gamma [\vec{\alpha }]);{\lambda }_H),$$

$$d_{CL}=\frac{1}{N}\sum _{i=1}^N\Vert {\Gamma }_R^i-{\Gamma }_T^i{\Vert }^2.$$

Full text

We introduce the closest-point proposal for the Metropolis-Hastings algorithm. Our proposal overcomes the limitation of slow convergence compared to a random-walk strategy. As the algorithm decouples inference from modeling the posterior using a propose-and-verify scheme, we show how to choose different distance measures for the likelihood model.

All presented results are fully reproducible using publicly available data and our open-source implementation of the registration framework.

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 for the analysis of data where the exact correspondence is unclear, such as smooth surfaces. Within the medical area, the uncertainty of a partial data reconstruction is needed to make informed surgical decisions [27]. Uncertainty estimates can also be used to build better generalizing Point Distribution Models (PDMs) by assigning an uncertainty measure to each landmark [20, 12].††Code available at https://github.com/unibas-gravis/icp-proposal

In this paper, we propose an efficient, fully probabilistic method for surface registration based on Metropolis-Hastings (MH) sampling. Efficiency is gained by introducing a specialized proposal based on finding the closest points between a template and a target mesh. With this, we benefit from the geometry-aware proposal, while at the same time obtaining an uncertainty estimate for the registration result. We formulate the non-rigid surface registration problem as an approximation of the posterior distribution over all possible instances of point-to-point correspondences, given a target surface. With this approach, the registration uncertainty is the remaining variance in the posterior distribution. We use the MH algorithm to sample surface registrations from the posterior distribution. Our method can escape local optima and aims to capture the full posterior distribution of registrations.

Our method improves on previous works in the literature in different ways. In [28, 17], the MH algorithm is used to estimate the uncertainty in non-rigid registration. These papers are working on the image domain, and are not transferable to the surface registration setting. MH has also been used in [21] to fit an Active Shape Model to images and in [30] 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 algorithm 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 [16], a Bayesian Neural Network is learned 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 [14], local random-walk is combined with an image-dependent global proposal distribution. This image dependent distribution is, however, not directly transferable to the problem of surface registration. We perform registration by warping a single template mesh to a target surface. No training data is required with our method, whereas most state-of-the-art neural-network-based non-rigid registration methods require thousands of meshes in correspondence for training [6, 10]. These methods work well for newer non-rigid registration challenges such as the FAUST dataset [4], where the focus is to learn shape articulation from a training dataset. This is not something we would advocate using our method for. Our method targets settings with no, or limited, available training data and adds to existing methods by providing an uncertainty measure, which is especially important within the medical domain.

Even though newer registration methods exist, Iterative Closest Point (ICP) and Coherent Point Drift (CPD) are still most commonly used for surface registration scenarios without available training data.222In this paper, we focus on non-rigid registration. We, therefore, refer to their non-rigid versions whenever ICP or CPD is mentioned. The ICP algorithm was originally developed for rigid alignment of point sets [5, 3]. ICP iteratively estimates the point-to-point correspondences between two surfaces and then computes a transformation based on the established correspondence. Finally, this transformation is applied to the reference surface. The algorithm has later been modified for non-rigid surface registration [7]. In general, standard ICP is very efficient and produces good results. However, a bad initialization 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 in scenarios such as missing data or large articulation deformation differences [1, 13, 25]. In [18], the ICP method is improved using Simulated Annealing and MCMC. Their method is a robust version of the ICP algorithm, which can find the global optimal rigid registration of point clouds. They do not, however, measure the registration uncertainty, nor are they able to perform non-rigid registrations. An extensive review of different ICP methods can be found in [26]. The CPD method [23, 22] is a probabilistic alternative to ICP. Yet, it does not provide an uncertainty estimate for the registration result.

As it is common in non-rigid registration, we do not infer the rigid alignment. Fortunately, MH allows for easy integration of proposal distributions of parameters other than the shape. The proposal distribution can therefore easily be extended to include translation, rotation, scaling as shown in [21], or texture, illumination, and camera position as in [30].

In this paper, we show how closest-point information can be incorporated into the MH algorithm to take geometric information into account. We introduce a novel closest-point-proposal (CP-proposal) to use within the MH algorithm. Our proposal can make informed updates while maintaining the theoretical 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 [8], 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. Several alternatives already exist to the random-walk proposal, such as MALA [9], Hamiltonian [24] or NUTS [11]. While the mentioned proposals work well in lower-dimensional spaces and for smoother posteriors, we experienced that they get computationally demanding in high dimensional spaces and have problems when the posterior is far from smooth.

In our experiments, we register femur bones, where the biggest challenge is establishing correspondence along the long smooth surface of the femur shaft. In another experiment, 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 or CPD, we can give an uncertainty estimate for each point in the surface reconstruction. Furthermore, we show how the standard non-rigid ICP algorithm can end up in local optima, while our method consistently provides good registrations and can also quantify the correspondence uncertainty. The three main contributions of this paper are:

• •

We introduce a theoretically sound, informed closest-point-proposal for the MH framework, with a well-defined transition ratio, section 3.2.

• •

We show the usefulness of the posterior distribution of registrations for completing partial shapes and for building better generalizing PDMs, section 4.2.

• •

We demonstrate that the MH algorithm with our proposal leads to better and more robust registration results than the standard non-rigid ICP and CPD algorithms, section 4.3.

2 Background

In this section, we formally introduce the GPMM as presented in [19], and we show how the analytic posterior, which we use in our CP-proposal, is computed.

2.1 Gaussian Process Morphable Model (GPMM)

GPMMs are a generalization of the classical point distribution models (PDMs) [19]. 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 $\Gamma_{R}\subset\mathbb{R}^{3}$ be a reference surface. We obtain a probabilistic model of possible target surfaces $\Gamma$ by setting

[TABLE]

where the deformation field $u$ is distributed according to a Gaussian process $u\sim GP(\mu,k)$ with mean function $\mu$ and covariance function $k:\Gamma_{R}\times\Gamma_{R}\to\mathbb{R}^{3\times 3}$. Let $\Gamma_{1},\ldots,\Gamma_{n}$ be a set of surfaces for which correspondence to the reference surface $\Gamma_{R}$ is known. From this, it follows that we can express any surface $\Gamma_{i}$ as:

[TABLE]

Classical PDMs define the mean and covariance function as:

[TABLE]

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:

[TABLE]

where $I_{3}$ is the identity matrix and $g(x,x^{\prime})$ is a Gaussian kernel.

The model, as stated above, is a possibly infinite-dimensional non-parametric model. In [19], they propose to use the truncated Karhunen-Loève expansion to obtain a low-rank approximation of the Gaussian process. In this representation, the Gaussian process $GP(\mu,k)$ is approximated as:

[TABLE]

where $r$ is the number of basis functions used in the approximation and $\lambda_{i},\phi_{i}$ are the i-th eigenvalue and eigenfunction of the covariance operator associated with the covariance function $k$. Consequently, any deformation $u$ is uniquely determined by a coefficient vector $\vec{\alpha}=(\alpha_{1},\ldots,\alpha_{r})$. Hence, we can easily express any surface $\Gamma$ as:

[TABLE]

with associated probability

[TABLE]

2.2 Analytical 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 taking the closest point. Indeed, the corresponding posterior model is again a Gaussian process, whose parameters are known in closed form.

Let $u\sim GP(\mu,k)$ be a GPMM and $\epsilon\sim\mathcal{N}(0,\Sigma)$, $\Sigma=\sigma_{noise}I_{3}$ be the certainty of each known landmark. Every landmark $l_{R}$ on the reference surface can then be matched with its corresponding landmark $l_{T}$ on the target. The set $L$ consists of the $n$ reference landmarks and its expected deformation to match the target

[TABLE]

with $\hat{u}$ being subject to Gaussian noise $\epsilon$. Using Gaussian process regression, we obtain the posterior model $u_{p}\sim GP(\mu_{p},k_{p})$, which models the possible surface deformations that are consistent with the given landmarks. Its mean and covariance function are given by:

[TABLE]

Here we defined $K_{X}(x)=\left(k(l_{R}^{i},x)\right)_{i=1,\ldots,n}$, a vector of the target deformation as $\hat{U}=\left(l_{T}^{i}-l_{R}^{i}\right)_{i=1,\ldots,n}$ and the kernel matrix $K_{XX}=\left(k(l_{R}^{i},l_{R}^{j})\right)_{i,j=1,\ldots,n}$.

3 Method

We formulate the registration problem as Bayesian inference, where we obtain the posterior distribution of parameters $\vec{\alpha}$ given the target surface as:

[TABLE]

The prior probability, computed with eq. 11, pushes the solution towards a more likely shape given the GPMM space, by penalizing unlikely shape deformations. The likelihood term can easily be customized with different distance measures and probability functions, depending on the application goal at hand. We are usually interested in modeling the L2 distance between two surfaces ($d_{l^{2}}$), for which we can use the independent point evaluator likelihood:

[TABLE]

as also used in [21]. The L2 distance between the $i$-th point $\Gamma[\vec{\alpha}]^{i}\in R^{3}$ and its closest point on the surface $\Gamma_{T}$ is rated using a zero-mean normal distribution with the expected standard deviation for a good registration. The variance $\sigma_{l^{2}}^{2}$ is the observation noise of the points of our target surface. We can register for a better Hausdorff distance [2] by changing the likelihood to:

[TABLE]

with $d_{H}$ being the Hausdorff distance between the two meshes and $Exp$ being the exponential distribution with pdf ${p(d)=\lambda_{H}e^{-\lambda_{H}d}}$.

3.1 Approximating the Posterior Distribution

The posterior distribution defined in eq. 14 can unfortunately not be obtained analytically. Yet, we can compute the unnormalized density value for any shape described by $\vec{\alpha}$. This allows us to use the Metropolis-Hastings algorithm [29] to generate samples from the posterior distribution in the form of a Markov-chain. The MH algorithm is summarized 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

[TABLE]

We usually combine differently scaled distributions, each with a specified $\sigma_{l}$, to allow for both local and global exploration of the posterior distribution. For each proposal, one distribution is chosen at random.

3.2 CP-proposal

A random-walk in the parameter space of the prior GPMM model $\mathcal{M}$ is time-consuming as it usually is high-dimensional. Instead, we propose to accelerate convergence by using an informed proposal. 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 CP-proposal, which takes the geometry of the model and the target into account to guide the proposed change. Internally, we use a posterior model $\mathcal{M}_{\alpha}$ based on estimated correspondences to propose randomized informed samples. From a current state $\vec{\alpha}$, we propose an update $\vec{\alpha}^{\prime}$ by executing the following steps (visualized in fig. 1):

Sample $m$ points $\{s_{i}\}$ on the current model instance $\Gamma[\vec{\alpha}]$. 2. 2.

For every point $s_{i}$, $i\in\left[0\ ...\ m\right]$ find the closest point $c_{i}$ on the target $\Gamma_{T}$. 3. 3.

Construct the set of observations $L$ based on corresponding landmark pairs ${(s_{i},c_{i})}$ according to eq. 12 and define the noise $\epsilon_{i}\sim\mathcal{N}(0,\Sigma_{s_{i}})$ using eq. 19. 4. 4.

Compute the analytic posterior $\mathcal{M}_{\alpha}$ (eq. 13) with $L$ and $\{\Sigma_{s_{i}}\}$. 5. 5.

Get $\vec{\alpha}_{o}$ by first drawing a random shape333Sampling all ${\alpha}_{i}$ independently from $\mathcal{N}(0,1)$ and constructing the shape with eq. 10. from the posterior model $\mathcal{M}_{\alpha}$ and then projecting it into the prior model $\mathcal{M}$. 6. 6.

We generate

[TABLE]

with $d\in[0.0\ ...\ 1.0]$ being a step-length.

The noise $\epsilon$ in step 3 is modeled with low variance along the normal direction and high variance along the surface. The variance at each point $s_{i}$ in $\Gamma[\vec{\alpha}]$ is computed by:

[TABLE]

where $\vec{n}$ is the surface normal at the position $s_{i}$ in the mesh and $\vec{v}_{1}$ and $\vec{v}_{2}$ are perpendicular vectors to the normal. The variances along the vectors are set to $\sigma_{n}^{2}=$3.0\text{,}{\mathrm{mm}}^{2} and $\sigma_{v}^{2}=$100.0\text{\,}{\mathrm{mm}}^{2}. This noise term ensures that the posterior model from step 4 takes the uncertain correspondence along the surface into account, which is not well defined in flat regions.

If a small step-length is chosen in step 6, the current proposal is only adjusted slightly in the direction of the target surface, resulting in a locally restricted step. With a step size of $1.0$, the proposed sample is an independent sample from the posterior in eq. 13.

In practice, closest-point based updates often find wrong correspondences if the sizes of $\Gamma$ and $\Gamma_{T}$ greatly differ. This is especially problematic in the case of elongated thin structures. It is, therefore, useful also to establish the correspondence from $\Gamma_{T}$ to $\Gamma$ from time to time.

3.2.1 Computing the Transition Probability

For each new proposal $\vec{\alpha}^{\prime}$ from the CP-proposal distribution, we need to compute the transition probability as part of the acceptance threshold (see algorithm 1 step 4). The transition probability $q(\vec{\alpha}^{\prime}|\vec{\alpha})$ is equal to the probability of sampling the shape corresponding to $\vec{\alpha}_{o}$ from the posterior model $\mathcal{M}_{\alpha}$, computed in step 4 of the CP-proposal. For $q(\vec{\alpha}|\vec{\alpha}^{\prime})$, the transition probability is computed in the same way. We solve eq. 18 for $\vec{\alpha}^{\prime}_{o}$ after swapping $\vec{\alpha}$ and $\vec{\alpha}^{\prime}$ and evaluate the corresponding shape likelihood under the posterior distribution $\mathcal{M}_{\alpha^{\prime}}$.

4 Experiments

In the following, we perform registration experiments on surfaces of femur bones as well as a reconstruction experiment of missing data on face scans. For the face experiment, the face template and the face GPMM from [8] are used together with 10 available face scans. For the femur experiments, we use 50 healthy femur meshes extracted from computed tomography (CT) images444Available via the SICAS Medical Image Repository [15].. Each surface is complete, i.e. no holes or artifacts. This setting is optimal for the standard ICP algorithm and therefore serves as a fair comparison to the CPD algorithm and our probabilistic implementation. The CPD experiments are performed with the MATLAB code from [22] and all other experiments with the Scalismo555https://scalismo.org library.

4.1 Convergence Comparison

We compare the convergence properties of our CP-proposal and random-walk. The CP-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. 17, with $\sigma_{l}$ being set to six different levels, from $1.0\text{\,}\mathrm{mm}$ to $0.01\text{\,}\mathrm{\SIUnitSymbolMicro m}$, and all six proposal distributions equally likely to be sampled from.

In fig. 2, the convergence time of the standard random-walk and the CP-proposal is shown. The experiment is performed with a GPMM with a low-rank approximation of 50 components, see section 2.1. We randomly initialize the model parameters and start 5 registrations in parallel. As expected, our proposal leads to much faster convergence. In fig. 3 we see a posterior plot comparison of the two proposals. Notice how less likely samples are often accepted, which makes it different from deterministic methods such as ICP and CPD.

4.2 Posterior Estimation of Missing Data

We use a face GPMM to estimate the posterior distribution of noses from partial face data, where the nose has been removed. In fig. 4, we see that there is no perfect overlap between the face model and the scan.

Therefore, we need to change our likelihood function to adjust for a possible changing number of correspondence points during the registration. To obtain a close fit on average, while also penalizing far away points, we use the collective average likelihood introduced in [30] and extend it with the Hausdorff likelihood,

[TABLE]

where

[TABLE]

Here the closest point of $\Gamma_{R}^{i}$ on the target is $\Gamma_{T}^{i}$ and $N$ is the number of landmarks in the reference mesh $\Gamma_{R}$.

On a technical note, using ICP to predict correspondence in a missing data scenario maps all points from the reference, which are not observed in the target to the closest border. To counter this effect, we filtered away all predicted correspondences, where the predicted target point is part of the target surface’s boundary. This is also done in e.g. [1].

In fig. 4, we show the correspondence uncertainty from the posterior registration distribution. Our method infers a larger correspondence uncertainty in the outer region. However, as the surface is observed, the uncertainty is low in the direction of the face surface normals but high within the surface. This is because there is no anatomical distinctive feature to predict the correspondence more precisely. High uncertainty is also inferred on the nose, where data is missing. In contrast to the outer region, the uncertainty in the direction of the normal of the reconstructed nose surface is large. This shows that uncertainty visualization can be used to detect missing areas or to build better PDMs by incorporating the uncertainty of the registration process as demonstrated next.

4.2.1 PDMs from Posterior Distributions

In this experiment, we show how to benefit from probabilistic registration when building new PDMs from a small dataset with partially missing data, which is common in medical imaging. If much more data is available, the influence of using the posterior is reduced as the variability gained from the probabilistic framework can be gained by including a lot more data. In section 4.2, we registered the 10 available face scans by computing the posterior distribution of reconstructions for each target. In fig. 5(b) we show samples from the posterior distribution of nose reconstructions.

In fig. 5(a), we compare two different PDMs’ generalization abilities, i.e. the capability to represent unseen data [31]. One PDM is built following the classical approach of only using the most likely registrations. The other PDM is built from 100 random samples from each of the 10 posterior distributions. We compute the average generalization of all 10 PDMs built using a leave-one-out scheme. The plot shows that the PDM built from the posteriors generalizes better.

4.3 Registration Accuracy - ICP vs CPD vs CP-proposal

In this experiment, we compare the best sample, also known as Maximum a posteriori (MAP) from our probabilistic method, with the ICP and CPD methods. We use a femur GPMM approximated with 200 basis functions. From the model, we sample 100 random meshes as starting points for the registration. We, therefore, end up with 100 registrations for each target. In fig. 7 we show the summary of all 5000 registrations (All) and some individual representative target meshes.

The naming scheme combines the registration method (ICP, CPD, or MH with CP) with the likelihood function which was used (L2). The box plot shows the variation of the average L2 surface distances from all 100 registrations of each target. As the ICP and CPD methods converge at maximum 100 iterations, we also restrict our sampling method to 100 samples. The chain usually converges within 100-300 samples as shown in fig. 2.

The CP-proposal is consistently better than the ICP and CPD registrations, and at the same time provides much less fluctuation in the quality. The few outliers for our method are cases where the chain has not converged. In fig. 6, we show the uncertainty of the established correspondence of a registration from 1000 samples (300 samples for the burn-in phase). Depicted is the uncertainty of individual points from the posterior distribution from a single registration. Note the high uncertainty values along the shaft, which indicate that the established correspondence is less reliable in that region. No variance along the normals indicates that the uncertainty is purely correspondence shift within the surface.

4.3.1 Alternative Hausdorff Distance Likelihood

In fig. 8(a), we compare the registration results based on their Hausdorff distance, and we compare results from sampling using the L2 likelihood, eq. 15, and the Hausdorff likelihood, eq. 16. As expected, we can focus the registration to avoid large Hausdorff distances. The equivalent L2 likelihood experiment is shown in fig. 8(b) and shows that while optimizing for the Hausdorff distance, the average L2 surface distance is increasing only slightly. This demonstrates the capability to change the likelihood in our framework based on the application’s needs.

4.3.2 Drawbacks of Deterministic Methods

The main problem with ICP and CPD is that they cannot recover from local optima. If the algorithm finds the closest point on the wrong part of the target, we end up with a bad registration. In fig. 9, we show a registration result of the 3 registration methods. The ICP method can get the overall length of the bone correct but ends up with a registration, where the structure is folding around itself. The CPD approach is more robust than ICP as it preserves the topological structure of the point sets. Our method additionally provides an uncertainty map.

4.3.3 Run-time Comparison

The number of components in the low-rank approximation can be seen as regularization of the deformations. More complex local deformations can be obtained using more components. The algorithm run-time scales linearly in the number of components, with the run-time being 2.5 times slower for each time the model rank doubles. For models with rank 50, 100 and 200, the CP-proposal takes: $46\text{\,}\mathrm{s}$, $110\text{\,}\mathrm{s}$ and $275\text{\,}\mathrm{s}$. In comparison, the ICP implementation takes $30\text{\,}\mathrm{s}$, $69\text{\,}\mathrm{s}$ and $155\text{\,}\mathrm{s}$. The CPD implementation does not make use of the GPMM model and takes $75\text{\,}\mathrm{s}$ for 100 samples with our setup. While our method is slower than the two others, we still get reasonable run-times while inferring more accurate results and estimate the full posterior instead of a single estimate.

5 Conclusion

In this paper, we presented a probabilistic registration framework. Our main contribution is the informed proposal for the MH algorithm, which makes it possible to work in a high-dimensional model space that would be difficult to explore with pure random-walk. Our informed proposal integrates geometry awareness in the update step, which results in faster convergence. In the case of missing data, our method provides an estimate of the posterior over possible reconstructions. Thus our framework can provide uncertainty measures for critical tasks such as surface reconstruction in the medical domain, as required for surgical decision making. Using our framework, different likelihood terms can be combined and used, while the choice is restricted to the L2 norm in standard ICP and CPD. Finally, we showed how to build PDMs that generalize better using the posterior distribution of registrations.

5.0.1 Acknowledgements.

This research is sponsored by the Gebert Rüf Foundation under the project GRS-029/17.