Representation of Surfaces with Normal Cycles and Application to Surface Registration

Abstract

In this paper, we present a framework for computing dissimilarities between surfaces which is based on the mathematical model of normal cycle from geometric measure theory. This model allows to take into account all the curvature information of the surface without explicitly computing it. By defining kernel metrics on normal cycles, we define explicit distances between surfaces that are sensitive to curvature. This mathematical framework also has the advantage of encompassing both continuous and discrete surfaces (triangulated surfaces). We then use this distance as a data attachment term for shape matching, using the large deformation diffeomorphic metric mapping for modelling deformations. We also present an efficient numerical implementation of this problem in PyTorch, using the KeOps library, which allows both the use of auto-differentiation tools and a parallelization of GPU calculations without memory overflow. We show that this method can be scalable on data up to a million points, and we present several examples on surfaces, comparing the results with those obtained with the similar varifold framework.

Appendices

Discrete Scalar Product with the Constant Kernel

1.1 Discrete Surfaces

For discrete surfaces, and with the constant normal kernel, it can be easily seen that the planar part is not involved. The scalar product of normal cycles above vertices (i.e. the spherical scalar product) involves terms as:

$$\begin{aligned} \begin{aligned}&k_p(x,y)\int _{S_1}\int _{S_2}k_n(u,v)\\&\quad \left\langle \tau _{{{\mathscr {N}}}_C}(x,u),\tau _{{{\mathscr {N}}}_C}(y,v)\right\rangle \mathrm{d}{\mathscr {H}}^2(u)\mathrm{d}{\mathscr {H}}^2(v) \\&\qquad = k_p(x,y)\left\langle \int _{S_1}\underbrace{\tau _{{{\mathscr {N}}}_C}(x,u)}_{= u \text { for the spherical part}}\mathrm{d}{\mathscr {H}}^2(u)\right. \\&\qquad \quad \left. {\int _{S_2}\tau _{{{\mathscr {N}}}_C}(y,v)\mathrm{d}{\mathscr {H}}^2(v)}\right\rangle \\&\qquad = k_p(x,y)\left\langle \int _{S_1}u \mathrm{d}{\mathscr {H}}^2(u),\int _{S_2}v \mathrm{d}{\mathscr {H}}^2(v)\right\rangle \end{aligned} \end{aligned}$$

If we focus on portion of sphere, one can show that if

$$\begin{aligned} S_1 = \left\{ \begin{pmatrix} s\theta _u c\varphi _u\\ s\theta _u s\varphi _u\\ c\theta _u \end{pmatrix} \bigg | \theta _u \in [0,\pi ], \, \varphi _u \in [0,\varphi _0]\right\} \end{aligned}$$

is a portion of sphere, then $\int _{S_1} u \mathrm{d}{\mathscr {H}}^2(u) = \pi \sin \big (\varphi _0/2\big ) \begin{pmatrix} \cos (\varphi _0/2)\\ \sin (\varphi _0/2) \\ 0 \end{pmatrix}$. Note that we retrieve the half sphere with $\varphi _0 = \pi$ and the total sphere with $\varphi _0 = 2\pi$. Now, taking into account the orientation, we have to compute for the spherical part:

$$\begin{aligned}&\left\langle {x \times \Big ([s.] -\sum [{\hbox {h.s}}] + \sum [{\hbox {p.s}}]\Big )}y\right. \\&\quad \left. \times \Big ([s.] - \sum [{\hbox {h.s}}] + \sum [{\hbox {p.s}}] \Big )\right\rangle _{W'} \end{aligned}$$

where s stands for sphere, ${\hbox {h.s}}$ for half sphere and ${\hbox {p.s}}$ for portion of spheres. In fact, one can show that for a given vertex, summing the contributions of the sphere, the half spheres (associated with the edges), and the portion of spheres (associated with the triangles), then the spherical part vanishes for a vertex that is not in the border. Moreover, we have the following equality:

$$\begin{aligned} \left\langle N(C)^\mathrm{sph},N(S)^\mathrm{sph}\right\rangle _{W'} = \left\langle N(\partial C),N(\partial S)\right\rangle _{W'} \end{aligned}$$

and thus, the spherical part is exactly the scalar product of the curves associated with the border, scalar product that has been computed in [43].

Fig. 20

Decomposition of the normal bundle for two triangles with a common edge. In this figure, the two normal bundles of the open triangles appear. Then, we add (only once) the normal bundle of the open edge (the red cylinder and the two green half spheres). Then, we add (only once) the normal bundle of the vertices of the edge (the two green spheres) (Color figure online)

Now let us focus on the cylindrical part. Using expression (8), with $k_n = 1$, one can see that the scalar product involving a full cylinder is null, and thus, only the half cylinders remains. Consider thus the scalar product between two half cylinders. If we denote $HCyl_1 = [a,b] \times S_{b-a}^{\perp }$, $HCyl_2 = [c,d] \times S_{d-c}^{\perp }$ two half cylinders (where $S_{b-a,\alpha }^{\perp +} =\Big \{ u \in {\mathbb {S}}^2 | \left\langle u,b-a\right\rangle = 0, \left\langle u,\alpha \right\rangle \ge 0 \Big \}$ is a half circle), we compute the scalar product in $W'$ between these two half cylinders. With the approximations of (8):

$$\begin{aligned}&\left\langle HCyl_1,HCyl _2\right\rangle _{W'}\simeq k_p\Big (\frac{a+b}{2},\frac{c+d}{2}\Big )\left\langle b-a,d-c\right\rangle \\&\qquad \times \int _{S_{b-a,\alpha }^{\perp ,+}}\int _{S_{d-c,\beta }^\perp ,+}\left\langle \frac{b-a}{|b-a|} \times u ,\frac{d-c}{|d-c|} \times v\right\rangle \\&\qquad \mathrm{d}{\mathscr {H}}^1(u)\mathrm{d}{\mathscr {H}}^1(v)\\&\quad \simeq k_p\Big (\frac{a+b}{2},\frac{c+d}{2}\Big )\left\langle b-a,d-c\right\rangle \\&\quad \quad \times \left\langle \frac{b-a}{|b-a|} \times \int _{S_{b-a,\alpha }^{\perp ,+}} u \mathrm{d}{\mathscr {H}}^1(u) \frac{d-c}{|d-c|}\right. \\&\left. \quad \times \int _{S_{d-c,\beta }^\perp ,+} v \mathrm{d}{\mathscr {H}}^1(v) \right\rangle \\&\quad \simeq \frac{\pi ^2}{4} k_p\Big (\frac{a+b}{2},\frac{c+d}{2}\Big )\left\langle {b-a}{d-c}\right\rangle \\&\quad \left\langle \frac{b-a}{|b-a|} \times \alpha ,\frac{d-c}{|d-c|} \times \beta \right\rangle \end{aligned}$$

In a triangle T, [a, b] corresponds to an edge and $\alpha$ corresponds to a unitary vector orthogonal to [a, b], in the plane defined by the triangle and oriented in the interior of the triangle. Finally, if we consider two triangulations ${{\mathscr {T}}}$ and ${{\mathscr {T}}}'$, with a decomposition of the unit normal bundle as in Fig. 20, we have

(21)

where $n_{T_i,f_i}$ is the normal vector of the triangle $T_i$ such that $n_{T_i, f_i} \times f_i$ is oriented inward for the triangle T.

For the boundary, $\left\langle N(\partial {{\mathscr {T}}}),N(\partial {{\mathscr {T}}}')\right\rangle _{W'}$, we can show similarly that:

$$\begin{aligned}&\left\langle N(\partial {{\mathscr {T}}}),N(\partial {{\mathscr {T}}}')\right\rangle _{W'}\nonumber \\&\quad {=} \frac{\pi ^2}{4}\sum _{x_k \in \partial {{\mathscr {T}}}} \sum _{y_l \in \partial {{\mathscr {T}}}'} k_p(x_k,y_l)\left\langle A_k,B_l\right\rangle \nonumber \\ \end{aligned}$$

(22)

where $A_k = \sum _{i}f_i^k/|f_i^k|$ is the sum of the normalized edges of the boundary with $x_k$ as vertex, and oriented outward from $x_k$.

Discrete Scalar Product with Linear Normal Kernel

We can show that only the planar and the spherical parts are involved in this scalar product. For the planar part, we have already seen that

$$\begin{aligned}&\left\langle N(T)^\mathrm{pln},N(T')^\mathrm{pln}\right\rangle _{W'}\\&\quad \simeq 4 \sum _{i = 1}^N\sum _{j = 1}^M|T_i| |T_j'|k_p(c_i,c'_{j}) \left\langle n_{T_i}, n_{T'_j}\right\rangle ^2. \end{aligned}$$

In this appendix, we want to compute explicitly the spherical scalar product for normal cycles, with the linear normal kernel. The generic expression involved is the following integral:

$$\begin{aligned} \int _{S_1}\int _{S_2}\left\langle u,v\right\rangle ^2 \mathrm{d}u\mathrm{d}v \end{aligned}$$

where $S_1$ and $S_2$ are two portions of spheres (that may be the whole sphere, half sphere) with no assumption on the relative position of one sphere compared to the other.

To compute this expression, without loss of generality, we can suppose that $S_1$ is parametrized as follows:

$$\begin{aligned} S_1 = \left\{ \begin{pmatrix} \sin \theta _u\cos \varphi _u \\ \sin \theta _u\sin \varphi _u \\ \cos \theta _u \end{pmatrix} \, \Big \vert \, \theta _u \in [0, \pi ], \varphi _u \in [0,\varphi _1] \right\} \end{aligned}$$

where $\varphi _1$ is the aperture angle of the portion of sphere ($\varphi _1 = 2\pi$ for a whole sphere, and $\pi$ for a half sphere). Suppose that $v = \begin{pmatrix} v_1 \\ v_2 \\ v_3 \end{pmatrix}$, then:

$$\begin{aligned}&\int _{S_1}\left\langle u,v\right\rangle ^2 \mathrm{d}u = \int _0^\pi \int _0^{\varphi _1}\left\langle \begin{pmatrix} \sin \theta _u\cos \varphi _u \\ \sin \theta _u\sin \varphi _u \\ \cos \theta _u \end{pmatrix},\begin{pmatrix} v_1 \\ v_2 \\ v_3 \end{pmatrix}\right\rangle ^2\\&\quad \sin \theta _u\mathrm{d}\theta _u\mathrm{d}\varphi _u \end{aligned}$$

which can be made explicit using the fact that

$$\begin{aligned}&\int _0^\pi \sin ^3\theta _u \mathrm{d}\theta _u = \frac{4}{3}, \, \int _0^\pi \sin ^2\theta _u\cos \theta _u\mathrm{d}\theta _u = 0,\\&\int _0^\pi \sin \theta _u \cos ^2\theta _u \mathrm{d}\theta _u =\frac{2}{3}. \end{aligned}$$

Integrating first with respect to $\theta _u$ and then to $\varphi _u$, we end up with:

$$\begin{aligned} \int _{S_1}\left\langle u,v\right\rangle ^2 \mathrm{d}u= & {} \frac{2}{3}\varphi _1 + \frac{1}{3}(v_1^2 - v_2^2)\sin (2\varphi _1)\\&+\, \frac{4}{3}v_1v_2\sin ^2(\varphi _1). \end{aligned}$$

The quantity of interest is now:

$$\begin{aligned} \begin{aligned} \int _{S_1}\int _{S_2}\left\langle u,v\right\rangle ^2 \mathrm{d}u\mathrm{d}v&= \int _{S_2} \Big [ \frac{2}{3}\varphi _1 + \frac{1}{3}(v_1^2 - v_2^2)\sin (2\varphi _1) \\&\quad +\, \frac{4}{3}v_1v_2\sin ^2(\varphi _1) \Big ]\mathrm{d}v \end{aligned} \end{aligned}$$

The main limitation is that we do not have an obvious parameterization of $S_2$ since there is no assumption on the relative disposition of $S_1$ and $S_2$. Suppose now that R is the rotation which brings $S_2$ to $S_2'$, where

$$\begin{aligned} S_2' = \left\{ \begin{pmatrix} \sin \theta _v\cos \varphi _v \\ \sin \theta _v\sin \varphi _v \\ \cos \theta _v \end{pmatrix} \, \Big \vert \, \theta _v \in [0, \pi ], \varphi _v \in [0,\varphi _2] \right\} \end{aligned}$$

We have:

$$\begin{aligned} \begin{aligned}&\int _{S_2}(v_1^2 - v_2^2)\mathrm{d}v = \int _{S_2} \left( \left\langle v,e_1\right\rangle ^2 - \left\langle v,e_2\right\rangle ^2 \right) \mathrm{d}v\\&= \int _{R S_2} \left( \left\langle R^{-1}v,e_1\right\rangle ^2 - \left\langle R^{-1}v,e_2\right\rangle ^2 \right) \mathrm{d}v\\&=\int _0^\pi \int _0^{\varphi _2} \left( \left\langle \begin{pmatrix} \sin \theta _v\cos \varphi _v \\ \sin \theta _v\sin \varphi _v \\ \cos \theta _v \end{pmatrix},Re_1\right\rangle ^2 - \left\langle \begin{pmatrix} \sin \theta _v\cos \varphi _v \\ \sin \theta _v\sin \varphi _v \\ \cos \theta _v \end{pmatrix},Re_2\right\rangle ^2 \right) \\&\quad \times \sin \theta _v\mathrm{d}\theta _v\mathrm{d}\varphi _v\\ \end{aligned} \end{aligned}$$

This computation is similar to the one of $\int _{S_1}\left\langle u,v\right\rangle ^2 \mathrm{d}u$ and we obtain:

$$\begin{aligned} \int _{S_2}(v_1^2 - v_2^2)\mathrm{d}v= & {} \frac{1}{3} \left[ r_{11}^2 - r_{12}^2 + r_{22}^2 - r_{21}^2 \right] \sin (2\varphi _2) \\&+\,\frac{4}{3} \left[ r_{11}r_{21} - r_{12}r_{22} \right] \sin ^2(\varphi _2) \end{aligned}$$

where $R = (r_{ij})_{1 \le i,j \le 3}$. The same reasoning for the term $v_1v_2$ leads to

$$\begin{aligned} \begin{aligned} \int _{S_2}v_1v_2 \mathrm{d}v&= \int _0^\pi \int _0^{\varphi _2}\left\langle v,Re_1\right\rangle \left\langle v,Re_2\right\rangle \mathrm{d}v. \\&= \frac{1}{3} \left[ r_{11}r_{12}-r_{21}r_{22} \right] \sin (2\varphi _2) \\&\quad + \frac{2}{3} \left[ r_{11}r_{22}+r_{12}r_{21} \right] \sin ^2(\varphi _2) \end{aligned} \end{aligned}$$

Finally, if we combine all the terms, we obtain:

$$\begin{aligned} \begin{aligned}&\int _{S_1}\int _{S_2}\left\langle u,v\right\rangle ^2\mathrm{d}u\mathrm{d}v \\&\quad = \frac{4}{3}\varphi _1\varphi _2 + \frac{1}{9} \left[ r_{11}^2 - r_{12}^2 + r_{22}^2 - r_{21}^2 \right] \sin (2\varphi _1) \sin (2\varphi _2) \\&\qquad +\,\frac{4}{9}\left[ r_{11}r_{21} + r_{12}r_{22} \right] \sin (2\varphi _1)\sin ^2(\varphi _2) \\&\qquad + \,\frac{4}{9} \left[ r_{11}r_{12}-r_{21}r_{22} \right] \sin ^2(\varphi _1)\sin (2\varphi _2) \\&\qquad +\,\frac{8}{9} \left[ r_{11}r_{22}+r_{12}r_{21} \right] \sin ^2(\varphi _1)\sin ^2(\varphi _2) \end{aligned} \end{aligned}$$

(23)

with

$$\begin{aligned} r_{11}= & {} \left\langle e_1,f_1\right\rangle \\ r_{21}= & {} \left\langle e_1,f_2\right\rangle = \left\langle e_1,\frac{1}{\sin \varphi _2}(f_{\varphi _2}-\cos \varphi _2 f_1)\right\rangle \\ r_{12}= & {} \left\langle e_2,f_1\right\rangle = \left\langle \frac{1}{\sin \varphi _1}(e_{\varphi _1}-\cos \varphi _1 e_1),f_1\right\rangle \\ r_{22}= & {} \left\langle e_2,f_2\right\rangle = \left\langle {\frac{1}{\sin \varphi _1}(e_{\varphi _1}-\cos \varphi _1 e_1)}\right. \\&\left. {\frac{1}{\sin \varphi _2}(f_{\varphi _2}-\cos \varphi _2 f_1)}\right\rangle \\ \end{aligned}$$

Note 8

Expression (23) simplifies greatly when one of the involved sphere is a half sphere or a sphere. Indeed, all the terms with a sinus vanish, and it remains: $\int _{S_1}\int _{S_2}\left\langle u,v\right\rangle ^2\mathrm{d}u\mathrm{d}v = \frac{4}{3}\varphi _1\varphi _2$.

Note 9

Note that this computation is useful to compute the scalar product between portions of sphere in the triangulation. For an implementation, remember that the portion of sphere is not defined by the edges of the triangles, say $e_1$, $e_2$, but by the orthogonals: $-e_2^{\perp }, e_1^\perp$.

Now, if we are given two triangulated meshes, we will express one part of the spherical scalar product: at two vertex x and y, as we have done for the constant normal kernel, we need to compute

$$\begin{aligned} \left\langle x \times \Big ([s.] -\sum [{\hbox {h.s}}] + \sum [{\hbox {p.s}}]\Big ) \right. \\ \left. y \times \Big ([s.] - \sum [{\hbox {h.s}}] + \sum [{\hbox {p.s}}] \Big )\right\rangle _{W'} \end{aligned}$$

If we use the previous expressions of $\int _{S_1}\int _{S_2}\left\langle u,v\right\rangle ^2 \mathrm{d}u\mathrm{d}v$, ad gathering all the terms $4/3\varphi _1\varphi _2$, we end up with:

$$\begin{aligned} \begin{aligned}&\left\langle N({{\mathscr {T}}})^\mathrm{sph},N({{\mathscr {T}}}')^\mathrm{sph}\right\rangle _{W'} \\&\quad = \frac{4}{3}\sum _{k = 1}^{N_v}\sum _{l = 1}^{M_v}k_p(x_k,y_l)G_{{\mathscr {T}}}(x_k)G_{{{\mathscr {T}}}'}(y_l)\\&\quad \quad + \text {second-order spherical terms} \end{aligned} \end{aligned}$$

where

$$\begin{aligned} \left\{ \begin{aligned} G_{{\mathscr {T}}}(x_k)&= \Big [\pi (2-n_{x_k}+N_{x_k}) - \sum _{i = 1}^{N_{x_k}}\varphi _{i,x_k}\Big ]\\ G_{{{\mathscr {T}}}'}(y_l)&= \Big [\pi (2-m_{y_l}+N_{y_l}) -\sum _{j = 1}^{M_{y_k}}\varphi _{j,y_l}\Big ] \end{aligned} \right. \end{aligned}$$

and where the second-order spherical terms appear when considering the crossed terms

$$\begin{aligned} \begin{aligned}&\frac{1}{9} \left[ r_{11}^2 - r_{12}^2 + r_{22}^2 - r_{21}^2 \right] \sin (2\varphi _1) \sin (2\varphi _2) \\&+\frac{4}{9}\left[ r_{11}r_{21} + r_{12}r_{22} \right] \sin (2\varphi _1)\sin ^2(\varphi _2) \\&+ \frac{4}{9} \left[ r_{11}r_{12}-r_{21}r_{22} \right] \sin ^2(\varphi _1)\sin (2\varphi _2) \\&+\frac{8}{9} \left[ r_{11}r_{22}+r_{12}r_{21} \right] \sin ^2(\varphi _1)\sin ^2(\varphi _2) \end{aligned} \end{aligned}$$

we do not explicit these second-order terms since we still not have a satisfying implementation of this metric.