Surface Networks

TL;DR

Surface Networks (SNs) enhance 3D mesh representations by leveraging extrinsic geometry, specifically the Dirac operator, to improve deformation stability and modeling power in shape analysis tasks.

Contribution

The paper introduces Surface Networks that incorporate the Dirac operator for extrinsic geometry, surpassing Laplacian-based models in stability and expressiveness.

Findings

SNs are stable to deformation and discretization.

SNs outperform existing models in mesh deformation prediction.

SNs enable effective generative modeling of 3D surfaces.

Abstract

We study data-driven representations for three-dimensional triangle meshes, which are one of the prevalent objects used to represent 3D geometry. Recent works have developed models that exploit the intrinsic geometry of manifolds and graphs, namely the Graph Neural Networks (GNNs) and its spectral variants, which learn from the local metric tensor via the Laplacian operator. Despite offering excellent sample complexity and built-in invariances, intrinsic geometry alone is invariant to isometric deformations, making it unsuitable for many applications. To overcome this limitation, we propose several upgrades to GNNs to leverage extrinsic differential geometry properties of three-dimensional surfaces, increasing its modeling power. In particular, we propose to exploit the Dirac operator, whose spectrum detects principal curvature directions --- this is in stark contrast with the classical Laplace operator, which directly measures mean curvature. We coin the resulting models \emph{Surface Networks (SN)}. We prove that these models define shape representations that are stable to deformation and to discretization, and we demonstrate the efficiency and versatility of SNs on two challenging tasks: temporal prediction of mesh deformations under non-linear dynamics and generative models using a variational autoencoder framework with encoders/decoders given by SNs.

Tables

Table 1. Table 1: Evaluation of different models on the temporal task

Model	Receptive field	Number of parameters	Smooth L1-loss (mean per sequence (std))
MLP	1	519672	64.56 (0.62)
PointCloud	-	1018872	23.64 (0.21)
Laplace	16	1018872	17.34 (0.52)
Dirac	8	1018872	16.84 (0.16)

Equations

$$x^{k+1}=\rho \left(A_k\Delta x^k+B_k{x}^k\right),A_k,B_k\in {\mathbb{R}}^{d_{k+1}\times d_k}.$$

$$\Delta x^1=-2H\mathbf{n},$$

$$D_{f,j}=\frac{-1}{2|{\mathcal{A}}_f|}{e}_j,f\in F,j\in V,$$

$$\text{mbox}Re{D}^{\ast }D=\Delta ,$$

$$y^{k+1}=\rho \left(C_{k}D{x}^k+E_k{y}^k\right),C_k,E_k\in {\mathbb{R}}^{d_{k+1}\times d_k},$$

$$x^{k+1}=\rho \left(A_k{D}^{\ast }{y}^{k+1}+B_k{x}^k\right),A_k,B_k\in {\mathbb{R}}^{d_{k+1}\times d_k},$$

$${\Delta }_{\mathcal{M}}=\sum _{k\le N}{\lambda }_k{e}_k{e}_k^T,e_k\in {\mathbb{R}}^N,0\le {\lambda }_1\le {\lambda }_2\cdots \le {\lambda }_N.$$

$$\Vert \hat{x}(k)\Vert \lesssim k^{-\beta },(k\to \infty ).$$

$$\Vert {\Phi }_{\Delta }(x_1;{\mathcal{M}}_N)-{\Phi }_{\Delta }(x_2;{\mathcal{M}}_N){\Vert }^2\le C(\beta )\Vert x_1-x_2{\Vert }^{h(\beta )},$$

$$\Vert \Phi (\mathcal{M})-\Phi (\tau (\mathcal{M}))\Vert \le \kappa \max (|\nabla \tau {|}_{\infty },\Vert {\nabla }^2\tau \Vert {)}^{h(\beta )}.$$

$$\underset{\theta ,\psi }{\min }\frac{1}{L}\sum _{l\le L}-{\mathbb{E}}_{h\sim q_{\psi }(h|{\mathcal{M}}_l)}\log p_{\theta }({\mathcal{M}}_l|h)+D_{KL}(q_{\psi }(h|{\mathcal{M}}_l)||p_0(h)).$$

$$p_{\theta }(f|h,\tilde{\mathcal{M}})=\mathcal{N}(\mu (h,\tilde{\mathcal{M}}),{\sigma }^2(h,\tilde{\mathcal{M}})1),$$

$$D_{f,j}=\frac{-1}{2|{\mathcal{A}}_f|}{e}_j,f\in F,j\in V,$$

$$\left[\begin{array}{cccc}a& -b& -c& -d\\ b& a& -d& c\\ c& d& a& -b\\ d& -c& b& a\end{array}\right]$$

$$\left[\begin{array}{cccc}a& b& c& d\\ -b& a& d& -c\\ -c& -d& a& b\\ -d& c& -b& a\end{array}\right].$$

$$\Delta =\text{diag}(\bar{\mathcal{A}}{)}^{-1}(U-W),$$

$$\Vert U-W\Vert$$

$$\Vert \Delta \Vert \le C\frac{\cot ({\alpha }_{\text{min}})S_{\text{max}}}{{\inf }_j{\bar{\mathcal{A}}}_j}:=L_{\mathcal{M}},$$

$$\Vert \rho \left(Ax+B\Delta x\right)-\rho \left(A{x}^{\prime }+B\Delta x^{\prime }\right)\Vert$$

$$\Vert \Phi (\mathcal{M};x)-\Phi (\mathcal{M};x^{\prime })\Vert \le \left(\prod _{k\le K}(\Vert A_k\Vert +\Vert B_k\Vert L_{\mathcal{M}})\right)\Vert x-x^{\prime }\Vert ,$$

$$\Vert D\Vert \le \sqrt{L_{\mathcal{M}}}.\square$$

$$\Vert p-q{\Vert }^2(1-|\nabla \tau {|}_{\infty }{)}^2$$

$$\mathcal{A}(p,q,r{)}^2(1-|\nabla \tau {|}_{\infty }C{\alpha }_{\text{min}}^{-2}-o({|\nabla \tau {|}_{\infty }}^2)$$

$$\mathcal{A}(p,q,r{)}^2=s(s-\Vert p-q\Vert )(s-\Vert p-r\Vert )(s-\Vert r-q\Vert ),$$

$$s_{\tau }-\Vert \tau (p)-\tau (q)\Vert \le s(1+|\nabla \tau {|}_{\infty })-\Vert p-q\Vert (1-|\nabla \tau {|}_{\infty })=s-\Vert p-q\Vert +|\nabla \tau {|}_{\infty }(s+\Vert p-q\Vert )\text{and }$$

$$s_{\tau }-\Vert \tau (p)-\tau (q)\Vert \ge s(1-|\nabla \tau {|}_{\infty })-\Vert p-q\Vert (1+|\nabla \tau {|}_{\infty })=s-\Vert p-q\Vert -|\nabla \tau {|}_{\infty }(s+\Vert p-q\Vert ),$$

$$\mathcal{A}(\tau (p),\tau (q),\tau (r){)}^2$$

$$\mathcal{A}(\tau (p),\tau (q),\tau (r){)}^2\le \mathcal{A}(p,q,r{)}^2\left[1+C|\nabla \tau {|}_{\infty }{\alpha }_{\text{min}}^{-2}-o({|\nabla \tau {|}_{\infty }}^2)\right].$$

$$w_{i,j}=\frac{-{\ell }_{ij}^2+{\ell }_{jk}^2+{\ell }_{ik}^2}{\mathcal{A}(i,j,k)}+\frac{-{\ell }_{ij}^2+{\ell }_{jh}^2+{\ell }_{ih}^2}{\mathcal{A}(i,j,h)},$$

$$\tau (w_{i,j})\le w_{i,j}{\left[1-C|\nabla \tau {|}_{\infty }{\alpha }_{\text{min}}^{-2}\right]}^{-1}+2|\nabla \tau {|}_{\infty }{\left[1-C|\nabla \tau {|}_{\infty }{\alpha }_{\text{min}}^{-2}\right]}^{-1}\left(\frac{{\ell }_{ij}^2+{\ell }_{jk}^2+{\ell }_{ik}^2}{\mathcal{A}(i,j,k)}+\frac{{\ell }_{ij}^2+{\ell }_{jh}^2+{\ell }_{ih}^2}{\mathcal{A}(i,j,h)}\right),$$

$$\tau (w_{i,j})\ge w_{i,j}{\left[1+C|\nabla \tau {|}_{\infty }{\alpha }_{\text{min}}^{-2}\right]}^{-1}-2|\nabla \tau {|}_{\infty }{\left[1+C|\nabla \tau {|}_{\infty }{\alpha }_{\text{min}}^{-2}\right]}^{-1}\left(\frac{{\ell }_{ij}^2+{\ell }_{jk}^2+{\ell }_{ik}^2}{\mathcal{A}(i,j,k)}+\frac{{\ell }_{ij}^2+{\ell }_{jh}^2+{\ell }_{ih}^2}{\mathcal{A}(i,j,h)}\right),$$

Full text

Despite offering excellent sample complexity and built-in invariances, intrinsic geometry alone is invariant to isometric deformations, making it unsuitable for many applications. To overcome this limitation, we propose several upgrades to GNNs to leverage extrinsic differential geometry properties of three-dimensional surfaces, increasing its modeling power. In particular, we propose to exploit the Dirac operator, whose spectrum detects principal curvature directions — this is in stark contrast with the classical Laplace operator, which directly measures mean curvature. We coin the resulting models Surface Networks (SN).

We prove that these models define shape representations that are stable to deformation and to discretization, and we demonstrate the efficiency and versatility of SNs on two challenging tasks: temporal prediction of mesh deformations under non-linear dynamics and generative models using a variational autoencoder framework with encoders/decoders given by SNs.

1 Introduction

3D geometry analysis, manipulation and synthesis plays an important role in a variety of applications from engineering to computer animation to medical imaging. Despite the vast amount of high-quality 3D geometric data available, data-driven approaches to problems involving complex geometry have yet to become mainstream, in part due to the lack of data representation regularity which is required for traditional convolutional neural network approaches. While in computer vision problems inputs are typically sampled on regular two or three-dimensional grids, surface geometry is represented in a more complex form and, in general, cannot be converted to an image-like format by parametrizing the shape using a single planar chart. Most commonly an irregular triangle mesh is used to represent shapes, capturing its main topological and geometrical properties.

Similarly to the regular grid case (used for images or videos), we are interested in data-driven representations that strike the right balance between expressive power and sample complexity. In the case of CNNs, this is achieved by exploiting the inductive bias that most computer vision tasks are locally stable to deformations, leading to localized, multiscale, stationary features. In the case of surfaces, we face a fundamental modeling choice between extrinsic versus intrinsic representations. Extrinsic representations rely on the specific embedding of surfaces within a three-dimensional ambient space, whereas intrinsic representations only capture geometric properties specific to the surface, irrespective of its parametrization. Whereas the former offer arbitrary representation power, they are unable to easily exploit inductive priors such as stability to local deformations and invariance to global transformations.

A particularly simple and popular extrinsic method [36, 37] represents shapes as point clouds in ${\mathbb{R}}^{3}$ of variable size, and leverages recent deep learning models that operate on input sets [44, 43]. Despite its advantages in terms of ease of data acquisition (they no longer require a mesh triangulation) and good empirical performance on shape classification and segmentation tasks, one may wonder whether this simplification comes at a loss of precision as one considers more challenging prediction tasks.

In this paper, we develop an alternative pipeline that applies neural networks directly on triangle meshes, building on geometric deep learning. These models provide data-driven intrinsic graph and manifold representations with inductive biases analogous to CNNs on natural images. Models based on Graph Neural Networks [40] and their spectral variants [6, 11, 26] have been successfully applied to geometry processing tasks such as shape correspondence [32]. In their basic form, these models learn a deep representation over the discretized surface by combining a latent representation at a given node with a local linear combination of its neighbors’ latent representations, and a point-wise nonlinearity. Different models vary in their choice of linear operator and point-wise nonlinearity, which notably includes the graph Laplacian, leading to spectral interpretations of those models.

Our contributions are three-fold. First, we extend the model to support extrinsic features. More specifically, we exploit the fact that surfaces in ${\mathbb{R}}^{3}$ admit a first-order differential operator, the Dirac operator, that is stable to discretization, provides a direct generalization of Laplacian-based propagation models, and is able to detect principal curvature directions [9, 18]. Next, we prove that the models resulting from either Laplace or Dirac operators are stable to deformations and to discretization, two major sources of variability in practical applications. Last, we introduce a generative model for surfaces based on the variational autoencoder framework [25, 39], that is able to exploit non-Euclidean geometric regularity.

By combining the Dirac operator with input coordinates, we obtain a fully differentiable, end-to-end feature representation that we apply to several challenging tasks. The resulting Surface Networks – using either the Dirac or the Laplacian, inherit the stability and invariance properties of these operators, thus providing data-driven representations with certified stability to deformations. We demonstrate the model efficiency on a temporal prediction task of complex dynamics, based on a physical simulation of elastic shells, which confirms that whenever geometric information (in the form of a mesh) is available, it can be leveraged to significantly outperform point-cloud based models.

Our main contributions are summarized as follows:

• •

We demonstrate that Surface Networks provide accurate temporal prediction of surfaces under complex non-linear dynamics, motivating the use of geometric shape information.

• •

We prove that Surface Networks define shape representations that are stable to deformation and to discretization.

• •

We introduce a generative model for 3D surfaces based on the variational autoencoder.

A reference implementation of our algorithm is available at https://github.com/jiangzhongshi/SurfaceNetworks.

2 Related Work

Learning end-to-end representations on irregular and non-Euclidean domains is an active and ongoing area of research. [40] introduced graph neural networks as recursive neural networks on graphs, whose stationary distributions could be trained by backpropagation. Subsequent works [27, 43] have relaxed the model by untying the recurrent layer weights and proposed several nonlinear updates through gating mechanisms. Graph neural networks are in fact natural generalizations of convolutional networks to non-Euclidean graphs. [6, 17] proposed to learn smooth spectral multipliers of the graph Laplacian, albeit with high computational cost, and [11, 26] resolved the computational bottleneck by learning polynomials of the graph Laplacian, thus avoiding the computation of eigenvectors and completing the connection with GNNs. We refer the reader to [5] for an exhaustive literature review on the topic. GNNs are finding application in many different domains. [2, 7] develop graph interaction networks that learn pairwise particle interactions and apply them to discrete particle physical dynamics. [12, 22] study molecular fingerprints using variants of the GNN architecture, and [14] further develop the model by combining it with set representations [44], showing state-of-the-art results on molecular prediction. The resulting models, so-called Message-Passing Neural Networks, also learn the diffusion operator, which can be seen as generalizations of the Dirac model on general graphs.

In the context of computer graphics, [30] developed the first CNN model on meshed surfaces using intrinsic patch representations, and further generalized in [4] and [32]. This last work allows for flexible representations via the so-called pseudo-coordinates and obtains state-of-the-art results on 3D shape correspondence, although it does not easily encode first-order differential information. These intrinsic models contrast with Euclidean models such as [48, 46], that have higher sample complexity, since they need to learn the underlying invariance of the surface embedding. Point-cloud based models are increasingly popular to model 3d objects due to their simplicity and versatility. [36, 37] use set-invariant representations from [44, 43] to solve shape segmentation and classification tasks. More recently, [29] proposes to learn surface convolutional network from a canonical representation of planar flat-torus, with excellent performance on shape segmentation and classification, although such canonical representations may introduce exponential scale changes that can introduce instabilities. Finally, [13] proposes a point-cloud generative model for 3D shapes, that incorporates invariance to point permutations, but does not encode geometrical information as our shape generative model. Learning variational deformations is an important problem for graphics applications, since it enables negligible and fixed per-frame cost [35], but it is currently limited to 2D deformations using point handles. In constrast, our method easily generalizes to 3D and learns dynamic behaviours.

3 Surface Networks

This section presents our surface neural network model and its basic properties. We start by introducing the problem setup and notations using the Laplacian formalism (Section 3.1), and then introduce our model based on the Dirac operator (Section 3.2).

3.1 Laplacian Surface Networks

Our first goal is to define a trainable representation of discrete surfaces. Let ${\mathcal{M}}=\{V,E,F\}$ be a triangular mesh, where $V=(v_{i}\in{\mathbb{R}}^{3})_{i\leq N}$ contains the node coordinates, $E=(e_{i,j})$ corresponds to edges, and $F$ is the set of triangular faces. We denote as $\Delta$ the discrete Laplace-Beltrami operator (we use the popular cotangent weights formulation, see [5] for details).

This operator can be interpreted as a local, linear high-pass filter in ${\mathcal{M}}$ that acts on signals $x\in{\mathbb{R}}^{d\times|V|}$ defined on the vertices as a simple matrix multiplication $\tilde{x}=\Delta x$. By combining $\Delta$ with an all-pass filter and learning generic linear combinations followed by a point-wise nonlinearity, we obtain a simple generalization of localized convolutional operators in ${\mathcal{M}}$ that update a feature map from layer $k$ to layer $k+1$ using trainable parameters $A_{k}$ and $B_{k}$:

[TABLE]

By observing that the Laplacian itself can be written in terms of the graph weight similarity by diagonal renormalization, this model is a specific instance of the graph neural network [40, 5, 26] and a generalization of the spectrum-free Laplacian networks from [11]. As shown in these previous works, convolutional-like layers (1) can be combined with graph coarsening or pooling layers.

In contrast to general graphs, meshes contain a low-dimensional Euclidean embedding that contains potentially useful information in many practical tasks, despite being extrinsic and thus not invariant to the global position of the surface. A simple strategy to strike a good balance between expressivity and invariance is to include the node canonical coordinates as input channels to the network: $x^{1}:=V\in{\mathbb{R}}^{|V|\times 3}$. The mean curvature can be computed by applying the Laplace operator to the coordinates of the vertices:

[TABLE]

where $H$ is the mean curvature function and ${\bf n}(u)$ is the normal vector of the surface at point $u$. As a result, the Laplacian neural model (1) has access to mean curvature and normal information. Feeding Euclidean embedding coordinates into graph neural network models is related to the use of generalized coordinates from [32]. By cascading $K$ layers of the form (1) we obtain a representation $\Phi_{\Delta}({\mathcal{M}})$ that contains generic features at each node location. When the number of layers $K$ is of the order of $\text{diam}({\mathcal{M}})$, the diameter of the graph determined by ${\mathcal{M}}$, then the network is able to propagate and aggregate information across the whole surface.

Equation (2) illustrates that a Laplacian layer is only able to extract isotropic high-frequency information, corresponding to the mean variations across all directions. Although in general graphs there is no well-defined procedure to recover anisotropic local variations, in the case of surfaces some authors ([4, 1, 32]) have considered anisotropic extensions. We describe next a particularly simple procedure to increase the expressive power of the network using a related operator from quantum mechanics: the Dirac operator, that has been previously used successfully in the context of surface deformation [9] and shape analysis [18].

3.2 Dirac Surface Networks

The Laplace-Beltrami operator $\Delta$ is a second-order differential operator, constructed as $\Delta=-\text{div}\nabla$ by combining the gradient (a first-order differential operator) with its adjoint, the divergence operator. In an Euclidean space, one has access to these first-order differential operators separately, enabling oriented high-pass filters.

For convenience, we embed ${\mathbb{R}}^{3}$ to the imaginary quaternion space $\text{Im}({\mathbb{H}})$ (see Appendix A in the Suppl. Material for details). The Dirac operator is then defined as a matrix $D\in{\mathbb{H}}^{|F|\times|V|}$ that maps (quaternion) signals on the nodes to signals on the faces. In coordinates,

[TABLE]

where $e_{j}$ is the opposing edge vector of node $j$ in the face $f$, and ${\mathcal{A}}_{f}$ is the area (see Appendix A) using counter-clockwise orientations on all faces.

To apply the Dirac operator defined in quaternions to signals in vertices and faces defined in real numbers, we write the feature vectors as quaternions by splitting them into chunks of 4 real numbers representing the real and imaginary parts of a quaternion; see Appendix A. Thus, we always work with feature vectors with dimensionalities that are multiples of $4$. The Dirac operator provides first-order differential information and is sensitive to local orientations. Moreover, one can verify [9] that

[TABLE]

where $D^{*}$ is the adjoint operator of $D$ in the quaternion space (see Appendix A). The adjoint matrix can be computed as $D^{*}=M^{-1}_{V}D^{H}M_{F}$ where $D^{H}$ is a conjugate transpose of $D$ and $M_{V}$, $M_{F}$ are diagonal mass matrices with one third of areas of triangles incident to a vertex and face areas respectively.

The Dirac operator can be used to define a new neural surface representation that alternates layers with signals defined over nodes with layers defined over faces. Given a $d$-dimensional feature representation over the nodes $x^{k}\in{\mathbb{R}}^{d\times|V|}$, and the faces of the mesh, $y^{k}\in{\mathbb{R}}^{d\times|F|}$, we define a $d^{\prime}$-dimensional mapping to a face representation as

[TABLE]

where $C_{k},E_{k}$ are trainable parameters. Similarly, we define the adjoint layer that maps back to a $\tilde{d}$-dimensional signal over nodes as

[TABLE]

where $A_{k},B_{k}$ are trainable parameters. A surface neural network layer is thus determined by parameters $\{A,B,C,E\}$ using equations (3) and (4) to define $x^{k+1}\in{\mathbb{R}}^{d_{k+1}\times|V|}$. We denote by $\Phi_{D}({\mathcal{M}})$ the mesh representation resulting from applying $K$ such layers (that we assume fixed for the purpose of exposition).

The Dirac-based surface network is related to edge feature transforms proposed on general graphs in [14], although these edge measurements cannot be associated with derivatives due to lack of proper orientation. In general graphs, there is no notion of square root of $\Delta$ that recovers oriented first-order derivatives.

4 Stability of Surface Networks

Here we describe how Surface Networks are geometrically stable, because surface deformations become additive noise under the model. Given a continuous surface $S\subset{\mathbb{R}}^{3}$ or a discrete mesh ${\mathcal{M}}$, and a smooth deformation field $\tau:{\mathbb{R}}^{3}\to{\mathbb{R}}^{3}$, we are particularly interested in two forms of stability:

• •

Given a discrete mesh ${\mathcal{M}}$ and a certain non-rigid deformation $\tau$ acting on ${\mathcal{M}}$, we want to certify that $\|\Phi({\mathcal{M}})-\Phi(\tau({\mathcal{M}}))\|$ is small if $\|\nabla\tau(\nabla\tau)^{*}-{\bf I}\|$ is small, i.e when the deformation is nearly rigid; see Theorem 4.1.

• •

Given two discretizations ${\mathcal{M}}_{1}$ and ${\mathcal{M}}_{2}$ of the same underlying surface $S$, we would like to control $\|\Phi({\mathcal{M}}_{1})-\Phi({\mathcal{M}}_{2})\|$ in terms of the resolution of the meshes; see Theorem 4.2.

These stability properties are important in applications, since most tasks we are interested in are stable to deformation and to discretization. We shall see that the first property is a simple consequence of the fact that the mesh Laplacian and Dirac operators are themselves stable to deformations. The second property will require us to specify under which conditions the discrete mesh Laplacian $\Delta_{\mathcal{M}}$ converges to the Laplace-Beltrami operator $\Delta_{S}$ on $S$. Unless it is clear from the context, in the following $\Delta$ will denote the discrete Laplacian.

Theorem 4.1

Let ${\mathcal{M}}$ be a $N$-node mesh and $x,\,x^{\prime}\in{\mathbb{R}}^{|V|\times d}$ be input signals defined on the nodes. Assume the nonlinearity $\rho(\,\cdot\,)$ is non-expansive ($|\rho(z)-\rho(z^{\prime})|\leq|z-z^{\prime}|$). Then

1. (a)

$\|\Phi_{\Delta}({\mathcal{M}};x)-\Phi_{\Delta}({\mathcal{M}};x^{\prime})\|\leq\alpha_{\Delta}\|x-x^{\prime}\|~{},$* where $\alpha_{\Delta}$ depends only on the trained weights and the mesh.* 2. (b)

$\|\Phi_{D}({\mathcal{M}};x)-\Phi_{D}({\mathcal{M}};x^{\prime})\|\leq\alpha_{D}\|x-x^{\prime}\|~{},$* where $\alpha_{D}$ depends only on the trained weights and the mesh.* 3. (c)

Let ${|\nabla\tau|_{\infty}}:=\sup_{u}\|\nabla\tau(u)(\nabla\tau(u))^{*}-{\bf 1}\|$, where $\nabla\tau(u)$ is the Jacobian matrix of $u\mapsto\tau(u)$. Then $\|\Phi_{\Delta}({\mathcal{M}};x)-\Phi_{\Delta}(\tau({\mathcal{M}});x)\|\leq\beta_{\Delta}{|\nabla\tau|_{\infty}}\|x\|~{},$ where $\beta_{\Delta}$ is independent of $\tau$ and $x$. 4. (d)

Denote by $\widetilde{|\nabla\tau|}_{\infty}:=\sup_{u}\|\nabla\tau(u)-{\bf 1}\|$. Then $\|\Phi_{D}({\mathcal{M}};x)-\Phi_{D}(\tau({\mathcal{M}});x)\|\leq\beta_{D}\widetilde{|\nabla\tau|}_{\infty}\|x\|~{},$ where $\beta_{D}$ is independent of $\tau$ and $x$.

Properties (a) and (b) are not specific to surface representations, and are a simple consequence of the non-expansive property of our chosen nonlinearities. The constant $\alpha$ is controlled by the product of $\ell_{2}$ norms of the network weights at each layer and the norm of the discrete Laplacian operator. Properties (c) and (d) are based on the fact that the Laplacian and Dirac operators are themselves stable to deformations, a property that depends on two key aspects: first, the Laplacian/Dirac is localized in space, and next, that it is a high-pass filter and therefore only depends on relative changes in position.

One caveat of Theorem 4.1 is that the constants appearing in the bounds depend upon a bandwidth parameter given by the reciprocal of triangle areas, which increases as the size of the mesh increases. This corresponds to the fact that the spectral radius of $\Delta_{\mathcal{M}}$ diverges as the mesh size $N$ increases.

In order to overcome this problematic asymptotic behavior, it is necessary to exploit the smoothness of the signals incoming to the surface network. This can be measured with Sobolev norms defined using the spectrum of the Laplacian operator. Given a mesh ${\mathcal{M}}$ of $N$ nodes approximating an underlying surface $S$, and its associated cotangent Laplacian $\Delta_{\mathcal{M}}$, consider the spectral decomposition of $\Delta_{\mathcal{M}}$ (a symmetric, positive definite operator):

[TABLE]

Under normal uniform convergence 111which controls how the normals of the mesh align with the surface normals; see [45]. [45], the spectrum of $\Delta_{\mathcal{M}}$ converges to the spectrum of the Laplace-Beltrami operator $\Delta_{S}$ of $S$. If $S$ is bounded, it is known from the Weyl law [47] that there exists $\gamma>0$ such that $k^{-\gamma(S)}\lesssim\lambda_{k}^{-1}$, so the eigenvalues $\lambda_{k}$ do not grow too fast. The smoothness of a signal $x\in{\mathbb{R}}^{|V|\times d}$ defined in ${\mathcal{M}}$ is captured by how fast its spectral decomposition $\hat{x}(k)=e_{k}^{T}x\in{\mathbb{R}}^{d}$ decays [42]. We define $\|x\|_{{\mathcal{H}}}^{2}:=\sum_{k}\lambda(k)^{2}\|\hat{x}(k)\|^{2}$ is Sobolev norm, and $\beta(x,S)>1$ as the largest rate such that its spectral decomposition coefficients satisfy

[TABLE]

If $x\in{\mathbb{R}}^{|V|\times d}$ is the input to the Laplace Surface Network of $R$ layers, we denote by $(\beta_{0},\beta_{1},\dots,\beta_{R-1})$ the smoothness rates of the feature maps $x^{(r)}$ defined at each layer $r\leq R$.

Theorem 4.2

Consider a surface $S$ and a finite-mesh approximation ${\mathcal{M}}_{N}$ of $N$ points, and $\Phi_{\Delta}$ a Laplace Surface Network with parameters $\{(A_{r},B_{r})\}_{r\leq R}$. Denote by $d(S,{\mathcal{M}}_{N})$ the uniform normal distance, and let $x_{1},x_{2}$ be piece-wise polyhedral approximations of $\bar{x}(t)$, $t\in S$ in ${\mathcal{M}}_{N}$, with $\|\bar{x}\|_{{\mathcal{H}}(S)}<\infty$. Assume $\|\bar{x}^{(r)}\|_{{\mathcal{H}}(S)}<\infty$ for $r\leq R$.

1. (a)

If $x_{1},x_{2}$ are two functions such that the $R$ feature maps $x_{l}^{(r)}$ have rates $(\beta_{0},\beta_{1},\dots,\beta_{R-1})$, then

[TABLE]

with $h(\beta)={\prod_{r=1}^{R}\frac{\beta_{r}-1}{\beta_{r}-1/2}}$, and where $C(\beta)$ does not depend upon $N$. 2. (b)

If $\tau$ is a smooth deformation field, then $\|\Phi_{\Delta}(x;{\mathcal{M}}_{N})-\Phi_{\Delta}(x;\tau({\mathcal{M}}_{N}))\|\leq C{|\nabla\tau|_{\infty}}^{h(\beta)}~{},$ where $C$ does not depend upon $N$. 3. (c)

Let ${\mathcal{M}}$ and ${\mathcal{M}}^{\prime}$ be $N$-point discretizations of $S$, If $\max(d({\mathcal{M}},S),d({\mathcal{M}}^{\prime},S))\leq\epsilon$, then $\|\Phi_{\Delta}({\mathcal{M}};x)-\Phi_{\Delta}({\mathcal{M}}^{\prime},x^{\prime})\|\leq C\epsilon^{h(\beta)}~{},$ where $C$ is independent of $N$.

This result ensures that if we use as generator of the SN an operator that is consistent as the mesh resolution increases, the resulting surface representation is also consistent. Although our present result only concerns the Laplacian, the Dirac operator also has a well-defined continuous counterpart [9] that generalizes the gradient operator in quaternion space. Also, our current bounds depend explicitly upon the smoothness of feature maps across different layers, which may be controlled in terms of the original signal if one considers nonlinearities that demodulate the signal, such as $\rho(x)=|x|$ or $\rho(x)=\text{ReLU}(x)$. These extensions are left for future work. Finally, a specific setup that we use in experiments is to use as input signal the canonical coordinates of the mesh ${\mathcal{M}}$. In that case, an immediate application of the previous theorem yields

Corollary 4.3

Denote $\Phi({\mathcal{M}}):=\Phi_{{\mathcal{M}}}(V)$, where $V$ are the node coordinates of ${\mathcal{M}}$. Then, if $A_{1}=0$,

[TABLE]

5 Generative Surface Models

State-of-the-art generative models for images, such as generative adversarial networks [38], pixel autoregressive networks [34], or variational autoencoders [25], exploit the locality and stationarity of natural images in their probabilistic models, in the sense that the model satisfies $p_{\theta}(x)\approx p_{\theta}(x_{\tau})$ by construction, where $x_{\tau}$ is a small deformation of a given input $x$. This property is obtained via encoders and decoders with a deep convolutional structure. We intend to exploit similar geometric stability priors with SNs, owing to their stability properties described in Section 4. A mesh generative model contains two distinct sources of randomness: on the one hand, the randomness associated with the underlying continuous surface, which corresponds to shape variability; on the other hand, the randomness of the discretization of the surface. Whereas the former contains the essential semantic information, the latter is not informative, and to some extent independent of the shape identity. We focus initially on meshes that can be represented as a depth map over an (irregular) 2D mesh, referred as height-field meshes in the literature. That is, a mesh ${\mathcal{M}}=(V,E,F)$ is expressed as $(\tilde{{\mathcal{M}}},f(\tilde{{\mathcal{M}}}))$, where $\tilde{{\mathcal{M}}}=(\tilde{V},\tilde{E},\tilde{F})$ is now a 2D mesh and $f~{}:~{}\tilde{V}\to{\mathbb{R}}$ is a depth-map encoding the original node locations $V$, as shown in Figure 1.

In this work, we consider the variational autoencoder framework [25, 39]. It considers a mixture model of the form $p({\mathcal{M}})=\int p_{\theta}({\mathcal{M}}~{}|~{}h)p_{0}(h)dh~{}$, where $h\in{\mathbb{R}}^{|}S|$ is a vector of latent variables. We train this model by optimizing the variational lower bound of the data log-likelihood:

[TABLE]

We thus need to specify a conditional generative model $p_{\theta}({\mathcal{M}}~{}|~{}h)$, a prior distribution $p_{0}(h)$ and a variational approximation to the posterior $q_{\psi}(h~{}|~{}{\mathcal{M}})$, where $\theta$ and $\psi$ denote respectively generative and variational trainable parameters. Based on the height-field representation, we choose for simplicity a separable model of the form $p_{\theta}({\mathcal{M}}~{}|~{}h)=p_{\theta}(f~{}|~{}h,\tilde{{\mathcal{M}}})\cdot p(\tilde{{\mathcal{M}}})~{},$ where $\tilde{{\mathcal{M}}}\sim p(\tilde{{\mathcal{M}}})$ is a homogeneous Poisson point process, and $f\sim p_{\theta}(f~{}|~{}h,\tilde{{\mathcal{M}}})$ is a normal distribution with mean and isotropic covariance parameters given by a SN:

[TABLE]

with $[\mu(h,\tilde{{\mathcal{M}}}),\sigma^{2}(h,\tilde{{\mathcal{M}}})]=\Phi_{D}(\tilde{M};h)~{}.$ The generation step thus proceeds as follows. We first sample a 2D mesh $\tilde{{\mathcal{M}}}$ independent of the latent variable $h$, and then sample a depth field over $\tilde{{\mathcal{M}}}$ conditioned on $h$ from the output of a decoder network $\Phi_{D}(\tilde{M};h)$. Finally, the variational family $q_{\psi}$ is also a Normal distribution whose parameters are obtained from an encoder Surface Neural Network whose last layer is a global pooling that removes the spatial localization: $q_{\psi}(h~{}|~{}{\mathcal{M}})={\cal N}(\bar{\mu},\bar{\sigma}^{2}{\bf 1})~{},~{}\text{ with }~{}[\bar{\mu},\bar{\sigma}]=\bar{\Phi}_{D}({\mathcal{M}})~{}.$

6 Experiments

For experimental evaluation, we compare models built using ResNet-v2 blocks [16], where convolutions are replaced with the appropriate operators (see Fig. 2): (i) a point cloud based model from [43] that aggregates global information by averaging features in the intermediate layers and distributing them to all nodes; (ii) a Laplacian Surface network with input canonical coordinates; (iii) a Dirac Surface Network model. We report experiments on generative models using an unstructured variant of MNIST digits (Section 6.1), and on temporal prediction under non-rigid deformation models (Section 6.2).

6.1 MeshMNIST

For this task, we construct a MeshMNIST database with only height-field meshes (Sec. 5). First, we sample points on a 2D plane $([0,27]\times[0,27])$ with Poisson disk sampling with $r=1.0$, which roughly generates $500$ points, and apply a Delaunay triangulation to these points. We then overlay the triangulation with the original MNIST images and assign to each point a $z$ coordinate bilinearly interpolating the grey-scale value. Thus, the procedure allows us to define a sampling process over 3D height-field meshes.

We used VAE models with decoders and encoders built using 10 ResNet-v2 blocks with 128 features. The encoder converts a mesh into a latent vector by averaging output of the last ResNet-v2 block and applying linear transformations to obtain mean and variance, while the decoder takes a latent vector and a 2D mesh as input (corresponding to a specific 3D mesh) and predicts offsets for the corresponding locations. We keep variance of the decoder as a trainable parameter that does not depend on input data. We trained the model for $75$ epochs using Adam optimizer [24] with learning rate $10^{-3}$, weight decay $10^{-5}$ and batch size $32$. Figures 3,4 illustrate samples from the model. The geometric encoder is able to leverage the local translation invariance of the data despite the irregular sampling, whereas the geometric decoder automatically adapts to the specific sampled grid, as opposed to set-based generative models.

6.2 Spatio-Temporal Predictions

One specific task we consider is temporal predictions of non-linear dynamics. Given a sequence of frames $X={X^{1},X^{2},\ldots,X^{n}}$, the task is to predict the following frames $Y={Y^{1}=X^{n+1},Y^{2},\ldots,Y^{m}=X^{n+m}}$. As in [31], we use a simple non-recurrent model that takes a concatenation of input frames $X$ and predicts a concatenation of frames $Y$. We condition on $n=2$ frames and predict the next $m=40$ frames. In order to generate data, we first extracted 10k patches from the MPI-Faust dataset[3], by selecting a random point and growing a topological sphere of radius $15$ edges (i.e. the 15-ring of the point). For each patch, we generate a sequence of $50$ frames by randomly rotating it and letting it fall to the ground. We consider the mesh a thin elastic shell, and we simulate it using the As-Rigid-As-Possible technique [41], with additional gravitational forces [20]. Libigl [21] has been used for the mesh processing tasks. Sequences with patches from the first 80 subjects were used in training, while the 20 last subjects were used for testing. The dataset and the code are available on request. We restrict our experiments to temporal prediction tasks that are deterministic when conditioned on several initial frames. Thus, we can train models by minimizing smooth-L1 loss [15] between target frames and output of our models.

We used models with 15 ResNet-v2 blocks with $128$ output features each. In order to cover larger context for Dirac and Laplace based models, we alternate these blocks with Average Pooling blocks. We predict offsets to the last conditioned frame and use the corresponding Laplace and Dirac operators. Thus, the models take 6-dimensional inputs and produce 120-dimensional outputs. We trained all models using the Adam optimizer [24] with learning rate $10^{-3}$, weight decay $10^{-5}$, and batch size 32. After 60k steps we decreased the learning rate by a factor of 2 every 10k steps. The models were trained for 110k steps in overall.

Table 1 reports quantitative prediction performance of different models, and Figure 5 displays samples from the prediction models at specific frames. The set-to-set model [44, 43], corresponding to a point-cloud representation used also in [36], already performs reasonably well on the task, even if the visual difference is noticeable. Nevertheless, the gap between this model and Laplace-/Dirac-based models is significant, both visually and quantitatively. Dirac-based model outperforms Laplace-based model despite the smaller receptive field. Videos comparing the performance of different models are available in the additional material.

Figure 6 illustrates the effect of replacing Laplace by Dirac in the formulation of the SN. Laplacian-based models, since they propagate information using an isotropic operator, have more difficulties at resolving corners and pointy structures than the Dirac operator, that is sensitive to principal curvature directions. However, the capacity of Laplace models to exploit the extrinsic information only via the input coordinates is remarkable and more computationally efficient than the Dirac counterpart. Figures 7 and 8 overlay the prediction error and compare Laplace against Dirac and PointCloud against Dirac respectively. They confirm first that SNs outperform the point-cloud based model, which often produce excessive flattening and large deformations, and next that first-order Dirac operators help resolve areas with high directional curvature. We refer to the supplementary material for additional qualitative results.

7 Conclusions

We have introduced Surface Networks, a deep neural network that is designed to naturally exploit the non-Euclidean geometry of surfaces. We have shown how a first-order differential operator (the Dirac operator) can detect and adapt to geometric features beyond the local mean curvature, the limit of what Laplacian-based methods can exploit. This distinction is important in practice, since areas with high directional curvature are perceptually important, as shown in the experiments. That said, the Dirac operator comes at increased computational cost due to the quaternion calculus, and it would be interesting to instead learn the operator, akin to recent Message-Passing NNs [14] and explore whether Dirac is recovered.

Whenever the data contains good-quality meshes, our experiments demonstrate that using intrinsic geometry offers vastly superior performance to point-cloud based models. While there are not many such datasets currently available, we expect them to become common in the next years, as scanning and reconstruction technology advances and 3D sensors are integrated in consumer devices. SNs provide efficient inference, with predictable runtime, which makes them appealing across many areas of computer graphics, where a fixed, per-frame cost is required to ensure a stable framerate, especially in VR applications. Our future plans include applying Surface Networks precisely to having automated, data-driven mesh processing, and generalizing the generative model to arbitrary meshes, which will require an appropriate multi-resolution pipeline.

Appendix A The Dirac Operator

The quaternions ${\mathbb{H}}$ is an extension of complex numbers. A quaternion $q\in{\mathbb{H}}$ can be represented in a form $q=a+bi+cj+dk$ where $a,b,c,d$ are real numbers and ${i,j,k}$ are quaternion units that satisfy the relationship $i^{2}=j^{2}=k^{2}=ijk=-1$.

As mentioned in Section 3.1, the Dirac operator used in the model can be conveniently represented as a quaternion matrix:

[TABLE]

where $e_{j}$ is the opposing edge vector of node $j$ in the face $f$, and ${\mathcal{A}}_{f}$ is the area, as illustrated in Fig. A, using counter-clockwise orientations on all faces.

The Deep Learning library PyTorch that we used to implement the models does not support quaternions. Nevertheless, quaternion-valued matrix multiplication can be replaced with real-valued matrix multiplication where each entry $q=a+bi+cj+dk$ is represented as a $4\times 4$ block

[TABLE]

and the conjugate $q^{*}=a-bi-cj-dk$ is a transpose of this real-valued matrix:

[TABLE]

Appendix B Theorem 4.1

B.1 Proof of (a)

We first show the result for the mapping $x\mapsto\rho\left(Ax+B\Delta x\right)$, corresponding to one layer of $\Phi_{\Delta}$. By definition, the Laplacian $\Delta$ of ${\mathcal{M}}$ is

[TABLE]

where $\bar{{\mathcal{A}}}_{j}$ is one third of the total area of triangles incident to node $j$, and $W=(w_{i,j})$ contains the cotangent weights [45], and $U=\text{diag}(W{\bf 1})$ contains the node aggregated weights in its diagonal.

From [10] we verify that

[TABLE]

where $d_{j}$ denotes the degree (number of neighbors) of node $j$, $\alpha_{\text{min}}$ is the smallest angle in the triangulation of ${\mathcal{M}}$ and $S_{\text{max}}$ the largest number of incident triangles. It results that

[TABLE]

which depends uniquely on the mesh ${\mathcal{M}}$ and is finite for non-degenerate meshes. Moreover, since $\rho(\,\cdot\,)$ is non-expansive, we have

[TABLE]

By cascading (B.1) across the $K$ layers of the network, we obtain

[TABLE]

which proves (a). $\square$

B.2 Proof of (b)

The proof is analogous, by observing that $\|D\|=\sqrt{\|\Delta\|}$ and therefore

[TABLE]

B.3 Proof of (c)

To establish (c) we first observe that given three points $p,q,r\in{\mathbb{R}}^{3}$ forming any of the triangles of ${\mathcal{M}}$,

[TABLE]

Indeed, (11) is a direct consequence of the lower and upper Lipschitz constants of $\tau(u)$, which are bounded respectively by $1-{|\nabla\tau|_{\infty}}$ and $1+{|\nabla\tau|_{\infty}}$. As for (12), we use the Heron formula

[TABLE]

with $s=\frac{1}{2}(\|p-q\|+\|p-r\|+\|r-q\|)$ being the half-perimeter. By denoting $s_{\tau}$ the corresponding half-perimeter determined by the deformed points $\tau(p),\tau(q),\tau(r)$, we have that

[TABLE]

[TABLE]

and similarly for the $\|r-q\|$ and $\|r-p\|$ terms. It results in

[TABLE]

and similarly

[TABLE]

By noting that the cotangent Laplacian weights can be written (see Fig. 9) as

[TABLE]

we have from the previous Bilipschitz bounds that

[TABLE]

[TABLE]

which proves that, up to second order terms, the cotangent weights are Lipschitz continuous to deformations.

Finally, since the mesh Laplacian operator is constructed as $\text{diag}(\bar{{\mathcal{A}}})^{-1}(U-W)$, with $\bar{{\mathcal{A}}}_{i,i}=\frac{1}{3}\sum_{j,k;(i,j,k)\in F}{\mathcal{A}}(i,j,k)$, and $U=\text{diag}(W{\bf 1})$, let us show how to bound $\|\Delta-\tau(\Delta)\|$ from

[TABLE]

and

[TABLE]

Using the fact that $\bar{{\mathcal{A}}}$, $\tau(\bar{{\mathcal{A}}})$ are diagonal, and using the spectral bound for $k\times m$ sparse matrices from [8], Lemma 5.12,

[TABLE]

the bounds (13) and (14) yield respectively

[TABLE]

It results that, up to second order terms,

[TABLE]

which shows that the Laplacian is stable to deformations in operator norm. Finally, by denoting $\tilde{x}_{\tau}$ a layer of the deformed Laplacian network

[TABLE]

it follows that

[TABLE]

Also,

[TABLE]

and therefore, by plugging (17) with $y=\tilde{x}_{\tau}$, $K$ layers of the Laplacian network satisfy

[TABLE]

B.4 Proof of (d)

The proof is also analogous to the proof of (c), with the difference that now the Dirac operator is no longer invariant to orthogonal transformations, only to translations. Given two points $p$, $q$, we verify that

[TABLE]

which, following the previous argument, leads to

[TABLE]

Appendix C Theorem 4.2

C.1 Proof of part (a)

The proof is based on the following lemma:

Lemma C.1

Let $x_{N},y_{N}\in\mathcal{H}({\mathcal{M}}_{N})$ such that $\forall~{}N$, $\|x_{N}\|_{{\mathcal{H}}}\leq c$,$\|y_{N}\|_{{\mathcal{H}}}\leq c$. Let $\hat{x_{N}}=\mathcal{E}_{N}(x_{N})$, where $\mathcal{E}_{N}$ is the eigendecomposition of the Laplacian operator $\Delta_{N}$ on ${\mathcal{M}}_{N}$, , with associated eigenvalues $\lambda_{1}\dots\lambda_{N}$ in increasing order. Let $\gamma>0$ and $\beta$ be defined as in (5) for $x_{N}$ and $y_{N}$. If $\beta>1$ and $\|x_{N}-y_{N}\|\leq\epsilon$ for all $N$,

[TABLE]

where $C$ is a constant independent of $\epsilon$ and $N$.

One layer of the network will transform the difference $x_{1}-x_{2}$ into $\rho(Ax_{1}+B\Delta x_{1})-\rho(Ax_{2}+B\Delta x_{2})$. We verify that

[TABLE]

We now apply Lemma C.1 to obtain

[TABLE]

where we redefine $C$ to account for the fact that $\|x_{1}-x_{2}\|^{(2\beta-1)^{-1}}$ is bounded. We have just showed that

[TABLE]

with $f_{r}=C(\|A_{r}\|+\|B_{r}\|)$ and $g_{r}=\frac{\beta_{r}-1}{\beta_{r}-1/2}$. By cascading (20) for each of the $R$ layers we thus obtain

[TABLE]

which proves (6) $\square$.

*Proof of (19): * Let $\{e_{1},\dots,e_{N}\}$ be the eigendecomposition of $\Delta_{N}$. For simplicity, we drop the subindex $N$ in the signals from now on. Let $\hat{x}(k)=\langle x,e_{k}\rangle$ and $\tilde{x}(k)=\lambda_{k}\hat{x}(k)$; and analogously for $y$. From the Parseval identity we have that $\|x\|^{2}=\|\hat{x}\|^{2}$. We express $\|\Delta(x-y)\|$ as

[TABLE]

The basic principle of the proof is to cut the spectral sum (22) in two parts, chosen to exploit the decay of $\tilde{x}(k)$. Let

[TABLE]

and analogously for $y$. For any cutoff $k_{*}\leq N$ we have

[TABLE]

where we denote for simplicity $F(k_{*})=\max(F(x)(k_{*}),~{}F(y)(k_{*}))$. By assumption, we have $\lambda_{k}^{2}\lesssim k^{2\gamma}$ and

[TABLE]

By denoting $\tilde{\beta}=\beta-\gamma-1/2$, it follows that

[TABLE]

Optimizing for $k_{*}$ yields

[TABLE]

[TABLE]

By plugging (25) back into (24) and dropping all constants independent of $N$ and $\epsilon$, this leads to

[TABLE]

which proves part (a) $\square$.

C.2 Proof of part (b)

We will use the following lemma:

Lemma C.2

Let ${\mathcal{M}}=(V,E,F)$ is a non-degenerate mesh, and define

[TABLE]

Then, given a smooth deformation $\tau$ and $x$ defined in ${\mathcal{M}}$, we have

[TABLE]

where $C$ depends only upon $\eta_{1}$, $\eta_{2}$ and $\eta_{3}$.

In that case, we need to control the difference $\rho(Ax+B\Delta x)-\rho(Ax+B\tau(\Delta)x)$. We verify that

[TABLE]

By Lemma C.2 it follows that $\|(\Delta-\tau(\Delta))x\|\leq C{|\nabla\tau|_{\infty}}\|\Delta x\|$ and therefore, by denoting $x_{1}^{(1)}=\rho(Ax+B\Delta x)$ and $x_{2}^{(1)}=\rho(Ax+B\tau(\Delta)x)$, we have

[TABLE]

By applying again Lemma C.1, we also have that

[TABLE]

which, by combining it with (28) and repeating through the $R$ layers yields

[TABLE]

which concludes the proof $\square$.

Proof of (27): The proof follows closely the proof of Theorem 4.1, part (c). From (13) and (14) we have that

[TABLE]

It follows that, up to second order $o({|\nabla\tau|_{\infty}}^{2})$ terms,

[TABLE]

By writing $\bar{{\mathcal{A}}}^{-1}H_{\tau}=\widetilde{H_{\tau}}\bar{{\mathcal{A}}}^{-1}$, and since $\bar{{\mathcal{A}}}$ is diagonal, we verify that

[TABLE]

$\frac{{\mathcal{A}}_{i,i}}{{\mathcal{A}}_{j,j}}\leq\eta_{1}$, and hence that

[TABLE]

We conclude by combining (30) and (31) into

[TABLE]

which proves (27) $\square$

C.3 Proof of part (c)

This result is a consequence of the consistency of the cotangent Laplacian to the Laplace-Beltrami operator on $S$ [45]:

Theorem C.3** ([45], Thm 3.4)**

Let ${\mathcal{M}}$ be a compact polyhedral surface which is a normal graph over a smooth surface $S$ with distortion tensor $\mathcal{T}$, and let $\bar{\mathcal{T}}=(\det\mathcal{T})^{1/2}\mathcal{T}^{-1}$. If the normal field uniform distance $d(\mathcal{T},{\bf 1})=\|\bar{\mathcal{T}}-{\bf 1}\|_{\infty}$ satisfies $d(\mathcal{T},{\bf 1})\leq\epsilon$, then

[TABLE]

If $\Delta_{\mathcal{M}}$ converges uniformly to $\Delta_{S}$, in particular we verify that

[TABLE]

Thus, given two meshes ${\mathcal{M}}$, ${\mathcal{M}}^{\prime}$ approximating a smooth surface $S$ in terms of uniform normal distance, and the corresponding irregular sampling $x$ and $x^{\prime}$ of an underlying function $\bar{x}:S\to{\mathbb{R}}$, we have

[TABLE]

Since ${\mathcal{M}}$ and ${\mathcal{M}}^{\prime}$ both converge uniformly normally to $S$ and $\bar{x}$ is Lipschitz on $S$, it results that

[TABLE]

thus $\|x-x^{\prime}\|\leq 2L\epsilon$. Also, thanks to the uniform normal convergence, we also have convergence in the Sobolev sense:

[TABLE]

which implies in particular that

[TABLE]

From (33) and (34) it follows that

[TABLE]

By applying again Lemma C.1 to $\tilde{x}=\rho(Ax+B\Delta_{{\mathcal{M}}}x)$, $\tilde{x}^{\prime}=\rho(Ax^{\prime}+B\Delta_{{\mathcal{M}}^{\prime}}x^{\prime})$, we have

[TABLE]

We conclude by retracing the same argument as before, reapplying Lemma C.1 at each layer to obtain

[TABLE]

Appendix D Proof of Corollary 4.3

We verify that

[TABLE]

The second term is $o({|\nabla\tau|_{\infty}})$ from Lemma C.2. The first term is

[TABLE]

where $\|\nabla^{2}\tau\|$ is the uniform Hessian norm of $\tau$. The result follows from applying the cascading argument from last section. $\square$

Appendix E Preliminary Study: Metric Learning for Dense Correspondence

As an interesting extension, we apply the architecture we built in Experiments 6.2 directly to a dense shape correspondence problem.

Similarly as the graph correspondence model from [33], we consider a Siamese Surface Network, consisting of two identical models with the same architecture and sharing parameters. For a pair of input surfaces $\mathcal{M}_{1},\mathcal{M}_{2}$ of $N_{1}$, $N_{2}$ points respectively, the network produces embeddings $E_{1}\in{\mathbb{R}}^{N_{1}\times d}$ and $E_{2}\in{\mathbb{R}}^{N_{2}\times d}$. These embeddings define a trainable similarity between points given by

[TABLE]

which can be trained by minimizing the cross-entropy relative to ground truth pairs. A diagram of the architecture is provided in Figure 10.

In general, dense shape correspondence is a task that requires a blend of intrinsic and extrinsic information, motivating the use of data-driven models that can obtain such tradeoffs automatically. Following the setup in Experiment 6.2, we use models with 15 ResNet-v2 blocks with 128 output features each, and alternate Laplace and Dirac based models with Average Pooling blocks to cover a larger context: The input to our network consists of vertex positions only.

We tested our architecture on a reconstructed (i.e. changing the mesh connectivity) version of the real scan of FAUST dataset[3]. The FAUST dataset contains 100 real scans and their corresponding ground truth registrations. The ground truth is based on a deformable template mesh with the same ordering and connectivity, which is fitted to the scans. In order to eliminate the bias of using the same template connectivity, as well as the need of a single connected component, the scans are reconstructed again with [19]. To foster replicability, we release the processed dataset in the additional material. In our experiment, we use 80 models for training and 20 models for testing.

Since the ground truth correspondence is implied only through the common template mesh, we compute the correspondence between our meshes with a nearest neighbor search between the point cloud and the reconstructed mesh. Consequently, due to the drastic change in vertex replacement after the remeshing, only 60-70 percent of labeled matches are used. Although making it more challenging, we believe this setup is close to a real case scenario, where acquisition noise and occlusions are unavoidable.

Our preliminary results are reported in Figure 11. For simplicity, we generate predicted correspondences by simply taking the mode of the softmax distribution for each reference node $i$: $\hat{j}(i)=\arg\max_{j}s_{i,j}$, thus avoiding a refinement step that is standard in other shape correspondence pipelines. The MLP model uses no context whatsoever and provides a baseline that captures the prior information from input coordinates alone. Using contextual information (even extrinsically as in point-cloud model) brings significative improvments, but these results may be substantially improved by encoding further prior knowledge. An example of the current failure of our model is depitcted in Figure 13, illustrating that our current architecture does not have sufficiently large spatial context to disambiguate between locally similar (but globally inconsistent) parts.

We postulate that the FAUST dataset [3] is not an ideal fit for our contribution for two reasons: (1) it is small (100 models), and (2) it contains only near-isometric deformations, which do not require the generality offered by our network. As demonstrated in [28], the correspondence performances can be dramatically improved by constructing basis that are invariant to the deformations. We look forward to the emergence of new geometric datasets, and we are currently developing a capture setup that will allow us to acquire a more challenging dataset for this task.

Appendix F Further Numerical Experiments