    # Computation of Circular Area and Spherical Volume Invariants via Boundary Integrals

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.

## Code Repositories

### Spherical-Volume-Invariant

Spherical Volume Invariant

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 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 with interior , at each point in the curve , the value of the (local) circular area invariant of radius at 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 centered at the point , denoted :

 AC,r(p)=A(Dr(p)∩Ω). (1.1)

The circular area is clearly invariant under Euclidean motions of the curve, of course assuming one relates the base points accordingly. The ability of the local circular area invariant to uniquely characterize the curve up to Euclidean motion is discussed in detail in . See Figure (a)a for an illustration. For sufficiently smooth curves, e.g. , the circular area invariant is related to the curvature at the point by the asymptotic expansion 

 AC,r(p)=πr22−13κ(p)r3+O(r4)   as r→0. (1.2)

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

 ˜AC,r=1L∮CAC,r(p(s))ds, (1.3)

where length denotes the length of .

Similarly, given a closed surface bounding a domain , we define the spherical volume invariant at each point to be the volume of the solid region given by intersecting the interior of the surface with a sphere of radius centered at the point :

 VS,r(p)=V(Ω∩Br(p)), (1.4)

as illustrated in Figure (b)b. Again, invariance under three-dimensional Euclidean motions is clear. For surfaces, the spherical volume invariant is related to the mean curvature of the surface via the expansion

 VS,r(p)=23πr3−14πH(p)r4+O(r5)   as r→0, (1.5)

where is the mean curvature of at , and is a

surface. The spherical volume invariant has proven useful for feature extraction

[39, 32, 33], and a further analysis of the shape of the region

provides a robust estimation of the second fundamental form — see

 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 . Our main result is the general formula (3.13) that expresses this integral invariant in terms of a hypersurface integral over . In the planar case, with , our general formula reduces to a useful formula (2.3) or (3.15) for the circular area invariant in terms of a suitable line integral over the curve . For surfaces in dimensional space, it reduces to the key formula (3.16) for the spherical volume invariant in terms of a surface integral over . Our results apply to Lipschitz codimension submanifolds, which allows

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

. 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, , has successfully applied the circular area integral invariant to planar jigsaw puzzle reassembly, following . 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 . 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  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, , and shape distributions, . 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), , and object-based query methods, .

### 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

. Finally, in Section 5 we discuss numerical implementations and present the results of numerical experiments on real data. We use the Euclidean norm on 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 is the oriented boundary of an open bounded domain with Lipschitz boundary. Consider a point with

. Consider the vector field

 V(x)=12(x−p)=12(x1−p1,x2−p2)

and notice that . By the Divergence Theorem, we can express the circular area invariant as

 =∬Ω∩Dr(p)(div V)dxdy
(2.1)

where denotes the unit outward normal to the curve in the first line integral and to the circular boundary in the second. Let us parametrize the circular boundary of the disk by , so that

 V⋅ν=12r2\rm on∂Dr(p).

Let be the angles at which the curve intersects the disk

, assuming for the moment there are only 2 intersections and that

lies inside the disk for . Now, the second term in (2.1) is

 ∮Ω∩∂Dr(p)V⋅νds=∫θ2θ1r22(sin2θ+cos2θ)dθ=r22(θ2−θ1). (2.2)

Therefore, our formula for the circular area invariant is

 AC,r(p)=∮C∩Dr(p)V⋅νds+r22(θ2−θ1). (2.3)

Notice this only involves integration along the curve . The contour integral is a correction from the flat setting where is a line and , since in this case and on .

It is straightforward to generalize (2.3) to more than two intersections of and . If the intersections occur at angles , and lies inside the disk111We ignore any intersection point where, nearby, remains on one side or the other of the boundary of the disk. for for . Then we have

 AC,r(p)=∮C∩Dr(p)V⋅νds+r22k∑i=1(θ2i−θ2i−1). (2.4)

## 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 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 , we proceed in general. We assume our hypersurface is the boundary of an open and bounded set with Lipschitz boundary. Without loss of generality, we take , and set to be the ball of radius centered at . The hyperspherical invariant at is thus

 VS,r:=VS,r(0)=∫Ω∩Brdx. (3.1)

Define the vector field

 V(x)=1nx,\rm and note thatdiv V=1. (3.2)

For any divergence free vector field , whereby , we can express via the Divergence Theorem as

 VS,r =∫Ω∩Brdiv(V+W)dx =∫S∩Br(V+W)⋅νdS+∫∂Ω∩Br(V+W)⋅νdS, (3.3)

where denotes the outward normal to in the first term, and to in the second. The first term is an integral over the surface , as we seek, while the second is an integral over , which is undesirable.

Now, the idea is to choose the vector field so that the second term vanishes, yielding our formula. Noting that on , we see that must satisfy

 W⋅ν+rn=0\rm on ∂Ω∩Br. (3.4)

We will construct as for a harmonic function . Then (3.4) is equivalent to the Poisson problem

 ⎧⎨⎩Δu=0in Br,∂u∂ν+rn=0on  ∂Ω∩Br. (3.5)

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

 ∫∂BR∂u∂νdS=0 (3.6)

must hold. This would require modifying the boundary condition away from , which is impractical, since the set could be arbitrarily small, and is dependent on the particular point chosen on the surface.

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

 ⎧⎨⎩Δu=0in Br∖{0},∂u∂ν+rn=0on ∂Br, (3.7)

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

 u(x)=αnrnΦ(x), (3.8)

where is the measure of the unit ball in , and

 Φ(x)=⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩−12πlog|x|,if n=21n(n−2)αn|x|n−2,if n≥3 (3.9)

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 , the argument leading to (3.3) is no longer valid, and we need to proceed more cautiously.

First, we note that, for any ,

 ∇u(x)=−rnnx|x|n\rm for x≠0. (3.10)

Let . By the Divergence Theorem and the boundary condition in (3.7) we have

 ∫S∩(Br∖Bε)(V+∇u)⋅νdS =∫∂(Ω∩(Br∖Bε))(V+∇u)⋅νdS+∫Ω∩∂Bε(V+∇u)⋅νdS =∫Ω∩(Br∖Bε)div(V+∇u)dx+∫Ω∩∂Bε(εn−rnnεn−1)dS =∫Ω∩(Br∖Bε)dx+(εn−rnnεn−1)Hn−1(Ω∩∂Bε) =VS,r−VS,ε+(εn−rnnεn−1)Hn−1(Ω∩∂Bε),

where denotes –dimensional Hausdorff measure. Therefore

 VS,r=VS,ε+1n∫S∩(Br∖Bε)(1−rn|x|n)(x⋅ν)dS+αn(rn−εn)Hn−1(Ω∩∂Bε)Hn−1(∂Bε). (3.11)

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

###### Theorem 1.

Let be open and bounded with Lipschitz boundary . Let and assume the limit

 Γ(p):=limε→0+Hn−1(Ω∩∂Bε(p))Hn−1(∂Bε(p)) (3.12)

exists. Then we have

 VS,r(p)=1n∫S∩Br(p)(1−rn|x−p|n)(x−p)⋅νdS+αnrnΓ(p). (3.13)

A few remarks are in order.

###### Remark 2.

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

 ∫S∩Br(p)(1−rn|x−p|n)(x−p)⋅νdS:=limε→0+∫S∩(Br(p)∖Bε(p))(1−rn|x−p|n)(x−p)⋅νdS,

which exists thanks to (3.11) and (3.12). If then we have

 (x−p)⋅ν=O(|x−p|1+α) as x→p

and so the kernel singularity is integrable on the dimensional surface.

###### Remark 3.

If the surface is differentiable at then , and thus

 VS,r(p)=1n∫S∩Br(p)(1−rn|x−p|n)(x−p)⋅νdS+12αnrn. (3.14)

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

###### Remark 4.

If the surface is a triangulated mesh and is a vertex of the mesh, then

 (x−p)⋅ν=0

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

###### Remark 5.

The limit (3.12) defining may fail to exist at a point of non-differentiability of a Lipschitz hypersurface . Consider, for example, and take the curve to be the graph of the Lipschitz function

 f(x)=|x|sin(log|x|).

Take the interior of to be the epigraph . Then the limit (3.12) does not exist at , since along the sequence we have .

###### Remark 6.

Finally, let us note that in dimension , the formula (3.13) reads

 AC,r(p)=12∮C∩Dr(p)(1−r2|x−p|2)(x−p)⋅νds+πr2Γ(p). (3.15)

In dimension , it becomes

 VS,r(p)=13∫S∩Br(p)(1−r3|x−p|3)(x−p)⋅νdS+43πr3Γ(p). (3.16)

### 3.1. An analytic expression for Γ(p) on a triangulated mesh

We give here an analytic expression for , defined in (3.12), when is a vertex of a triangulated mesh surface in . Let us assume we have made a translation and rotation so that the vertex under consideration is and the unit outward normal vector at the origin is . Of course, there is no well-defined normal at the vertex itself, and so 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 .

The computation of involves only the vertex triangles that are adjacent to the vertex . See Figure 2 for a depiction of these triangles and the area of the sphere we wish to compute. Since the outward normal at is , we will also assume that the outward unit normal vector to each vertex triangle satisfies222We will exclude “bizarre” vertices where this assumption does not hold under any reasonable choice of the normal at . . Finally, in view of the definition (3.12) of , we may extend the vertex triangles to in the radial direction, and compute

 Γ:=14π∫Ω∩∂B1dS, (3.17)

where is the region above the (extended) vertex triangles in the -direction.

We work in spherical coordinates

 x1=rsinφcosθ,x2=rsinφsinθ,x3=rcosφ. (3.18)

The edges of the vertex triangles containing the origin will be called vertex edges, and we let be their corresponding spherical angles. We order the vertex edges and triangles so that

 0≤θ1<θ2<θ3<⋯<θk<2π,

and, for convenience, set with azimuthal angle . The vertex triangles are similarly ordered, so that has vertex edges and . Observe that the vertex edges of are ordered so that form a left-handed frame, keeping in mind that, by the preceding assumption, the normal points downwards.

Each vertex triangle intersects the unit sphere along a curve

 Ci={φ=hi(θ)}=Ti∩S1 (3.19)

connecting to . In terms of the intersection curves we can compute

 Γ=14πk∑i=1∫θi+1θi∫hi(φ)0sinφdφdθ=12−14πk∑i=1∫θi+1θigi(θ)dθ, (3.20)

where

 gi(θ)=coshi(θ).

We can find an explicit formula for . Indeed, the face is described by the plane

 x3=aix1+bix2,\rm whereai=−νi1νi3\rm andbi=−νi2νi3. (3.21)

Therefore, the intersection curve (3.19) satisfies

 gi(θ)=coshi(θ)=(aicosθ+bisinθ)sinhi(θ),

and hence

 hi(θ)=cot−1(aicosθ+bisinθ). (3.22)

Noting the identity

 cos(cot−1x)=x√1+x2,

we have

 gi(θ)=aicosθ+bisinθ√1+(aicosθ+bisinθ)2.

Since

 aicosθ+bisinθ=cicos(θ−δi),\rm where% δi=atan2(bi,ai),ci=√a2i+b2i,

we can simplify the preceding formula to read

 gi(θ)=cicos(θ−δi)√1+c2icos2(θ−δi). (3.23)

We note that is the two-argument function, which gives the angle in radians between the positive -axis and the ray from the origin to the point , returning values in the interaval . Integrating yields

 ∫gi(θ)dθ=arcsin(disin(θ−δi))+Constant,\rm where di=ci√1+c2i. (3.24)

This yields the following explicit formula:

 Γ=12−14πk∑i=1[arcsin(disin(θi+1−δi))−arcsin(disin(θi−δi))]. (3.25)

Unwrapping the definitions we have

 1+c2i=1+a2i+b2i=1+(νi1)2(νi3)2+(νi2)2(νi3)2=1(νi3)2,\rm soc2i=(νi1)2+(νi2)2(νi3)2.

It follows that

 di=√(νi1)2+(νi2)2. (3.26)

We also note that

 δi=atan2(νi2,νi1),θi=atan2(yi,xi), (3.27)

where is any point along the edge .

### 3.2. Principal component analysis on local neighborhoods

The spherical volume invariant of a surface in 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  and analyze the shape of the region . In particular, it is suggested in 

to perform principal component analysis (PCA) on this region, that is, we compute the eigenvalues

of the symmetric matrix333Here we take to be a column vector.

 MS,r(p):=∫Ω∩Br(p)(x−¯¯¯x(p))(x−¯¯¯x(p))Tdx, (3.28)

where

 ¯¯¯x(p):=1VS,r(p)∫Ω∩Br(p)xdx (3.29)

is the centroid of , cf. (1.4). Assuming is sufficiently smooth, it was shown in  that the eigenvalues of have the asymptotic expansions

 (3.30)

where are the principal curvatures of the surface at the point , and, in the last formula, the term gives the mean curvature

 H(p)=12[κ1(p)+κ2(p)].

Moreover, the first two corresponding eigenvectors

are approximately tangent to the surface, and, assuming we are at a non-umbilic point, offer an approximation of the directions of principal curvatures, while is approximately normal to the surface and is an approximation of the unit normal. Thus, the matrix provides a robust estimation of the second fundamental form of at a non-umbilic point .

Let us now show how to compute the matrix via surface integrals, as we did for the spherical volume invariant in Theorem 1. While these results are mainly of interest in dimension , we carry out the derivation for an arbitrary dimension . Noting that

 MS,r(p)=∫Ω∩Br(p)(x−p)(x−p)Tdx−VS,r(p)(¯¯¯x(p)−p)(¯¯¯x(p)−p)T. (3.31)

it suffices to compute the first two moments

 mi(p):=∫Ω∩Br(p)(xi−pi)dx,cij(p):=∫Ω∩Br(p)(xi−pi)(xj−pj)dx, (3.32)

in terms of which the entry of is given by

 [MS,r(p)]i,j=cij(p)−1VS,r(p)mi(p)mj(p). (3.33)

The computation of and in terms of surface integrals is relatively straightforward, compared to the computation of . In what follows, denote the standard basis vectors in , and is the Kronecker delta.

###### Lemma 7.

Let us abbreviate . Then, for any , we have

 mi(p)=1n+1∫S∩Br(p)(yiy−r2ei)⋅νdS(x). (3.34)

and

 cij(p)=r2n+2VS,r(p)δij+12n+4∫S∩Br(p)(2yiyjy−r2(yjei+yiej))⋅νdS(x). (3.35)
###### Proof.

Without loss of generality, we may assume . Then and we write . We first prove (3.34). Define the vector field

 V(x)=xix−r2ein+1\rm so that% div V=xi.

By the Divergence Theorem,

 mi=∫Ω∩Brdiv Vdx=∫S∩BrV(x)⋅νdS+∫∂Ω∩BrV(x)⋅νdS.

On the spherical portion of the boundary , we have and so

 V(x)⋅ν=1rV(x)⋅x=xi(|x|2−r2)(n+1)r=0

since on . This completes the proof of (3.34).

We now prove (3.35). Define the vector field

 W(x)=2xixjx−r2(xjei+xiej)2n+4,% \rm wherebydiv W=xixj−1n+2r2δij.

By the Divergence Theorem, we have

 cij =∫Ω∩Br(1n+2r2δij+div W)dx =1n+2r2δijVS,r+∫S∩BrW(x)⋅νdS+∫∂Ω∩BrW(x)⋅νdS.

On the portion of the boundary

 W(x)⋅ν=1rW(x)⋅x=2xixjr2−r2(xjxi+xixj)(2n+4)r=0,

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

 ∫S∩Br(p)f(x)dS (4.1)

for various choices of kernel function . We adopt the convention that if , and hence rewrite (4.1) as simply

 ∫Sf(x)dS. (4.2)

Let denote the triangles in the triangulated surface . Then we can write

 ∫Sf(x)dS=M∑m=1 ∫Tmf(x)dS. (4.3)

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 used in this paper. Let us note that on the right hand side of (4.3), we need only sum over triangles that have non-empty intersection with . 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 , and terminating when all triangles in the connected component of containing 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 that do not contain . 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 . For simplicity, we take , write , and consider a triangle .

For the spherical volume invariant, for any triangle with the surface integral (3.13) from Theorem 1 requires us to compute

 A:=13∫T(1−r3|x|3)x⋅νdS.

Since is constant over the triangle , we have

 A=13z⋅ν(|T|−r3∫T1|x|3dS),

where is any point belonging to , such as its centroid or one of its vertices, while denotes the surface area of . The remaining integrand is known as a hypersingular

kernel, and arises, for instance, in the boundary element method for solving partial differential equations

. The integral of this hypersingular kernel over any planar triangle can be computed analytically  provided , which we may freely assume since when . For convenience, we recall the analytic formula, which is rather tedious and derived in , in Appendix A.

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

 14∫T(xix−r2ei)⋅νdS,\rm and110∫T[2xixjx−r2(xjei+xiej)]⋅νdS. (4.4)

Since and are constant over , we just need to compute the quantities

 ai:=∫TxidS,bij:=∫TxixjdS. (4.5)

Let us denote the vertices of by . The first integrand in (4.5) is linear, and so the integral can be computed analytically with the three point stencil

 ai=13|T|(xi+yi+zi). (4.6)

For , we compute the integral in barycentric coordinates

 bij =2|T|∫10∫1−t0((1−s−t)xi+syi+tzi)((1−s−t)xj+syj+tzj)dsdt, =2|T|[∫10∫1−t0(1−s−t)2xixj+s2yiyj+t2zizj+st(yizj+yjzi) +(1−s−t)s(xiyj+yjyi)+(1−s−t)t(xizj+zjzi)dsdt]

Computing

 ∫10∫1−t0(1−s−t)2dsdt=∫10∫1−t0s2dsdt=∫10∫1−t0t2dsdt=112,

and

 ∫10∫1−t0(1−s−t)sdsdt=∫10∫1−t0(1−s−t)tdsdt=∫10∫1−t0stdsdt=124,

we have