Inexact elastic shape matching in the square root normal field framework

TL;DR

This paper introduces a new elastic shape matching method for unparametrized curves and surfaces that leverages square root normal fields and varifold metrics, enabling flexible, efficient, and topology-agnostic matching with texture incorporation.

Contribution

It presents a novel variational framework that minimizes over reparametrizations without discretizing the reparametrization group, allowing for efficient and flexible shape matching.

Findings

Effective matching of curves and surfaces with different samplings and topologies.

Incorporation of texture as additional matching information.

Efficient numerical implementation of the objective and gradient.

Abstract

This paper puts forth a new formulation and algorithm for the elastic matching problem on unparametrized curves and surfaces. Our approach combines the frameworks of square root normal fields and varifold fidelity metrics into a novel framework, which has several potential advantages over previous works. First, our variational formulation allows us to minimize over reparametrizations without discretizing the reparametrization group. Second, the objective function and gradient are easy to implement and efficient to evaluate numerically. Third, the initial and target surface may have different samplings and even different topologies. Fourth, texture can be incorporated as additional information in the matching term similarly to the fshape framework. We demonstrate the usefulness of this approach with several numerical examples of curves and surfaces.

Equations

$${\mathrm{dist}}_{\mathrm{Imm}}(f_0,f_1{)}^2=\underset{\begin{array}{c}f\in C^{\infty }([0,1],\mathrm{Imm}(M,{\mathbb{R}}^d))\\ f(0)=f_0,f(1)=f_1\end{array}}{\inf }{\int }_0^1{G}_f({\partial }_{t}f,{\partial }_{t}f)dt,$$

$${\mathrm{dist}}_{B_i}([f_0],[f_1])=\underset{\phi \in \mathrm{Diff}(M)}{\inf }{\mathrm{dist}}_{\mathrm{Imm}}(f_0,f_1\circ \phi ).$$

$$G_f(h,k)={\int }_M⟨D_{(f,h)}{\tilde{n}}_f,D_{(f,k)}{\tilde{n}}_f{⟩}_{{\mathbb{R}}^d},$$

$${\mathrm{dist}}_{\mathrm{Imm}}(f_0,f_1)\approx \Vert {\tilde{n}}_{f_0}-{\tilde{n}}_{f_1}{\Vert }_{L^2}.$$

$$\forall w\in W:({\mu }_f|w{)}_{W^{\ast },W}={\int }_{M}w(f(x),n(x)){\mathrm{vol}}_f(dx).$$

$${\mathrm{dist}}_{\mathrm{Var}}([f_0],[f_1])=\Vert {\mu }_{f_0}-{\mu }_{f_1}{\Vert }_{W^{\ast }}.$$

$$\underset{f}{\mathrm{minimize}}\mathrm{dist}(f_0,f)\text{subject to}{\mathrm{dist}}_{\mathrm{Var}}([f],[f_1])=0.$$

$$\underset{f}{\mathrm{minimize}}\Vert {\tilde{n}}_{f_0}-{\tilde{n}}_f{\Vert }_{L^2}^2+\lambda {\mathrm{dist}}_{\mathrm{Var}}([f],[f_1]{)}^2.$$

Full text

P. Harms is supported by the Freiburg Institute of Advanced Studies in the form of a Junior Fellowship. N. Charon is supported by the NSF grant 1819131.

1. Introduction

1.1. Context.

The statistical analysis of datasets of curves and surfaces is an active research field with many applications in e.g. computer vision, robotics, and medical imaging; see [24, 4, 20] and references therein. A recurring and fundamental task is finding optimal point correspondences between given shapes (i.e., the matching or registration problem), where optimality is typically expressed in terms of an elastic deformation energy. Solving the elastic matching problem in a numerically efficient way, which scales well to high-dimensional data encountered in real-world applications, remains a major challenge to date.

1.2. Relation to previous work.

This paper draws on two lines of work: square root normal fields (SRNFs) [19, 14, 11, 12], which allow one to efficiently calculate elastic distances between parametrized shapes, and varifold distances [22, 9, 7, 18, 13], which are distances between unparametrized shapes without any elastic interpretation. For each of these frameworks, efficient numerical implementations have been developed.

1.3. Contribution.

We propose a new algorithm which combines SRNFs with varifold distances and inherits many advantages of both approaches. The key idea is to use varifold distances to relax the terminal constraint in the elastic matching problem. This bypasses the discretization of the reparametrization group, thereby eliminating the main computational burden in previous implementations of SRNF-based elastic shape matching. The resulting optimization problem is easy to implement and yields good results on some preliminary experiments on curves and surfaces. Moreover, the varifold distances allow one to match shapes with different meshes and even different topologies and to use texture information as in the fshape framework.

2. Shape analysis of curves and surfaces

2.1. Elastic shape analysis.

Elastic shape analysis operates in a Riemannian framework where infinitesimal shape deformations are measured by a Riemannian metric, which is often related to an elastic (or plastic) deformation energy; see the surveys [4, 12]. We consider parameterized shapes as elements of the Fréchet manifold $\operatorname{Imm}(M,\mathbb{R}^{d})$ of immersed hypersurfaces of a $(d-1)$-dimensional compact manifold $M$ into $\mathbb{R}^{d}$. The corresponding space of unparameterized shapes is the quotient space $B_{i}(M,\mathbb{R}^{d})=\operatorname{Imm}(M,\mathbb{R}^{d})/\operatorname{Diff}(M)$, whose elements are denoted by $[f]=\{f\circ\varphi;\varphi\in\operatorname{Diff}(M)\}$. Given a $\operatorname{Diff}(M)$-invariant weak Riemannian metric $G$ on $\operatorname{Imm}(M,\mathbb{R}^{d})$, one defines a pseudo-distance between any two immersions $f_{0},f_{1}\in\operatorname{Imm}(M,\mathbb{R}^{d})$ and their equivalence classes $[f_{0}],[f_{1}]\in B_{i}(M,\mathbb{R}^{d})$ by

[TABLE]

Symmetry of the pseudo-distance on $B_{i}(M,\mathbb{R}^{d})$ follows from the invariance of the metric $G$ with respect to reparametrizations. Under suitable conditions on $G$ the pseudo-distance is a distance, i.e., it separates points in the shape space of curves [17, 16] or surfaces [5].

From a numerical perspective, the challenge is to calculate the above distances and the corresponding optimizers efficiently. This minimization can be solved numerically by path straightening and geodesic shooting methods (see e.g. [10, 2]) or as in the next section by exploiting isometries to simpler spaces.

2.2. Square root normal fields.

Problem (1) simplifies considerably for certain first order Sobolev metrics [23, 25, 19, 21, 14, 11, 3]. One class of such metrics is defined using square root normal fields (SRNFs), which were introduced by Srivastava e.a. [19, 14] for planar curves and later generalized to surfaces by Jermyn e.a. [11]. The SRNF of an oriented immersed hypersurface $f\in\operatorname{Imm}(M,\mathbb{R}^{d})$ is defined as $\smash{\tilde{n}_{f}=n_{f}\operatorname{vol}_{f}^{1/2}}$, where $n_{f}$ is the unit normal field and $\smash{\operatorname{vol}_{f}^{1/2}}$ the half density of $f$. For example, the SRNF of a planar curve $f\in\operatorname{Imm}(S^{1},\mathbb{R}^{2})$ is given in coordinates $\theta\in S^{1}$ by $\smash{\tilde{n}_{f}=if_{\theta}\|f_{\theta}\|_{\mathbb{R}^{2}}^{-1/2}}$, where $i$ denotes rotation by 90 degrees and coordinates in subscripts denote derivatives. Similarly, the SRNF of a surface $f\in\operatorname{Imm}(S^{2},\mathbb{R}^{3})$ is given in coordinates $(u,v)\in S^{2}$ as $\smash{\tilde{n}_{f}=(f_{u}\times f_{v})\|f_{u}\times f_{v}\|_{\mathbb{R}^{3}}^{-1/2}}$. In general, one obtains an elastic pseudo-Riemannian metric on $\operatorname{Imm}(M,\mathbb{R}^{d})$ by setting

[TABLE]

where $D_{(f,h)}\tilde{n}_{f}$ denotes the directional derivative of $\tilde{n}_{f}$ at $f$ in the direction $h$. This pseudo-Riemannian metric $G$ is $\operatorname{Diff}(M)$-invariant, and by construction the map $f\mapsto\tilde{n}_{f}$ is a Riemannian isometry into the flat space of square integrable vector-valued half densities. For curves one obtains a Riemannian metric by modding out translations. For surfaces the situation is more complicated, as described in [11], and the kernel of the pseudo-metric may be larger than only translations. The metric belongs to the class of first order Sobolev metrics, which have been studied in great detail [17, 16, 5].

The advantage of this construction is that the Riemannian distance of $G$ on $\operatorname{Imm}(M,\mathbb{R}^{d})$ can be approximated efficiently as follows:

[TABLE]

Equality holds whenever the straight line between $\tilde{n}_{f_{0}}$ and $\tilde{n}_{f_{1}}$ is contained in the range of the map $f\mapsto\tilde{n}_{f}$. In general, equality holds up to first order for $f_{0}$ close to $f_{1}$ because the map $f\mapsto\tilde{n}_{f}$ is a Riemannian isometry.

The approximate distance (3) descends to the quotient space $B_{i}(M,\mathbb{R}^{3})$ as described in (2). However, (2) involves a minimization over the reparametrization group, which is computationally costly. For curves this can be solved by dynamic programming [19] or using an explicit formula [15]. For surfaces in spherical coordinates, Jermyn e.a.  [11] proposed to discretize the diffeomorphism group of the two-dimensional sphere using spherical harmonics. This article puts forth an alternative method for minimization over the reparametrization group, which is based on varifold distances.

2.3. Varifold distances.

Geometric measure theory provides several embeddings of shape spaces into Banach spaces of distributions [22, 9, 7, 18, 13] with corresponding metrics. Varifold embeddings are one instance of this construction and are defined as follows (cf. [13] for details). Given a reproducing kernel Hilbert space $W$ of real-valued functions on $\mathbb{R}^{d}\times S^{d-1}$, one associates to any immersion $f\in\operatorname{Imm}(M,\mathbb{R}^{d})$ the varifold $\mu_{f}\in W^{*}$ which satisfies

[TABLE]

The map $f\mapsto\mu_{f}$ is reparametrization-invariant and, under suitable assumptions on the kernel of $W$, injective [13]. Thus, one obtains a well-defined distance on the quotient space $B_{i}(M,\mathbb{R}^{d})$ by defining for any two immersions $f_{0},f_{1}\in\operatorname{Imm}(M,\mathbb{R}^{d})$:

[TABLE]

From a computational point of view, these distances have explicit expressions in terms of the kernel function of $W$ and are easy to implement for discrete curves or surfaces. We will use such distances to relax the terminal constraint in the boundary value problem for geodesics on shape space, as described next.

2.4. Combining SRNFs and varifold distances.

Square root normal fields and varifold distances can be combined in an efficient matching algorithm for unparametrized shapes. This idea has been previously used in combination with large deformation models in [7, 13] and with $H^{2}$ metrics on the space of curves in [1]. The boundary value problem (2) for geodesics on $B_{i}(M,\mathbb{R}^{d})$ can be formulated as the program

[TABLE]

Relaxation using a (large) Lagrange multiplier $\lambda$ and approximation of the elastic distance as in (3) yields

[TABLE]

This program has several advantages over previous alternative formulations of the SRNF matching problem [19, 14, 11]. First, the objective function and its gradient are easy to implement and can be computed efficiently. Second, the initial and target surface may have different discretizations and even different topologies. Third, texture information can be incorporated into the varifold matching term similarly to the fshape framework [8, 6].

3. Numerical implementation and results

3.1. Algorithm.

Given a pair $(f_{0},f_{1})$ of curves or surfaces, the program (5) looks for a minimizer $f$ with $\operatorname{dist}_{B_{i}}([f_{0}],[f_{1}])=\operatorname{dist}_{\operatorname{Imm}}(f_{0},f)$ and $[f]=[f_{1}]$. Thus, the algorithm solves the registration problem and calculates the distance between the unparametrized shapes $[f_{0}]$ and $[f_{1}]$. Note that it does, however, not provide a geodesic homotopy between these shapes. Such a homotopy can be obtained from the linear homotopy between $\tilde{n}_{f_{0}}$ and $\tilde{n}_{f_{1}}$ by (approximate) inversion of the SRNF map $f\mapsto\tilde{n}_{f}$. For open curves this inversion is exact and easy to implement. For closed curves, the range of the SRNF map is not convex, and an approximate inverse has to be used [19]. For surfaces, this is a delicate issue [11], and to the best of our knowledge there exists no publicly available implementation for general triangulated surfaces.

3.2. Implementation.

To implement the program (5) numerically, one has to discretize the space of parametrized shapes. An advantage over [11] is that the reparametrization group does not need to be discretized. Piecewise linear curves and triangular meshes are suitable discretizations in our context, the reason being that square root normal fields and kernel-based varifold distances extend naturally to these spaces. The minimization is performed using an L-BFGS method. The gradient of the discretized energy functional (5), which is needed by the L-BFGS method, has an explicit form and can be implemented efficiently.

3.3. Curves.

For curves, our algorithm is comparable to state of the art methods. On the Kimia dataset111Computer Vision Group at LEMS at Brown University: Database of 99 binary shapes. https://vision.lems.brown.edu/content/available-software-and-databases it produces distances and clusters which are similar to those based on dynamic programming, as shown in Fig 1. A nice feature of our algorithm, which stems from the use of varifold distances, is that the initial and target shapes are allowed to have different topologies. For example, one can match a single circle to a pair of circles, as demonstrated in Fig. 2. This is not possible using previous methods for shape matching using SRNFs or SRVTs. There are potential applications in cell division and removal of topological noise. Another feature of our algorithm is that it can account for functional data on the given shapes, as demonstrated in Fig. 3. To this aim, the varifold distance in (5) is replaced by a functional shape distance, as developed in [8, 6]. This has several applications. The functional data may be dictated by the application at hand, as e.g. in the case of texture information. An interesting alternative to be explored in future work is to use shape descriptors as functional data to guide the matching algorithm.

3.4. Surfaces.

For surfaces, we obtain some promising first results and see a high potential of improvement over alternative methods. An example is presented in Fig. 4, where the optimal point correspondences between two hand postures were calculated. As the two triangulated surfaces in this experiment had different mesh connectivities, and no point-to-point correspondences were initially available, we had to initialize the optimization procedure with the template surface. After optimization using an adaptive choice of Lagrange multiplier $\lambda$ in (5), we obtained an excellent fit of the deformed template onto the target with anatomically correct point correspondences.