Computation of Circular Area and Spherical Volume Invariants via Boundary Integrals

TL;DR

This paper presents a method to compute circular area and spherical volume invariants of surfaces using boundary integrals, enabling efficient analysis of 3D shapes without ambient space discretization.

Contribution

It introduces a boundary integral approach for computing surface invariants, extending to higher dimensions and applicable to triangulated meshes, with potential applications in anthropology.

Findings

Provides analytical formulas for invariants on triangulated meshes

Enables computation without ambient space discretization

Applicable to feature detection on 3D surfaces

Abstract

We show how to compute the circular area invariant of planar curves, and the spherical volume invariant of surfaces, in terms of line and surface integrals, respectively. We use the Divergence Theorem to express the area and volume integrals as line and surface integrals, respectively, against particular kernels; our results also extend to higher dimensional hypersurfaces. The resulting surface integrals are computable analytically on a triangulated mesh. This gives a simple computational algorithm for computing the spherical volume invariant for triangulated surfaces that does not involve discretizing the ambient space. We discuss potential applications to feature detection on broken bone fragments of interest in anthropology.

Tables

Table 1. Table 1. Wall-clock times for computing the spherical volume invariant on the Stanford dragon [ 14 ] with ε = 1 𝜀 1 \varepsilon=1 for boundary triangle refinement. Computations were performed on a standard laptop computer using a single 3.2 3.2 3.2 GHz core, and CPU times were found to be very similar with and without boundary triangle refinement. For reference, on the finest mesh (181,444 triangles) with r = 5 𝑟 5 r=5 , each ball B ​ ( p , r ) 𝐵 𝑝 𝑟 B(p,r) contains on average 7,384.3 triangles.

Mesh size	Radius
(# triangles/#vertices)	$r=0.5$	$r=1$	$r=2$	$r=3$	$r=4$	$r=5$
45,360/22,678	0.19s	0.69s	2.5s	6.1s	10.3s	16.8s
90,722/45,359	0.67s	2.1s	8.9s	26.2s	40.7s	66.7s
181,444/90,720	2.0s	7.8s	32.8s	83.3s	151.4s	268.4s

Equations

$$A_{C,r}(p)=A(D_r(p)\cap \Omega ).$$

$$A_{C,r}(p)=\frac{\pi r^2}{2}-\frac{1}{3}\kappa (p)r^3+\text{O}(r^4)\text{ as }r\to 0.$$

$${\tilde{A}}_{C,r}=\frac{1}{L}{\oint }_C{A}_{C,r}(p(s))ds,$$

$$V_{S,r}(p)=V(\Omega \cap B_r(p)),$$

$$V_{S,r}(p)=\frac{2}{3}\pi r^3-\frac{1}{4}\pi H(p)r^4+\text{O}(r^5)\text{ as }r\to 0,$$

$$\mathbf{V}(x)=\frac{1}{2}(x-p)=\frac{1}{2}(x_1-p_1,x_2-p_2)$$

$$\text{hfil}{A}_{C,r}(p)={\iint }_{\Omega \cap D_r(p)}dxdy$$

$$\mathbf{V}\cdot \nu =\frac{1}{2}{r}^2\mathrm{on}\partial D_r(p).$$

$${\oint }_{\Omega \cap \partial D_r(p)}\mathbf{V}\cdot \nu ds={\int }_{{\theta }_1}^{{\theta }_2}\frac{r^2}{2}({\sin }^2\theta +{\cos }^2\theta )d\theta =\frac{r^2}{2}({\theta }_2-{\theta }_1).$$

$$A_{C,r}(p)={\oint }_{C\cap D_r(p)}\mathbf{V}\cdot \nu ds+\frac{r^2}{2}({\theta }_2-{\theta }_1).$$

$$A_{C,r}(p)={\oint }_{C\cap D_r(p)}\mathbf{V}\cdot \nu ds+\frac{r^2}{2}\sum _{i=1}^k({\theta }_{2i}-{\theta }_{2i-1}).$$

$$V_{S,r}:=V_{S,r}(0)={\int }_{\Omega \cap B_r}dx.$$

$$\mathbf{V}(x)=\frac{1}{n}x,and\; note\; that\text{div}\mathbf{V}=1.$$

$$V_{S,r}$$

$$\mathbf{W}\cdot \nu +\frac{r}{n}=0\mathrm{on}\partial \Omega \cap B_r.$$

$$\{\begin{array}{rlrl}\Delta u& =0& & \text{in }{B}_r,\\ \frac{\partial u}{\partial \nu }+\frac{r}{n}& =0& & \text{on }\partial \Omega \cap B_r.\end{array}$$

$${\int }_{\partial B_R}\frac{\partial u}{\partial \nu }dS=0$$

$$\{\begin{array}{rlrl}\Delta u& =0& & \text{in }{B}_r\setminus \{0\},\\ \frac{\partial u}{\partial \nu }+\frac{r}{n}& =0& & \text{on }\partial B_r,\end{array}$$

$$u(x)={\alpha }_n{r}^n\Phi (x),$$

$$\Phi (x)=\{\begin{array}{ll}-\frac{1}{2\pi }\log |x|,& \text{if }n=2\\ \frac{1}{n(n-2){\alpha }_n|x{|}^{n-2}},& \text{if }n\ge 3\end{array}$$

$$\nabla u(x)=-\frac{r^n}{n}\frac{x}{|x{|}^n}\mathrm{for}x\ne 0.$$

$${\int }_{S\cap (B_r\setminus B_{\epsilon })}(\mathbf{V}+\nabla u)\cdot \nu dS$$

$$={\int }_{\Omega \cap (B_r\setminus B_{\epsilon })}\text{div}(\mathbf{V}+\nabla u)dx+{\int }_{\Omega \cap \partial B_{\epsilon }}\left(\frac{\epsilon }{n}-\frac{r^n}{n{\epsilon }^{n-1}}\right)dS$$

$$={\int }_{\Omega \cap (B_r\setminus B_{\epsilon })}dx+\left(\frac{\epsilon }{n}-\frac{r^n}{n{\epsilon }^{n-1}}\right){\mathcal{H}}^{n-1}(\Omega \cap \partial B_{\epsilon })$$

$$=V_{S,r}-V_{S,\epsilon }+\left(\frac{\epsilon }{n}-\frac{r^n}{n{\epsilon }^{n-1}}\right){\mathcal{H}}^{n-1}(\Omega \cap \partial B_{\epsilon }),$$

$$V_{S,r}=V_{S,\epsilon }+\frac{1}{n}{\int }_{S\cap (B_r\setminus B_{\epsilon })}\left(1-\frac{r^n}{|x{|}^n}\right)(x\cdot \nu )dS+{\alpha }_n(r^n-{\epsilon }^n)\frac{{\mathcal{H}}^{n-1}(\Omega \cap \partial B_{\epsilon })}{{\mathcal{H}}^{n-1}(\partial B_{\epsilon })}.$$

$$\Gamma (p):=\underset{\epsilon \to 0^{+}}{\lim }\frac{{\mathcal{H}}^{n-1}(\Omega \cap \partial B_{\epsilon }(p))}{{\mathcal{H}}^{n-1}(\partial B_{\epsilon }(p))}$$

$$V_{S,r}(p)=\frac{1}{n}{\int }_{S\cap B_r(p)}\left(1-\frac{r^n}{|x-p{|}^n}\right)(x-p)\cdot \nu dS+{\alpha }_n{r}^n\Gamma (p).$$

$${\int }_{S\cap B_r(p)}\left(1-\frac{r^n}{|x-p{|}^n}\right)(x-p)\cdot \nu dS:=\underset{\epsilon \to 0^{+}}{\lim }{\int }_{S\cap (B_r(p)\setminus B_{\epsilon }(p))}\left(1-\frac{r^n}{|x-p{|}^n}\right)(x-p)\cdot \nu dS,$$

$$(x-p)\cdot \nu =\text{O}(|x-p{|}^{1+\alpha })\text{ as }x\to p$$

$$V_{S,r}(p)=\frac{1}{n}{\int }_{S\cap B_r(p)}\left(1-\frac{r^n}{|x-p{|}^n}\right)(x-p)\cdot \nu dS+\frac{1}{2}{\alpha }_n{r}^n.$$

$$(x-p)\cdot \nu =0$$

Full text

Funding: The authors are grateful to the National Science Foundation for support under grant NSF-DMS:1816917, the University of St. Thomas Center for Applied Math, and a University of Minnesota Grant in Aid award.   Source Code: https://github.com/jwcalder/Spherical-Volume-Invariant

1. Introduction

The aim of this paper is to facilitate the computation of certain integral invariants that have been proposed for applications in digital image processing, namely, the circular area and spherical volume invariants, as defined below. We show that both can be efficiently evaluated by reducing them to boundary integrals — line or surface integrals, respectively, — plus an additional term that depends only on the local surface geometry, thus enabling them to be computed directly from the curve or surface image data.

More specifically, given a Jordan plane curve $C\subset{\mathbb{R}}\mkern 1.0mu^{2}$ with interior $\Omega={\rm int}\>C$, at each point $p$ in the curve $C$, the value of the (local) circular area invariant of radius $r>0$ at $p$ is defined as the area (Lebesgue measure) of the region given by the intersection of the interior of the curve with a disk of radius $r$ centered at the point $p$, denoted $D_{r}(p)$:

[TABLE]

The circular area is clearly invariant under Euclidean motions of the curve, of course assuming one relates the base points $p$ accordingly. The ability of the local circular area invariant to uniquely characterize the curve up to Euclidean motion is discussed in detail in [9]. See Figure 1a for an illustration. For sufficiently smooth curves, e.g. $\textrm{C}^{3}$, the circular area invariant is related to the curvature $\kappa(p)$ at the point $p\in C$ by the asymptotic expansion [9]

[TABLE]

A global invariant can be obtained by averaging over the curve:

[TABLE]

where length $l(C)$ denotes the length of $C$.

Similarly, given a closed surface $S\subset{\mathbb{R}}\mkern 1.0mu^{3}$ bounding a domain $\Omega={\rm int}\>S$, we define the spherical volume invariant at each point $p\in S$ to be the volume of the solid region given by intersecting the interior of the surface with a sphere of radius $r>0$ centered at the point $p$:

[TABLE]

as illustrated in Figure 1b. Again, invariance under three-dimensional Euclidean motions is clear. For $\textrm{C}^{3}$ surfaces, the spherical volume invariant is related to the mean curvature of the surface via the expansion

[TABLE]

where $H(p)$ is the mean curvature of $S$ at $p$, and $S$ is a $\textrm{C}^{3}$ surface. The spherical volume invariant has proven useful for feature extraction [39, 32, 33], and a further analysis of the shape of the region ${\rm int}\>S\,\cap\,B_{r}(p)$ provides a robust estimation of the second fundamental form — see [33] and Subsection 3.2. Again, one can produce the corresponding global spherical volume invariant by integrating the local invariant over the entire surface.

These quantities clearly extend to the corresponding hyperspherical volume invariant of closed hypersurfaces in $S\subset{\mathbb{R}}^{n}$. Our main result is the general formula (3.13) that expresses this integral invariant in terms of a hypersurface integral over $S$. In the planar case, with $n=2$, our general formula reduces to a useful formula (2.3) or (3.15) for the circular area invariant $A_{C,r}(p)$ in terms of a suitable line integral over the curve $C$. For surfaces in $n=3$ dimensional space, it reduces to the key formula (3.16) for the spherical volume invariant $V_{S,r}(p)$ in terms of a surface integral over $S$. Our results apply to Lipschitz codimension $1$ submanifolds, which allows $S$ to be a triangulated mesh, as is often used to approximate surfaces in practice. These new formulas are simple and fast to implement on a triangulated mesh. In particular, our method does not require discretizing the ambient three dimensional space off the surface, as was done using octrees and the Fast Fourier Transform in [32]. Similar ideas can be used to evaluate other integral invariants, although a number of them are already expressed in terms of integrals of the type sought after here.

This paper was motivated by an ongoing project to analyze and reassemble broken bone fragments, a problem of significant interest in anthropology, paleontology, and surgery, building on earlier work of two of the authors on planar and surface jigsaw puzzle reassembly, [20, 18]. A recent undergraduate REU project, [38], has successfully applied the circular area integral invariant to planar jigsaw puzzle reassembly, following [20]. Indeed, one can easily envision modifying the circular area invariant in order to incorporate designs (writing, pictures, texture) that may appear on the puzzle pieces, potentially relying on some form of digital inpainting algorithm, [6, 7, 11, 12], to extend the design in the circular region on one side of the curve to the other, after which it could be compared to other potential matches, or, alternatively use of texture information to effect the reconstruction, as advocated in [35, 34].

Another potential application of these invariants is the detection of fracture edges, meaning ridges delineating the boundaries between original surfaces of the bone and break surfaces. Paleoanthropologists and zooarchaeologists study human biological and behavioral evolution and are interested in fracture edges because they provide valuable information about the agent of fragmentation [13, 17, 30], which may be, for example, humans, large carnivores, trampling, geological processes, or hydraulic action [26, 21]. Determining the agent of fragmentation is essential for reconstructing how archaeological sites were formed. Fracture edges can also be used to find bones that refit, which aids in the identification of taxa and skeletal elements of vertebrates found at sites [4]. We propose to detect fracture edges by thresholding the spherical volume invariant, and demonstrate by showing results of detecting fracture edges on bone fragments in Section 5.

The circular area and spherical volume invariants are particular cases of the general theory of integral invariants, [19, 23, 31], which have also been successfully applied to a variety of image processing problems. See [16] for applications of the moving frame method to their classification and signature construction under basic group actions, e.g., Euclidean and equi-affine geometries. Distance histograms underly the widely used methods of shape contexts, [5], and shape distributions, [28]. Histograms based on various geometric invariants (lengths, areas, etc.) play a fundamental role throughout a broad range of modern image processing algorithms, including shape representation and classification, [2, 37], image enhancement, [37, 36], the scale-invariant feature transform (SIFT) [22, 29], its affine-invariant counterpart (ASIFT), [40], and object-based query methods, [8].

1.1. Outline

In Section 2 we give a simple formula for the circular area invariant in terms of a line integral. In Section 3 we study the spherical volume invariant, and show how to use the Divergence Theorem to convert the volume integral into a surface integral, yielding a new formula for the invariant. Furthermore, in Subsection 3.2, we show how to extend our methods to estimate the principal curvatures of the surface by adapting the methods based on Principal Component Analysis (PCA) on local neighborhoods developed in [33]. Finally, in Section 5 we discuss numerical implementations and present the results of numerical experiments on real data. We use the Euclidean norm on ${\mathbb{R}}^{n}$ throughout, leaving the investigation of more general norms to a future project.

2. The Circular Area Invariant

As a warmup, we consider the local circular area invariant (1.1). We assume $C$ is the oriented boundary of an open bounded domain $\Omega\subset{\mathbb{R}}^{2}$ with Lipschitz boundary. Consider a point $p\in C$ with $p=(p_{1},p_{2})$. Consider the vector field

[TABLE]

and notice that $\text{div}\ {\bf V}=1$. By the Divergence Theorem, we can express the circular area invariant as

[TABLE]

where $\nu$ denotes the unit outward normal to the curve $C$ in the first line integral and to the circular boundary $\partial D_{r}(p)$ in the second. Let us parametrize the circular boundary of the disk $D_{r}(p)$ by $c(\theta)=p+r(\cos\theta,\sin\theta)$, so that

[TABLE]

Let $\theta_{1}<\theta_{2}$ be the angles at which the curve $C$ intersects the disk $D_{r}(p)$, assuming for the moment there are only 2 intersections and that $C$ lies inside the disk $D_{r}(p)$ for $\theta_{1}<\theta<\theta_{2}$. Now, the second term in (2.1) is

[TABLE]

Therefore, our formula for the circular area invariant is

[TABLE]

Notice this only involves integration along the curve $C$. The contour integral is a correction from the flat setting where $C$ is a line and $A_{C,r}(p)=\pi r^{2}/2$, since in this case $\theta_{2}-\theta_{1}=\pi$ and ${\bf V}\cdot\nu=0$ on $C\,\cap\,D_{r}(p)$.

It is straightforward to generalize (2.3) to more than two intersections of $C$ and $\partial D_{r}(p)$. If the intersections occur at angles $\theta_{1}<\theta_{2}<\cdots<\theta_{2k}$, and $C$ lies inside the disk111We ignore any intersection point where, nearby, $C$ remains on one side or the other of the boundary of the disk. $D_{r}(p)$ for $\theta_{2i-1}<\theta<\theta_{2i}$ for $i=1,\ldots,k$. Then we have

[TABLE]

3. The Spherical Volume Invariant

Having established a formula in the simple case of the two dimensional circular area invariant, we now turn to the spherical volume invariant (1.4). The argument used in Section 2 is not practical in three dimensions, since the integration over $\Omega\cap\partial D_{r}(p)$ in (2.1) becomes a surface integral, which defeats the point of reducing the calculation to an integral on the boundary surface.

We thus take a slightly different approach. Since the resulting formula will be applicable in all dimensions $n\geq 2$, we proceed in general. We assume our hypersurface $S\subset{\mathbb{R}}^{n}$ is the boundary of an open and bounded set $\Omega\subset{\mathbb{R}}^{n}$ with Lipschitz boundary. Without loss of generality, we take $p=0\in S$, and set $B_{r}=B_{r}(p)=B_{r}(0)$ to be the ball of radius $r$ centered at $p=0$. The hyperspherical invariant at $p=0$ is thus

[TABLE]

Define the vector field

[TABLE]

For any divergence free vector field ${\bf W}\colon{\mathbb{R}}^{n}\to{\mathbb{R}}^{n}$, whereby $\text{div}{\bf W}=0$, we can express $V_{S,r}$ via the Divergence Theorem as

[TABLE]

where $\nu$ denotes the outward normal to $S\,\cap\,B_{r}$ in the first term, and to $\partial\Omega\,\cap\,B_{r}$ in the second. The first term is an integral over the surface $S$, as we seek, while the second is an integral over $\partial B_{r}$, which is undesirable.

Now, the idea is to choose the vector field ${\bf W}$ so that the second term vanishes, yielding our formula. Noting that ${\bf V}\cdot\nu=r/n$ on $\partial B_{r}$, we see that ${\bf W}$ must satisfy

[TABLE]

We will construct ${\bf W}$ as ${\bf W}=\nabla u$ for a harmonic function $u$. Then (3.4) is equivalent to the Poisson problem

[TABLE]

If we look for a smooth solution of (3.5) then the compatibility condition

[TABLE]

must hold. This would require modifying the boundary condition away from $\partial\Omega\,\cap\,B_{r}$, which is impractical, since the set $\partial B_{r}\setminus\Omega$ could be arbitrarily small, and is dependent on the particular point $p$ chosen on the surface.

Instead of seeking to satisfy the compatibility condition (3.6), we relax the requirement that $u$ is smooth but continue to impose the boundary condition in (3.5). We allow $u$ to have a singularity at the origin, and thus consider the Poisson problem

[TABLE]

on the punctured ball. A solution to the latter boundary value problem is given by

[TABLE]

where $\alpha_{n}$ is the measure of the unit ball in ${\mathbb{R}}^{n}$, and

[TABLE]

is the fundamental solution of Laplace’s equation. Thus, we are effectively circumventing the compatibility condition (3.6) by placing a point source at the origin. Due to the singularity of $u$, the argument leading to (3.3) is no longer valid, and we need to proceed more cautiously.

First, we note that, for any $n$,

[TABLE]

Let $0<\varepsilon<r$. By the Divergence Theorem and the boundary condition in (3.7) we have

[TABLE]

where ${\mathcal{H}}^{n-1}$ denotes $(n-1)$–dimensional Hausdorff measure. Therefore

[TABLE]

All that is left is to send $\varepsilon\to 0$, and we state the consequence as a theorem.

Theorem 1**.**

Let $\Omega\subset{\mathbb{R}}^{n}$ be open and bounded with Lipschitz boundary $S:=\partial\Omega$. Let $p\in S$ and assume the limit

[TABLE]

exists. Then we have

[TABLE]

A few remarks are in order.

Remark 2**.**

Notice the integrand in (3.13) has a singularity at $x=p$. Since $S$ is only assumed to be Lipschitz, the singularity may not be integrable, and so we define the integral via its principal value

[TABLE]

which exists thanks to (3.11) and (3.12). If $S\in\textrm{C}^{1,\alpha}$ then we have

[TABLE]

and so the kernel singularity $|x-p|^{1-n+\alpha}$ is integrable on the $n-1$ dimensional surface.

Remark 3**.**

If the surface $S$ is differentiable at $p$ then $\Gamma(p)=\frac{1}{2}$, and thus

[TABLE]

Since a Lipschitz surface is differentiable almost everywhere, the formula (3.14) holds at almost every point of $S$.

Remark 4**.**

If the surface $S\subset{\mathbb{R}}\mkern 1.0mu^{3}$ is a triangulated mesh and $p\in S$ is a vertex of the mesh, then

[TABLE]

at all points $x$ in the vertex polygon associated to $p$ (i.e., the triangles adjacent to $p$), and where $\nu$ denotes the unit normal to the triangle containing $x$. Thus, the kernel is integrable on triangulated meshes. Moreover, $\Gamma(p)$ exists, and (3.13) holds, for every $p\in S$. In Subsection 3.1, we derive an explicit formula for $\Gamma(p)$ on a triangulated mesh in terms of the vertex polygon of $p$.

Remark 5**.**

The limit (3.12) defining $\Gamma(p)$ may fail to exist at a point of non-differentiability of a Lipschitz hypersurface $S$. Consider, for example, $n=2$ and take the curve $C$ to be the graph of the Lipschitz function

[TABLE]

Take the interior of $C$ to be the epigraph $\{x\in{\mathbb{R}}^{2}\,:\,f(x)>0\}$. Then the limit (3.12) does not exist at $p=0$, since along the sequence $x_{k}=e^{\frac{\pi}{2}-k\pi}$ we have $f(x_{k})=f(x_{k})=(-1)^{k}x_{k}$.

Remark 6**.**

Finally, let us note that in dimension $n=2$, the formula (3.13) reads

[TABLE]

In dimension $n=3$, it becomes

[TABLE]

3.1. An analytic expression for $\Gamma(p)$ on a triangulated mesh

We give here an analytic expression for $\Gamma(p)$, defined in (3.12), when $p$ is a vertex of a triangulated mesh surface in ${\mathbb{R}}^{3}$. Let us assume we have made a translation and rotation so that the vertex under consideration is $p=0$ and the unit outward normal vector at the origin is $\nu=(0,0,-1)$. Of course, there is no well-defined normal at the vertex $p=0$ itself, and so $\nu$ should chosen to be “close” to the nearby unit normals, in that it approximates the normal to the smooth surface represented by the mesh. For example, it could be the normalized average of the normals to the triangles in the vertex polygon; another possibility is that it is the normal to the least squares approximating plane to the vertices adjacent to $p$.

The computation of $\Gamma:=\Gamma(0)$ involves only the vertex triangles $T_{1},\dots,T_{k}$ that are adjacent to the vertex $p=0$. See Figure 2 for a depiction of these triangles and the area of the sphere we wish to compute. Since the outward normal at $p=0$ is $(0,0,-1)$, we will also assume that the outward unit normal vector $\nu^{i}=(\nu^{i}_{1},\nu^{i}_{2},\nu^{i}_{3})$ to each vertex triangle $T_{i}$ satisfies222We will exclude “bizarre” vertices where this assumption does not hold under any reasonable choice of the normal $\nu$ at $p$. $\nu^{i}_{3}<0$. Finally, in view of the definition (3.12) of $\Gamma$, we may extend the vertex triangles to $\infty$ in the radial direction, and compute

[TABLE]

where $\Omega$ is the region above the (extended) vertex triangles in the $x_{3}$-direction.

We work in spherical coordinates

[TABLE]

The edges $E_{1},\ldots,E_{k}$ of the vertex triangles $T_{1},\ldots,T_{k}$ containing the origin will be called vertex edges, and we let $(\theta_{i},\varphi_{i})$ be their corresponding spherical angles. We order the vertex edges and triangles so that

[TABLE]

and, for convenience, set $E_{k+1}=E_{1}$ with azimuthal angle $\theta_{k+1}=\theta_{1}+2\pi\geq 2\pi$. The vertex triangles are similarly ordered, so that $T_{i}$ has vertex edges $E_{i}$ and $E_{i+1}$. Observe that the vertex edges of $T_{i}$ are ordered so that $E_{i},E_{i+1},\nu_{i}$ form a left-handed frame, keeping in mind that, by the preceding assumption, the normal $\nu_{i}$ points downwards.

Each vertex triangle $T_{i}$ intersects the unit sphere $S_{1}=\partial B_{1}$ along a curve

[TABLE]

connecting $(\theta_{i},\varphi_{i})$ to $(\theta_{i+1},\varphi_{i+1})$. In terms of the intersection curves $C_{i}$ we can compute

[TABLE]

where

[TABLE]

We can find an explicit formula for $g_{i}(\theta)$. Indeed, the face $T_{i}$ is described by the plane

[TABLE]

Therefore, the intersection curve (3.19) satisfies

[TABLE]

and hence

[TABLE]

Noting the identity

[TABLE]

we have

[TABLE]

Since

[TABLE]

we can simplify the preceding formula to read

[TABLE]

We note that $\text{atan2}(y,x)$ is the two-argument $\arctan$ function, which gives the angle in radians between the positive $x$-axis and the ray from the origin to the point $(x,y)$, returning values in the interaval $[0,2\pi)$. Integrating $g_{i}$ yields

[TABLE]

This yields the following explicit formula:

[TABLE]

Unwrapping the definitions we have

[TABLE]

It follows that

[TABLE]

We also note that

[TABLE]

where $(x_{i},y_{i},z_{i})$ is any point along the edge $E_{i}$.

3.2. Principal component analysis on local neighborhoods

The spherical volume invariant of a surface in ${\mathbb{R}}^{3}$ is a robust estimator of its mean curvature, due to the asymptotic expansion given in (1.5). However, it gives no information about other differential geometric quantities of interest, such as the second fundamental form, the individual principal curvatures, the Gauss curvature, or the directions of principal curvature.

To capture additional geometric information, we follow [33] and analyze the shape of the region $\Omega\,\cap\,B_{r}(p)$. In particular, it is suggested in [33] to perform principal component analysis (PCA) on this region, that is, we compute the eigenvalues $\lambda_{1}(p)\geq\lambda_{2}(p)\geq\lambda_{3}(p)$ of the $3\times 3$ symmetric matrix333Here we take $x$ to be a column vector.

[TABLE]

where

[TABLE]

is the centroid of $\Omega\,\cap\,B_{r}(p)$, cf. (1.4). Assuming $S$ is sufficiently smooth, it was shown in [33] that the eigenvalues of $M_{S,r}(p)$ have the asymptotic expansions

[TABLE]

where $\kappa_{1}(p),\kappa_{2}(p)$ are the principal curvatures of the surface $S$ at the point $p\in S$, and, in the last formula, the $\textrm{O}(r^{6})$ term gives the mean curvature

[TABLE]

Moreover, the first two corresponding eigenvectors ${\bf v}_{1},{\bf v}_{2}$ are approximately tangent to the surface, and, assuming we are at a non-umbilic point, offer an $\textrm{O}(r/|\kappa_{1}-\kappa_{2}|)$ approximation of the directions of principal curvatures, while ${\bf v}_{3}$ is approximately normal to the surface and is an $\textrm{O}(r^{2})$ approximation of the unit normal. Thus, the matrix $M_{S,r}(p)$ provides a robust estimation of the second fundamental form of $S$ at a non-umbilic point $p$.

Let us now show how to compute the matrix $M_{S,r}(p)$ via surface integrals, as we did for the spherical volume invariant $V_{S,r}(p)$ in Theorem 1. While these results are mainly of interest in dimension $n=3$, we carry out the derivation for an arbitrary dimension $n$. Noting that

[TABLE]

it suffices to compute the first two moments

[TABLE]

in terms of which the $(i,j)$ entry of $M_{S,r}(p)$ is given by

[TABLE]

The computation of $m_{i}(p)$ and $c_{ij}(p)$ in terms of surface integrals is relatively straightforward, compared to the computation of $V_{S,r}$. In what follows, $e_{1},e_{2},\dots,e_{n}$ denote the standard basis vectors in ${\mathbb{R}}^{n}$, and $\delta_{ij}$ is the Kronecker delta.

Lemma 7**.**

Let us abbreviate $y=x-p$. Then, for any $1\leq i,j\leq n$, we have

[TABLE]

and

[TABLE]

Proof.

Without loss of generality, we may assume $p=0$. Then $y=x$ and we write $B_{r}=B_{r}(p)=B_{r}(0)$. We first prove (3.34). Define the vector field

[TABLE]

By the Divergence Theorem,

[TABLE]

On the spherical portion of the boundary $\partial\Omega\,\cap\,B_{r}$, we have $\nu=x/r$ and so

[TABLE]

since $|x|^{2}=r^{2}$ on $\partial B_{r}$. This completes the proof of (3.34).

We now prove (3.35). Define the vector field

[TABLE]

By the Divergence Theorem, we have

[TABLE]

On the portion of the boundary $x\in\partial\Omega\,\cap\,B_{r}$

[TABLE]

which completes the proof. ∎

4. Implementation

Let us next discuss how to compute the surface integrals from Theorem 1 and Lemma 7 on a surface given as a triangulated mesh, which is often the case in practice. The integrals we wish to compute all have the form

[TABLE]

for various choices of kernel function $f(x)$. We adopt the convention that $f(x)=0$ if $|x-p|>r$, and hence rewrite (4.1) as simply

[TABLE]

Let $T_{1},\dots,T_{M}$ denote the triangles in the triangulated surface $S$. Then we can write

[TABLE]

We show in Sections 4.1 and 4.2 that the triangular integrals appearing in the summation can be computed analytically for all of the kernels $f$ used in this paper. Let us note that on the right hand side of (4.3), we need only sum over triangles $T_{m}$ that have non-empty intersection with $B_{r}(p)$. However, it is computationally expensive to perform a range search to find all such triangles, especially for large meshes. In our implementation, we instead perform a depth first search on the triangle graph of the mesh, starting at any triangle adjacent to $p$, and terminating when all triangles in the connected component of $S\cap B_{r}(p)$ containing $p$ are found. While the depth first search has linear complexity and is very fast in practice, it will fail to find any additional connected components of $S\cap B_{r}(p)$ that do not contain $p$. On the other hand, this may be a desirable property of the algorithm, especially if one is primarily interested in the local geometry of the mesh.

4.1. Analytic integration over triangles

Let us show how all the integrals considered in this paper can be computed analytically over triangles $T_{m}\subset B_{r}(p)$. For simplicity, we take $p=0$, write $B_{r}=B_{r}(0)$, and consider a triangle $T$.

For the spherical volume invariant, for any triangle $T$ with $T\subset B_{r}$ the surface integral (3.13) from Theorem 1 requires us to compute

[TABLE]

Since $x\cdot\nu$ is constant over the triangle $T$, we have

[TABLE]

where $z$ is any point belonging to $T$, such as its centroid or one of its vertices, while $|T|$ denotes the surface area of $T$. The remaining integrand $|x|^{-3}$ is known as a hypersingular kernel, and arises, for instance, in the boundary element method for solving partial differential equations [3]. The integral of this hypersingular kernel over any planar triangle can be computed analytically [27] provided $p=0\not\in T$, which we may freely assume since $z\cdot\nu=0$ when $0\in T$. For convenience, we recall the analytic formula, which is rather tedious and derived in [27], in Appendix A.

For PCA on local neighborhoods, the integrals we need to compute from Lemma 7 correspond to

[TABLE]

Since $x\cdot\nu$ and $e_{i}\cdot\nu$ are constant over $T$, we just need to compute the quantities

[TABLE]

Let us denote the vertices of $T$ by $x,y,z\in{\mathbb{R}}^{3}$. The first integrand in (4.5) is linear, and so the integral can be computed analytically with the three point stencil

[TABLE]

For $b_{ij}$, we compute the integral in barycentric coordinates

[TABLE]

Computing

[TABLE]

and

[TABLE]

we have

[TABLE]

If we were to denote the vertices by $v_{i}=(v_{i}^{1},v_{i}^{2},v_{i}^{3})$, say, then (4.7) would have the simple form

[TABLE]

and similarly for (4.6).

4.2. Boundary triangles

For triangles $T$ that have a non-empty intersection with the boundary $\partial B_{r}$ of the ball $B_{r}$, the integral over $T$ cannot be computed analytically. To determine whether a triangle $T$ intersects $\partial B_{r}$, we compute

[TABLE]

and check whether $r_{1}\leq r\leq r_{2}$. To compute $r_{2}$, it is sufficient to check the vertices of the triangle, since $x\mapsto|x|$ is convex. The computation of $r_{1}$ is more tedious, since the minimum distance may occur interior to $T$. To compute $r_{1}$ we orthogonally project the origin $p=0$ onto the plane containing the triangle $T$, calling the projection $x_{P}$. If $x_{P}\in T$, then $r_{1}=|x_{P}|$. If $x_{P}\not\in T$, then we find the closest point $x_{T}\in T$ to the projection $x_{P}$, and therefore $r_{1}^{2}=|x_{T}|^{2}+|x_{T}-x_{P}|^{2}$ by the Pythagorean Theorem.

To compute the integral over such boundary triangles, we fix a maximum desired side length $\ell>0$ and recursively bisect the triangle along the line segment connecting the midpoint of its longest side with the opposing vertex. We stop the bisection procedure on a given subtriangle $T_{s}$ if $T_{s}\cap\partial B_{r}=\varnothing$, or the maximum side length of $T_{s}$, denoted $L(T_{s})$, falls below $\ell$. See Figure 3 for an illustration of the bisection process. We compute the integration over $T_{s}$ analytically if $T_{s}\subset B_{r}$, or with the approximation

[TABLE]

if $T_{s}\cap\partial B_{r}\not=\varnothing$, where $x,y,z$ are the vertices of $T_{s}$. We note the approximation error is bounded by

[TABLE]

where

[TABLE]

denotes the oscillation of the function $f$ over the triangle $T_{s}$. Now, let $\delta>0$ so that

[TABLE]

Note that $\delta$ is reduced throughout the bisection procedure, since the diameter of triangles that intersect $\partial B_{r}$ is decreasing. Since the triangles $T_{s}$ in (4.10) belong to $S\cap B_{r+\ell}\setminus B_{r-\ell}$, the error in computing (4.2) is bounded by

[TABLE]

where $|S\cap A|$ denotes the surface area of $S\cap A$. We assume that for the mesh $S$, there exists a constant $C>0$, independent of $r$ and $\ell$, such that

[TABLE]

Therefore, our approximation error is at most

[TABLE]

where $\delta$ is defined in (4.10). We note the volume growth assumption (4.12) is convenient, in that it leads to a simple form for the integration error (4.13). However, the analysis below can be easily carried out with other assumptions in place of (4.12), if needed. The volume growth assumption is true for smooth surfaces, with $r>0$ small, and hence for any triangulated mesh that well-approximates a smooth surface.

The application of (4.13) depends on the context. For the spherical volume invariant, we have

[TABLE]

where $a_{+}:=\max\{a,0\}$. For any triangle $T_{s}$ with maximum side length less than $\ell$ and satisfying $T_{s}\cap\partial B_{r}\not=\varnothing$, we have

[TABLE]

Thus, $\delta$ from (4.10) can be chosen as $\delta=r^{4}\ell/(r-\ell)^{4}$. Since the spherical volume invariant scales with $r^{3}$, it is reasonable to select an error tolerance $\varepsilon>0$ and ask that the integration error is bounded by $\varepsilon\,r^{3}$. Thus, invoking (4.13) we find that $\ell$ should be selected so that $\ell<r$ and

[TABLE]

Note that we are discarding the constant $2C\pi$ in (4.12), since we are only interested in how $\ell$ should scale with $r$ and $\varepsilon$. If $\varepsilon\ll 1$ so that $\ell\ll r$, this condition can be approximated by $\ell^{\,2}\leq\varepsilon\,r^{2}$. In particular, the triangle refinement is more important for small radii $r>0$, and for sufficiently large $r>0$, no refinement is needed.

For PCA on local neighborhoods, we have two integrals to compute. The first (3.34) corresponds to

[TABLE]

Since $f$ is not Lipschitz, the oscillation bound is at best

[TABLE]

provided $\ell\leq r$. Thus, $\delta$ from (4.10) can be chosen as $\delta=2r^{2}$. Inspecting (3.31), we see that it is reasonable to ask that the integration error is bounded by $\varepsilon\,r^{4}$, for an error tolerance parameter $\varepsilon>0$. Combining this with (4.13) the restriction on $\ell$ becomes $\ell\leq\min\{\varepsilon\,r,r\}$.

The second integral (3.35) required by PCA on local neighborhoods corresponds to

[TABLE]

As before, we bound the oscillation by

[TABLE]

provided $\ell\leq r$, and so $\delta=4\,r^{3}$. By (3.31) we see that it is natural to bound the integration error by $\varepsilon\,r^{5}$, yielding again the condition $\ell\leq\min\{\varepsilon\,r,r\}$.

5. Numerical experiments

We now present the results of numerical experiments using our method to compute the spherical volume invariant for triangulated surfaces arising from standard images, and for real experimental data arising from a project to classify and reassemble broken bone fragments in an archaeological context. For brevity, we will not discuss the much simpler case of curves and the circular area invariant. Our code is written in C and can be run from Matlab via the MEX interface, and from Python via an extension module. The code is available for download on GitHub:                             https://github.com/jwcalder/Spherical-Volume-Invariant

We first consider the standard test case of the Stanford dragon [14]. Figure 4 shows the spherical volume invariant for radii $r=1,2,5$ computed on the dragon. In Figure 4, and in all other plots below (unless otherwise specified), the colors indicate the values of the spherical volume invariant, with red indicating the lowest value and blue corresponding to the highest. For the dragon, and all other experiments, we used an error tolerance of $\varepsilon=1$ for bisecting boundary triangles. In the case of the dragon, the maximum triangle bisection depth was 8 and the maximum number of sub-triangles in any refinement was 57.

We mention that the original version of the Stanford dragon exhibits some non-manifold geometry. In particular, there are stray vertices not connected to triangles, and some edges are shared by more than 2 triangles. On such meshes, our method can produce unpredictable results, such as negative values for volumes, since Theorem 1 no longer holds. This can be easily remedied by cleaning the mesh with any standard mesh software package before running our code, or obtaining the mesh from a reliable algorithm, such as isosurfacing. For our experiment with the Stanford dragon reported in Figure 4, we obtained a version of the Stanford dragon in PWN (Points with Normals) format from [1], and converted to a clean triangulated mesh using the code provided in [1].

Our method is computationally efficient for large meshes. Table 1 shows the wall-clock times444“Wall-clock” time refers to the actual amount of time taken to perform the operation, as opposed to CPU time, which is often used to refer to how much time the processor spent on the job. for computing the spherical volume invariant with our method on the dragon for various radii. We also include results for lower resolution versions of the dragon for comparison. We can see the complexity of our method scales quadratically with the radius $r$, as expected. These CPU times are comparable to the FFT methods reported in [33]. We note FFT methods require coarsely discretizing the ambient space, resulting in larger numerical errors.

Let us next compute the spherical volume invariant on broken bone fragments that have been scanned and digitized for anthropological applications, as outlined in the introduction. Figure 5 shows the spherical volume invariant plotted over bone fragments at different radii, demonstrating how varying the radius allows one to change the scale of detected features. Here, the values of the spherical volume invariant are normalized with a power-law correction $v\mapsto v^{p}$, with $p=.5$ unless otherwise stated, to maximize contrast for visualization. We note that for larger radii in Figure 5, the spherical volume appears to be discontinuous at distance $r$ from a sharp fracture edge, which may be a desirable feature, depending on the application. This is due to our use of the connected component of $B(p,r)\cap S$ in computations, which fails to explore the opposite side of the fragment if $B(p,r)$ does not intersect the fracture edge. We also computed the principal curvatures via PCA on local neighborhoods. Figure 6 shows the Gauss curvature and Figures 7 and 8 show the two principal curvatures for some of the fragments. We did not include figures for mean curvature, since they are identical to those in Figure 5 for the spherical volume invariant, except with the colors reversed.

We can detect fracture edges by thresholding the spherical volume invariant. Edge points are taken as those with spherical volumes less than one standard deviation below the mean spherical volumes for the whole fragment. Figure 9 shows the results of fracture edge detection on several bone fragments and the Stanford Dragon. This simple approach gives a rough outline of most fracture edges. In future work, we plan to investigate automated algorithms for choosing the thresholds as well as the prospect of using the spherical volume invariant with more sophisticated edge detection methods, such as active contours [10] on surfaces, or graph-cut segmentation algorithms [24, 25, 15].

6. Conclusion

In this paper, we showed how to compute a class of integral invariants, including the circular area invariant and the spherical volume invariant, in terms of line and surface integrals over the bounding curve and surface. The method is computationally efficient to implement on a triangulated mesh, since it involves simply integrating a function over the mesh triangles, which, when the triangle lies inside the ball, can be done explicitly. In particular, it does not require discretizing the ambient three dimensional space. We showed how to numerically implement the integration accurately and efficiently, and presented the results of some numerical experiments with real data.

Acknowledgements: The bone fragments depicted are from an adult elk (Cervus canadensis) femur that was broken by adult male spotted hyena (Crocuta crocuta) named Scruffy who resides at the Milwaukee County Zoo. The femur was disarticulated and defleshed prior to being fed to the hyena. All fragments were scanned with a DAVID white light scanner that was made available by the Evolutionary Anthropology Labs at the University of Minnesota. The authors also gratefully acknowledge discussions with Martha Tappen, Jacob Elafandi, and Jacob Theis.

Appendix A Analytic formula for a hypersingular integral

Here, for the reader’s convenience, we recall the analytic formula from [27] for the hypersingular integral

[TABLE]

where $T$ is a planar triangle in ${\mathbb{R}}^{3}$, such that $0\not\in T$. (We note that [27] includes analytic formulas for several such triangular hypersingular integrals involving other negative integer powers of $|x|$.) Let $P$ denote the plane containing $T$ and $\nu$ the unit outward normal vector to $T$ and $P$. In what follows, we take the triangle $T$ to be an open subset of $P$, i.e., $T\cap\partial T=\mathchar 0\relax$.

Let $x^{*}\in{\mathbb{R}}^{3}$ denote the orthogonal projection of the origin onto the plane $P$. Let $x^{1},x^{2},x^{3}\in{\mathbb{R}}^{3}$ be the vertices of $T$, given with positive orientation; for convenience of notation we write $x^{4}=x^{1}$. Define

[TABLE]

where $\theta_{i}$ is the interior angle of $T$ at the vertex $x^{i}$.

Let $L^{i}$ denote the oriented edge of the triangle $T$ from $x^{i}$ to $x^{i+1}$. Associated with each edge $L^{i}$, we construct an orthonormal basis $(\mathbf{e}^{i}_{1},\mathbf{e}^{i}_{2})$ for the plane $P$ with origin $x^{*}$, $\mathbf{e}^{i}_{1}$ taken in the direction of the edge $L^{i}$, and $\mathbf{e}^{i}_{2}=\nu\times\mathbf{e}^{i}_{1}$ chosen so that $(\mathbf{e}^{i}_{1},\mathbf{e}^{i}_{2},\nu)$ is an orthonormal basis for ${\mathbb{R}}^{3}$. Let

[TABLE]

be the planar coordinates of the vertex $x^{j}$ in the basis $(\mathbf{e}^{i}_{1},\mathbf{e}^{i}_{2})$. By definition, $q^{1}_{1}=q^{1}_{2}$, $q^{2}_{2}=q^{3}_{2}$, and $q^{3}_{3}=q^{4}_{3}$, since the vertices $x^{j}$ and $x^{j+1}$ lie along the line spanned by $\mathbf{e}^{j}_{1}$. We denote the common values as

[TABLE]

Finally, set $\eta=x^{1}\cdot\nu$, noting that $\eta\neq 0$, since $0\not\in T$. We then define

[TABLE]

using the branch of $\arctan$ with values in $(-\pi/2,\pi/2)$. Finally, the hypersingular integral (A.1) is given by

[TABLE]