Eyes on the Prize: Improved Biological Surface Registration via Forward Propagation
Robert J. Ravier











Robert J. Ravier Department of Electrical and Computer Engineering 

              Duke University 

              robert.ravier@duke.edu           







Eyes on the Prize: Improving Biological Surface Registration via Forward Propagation



    Robert J. Ravier

    Received: date / Accepted: date
=======================================================================================




Many algorithms for surface registration risk producing significant errors if surfaces are significantly nonisometric. Manifold learning has been shown to be effective at improving registration quality, using information from an entire collection of surfaces to correct issues present in pairwise registrations. These methods, however, are not robust to changes in the collection of surfaces, or do not produce accurate registrations at a resolution high enough for subsequent downstream analysis. We propose a novel algorithm for efficiently registering such collections given initial correspondences with varying degrees of accuracy. By combining the initial information with recent developments in manifold learning, we employ a simple metric condition to construct a measure on the space of correspondences between any pair of shapes in our collection, which we then use to distill soft correspondences. We demonstrate that this measure can improve correspondence accuracy between feature points compared to currently employed, less robust methods on a diverse dataset of surfaces from evolutionary biology. We then show how our methods can be used, in combination with recent sampling and interpolation methods, to compute accurate and consistent homeomorphisms between surfaces.







§ INTRODUCTION

Shape registration is a fundamental problem in computer vision. Accurate registrations are required for many applications, such as object retrieval, statistical shape analysis, medical imaging, among others <cit.>. A plethora of algorithms for registering pairs of shapes have been proposed, many of which can be directly formulated as minimizing some well-defined objective function (see <cit.>) for a small sample of such works. For practical applications involving collections of shapes, it is often necessary to have a consistent set of registrations: given shapes {S_i}_i ∈ℐ, the collection of registrations {ϕ_ji}_i,j ∈ℐ is consistent if ϕ_ik∘ϕ_kj∘ϕ_ji = Id_S_i for all i,j,k ∈ℐ, where Id_S_i is the identity map on S_i. Consistency is a requirement for posthoc statistical analysis <cit.>, and is generally obtained by a refinement of an initially computed collection of registrations, either through optimization <cit.> or by directly defining a consistent collection based on composition <cit.>.



Anatomical surfaces studied in evolutionary biology, the shapes that serve as the main motivation for this work, have historically proven to be difficult for algorithms to accurately register, especially if one is interested in homeomorphism. In contrast to those considered in medical biology, surfaces considered in evolutionary biology often consist of multiple species <cit.>, which in turn often leads to high amounts of nonisometry. Though registration of nonisometric shapes is of increasing research interest <cit.>, the mechanisms of nonisometry studied in recent works are generally distinct from those in evolutionary biology. In this setting, the main difficulty stems from both the lack of complexity of shape compared to those often considered, potential differences in the amount and location of regions with high curvature, as well as the potential type of curvature itself.  Figure <ref> shows an example of a failure of the algorithm in <cit.> to register two primate molars from <cit.>, for which ground truth registrations are available. In this case, the nonisometry between the two surfaces, occurring in the yellow regions of both texture maps, has resulted in a failure to find a good registration, as represented by the extreme stretching on the left textured surface. 

The failure in Figure <ref> is perhaps unsurprising to the reader: the geometries of the top portions of both surfaces are quite distinct in the sense that some of the features on the left surface do not have obvious counterparts on the right surface. Though this stark difference naturally begs the question as to whether a meaningful correspondence is feasible to compute without significant uncertainty, we again recall from the previous paragraph that well-accepted ground truth correspondences exist. <cit.>. Correspondences for this example and other datasets within evolutionary biology are determined by experts manually selecting such features (termed landmarks <cit.>) based on both geometric and non geometric biological information from the collection as a whole. Corresponding landmarks on pairs of disparate shapes are often determined by composing correspondences along paths of intermediates, where the degree of nonisometry between any two intermediates is smaller than those of the initial shapes of interest. This procedure matches recent findings within the registration literature <cit.>, which suggest that the accuracy of registrations (and, in the case of <cit.>, rigid alignments) between disparate pairs of shapes can be improved by using compositions of registrations between similar shapes. Here, the methods in <cit.> define similar by an explicit notion of distance, while the methods in <cit.> use time difference as a similarity measure, appropriate as the data of interest consists of motion capture of humans. 

All methods mentioned in the previous paragraph yield consistent registrations amongst all shapes in a collection; this is due to explicit choices of particular sequences of intermediate for each pair of shapes that result in consistency. As noted in the beginning of the section, this is not strictly necessary: the method proposed in <cit.> allows one to obtain consistent registrations from an initial set of pairwise ones by means of solving an optimization problem. We note that the methodology proposed in <cit.> makes the intrinsic and completely reasonable assumption that the initial pairwise registrations are of sufficient quality. For disparate pairs of shapes in our setting, this requires a sequence of intermediates over which to compute and compose registrations per experimental evidence. Though one could use additional information not immediately present within the collection of shapes of interest, it is likely that such information already results in consistent registrations, defeating the point of the optimization method in this setting.



The distance-based methods for achieving consistent registrations  <cit.> may still not yield ones that are of sufficient quality between disparate shapes. Per the analysis of <cit.>, accuracy of registration was observed to need not perfectly correlate with composition over small distances. Furthermore, there are few stability guarantees in this setting: for example, minimum spanning trees, which are used in <cit.> as a means of determining consistent registrations, may wildly change topology if an additional sample is added or removed. An example of this can be found in Figure <ref>, meaning that incorporation of additional samples for analysis can cause monumental issues in registration. Addressing this issue, namely obtaining consistent, high-quality registrations in the nonisometric scenario, is the focus of our work.

An obvious way to solve this would be to improve the quality of all registrations, or at enough to generate a consistent collection via some heuristic. Based on the aforementioned experimental evidence, it is thus natural to define an optimization problem for each pair of shapes in the collection, where the domain of the objective is over the space of all possible sequences of intermediates. This is naturally posed within the framework of combinatorial optimization, for which many solution techniques exist. The main problem is the choice of objective function, which is extremely nontrivial. Though one could minimize an energy functional that depends on a given pair of shapes and the registration between them, such as the continuous Procrustes distance <cit.>, as previously mentioned, smaller values of such energies do not perfectly correlate with increased accuracy of registration <cit.>. At the same time, it is also infeasible to use direct accuracy measures as an objective function, as this implies it is known a priori, i.e. the problem is already solved. Thus, even in a quite general framework, we apparently stuck in a catch-22. The goal of this paper is to propose a mathematically principled solution to this problem by eschewing the language of optimization in favor of that of differential geometry and manifold learning.

In this paper, we propose a novel solution to the problem of registering disparate shapes in a collection. Rather than directly optimizing an objective function, we instead consider a surrogate problem of computing parallel transports along geodesics of a manifold. Geodesic parallel transports have been previously considered in the realm of registration, most notably for medical shape analysis as a means of transforming data to a common template on which to perform statistical analysis <cit.>. We will instead based our approach off recent advances within the realms of manifold learning <cit.>, in which registrations themselves are modeled as parallel transports on a fiber bundle. We will use this model to significantly refine our search space to registrations induced by eyes on the prize flows (EOP flows) within the shape space. This allows us to naturally arrive at a robust solution to our registration that is surprisingly simple to employ in practice. Namely, we will detail an algorithm, EOP-PM, that uses all information in the collection to conservatively determine which points on any two shapes are likely in correspondence, and show how one can adapt this to naturally get consistent registration of a wide class of surfaces. This is done by employing what we call the approximately cycle-consistent condition, or ACC condition, a constraint in which registrations are only made if they are very likely with respect to those corresponding to EOP flows.

The remainder of the paper is organized as follows. In Section <ref>, we go over necessary preliminaries as well as review related work. In Section <ref>, we investigate our parallel transport model in order to derive simple constraints on which to restrict our search space. In Section <ref>, we use the discussion from Section <ref> to detail a novel algorithm for registering disparate shapes given information from an entire collection, which we then adapt to compute consistent registrations in Section <ref>. We make concluding remarks in Section <ref>



§ PRELIMINARIES



 §.§ Registration

We denote by 𝒮 a collection of shapes S_1,...,S_n. For the purposes of our paper, we will be concerned with homeomorphic triangular meshes in ℝ^3, i.e. S_i = (V_i,F_i), where V_i and F_i are the set of vertices and faces respectively. Though all shapes we work with are homeomorphic to either the unit disk or the unit sphere, the ideas naturally generalize to other topologies. The astute reader will observe that many of our proposed algorithms and their implementation depend only on the point cloud structure; the mesh structure is used for initial computations as needed as well as for reparametrization as discussed in Section <ref>.

We denote by 𝒞_𝒮 a collection of precomputed pairwise registrations f_ji: S_i→ S_j, with f_ii = Id_S_i. We do not make any assumptions a priori as to the consistency or invertibility of the registrations in 𝒞_𝒮. We only assume that correspondences are everywhere defined so that an arbitrary composition is meaningful. The collection 𝒞_𝒮 can be used to define notions of dissimilarities between two shapes. For the purposes of this work, we focus on dissimilarities of the following form:


    d_𝒞_𝒮(S_i,S_j)    = 1/2( ∫_S_i g(x,f_ji(x))^2  dS_i)^1/2
       + 1/2( ∫_S_j g(x,f_ij(x))^2  dS_j)^1/2


where g ≥ 0 is a real-valued dissimilarity, and both integrals are with respect to the usual area measures on each shape. Note that d_𝒞_𝒮 is nonnegative and symmetric by definition. In many cases, functionals of forms similar to that of Equation (<ref>) can be used to obtain the initial registrations in 𝒞_𝒮 by an optimization procedure; see, for example, <cit.>. To simplify notation, we will often write d_ij in place of d_𝒞_𝒮(S_i,S_j). For the remainder of the work, we will assume that the d_𝒞_𝒮 is an actual distance function, i.e. satisfies a triangle inequality. We do not view this as a huge restriction, as one can always obtain a distance from a general class of similarity scores via diffusion maps <cit.>. 

With a notion of distance, we can naturally represent the collection of shapes 𝒮 with registrations 𝒞_𝒮 as a weighted complete directed graph Γ(𝒮,𝒞_𝒮), or Γ for short. Each vertex v_i∈Γ corresponds to the shape S_i in our collection, and each directed edge from v_i to  v_j has as a weight the formal pair (w_ij,f_ji) = (K_σ(d_ij),f_ji), where K_σ is a nonnegative kernel, such as the Gaussian kernel K_σ(d_ij) = e^- d_ij^2/σ. This particular representation has been extensively studied in the bundle diffusion maps literature (e.g. <cit.>) as a way to analyze datasets in which any two datapoints can be related to one another by means of an G-action, where G is usually a group (O(2) in <cit.>) or groupoid (registrations in <cit.>). We observe that our choice of weights on directed edges naturally yields an analogous formal pair of weights for all paths between any two vertices on Γ. More precisely, an arbitrary sequence of shapes S_i_0,...,S_i_n naturally corresponds to the formal pair 


    ( ∏_k=0^n-1 w_(i_k+1)i_k, f_i_ni_n-1∘ ... ∘ f_i_1i_0)


where ∘ denotes composition. We will often shorten this to f_γ, for γ a path in Γ.



 §.§ Differential Geometry and Fiber Bundles


To properly develop the methods proposed in this paper, we require some language from differential geometry that has not traditionally been used in applications, but has been recently proposed in <cit.>. We will model our collection of data as a fiber bundle, whose definition we take from <cit.>.


A fiber bundle is a formal 4-tuple (E,B,F,π), where E, B, and F are smooth manifolds, and π: E → B is a smooth mapping such that each b ∈ B has an open neighborhood U_b such that π^-1(U_b) is diffeomorphic to U_b× F. Here, E is known as the total manifold, B is known as the base manifold, and F is known as the fiber.


In other words, a fiber bundle is locally a product manifold. Every b ∈ B has a corresponding fiber F_b that is diffeomorphic to F. In our setting, we model our collection of shapes as fibers in a fiber bundle , with shape S_i having corresponding point b_i in the base manifold (analogous to their corresponding vertices v_i in Γ.) The distance between b_i and b_j is thus given by d_ij. 

The fiber bundle framework allows us to model registrations as parallel transports, which lie at the heart of our method. In Riemannian geometry, a smooth curve γ between two points p and q on a Riemannian manifold M naturally induces a parallel transport operator P_γ: T_pM → T_qM. Parallel transports in this setting have been used in medical imaging <cit.> to align vector field deformations of images to a common template. Our use of meshes does not allow for any direct adaptation of such methods to our setting, e.g. fiber bundles. We reserve the discussion of full details for the Appendix. The key fact that we require is a basic consequence of the defintion, which we summarize here: parallel transports on the base manifold naturally induce diffeomorphisms of the fibers <cit.>. This leads us to our ability to mathematically model computed correspondences as parallel transports. 



The utility of this modeling mechanism can be seen through holonomy. The precise definitions are out of scope for this work, but we illustrate the intuitive idea in Figure <ref>, which shows an example of the Riemannian manifold setting; we refer the interested reader to recent work on synchronization for more thorough discussions<cit.>. Here, we see that different paths result in different parallel transports of tangent vectors. Since parallel transports are registrations in our fiber bundle setting, this corresponds to inconsistency of a collection of registrations. A natural question one can ask is whether it is possible to determine a collection of paths for which the parallel transports are reasonably close to one another; this is central to our methodology. 





§ FINDING MEANINGFUL PATHS FOR PARALLEL TRANSPORT

Though our ultimate interest is in obtaining consistent registrations, we first concern ourself with registration quality. For each pair of shapes (S_i,S_j) our assumed collection of registrations naturally yield an infinite number of registrations between S_i and S_j. Indeed, every composition of the form f_jk_n∘ ... ∘ f_k_0i is a valid registration. It is thus natural to ask if there is any proxy to estimate the quality of such registrations without requiring ground truth information. In this section, we detail such a condition motivated by parallel transport.

We first immediately note that the infinite number of registrations can be trimmed to a finite (but possibly large) subset without issue. Indeed, since each such registration corresponds to a path γ in Γ, we can consider two different classes of paths: ones which contain loops, and ones which do not. The collection of paths that do not contain loops is finite as Γ only has a finite number of vertices, meaning that the collection containing loops is infinite. Recall that we assumed that the collection of registrations is not consistent. Thus, if γ̃ is a loop in Γ, the corresponding registration f_γ̃ is not necessarily the identity map; if the collection was consistent, then f_γ̃ would necessarily be the identity. This means that loops function only as error perturbations; each loop only adds error to a registration corresponding to a simple path without loops. 

We can thus safely remove all registrations with loops from consideration, leaving only a finite set left. Nevertheless, this finite set is potentially large as the number of simple paths in a graph grows factorially with the number of vertices, potentially making it prohibitively large to work with. Though one could proceed via sampling, it is natural to ask if we can further reduce the collection under consideration. A common assumption within the general area of shape analysis is that good registrations can be obtained from geodesics in shape space <cit.>. If we keep this assumption, then a reasonable way forward is to restrict the collection of paths of interest to those whose parallel transports are close to those of geodesics. A particular way of measuring this is given by the following result, which is an exercise in <cit.>.

  Let H(s,t) be an arbitrary homotopy between two smooth curves γ_0(t):=H(0,t) and γ_1(t):=H(1,t) that do not intersect except at endpoints p:=H(s,0) and q:=H(s,1). Let X ∈ T_pM be a unit vector, and let X_s be the vector field defined by parallel transporting X along γ_s(t):=H(s,t) for s fixed. Then


    ||X_1(1)-X_0(1)|| ≤4/3 K_maxA(H)


where A(H) is the area of the graph of the homotopy and K_max is the maximum of the absolute Gaussian curvature of the graph of the homotopy.


In other words, the difference between the parallel transport of a geodesic and a curve homotopic to a geodesic is a function of the curvature of the manifold and the area of any homotopy. Though we do not necessarily have control over the curvature of a given shape space, it is possible using metric information to control the area term. To this end, we make a definition.

 
For vertices v_i and v_j in a connected weighted graph Γ = (V,E) with edge weights given by a matrix D = (d_ij) with d_ij = d(v_i,v_j) for some metric d, we define the Eyes on the Prize (EOP) matrix between v_i and v_j, denoted by F_i → j^Γ, is the |V| × |V| binary matrix defined by


    (F_i → j^Γ)_mn:=   1    if  d_ij < d_in and  d_jm > d_jn
    0    otherwise

 

The EOP name comes purely from intuition: a path on this adjacency matrix is simultaneously moving towards its target and can never move back towards its source, e.g. ”its eyes are always on the prize." The directed paths corresponding to the adjacency matrix F_i → j^Γ are those that both always move away from v_i and always move towards v_j. We will refer to these paths as EOP directed flows.From the perspective of Theorem <ref>, the efficient movement of the paths of F_i → j^Γ limits the area of the minimizing homotopy of Equation (<ref>).

We conclude this section by briefly remarking as to potential theoretical guarantees of paths from the F_i → j^Γ. Given the minimal information assumed, it is difficult to derive explicit estimates of Equation (<ref>), as one requires information about the underlying shape space that is not necessarily available a priori in order to derive such estimates. Nevertheless, it is possible to conclude the following:


Let x and y be two points in M with minimizing geodesic γ. Let C_x,y^ε be the collection of unit speed C^1 curves s in M that satisfy the following conditions


    g (s'(t),v_x,s(t))  > ε
    
    g (s'(t),v_s(t),y) > ε


where 1 > ε > 0 and v_x,s(t) and v_s(t),y are the initial condition vectors for the unit speed geodesic between x and s(t), and s(t) and y respectively.  Then, given a Fermi coordinate neighborhood of γ, there is a subneighborhood such that every path in C_x,y^ε whose image is contained in the subneighborhood is necessarily a function of the coordinate in the neighborhood corresponding to the geodesic.


All details are deferred to the Appendix. To summarize, there is a choice of coordinates for which the more efficient paths in F_i → j^Γ can necessarily be written as functions of one coordinate. Given a specific model of a shape space, one can then derive estimates of the bounds in Equation (<ref>). As we are explicitly attempting to not reference a particular model of shape space, we defer such discussions for future work. 



§ IMPROVING CORRESPONDENCE QUALITY VIA EOP FLOWS


Algorithm <ref> details EOP-SR, one of the main contributions of our paper. The algorithm seeks to distill soft registrations from an initial collection of registrations and distances <cit.>. For each pair of shapes, the EOP matrix from Definition <ref> is first computed. From that, we construct the corresponding a weighted version of the EOP matrix, done by taking the Hadamard (pointwise) product with the weight matrix W = (w_ij), where w_ij = K(d_ij) for some kernel function K. The soft registrations are then constructed via a depth first search approach that keeps track of both the composed correspondences as well as the weights. Recall from Equationẽqrefeqn:pathWt that the weight of a path is the product of the weights of the edges. The hard threshold λ is user-defined: such thresholds could include the number of edges used in a given path (where paths with too many edges are omitted), or path weight (where paths with low weight are omitted). In practice, the depth first search procedure has an iterative relaxation procedure. If the minimum path weight specified is too high, i.e. no paths have the minimum weight, the search is repeated with a gradually decreasing threshold until a specified number of paths are admitted.

The resulting output is thus a collection of probability distributions: each point in S_i has a corresponding probability distribution of points in S_j. Our use of soft correspondences naturally leads to robustness; small changes in the number of surfaces considered will naturally lead to small changes in the soft correspondences. In order to obtain an actual correspondence from the soft correspondences computed above, one need only to take a statistic of the soft correspondences. We restrict our attention to the maximum likelihood estimate and Frechet mean with respect to geodesic distance. 



 §.§ Accuracy of Ground Truth Propagation




To evaluate the quality of the soft registrations computed by Algorithm <ref>, we test our ability to propagate given ground truth for a collection of anatomical surfaces. We use as data the 116 tooth crowns used in <cit.>. These teeth are equipped with a consistent set of eighteen observer landmarks divided into two categories: type 2, which correspond to nongeometric features relative to the position of the tooth in the jaw, and type 3 landmarks, which corresponding to more familiar geometric features such as cusps and saddles. Two such landmarked surfaces are given in Figure <ref>. 

For each pair of surfaces, we compute two different initial collections of registrations and induced distances via two unsupervised methods: those via the continuous Procrustes distance <cit.>, and those via the GP-BD procedure defined in <cit.>. Registrations from the continuous Procrustes method are constructed by first computing local curvature extrema, finding an optimal Möbius transform via exhaustive search that minimizes distances between these extrema after the surfaces are embedded in the hyperbolic disk, and finally interpolates the correspondences via a thin plate spline. The GP-BD methodology focuses on sampling a number of Gaussian Process landmarks, for which initial putative matches are then made via the wave kernel signature <cit.> and then subsequently refined and interpolated via methods optimizing bounded distortion <cit.>. Note that both methods neither require knowledge of ground truth nor significant preprocessing such as alignment; all correspondences are learned. For both collections, we consider the following five methods for landmark propagation



  * Direct pairwise propagation

  * Propagation along a minimum spanning tree

  * Propagation along the shortest path in the 5-NN graph

  * Maximal likelihood estimate of the output of EOP-SR (EOP Mode)

  * Frechet mean of the output of EOP-SR (EOP Mean)


The first two methods listed are currently utilized and are present to serve as a benchmark as well as to illustrate the general inaccuracies that can occur by solely relying on direct propagation. The third method is natural to consider given the context of the paper, as such a path is an approximation of the geodesic in shape space. The last methods are statistics derived from EOP-SR. 

The kernel we use in defining our weight matrices is the Gaussian kernel K(d) = e^-d^2/σ. Letting d̅ denote average of the distances between each pair of surfaces, the default parameters used for the correspondences derived from Algorithm <ref> are as follows: we let σ = 0.25 d̅, the maximum number of edges in any path is 4, and the minimum weight of any path is e^-d̅^2.




To visualize the accuracy of our method, we utilize the Princeton benchmark introduced in <cit.>. The axes of this graph are geodesic distance and percentage of computed landmarks that are at most the corresponding geodesic distance away from the ground truth. Each curve on this graph represents the resulting accuracy of the correspondences. We see the results for all pairwise correspondences for all pairs of surfaces in Figure <ref>. We see that for both collections of registrations, propagation via the mode of our soft correspondences is more accurate than the other methods utilized. That the mode is generally more accurate than both the direct and MST-based propagation is not surprising; the former tends to fail for highly nonisometric pairs, and our soft correspondences do not have any consistency requirements as imposed by the MST-based method. It is perhaps more surprising, however, that the soft correspondence mode is generally superior to the shortest path method, and that the Frechet mean does as poorly as it does. To the former, we remark that any shortest path computed is one particular approximation to a geodesic, while our soft correspondences essentially incorporate multiple, allowing for more robust decision making. Given the success of the soft correspondence mode, the failure of the soft correspondence Frechet mean is likely due to the shape of the support of the distribution.



 §.§ Evaluation of Parameter Choices




In the defintion of Algorithm <ref>, a number of different parameters choices are mentioned. It is naturally interesting to see what, if any effect, they have on the accuracy results. Figure <ref> shows the accuracy results for the soft correspondence mode propagation when changing three separate parameters. The top of the figure considers three different values of σ = a σ̅, where a = 0.1, 0.25, 0.5. The middle considers changes in admissible path weight
of the form b e^-d̅^2, where b = 0.5, 1, 2. Lastly, the bottom of the figure considers changes in the maximum number of edges in any admissible path, ranging between 3 and 5. The figure shows negligible change in all cases, suggesting robustness of the parameters chosen.









§ APPLICATION: ROBUST AND CONSISTENT SURFACE REGISTRATION

The results from the previous section suggest that our forward propagation methodology can significantly improve registration quality. Those results, however, are alone insufficient for real-world application for two reasons. First, any rigorous analysis, such as one using standard shape space methods <cit.>, ultimately requires consistent registration of all surfaces in a collection. Both Algorithm <ref> and all experiments related to it conducted in Section <ref> make no assumption about end result concerning consistency; the algorithm makes a choice based on the computed soft correspondence, but has no requirements on consistency, nor is there any requirement that the points used in each surface belong to some canonical subset such as the landmarks from <cit.>. The algorithm could be easily modified if such a subset were known: for example, a projection operation could have been added to the experiments in Section <ref> that would require any decision be in terms of one of the eighteen landmarks used. However, the second issue we bring up is that such an expectation of always having a canonical set of expert landmarks is not always feasible in practice, especially for large collections of surfaces. Even if they were to exist, we would be immediately presented with a paradox: from a practical perspective, if we knew such landmarks a priori, we would also likely know how they are registered! With such registrations, we could then easily obtain registrations of the entire surface, assuming that such a surface is a topological disk or sphere, using established methods for interpolating sparse correspondences <cit.>. Thus, the main practical bottlenecks are to first determine sparse correspondences between a given pair of surfaces which can be interpolated to yield a high-quality registration, and then enforce that such registrations are consistent.



 §.§ Robust Partial Matching of Gaussian Process Landmarks




Maps generated by the continuous Procrustes (cP) <cit.> and Gaussian Process-Bounded Distortion (GP-BD) <cit.> procedure operate at a high level in the same way. First, candidate points on each surface of interest are generated, matches between said candidate points (if they exist) are obtained, and the resulting sparse correspondences are then interpolated in order to complete the registration procedure. In particular, not all of the generated points are placed in correspondence; the output of each method prior to interpolation is a partial matching with some of the generated candidate points omitted. Given the success of our soft correspondence mode in improving registration accuracy, as shown in Figure <ref>, we would like to adapt our procedure in order to generate robust sparse correspondences of surfaces using only the prior information given, i.e. namely the surfaces and the initial collection of registrations and distances.



Our proposed partial matching algorithm can be found in Algorithm <ref>, for which a demonstrated workflow can be found in Figure <ref> using registrations from the cP procedure on a pair of surfaces from <cit.> using 50 Gaussian process landmarks<cit.>, which have been empirically shown to generate a collection of landmarks used that outlines regions of significant curvature on a surface. First, for a given pair of shapes S_a and S_b,  we follow Algorithm <ref> to generate soft correspondences for each of the possible landmarks in A and B. We then employ a conservative mutual matching procedure. A pair of landmarks (v_i^a,v_j^b) are matched if most of the probability mass of f_ba^s(v_i^a) is concentrated around a small geodesic ball B_ε(v_j^b)) around v_j^b, where ε>0 is a user-input quantity, and vice versa. Given that the soft correspondences are generated from a family of efficient paths between S_a and S_b, we are in essence declaring that a pair of landmarks are registered if their soft correspondences are approximately cycle-consistent (ACC) with respect to the family of paths chosen. Note in particular that, for ε sufficiently small, the ACC constraint guarantees that each landmark on S_a can be paired with at most one landmark on S_b and vice versa, removing any possibility of tie-breaking. Note that Figure <ref> shows the effect of this conservative constraint: the back cusp on the top specimen is not registered to any point on the cuspless specimen on the bottom. We do not, however, view this as a flaw in light of the expert registered landmarks given in <cit.>. Algorithm <ref> is conservative by design, and only registers pairs based on their satisfaction of a relatively stringent condition. The lack of registration can only indicate sufficient uncertainty, given the initial collection of distances and registrations, to register areas of distinct nonisometry.

As a brief aside, it is perhaps of mathematical interest to note that our sparse registration procedure does not involve the direct optimization of a particular objective function. Rather, our procedure is one of determining which points satisfy the ACC constraint given by the two if-statements, and not just choosing the most cycle-consistent pairs via minimization of an objective function as in <cit.>. We suspect that this distinction is important and is worth studying of its own accord, especially given the recently established importance of cycle consistency outside of registration, e.g. CycleGAN <cit.>. Such a study is, however, out of scope for this paper and is deferred to future, much more general work.



Figure <ref> shows two examples, one in each row, of scenarios commonly encountered in collections of highly nonisometric surfaces <cit.>, again using surfaces from <cit.> and Gaussian process landmarks. The left shows a surface that has been textured based on its conformal parametrization to the plane via the mid-edge uniformization <cit.>. The middle shows another surface with transferred texture from the left based on its cP registration. The right is the resulting texture transfer based on a combination of the output of Algorithm <ref> with interpolation based on thin plate splines, as in the cP registration method. In the top row, we note that the texture based on the cP registration is not aligned with that of the base texture; this can be seen by looking at the change in the texture along the bottom ridge. The texture obtained from the Gaussian Process partial matching procedure does not suffer from this issue. In the bottom row, we see that the pink portion of the texture on the Gaussian Process partial matching registration (around the top right of each surface) exhibits significantly less warping than that of the cP registration. Thus, our partial matching procedure using an entire collection of registrations along with the ACC constraint can allow for higher quality registration than those obtained by directly minimizing an objective, making it easy to fix these common errors.



 §.§ An Algorithm to Consistently Register Aligned Surfaces

We are now ready to define our method for obtaining consistent registration of surfaces, but first add one additional motivating remark. Many methods used in end-to-end registration of collections of surfaces, either in their proposal or their implementation, make explicit or implicit assumptions as to the topology of surfaces in question. For example, the implementation of both the cP and GP-BD registration methods require surfaces to be topological disks. This is a major bottleneck, as many surfaces of practical interest need not be simply connected. Some methods, however, do not: the hyperbolic orbifold procedure for interpolating sparse correspondences proposed in <cit.>, for example, is theoretically valid for surfaces of higher genus (though to the best of our knowledge, no such implementation exists). Given that the methods by design depend on an initial collection of registrations, it is thus crucial that any pipeline we propose be theoretically amenable to surfaces regardless of topology; merely relying on the previously used methods depending on disk topology is thus insufficient.

Note, however, that none of our algorithms make any assumptions on topology. Indeed, the methods we have developed readily adapt to any particular choice of registration. Given our restricted interest in biological surfaces, which are generally rigid in nature, a natural candidate input to our methods appears: projection after alignment. Indeed, though biological surfaces such as the ones considered here have proven difficult to fully register, there are well established algorithms (such as Auto3DGM <cit.>) that will accurately align large collections of such surfaces independent of assumed topology. Thus, distances and registrations from these algorithms are natural to feed in to the procedures we've outlined.


The formal pipeline we propose is titled EOP-Proj and is given by Algorithm <ref>. Though the pipeline does not explicitly depend on aligned surfaces (and hence, projections after alignment), we restrict our attention to this setting. The algorithm works by first extracting two sets of points of interest from each shape: local maxima of some curvature (with radius of local maxima user-defined), and Gaussian Process landmarks. We pick as a template shape the Frechet mean with respect to the input distance matrix so as to have a canonical, data-driven template for which to base all registration off. We then run the EOP-PM on both sets of landmarks (the curvature extrema and the Gaussian process landmarks) to obtain two sets of matched landmarks, though using different mismatch thresholds for each set of landmarks. We are then left with two sets of sparse correspondences. 

Since both sets of landmarks are generated independently of one another, there is a significant possibility that both collections of registered pairs overlap. In order to deal with this overlap, we employ the following procedure. Every registered pair of curvature extrema is admitted into our collection of registered landmarks. Additional registered pairs are added iteratively based on maximizing the following functional, for A_i = {(a_1^k,a_2^k)}⊆ S_i× S_j_Fr: 


    E(p,A_i) =  min_(a_1,a_2) ∈ A_i d_S_i(p,a_1) +  min_(a'_1,a'_2) ∈ A_i d_S_j_Fr(p,a'_2)


where d_S_i and d_j_Fr are the geodesic distances on S_i and S_j_Fr respectively. In other words, we greedily add more registered pairs based on how far each pair is away from previously admitted pairs. Our choice of starting with registered pairs stemming from curvature extrema is motivated by empirical observations for small collections of nonisometric surfaces, where homologous curvature maxima may be located farther away from one another, which may not be properly accounted for in Gaussian Process landmark registration. Thus, we can think of the addition of Gaussian Process landmark pairs as filling in details unaccounted for by the curvature extrema.

After obtaining the sparse registration, we proceed by interpolating these pairs to obtain full registrations of the surfaces with the Frechet mean. Once these are obtained, we then reparametrize each of the surfaces so that the mesh structure is the same as that of the Frechet mean. This is possible if the registrations computed are homeomorphisms; though homeomorphisms are not guaranteed in general, some methods for computing full registrations are guaranteed to be locally homeomorphic (e.g. <cit.>). We do not, however, propose a specific procedure for this that we always follow, and leave that for the user to determine. At this point, all surfaces are consistently registered as all registrations follow naturally from the mesh structure of the Frechet mean.



 §.§ Experimental Results




  §.§.§ Tooth Crowns



Continuing our previous experiments, we apply Algorithm <ref> to the entire collection of tooth crowns from <cit.> to generate a consistently registered collection of surfaces. All surfaces are first aligned via Auto3DGM <cit.>. Here, we use as curvature the conformal factor corresponding to the mid-edge uniformization <cit.> along with 400 Gaussian process landmarks, and restrict the total collection of matches used in registration to at most 14. As the meshes have reasonably uniform triangulation, we define our threshold radii in terms of the discrete distance based on the triangulation (i.e. the number of edges between one vertex and another): we use a radius of 8 for the conformal factor maxima, and 1 for the Gaussian Process landmarks. All parameters for the soft correspondence are the same default parameters used in Section <ref>. The surfaces were all embedded into the plane using the same landmark constraint methodology employed in <cit.>. We show texture maps on thirty of the consistently  registered surfaces in Figure <ref>.



As a test of the quality of our registrations, we repeat the leave one out classification benchmark introduced in <cit.> and report the results in Table <ref>. Each method in the table naturally yields a distance matrix, making leave one out classification a matter of finding the smallest nonzero element of each row. Our methods are given by EOP-Proj-FPS-200 and EOP-Proj-GP-200. To generate our distance matrix for EOP-Proj-FPS-200, we first compute the conformal factor extrema of the Frechet mean surface used in Algorithm <ref>, and then feed those points into either the usual farthest point sampling algorithm until we have a collection of 200 points. We then take, from each surface in our colletion, the corresponding 150 points, and then compute the usual Procrustes distance (after normalization and alignment to the point cloud of the Frechet mean). The sample points used in EOP-Proj-GP-200 are the corresponding points after selecting 200 Gaussian Process landmarks on the Frechet mean surface.

We compare our method to a number previously studied on this example: expert landmarks <cit.>, continuous Procrustes registration <cit.>, Auto3DGM <cit.>, and PM-SDP <cit.>, a more recent method based on semi-definite programming that, to the best of our knowledge, achieved the best prior results on this dataset. We see in the table that EOP-Proj-GP-200 achieves state of the art results for both classification by genus and family, and is better than expert landmarks for classification by order. This heavily suggests that our method is able to better capture the fine level details of shape space than other automatic and expert-based methods. That the Gaussian Process method outperforms the farthest point sampling method gives further credence to the idea that the way in which the points are subsampled over the surface matters



  §.§.§ Talus bones








As previously mentioned, there is no prior restriction of our methods on surface topology. To illustrate this, we consider a small sample of talus bones from three genera (Adapis, Eulemur, and Perodicticus) that were obtained from Morphosource <cit.>. These surfaces were registered using Gaussian curvature maxima and 200 Gaussian process landmarks. Registered landmarks were interpolated by the hyperbolic orbifold method proposed in <cit.>. Registrations are shown in Figure <ref>, showing homologous curvature extrema in correspondence. We can also see in Figure <ref> the result of a two dimensional PCA embedding of these surfaces with respect to the Procrustes distances. The distances generated by EOP-Proj are capable of clustering each genus, and showing that both Adapis and Eulemur are more closely related to each other than either are to Perodicticus despite using a relatively low dimensional embedding.



 §.§ Limitations

The main limitation of EOP-Proj and variants that use some other class of registrations, as well as the other algorithms proposed in the paper, is the implicit requirements on the datasets at hand. Namely, the implicit assumption that the collection of shapes is sufficiently densely sampled is not a ubiquitous property of many datasets of interest within the realm of registration (though many in addition to the ones considered above, such as DynaFAUST <cit.> and faces <cit.>, at least somewhat satisfy this property). Without this assumption, the differential geometric ideas on which we built our algorithms breaks down. We do not, however, see this as a total negative, especially in light of the promising state of the art results achieved by the above methods. The methods proposed are a natural way to deal with registration of nonisometric shapes, in which pure optimization can fail. We thus see our methodology as addressing an important use case that had previously gone unaddressed within the realm of registration.

Another limitation of EOP-Proj is the particular way in which the registrations are made consistent. Our use of the Frechet mean, though a reasonably heuristic, is by no means guaranteed to be optimal. An easy thought experiment suggests potential problems can arise if the Frechet mean itself has radically different geometry than other shapes in the collection. Alternative methods to guarantee consistency are thus subjects for future work.


§ CONCLUSION

In this paper, we investigated the problem of consistent registrations of nonisometric surfaces. Motivated by biology, manifold learning and the previous role of parallel transport within registration, we proposed soft registration via EOP flows (EOP-SR), which have appealing properties both theoretically and computationally. We showed that employing EOP-SR can drastically improve accuracy over registration based purely on optimization methods while being robust by registering a point on one shape to a distribution of points on another shape. We adapted our methodology to build a partial matching method, EOP-PM, that can fix even subtle issues in registration quality. We then showed how EOP-PM can be used to accurately and consistently register large collections of shapes, achieving state of the art classification results in the process.

There are many directions of interest in which to adapt our work. In the course of our paper, interesting theoretical questions concerning parallel transport arose; progress on these questions could yield insight for guarantees about registration quality. From a computational perspective, our focus on surfaces from evolutionary biology was purely of interest and not a restriction. In particular, other datasets, such as face scans and certain collections of images, are natural candidate in which to adapt the algorithms of our paper, which we leave for future work. Finally, as mentioned in the Section <ref>, the ACC condition proposed in the definition of EOP-PM has potential merit for further investigation given the recent increase that cycle consistency has received within the general machine learning literature.



 §.§ Declarations

The author was partially supported by AFOSR through an NDSEG Research Fellowship. The author declares no conflicts of interest. Availability information of data and code relevant to the above discussions is detailed in the Appendix. The author would like to thank Ingrid Daubechies, Tingran Gao, Shahar Kovalsky, Chen-Yun Lin, Sayan Mukherjee, Shan Shan, and Barak Sober for many helpful discussions in compiling this work.

plainnat






§ APPENDIX



 §.§ Further Review of Differential Geometry

We briefly elaborate on details relevant to discussions in Sections <ref> and <ref>. For brevity, we refer to readers largely unfamiliar with smooth manifold theory and Riemannian geometry to standard texts previously mentioned, namely <cit.>, and review only details most relevant for our work.



  §.§.§ Riemannian Geometry

We let (M^m,g) be an m-dimensional Riemannian manifold with metric g and tangent bundle TM. We also denote by 𝒱(M) the space of smooth vector fields on M. Recall that vector fields naturally as differential operators on smooth functions; e.g. if a vector field X = ∑_i=1^m X_i∂/∂ x_i is an expansion of X in local coordinates, then for a smooth function f on M, we have Xf = ∑_i=1^m X_i∂ f/∂ x_i.


A connection on a Riemannian manifold (M,g) is a map ∇: 𝒱(M) ×𝒱(M) →𝒱(M), written (X,Y) →∇_X Y, that satisfies the following three properties:



  * For f, g smooth functions on M and X, Y, Z smooth vector fields, we have


    ∇_fX + gY Z = f ∇_X Z + g ∇_Y Z



  * For a, b ∈ℝ and X,Y,Z smooth vector fields on M, we have


    ∇_X (aY+bZ) = a ∇_XY + b ∇_XZ



  * For f a smooth function on M and X,Y smooth vector fields, we have a Leibniz rule


    ∇X(fY) = f ∇_XY + (Xf)Y


where X acts a partial derivative on f



Connections are the appropriate analogue of a directional derivative for vector fields. Specifically, for a vector field X along a curve γ:[a,b] → M, X'(t) = d/dtX := ∇_γ̇(t), sometimes written as ∇_d/dt X. It is well-known in Riemannian geometry that every manifold has a unique connection called the Levi-Civita connection that satisfies the above as well as two additional properties: it is symmetric, i.e.


    ∇_XY -∇_YX = [X,Y]


where [·, ·] is the Lie bracket on TM, the tangent bundle of M, and it is compatible with the metric, i.e.


    ∇_Xg(Y,Z) = g ( ∇_XY,Z ) + g (X, ∇_Y Z )


The Levi-Civita connection satisfies the following property, which we need in proving Theorem <ref>. Note, however, that we only use the compatibility property.

 
Let X(t) be a smooth vector field defined along a curve γ: [a,b] → M nowhere equal to zero. Then
we have


    d/dt||X|| := d/dt(g(X,X))^1/2≤( g (X,X ) )^1/2 = ||∇_d/dtX||


where all derivatives are taken with respect to the Levi-Civita connection.


Since the Levi-Civita connection is symmetric, we have


    d/dt(g(X,X))^1/2   = 1/2g(X,X)^- 1/2d/dt g(X,X) 
       =1/2g(X,X)^- 1/2 2 g(∇_d/dtX,X) 
       =g(X,X)^- 1/2g(∇_d/dtX,X)

by symmetry of the Riemannian metric and compatibility of the Levi-Civita connection. A direct application of the Cauchy-Schwarz inequality on g(X',X) completes the proof.


We assume for the remainder that we are only working with the Levi-Civita connection.



  §.§.§ Parallel Transport


The heart of our methodology relies on parallel transport. The usual definition pertaining to vector fields is easy to state given for a Levi-Civita connection; given a smooth curve γ:[a,b] → M, there is a natural vector field ·γ defined on the image of γ given by the usual time derivative. Then, the parallel transport of a vector v_γ(a)∈ T_γ(a)M along γ is the unique vector field V defined on γ(t) such that V(a) = v and


    ∇_·γ(t)V(t) = 0.


For fiber bundles, the standard definition of parallel transportis similar but quite technical, requiring a number of definitions that we would not directly reference again. In light of this, we adopt an equivalent definition, in which the notion of connection is itself abstracted away from the usual language. Despite this abstraction, the concepts are much more intuitive and clearly apply to our setting. Let H(F) denote the space of homeomorphisms (or diffeomorphisms, as appropriate) of F.


Let (E,B,F,π) be a fiber bundle with base manifold B and fiber F. A connection is a function that maps every curve γ:[0,1] → B, to a family of homeomorphisms of F (or diffeomorphisms, if the fiber bundle is smooth) P_γ: [0,1] × [0,1] → H(F) such that P_γ is continuous on [0,1] × [0,1], P_γ(t,t) = Id_F, and satisfies the following transitivity property for all a,b,c ∈ [0,1]:


    P_γ(c,b) ∘ P_γ(b,a) = P_γ(c,a)




Let (E,B,F,π) be a fiber bundle with base manifold B and fiber F. A connection is a function that maps every curve γ:[0,1] → B, to a family of homeomorphisms of F (or diffeomorphisms, if the fiber bundle is smooth) P_γ: [0,1] × [0,1] → H(F) such that P_γ is continuous on [0,1] × [0,1], P_γ(t,t) = Id_F, and satisfies the following transitivity property for all a,b,c ∈ [0,1]:


    P_γ(c,b) ∘ P_γ(b,a) = P_γ(c,a)


Given such a connection, the parallel transport of γ is given by P_γ(1,0).


This property can be derived from the usual definition of parallel transport on fiber bundles as present in <cit.>. Of particular interest, however, is that our definition works for more general cases when smooth structures don't exist, and would readily apply to spaces of meshes.



  §.§.§ Fermi normal coordinates


The last piece of exposition we require is Fermi normal coordinates, which were mentioned briefly in Theorem <ref>. Fermi normal coordinates are the natural analogue to Cartesian coordinates in Euclidean space as exponential coordinates are the manifold analogue of polar coordinates in Euclidean space.


Let P^q⊂ M^m be a regular submanifold, let x_1,...,x_q be a coordinate chart on P, and let E_q+1,...,E_m be a section of the tangent bundle along P such that the E_j|_p are an orthonormal completion of T_pP in T_pM. Then the functions y_i, 1 ≤ i ≤ m defined by


    y_i( exp_p( ∑_j=q+1^m a_jE_j|_p) ) = x_i(p), 1 ≤ i ≤ q 
    
    y_i( exp_p( ∑_j=q+1^m a_jE_j|_p) ) = a_i, q+1 ≤ i ≤ m


form a Fermi normal coordinate chart with respect to P.
 

Though the definition we state is for the general setting, we are only interested in the case where q = 1 in the above definition. Thus, in our work, we let P be a geodesic between two points in a manifold. Theorem <ref> that there is a family of curves that can be written as functions along the geodesic axis of such a Fermi coordinate chart.



 §.§ Proofs



  §.§.§ Proof of Theorem <ref>


By the fundamental theorem of calculus, we have


    ‖ X_1(1)-X_0(1) ‖   = ‖∫_0^1∇_d/dsX_s(1) ds ‖
       ≤∫_0^1‖∇_d/ds(X_s(1)) ‖ ds   
       = ∫_0^1( ‖∇_d/ds(X_s(1)) ‖- ‖(∇_d/dsX_s(0)) ‖)ds


where the last equality follows from the assumption that X_s(0) is constant. Another application of the fundamental theorem gives


    ‖ X_1(1)-X_0(1) ‖   ≤∫_0^1( ‖∇_d/ds(X_s(1)) ‖- ‖(∇_d/dsX_s(0)) ‖)ds 
       = ∫_0^1( ∫_0^1d/dt‖∇_d/ds(X_s(t)) ‖ dt )ds 
       ≤∫_0^1( ∫_0^1‖∇_d/dt∇_d/ds(X_s(t)) ‖ dt )ds 
       = ∫_0^1( ∫_0^1‖ (∇_d/dt∇_d/ds - ∇_d/ds∇_d/dt) (X_s(t)) ‖ dt )ds



where the second inequality follows from Lemma <ref> and the last equality follows from the definition of parallel transport. Since d/ds and d/dt are coordinate vector fields, we have from Definition <ref>


    ‖ X_1(1)-X_0(1) ‖≤∫_0^1( ∫_0^1‖ (R(H_t,H_s)(X_s(t)) ‖ dt )ds


where H_t and H_s are the corresponding partial derivative vector fields of the homotopy H. The Riemann curvature term in the above integral is bounded above by 4/3 K_max‖ H_s∧ H_t‖, where ‖ H_s∧ H_t‖ is the induced area form of the homotopy; this is a direct consequence of Lemma 3.7 in <cit.>. The desired result immediately follows.




  §.§.§ Proof of Theorem <ref>

Assume that the inner products of s'(t) with v_s(t),x and v_s(t),y are strictly bounded above and below respectively by some uniform numbers M_1 and M_2 whenever the inner products are well defined, i.e. when s(t) ≠ x, y. We now define two sets: let U_x,ε be the subset of M such that every point u ∈ U_x,ε satisfies g ( v_x,u, ∂/∂γ) ≥ 1 - ε, and analogously  U_y,ε. If M_1 and M_2 are at least ε, then any such curve s(t) contained in the intersection of the Fermi coordinate neighborhood of γ with the union of U_x,ε and U_y,ε necessarily has to satisfy g ( s'(t), ∂/∂γ) > 0. An application of the implicit function theorem then completes the proof.



 §.§ Experiment Details



  §.§.§ Meshes

All meshes used throughout this paper are publicly available on the Morphosource Repository <cit.>. All specimens were cleaned to be either homeomorphic to a disk or sphere as appropriate and downsampled to approximately 5000 vertices.


  §.§.§ Code

We made use of the following external repositories for methods referenced in the text.



  * Continuous Procrustes Maps: https://github.com/trgao10/cPdist  
  * GP-BD Registration: https://github.com/shaharkov/GPLmkBDMatch 
  * Auto3DGM: https://github.com/trgao10/PuenteAlignment

Default parameters were used. 

The code used to run the experiments of the paper was derived from the SAMS repository found at https://github.com/RRavier/SAMS, in which the registration method is an implementation of EOP-Proj written by the author of this work. A separate repository, which https://github.com/RRavier/EOP-Experiments, will contain a more concise version focused on replicating the experiments found in the paper.
