Categories
Math Research

Reduced deformation collision detection

Collision detection is an important problem in interactive computer graphics and physics-based simulation that seeks to determine if, when and where two or more objects come into contact. [4] In this project, we implement bounded deformation trees (BD-Trees) and adapt this method to represent complex deformations of any geometry as linear superpositions of displacement fields.

When an object collides with a surface, we should expect the object to deform in some way, e.g. if a bouncing ball is thrown against a wall or dropped from a building, it should momentarily be “squished” or flattened at the site of collision. This is the effect we aim to accomplish using modal analysis.

We start with a manifold triangular mesh e.g. Spot the cow, and tetrahedralize it using the python library of TetGen, a Delaunay-based tetrahedral mesh generator. [1] The resulting mesh is given as a \((V, C)\), where \(C\) is a set of tetrahedral cells whose vertices are in \(V\), as shown in Figure 1 below.

We use the Physics Based Animation Toolkit (PBAT) to compute the free vibrational modes of our model. Physically, one can describe vibration as the oscillatory motion of a physical structure, induced by energy exchanges of the potential (elastic deformation) and the kinetic (moving mass) energies. Vibrations are typically classified as either free or forced. In free vibrations, there are no continuous external forces acting on the structure, e.g. when a guitar string is plucked, while forced vibrations result from ongoing external forces. By looking at these free vibrations, we can determine the natural frequencies and normal modes of the structure.

First, we convert our geometric mesh into a FEM mesh and compute its Jacobian determinants and gradients of its shape function. You can check the documentation to learn more about FEM meshes.

import meshio
import numpy as np
import pbatoolkit as pbat

imesh = meshio.read("./spot.mesh")
V, C = imesh.points, imesh.cells_dict["tetra"]

mesh = pbat.fem.Mesh(V.T, C.T, element=pbat.fem.Element.Tetrahedron, order=1)

GNeU = pbat.fem.shape_function_gradients(mesh, quadrature_order=1)
detJeM = pbat.fem.jacobian_determinants(mesh, quadrature_order=2)
detJeU = pbat.fem.jacobian_determinants(mesh, quadrature_order=1)

Using these FEM quantities, we can model a hyperelastic material given its Young’s modulus \(Y\), Poisson’s ratio \(\nu\) and mass density \(\rho\).

rho = 1000.
Y = np.full(mesh.E.shape[1], 1e6)
nu = np.full(mesh.E.shape[1], 0.45)

# Compute mass matrix
M = pbat.fem.MassMatrix(mesh, detJeM, rho=rho, dims=3, quadrature_order=2).to_matrix()

# Define hyperelastic potential
hep = pbat.fem.HyperElasticPotential(mesh, detJeU, GNeU, Y, nu, energy=pbat.fem.HyperElasticEnergy.StableNeoHookean, quadrature_order=1)

Now we compute the Hessian matrix of the hyperelastic potential, and solve the generalized eigenvalue problem \(Av = \lambda M v\) using SciPy, where \(A\) denotes the Hessian matrix (a real symmetric matrix) and \(M\) denotes the mass matrix.

import scipy as sp

# Reshape matrix of vertices into a one-dimensional array
vs = mesh.X.reshape(mesh.X.shape[0]*mesh.X.shape[1], order="f")

hep.precompute_hessian_sparsity()
hep.compute_element_elasticity(vs)

HU = hep.hessian()
leigs, Veigs = sp.sparse.linalg.eigsh(HU, k=30, M=M, sigma=-1e-5, which="LM")

The resulting eigenvectors represent different deformation modes of the mesh. They can be animated as time continuous signals, as shown in Figure 2 below.

Reduced Deformation Models

The BD-Tree paper [2] introduced the bounded deformation tree, which can perform collision detection for reduced deformable models at similar costs to standard algorithms for rigid bodies. But what do we mean exactly by reduced deformable models? First, unlike rigid bodies, where collisions affect only the position or movement of the object, deformable bodies can dynamically change their shape when forces are applied. Naturally, collision detection is simpler for rigid bodies than for deformable ones. Second, instead of explicitly tracking every individual triangle in a mesh, reduced deformable models represent complex deformations efficiently by a smaller set of parameters. This is achieved by using a linear superposition of pre-computed displacement fields that capture the essential ways a model can deform.

Suppose we have a triangular mesh with \(|V| = n\). Let \(\boldsymbol{p} \in \mathbb{R}^{3n}\) denote the undeformed vertices locations, and let \(U \in \mathbb{R}^{3n \times r}\) be a matrix with \(r \ll n\). Then the new deformed vertices location \(\boldsymbol{p’}\) are approximated by a linear superposition of \(r\) displacement fields given by the columns of \(U\) such that

\(\displaystyle \boldsymbol{p’} = \boldsymbol{p} + U\boldsymbol{q},\)

where the amplitude of each displacement field is determined by the reduced coordinates \(\boldsymbol{q} \in \mathbb{R}^{r}\). Both \(U\) and \(\boldsymbol{q}\) must already be known in advance. In our case, the columns of \(U\) are the eigenvectors obtained from modal analysis described earlier, although they could also result from methods, e.g. an interpolation process. The reduced coordinates \(\boldsymbol{q}\) could also be determined by some possibly non-linear black box process. This is important to note: although the shape model is linear, the deformation process itself can be arbitrary!

Bounded deformation trees

Welzl’s algorithm

The BD-Tree works by constructing a hierarchy of minimum bounding spheres. As a first step, we need a method to construct the smallest enclosing sphere for some set of points. Fortunately, this problem has been well studied in the field of computational geometry, and we can use the randomized recursive algorithm of Welzl [3] that runs in expected linear time.

The Welzl’s algorithm is based on a simple observation: assume a minimum bounding sphere \(S\) has been computed a set of points \(P\). If a new point \(p\) is added to \(P\), then \(S\) needs to be recomputed only if \(p\) lies outside of \(S\), and the new point \(p\) must lie on the boundary of the new minimum bounding sphere for the points \(P \cup \{p\}\). So the algorithm keeps track of the set of input points and a set of support, which contains the points from the input set that must lie on the boundary of the minimum bounding sphere.

The pseudocode below [4] outlines the algorithm:

Sphere WelzlSphere(Point pt[], unsigned int numPts, Point sos[], unsigned int numSos)
{ 
	// if no input points, the recursion has bottomed out.
    // Now compute an exact sphere based on points in set of support (zero through four points)
	if (numPts == 0) {
		switch (numSos) {
			case 0: return Sphere();
			case 1: return Sphere(sos[0]);
			case 2: return Sphere(sos[0], sos[1]);
			case 3: return Sphere(sos[0], sos[1], sos[2]);
			case 4: return Sphere(sos[0], sos[1], sos[2], sos[3]); 
		}
	}
	// Pick a point at "random" (here just the last point of the input set)
	int index = numPts - 1;
	// Recursively compute the smallest bounding sphere of the remaining points
	Sphere smallestSphere = WelzlSphere(pt, numPts - 1, sos, numSos); 
	// If the selected point lies inside this sphere, it is indeed the smallest
	if(PointInsideSphere(pt[index], smallestSphere))
		return smallestSphere;
	// Otherwise, update set of support to additionally contain the new point
	sos[numSos] = pt[index];
	// Recursively compute the smallest sphere of remaining points with new s.o.s. 
	return WelzlSphere(pt, numPts - 1, sos, numSos + 1);
}
The BD-Tree Method
Figure 3: Wrapped BD-Tree for Spot at increasing recursion levels.

Now that we can compute the minimum bounding spheres for any set of points, we are ready to construct a hierarchical sphere tree on the undeformed model, after which it can be updated following deformation. First we note that the BD-Tree is a wrapped hierarchy, wherein the bounding spheres tightly enclose the underlying geometry but any bounding sphere at one level need not contain its child spheres. This is different from a layered hierarchy in which spheres must enclose their child spheres, but can fit the underlying geometry more loosely (see Figure 4).

As shown in Figure 5, there are many possible approaches to building a binary tree. In our case, we use a simple top-down approach while partitioning the underlying geometry, i.e. we recursively split Spot at its median into two parts with respect to its local coordination axes, such that a leaf node (the lowest level) contains only one triangle.

Figure 5: Hierarchical binary tree construction with four objects using a top-down (top), bottom-up (middle) and insertion approach (bottom).

As the object deforms, how do we compute the new bounding spheres quickly and efficiently?

Let \(S\) denote a sphere in the hierarchical tree with center \(c\) and radius \(R\) containing the \(k\) points of the geometry \(\{p_{i}\}_{1 \leq i \leq k}\). After the deformation, the center of the sphere \(c\) is displaced by a weighted average of the contained points’ displacements \(u_{i}\) with weights \(\beta_i\) e.g. \(\beta_i := 1/k\) for \(1 \leq i \leq k\). So the new center can be expressed as

\(\displaystyle c’ = c + \sum_{i = 1}^{k} \beta_i u_i.\)

Using the displacement equation above, we can write \(u_i\) as the sum \(\sum_{j = 1}^{r} U_{ij} q_{j}\) and substitute this into the previous equation to obtain:

\(\displaystyle c’ = c + \sum_{j = 1}^{r} \left(\sum_{i = 1}^{k} \beta_i U_{ij}\right) q_j = c + \bar{U} \boldsymbol{q}.\)

To compute the new radius \(R’\), we make use of the triangle inequality

\(\displaystyle \max_{i = 1, \dots, k} || p’_{i} – c’_{i} ||_{2} \leq \max_{i = 1, \dots, k} \left(||p’_{i}||_{2} + ||c’_{i}||_{2} \right).\)

Expanding both \(p’_{i}\) and \(c’_{i}\) using their displacement equations and re-arranging the terms, we get that

\(\displaystyle R’ = R + \sum_{j = 1}^{r} \left( \max_{i = 1, \dots, k} ||U_{ij} – \bar{U}_{j}||_{2} \right) |q_j|.\)

Thus, we have an updated bounding sphere \(S’\) with its center \(c’\) and radius \(R’\) computed as functions of the reduced coordinates \(\boldsymbol{q}\).

References

[1] Hang Si. TetGen, a Delaunay-based quality tetrahedral mesh generator. ACM Transactions on Mathematical Software, 41(2), 2015.

[2] Doug L. James and Dinesh K. Pai. BD-Tree: Output-Sensitive Collision Detection for Reduced Deformable Models. ACM Transactions on Graphics, 23(3), 2004.

[3] Emo Welzl. Smallest enclosing disks (balls and ellipsoids). New Results and New Trends in Computer Science, 555, 1991.

[4] Christer Ericson. Real-Time Collision Detection. CRC Press, Taylor & Francis Group. Chapter 4, p. 99-100. 2005.

Categories
Research

Higher-Order Discretization of Mean Curvature Flow

Author: Ehsan Shams (Alexandria University, EG)

Project Mentor: Justin Solomon (MIT, USA)
Volunteers: Biruk Ambaw (Université Paris-Saclay, Fr), Andrew Rodriguez (Georgia Institute of Technology, USA), and Josh Vekhter ( UT Austin, USA)

Acknowledgments. Sincere thanks to Professor Justin Solomon for his invaluable guidance throughout this project. His carefully designed questions and coding tasks not only deepened my understanding of core topics in geometry processing in a short amount of time, but also sharpened my coding skills—and ensured my lunch breaks were notably shorter than they might have been otherwise :). I would also like to thank Josh Vekhter, Andrew, and Biruk for their support and valuable feedback to my teammates and me. In addition to the math and codes, Professor Justin and Josh taught me the value of a well-timed joke to lighten the load!

Introduction

Mean Curvature Flow (MCF) is a fundamental concept in differential geometry, that describes the process by which a surface evolves over time under the influence of its mean curvature. The concept is of great utility to many applications in geometry processing (GP), including surface smoothing and mesh denoising. MCF is described by a parabolic partial differential equation (PDE). While spatial discretization of this PDE on triangle meshes is well-established, temporal discretization presents significant challenges.

Three primary approaches exist for temporal discretization: explicit, implicit, and semi-implicit methods. Explicit methods, like forward Euler, often become unstable with larger time steps. Implicit methods, such as backward Euler, offer stability but at a high computational cost. Semi-implicit methods, like the one introduced by Desbrun et al. (1999), strike a compromise between these two extremes but may still fall short in terms of accuracy and stability. In practice, the temporal discretizations that are commonly employed in the literature are only first-order accurate thus, may not provide the desired level of accuracy, and stability we ideally wish for.

The central focus of this project was to derive and implement higher-order temporal discretization methods, both explicit and semi-implicit, for MCF on triangle meshes, with the goal of improving the accuracy of the numerical solution to the PDE while also ensuring stability throughout the flow.

In this article, I externalize my internal exploratory journey and insights gained during my last SGI research week under the guidance of my mentors in details. The article begins with an introduction to MCF and its significance in geometry processing, followed by a brief discussion on the process of discretization followed by the derivation of spatial discretization via finite element methods and first-order accurate temporal discretization via finite difference methods for MCF. A theoretical comparison of the methods in question is provided, highlighting their pros and cons. Finally, we derive the second-order accurate temporal discretization of MCF, in both, the explicit, and semi-implicit schemes. Empirical validations of theoretical results are presented, along with a humorous fail. The goal here is on exposition, rather than taking the shortest path.

Understanding Mean Curvature Flow in \(\mathbb{R}^3 \)

Mean curvature flow (MCF) of a surface in \(\mathbb{R}^3 \) is a geometric evolutionary process in which each point on the surface moves in the direction of the surface’s normal vector, with a velocity proportional to the mean curvature at that point which results in smoothing the irregularities in the surface. For strictly convex surfaces, the flow causes the surface to shrink uniformly in volume until it collapses into a point over a finite time while preserving convexity (Huisken, 1984). MCF is a key tool in differential geometry and geometric analysis due to its ability to drive the evolution of a surface purely based on its intrinsic geometry.

Formal Definition

Let \(M_t \subset \mathbb{R}^3\) represent a family of smoothly embedded surfaces parameterized by time \(t\). The surface at time \(t\) can be described by a smooth mapping:

\[ X(\mathbf{u}, t) : U \times [0, T) \rightarrow \mathbb{R}^3 \]

Here:

  • \(U\) is an open set in \(\mathbb{R}^2\) representing the parameter space of the surface, with local coordinates \( \mathbf{u} = (u_1, u_2)\).
  • \([0,T) \) represents time, where \(T\) is the time until which the flow is considered.
  • \(X(\mathbf{u}, t)\) is the position vector of a point on the surface in \(\mathbb{R}^3\) at time \(t\), and can be explicitly written as: \[ X(\mathbf{u},t) = \begin{pmatrix} X_1(\mathbf{u},t) \\ X_2(\mathbf{u},t) \\ X_3(\mathbf{u},t) \end{pmatrix} \] and \(X_1(\mathbf{u},t), X_2(\mathbf{u},t), X_3(\mathbf{u},t)\) are the coordinate functions that determine the \(x, y,\) and \(z\) components of the position vector in \( \mathbb{R}^3\)

The mean curvature flow is then dictated by the following partial differential equation:

\[ \frac{\partial X}{\partial t}(\mathbf{u}, t) = – H(\mathbf{u}, t) \mathbf{n}(\mathbf{u}, t) =\Delta X(\mathbf{u}, t)\]

Where:

  • \(\frac{\partial X}{\partial t}(\mathbf{u}, t)\) is the velocity of the surface at the point \(X(\mathbf{u}, t)\).
  • \(H(\mathbf{u}, t)\) is the mean curvature at the point \(X(\mathbf{u}, t)\).
  • \(\mathbf{n}(\mathbf{u}, t)\) is the unit normal vector at the point \(X(\mathbf{u}, t)\), pointing outward or inward. It indicates the direction in which the surface will move during the flow.
  • \( \Delta \) is the Laplace-Beltrami operator associated with the surface \(M_t\).

Abuse of Notation. For simplicity, any operator \( \phi (\mathbf{u},t) \) at a point parameterized by \( (\mathbf{u},t))\) will be referred to simply as \( \phi \) in this article. The Laplace-Beltrami operator, might sometimes be referred to as the Laplacian depending on my mood.

Mean Curvature \(H\):

The mean curvature \( H\) at a point is defined as the average of the principal curvatures \( k_1 \) and \( k_2 \) at that point on the surface: \[ H= \frac{1}{2} (k_1 + k_2) \]

Where \(k_1\) and \(k_2\) are the eigenvalues of the second fundamental form of the surface. These principal curvatures measure the maximum and minimum bending of the surface in orthogonal directions at this particular point in question.

Alternative Expressions for \( H \):

  • Divergence of the Normal Vector Field: \(H\) can be expressed as the negative half of the divergence of the unit normal vector field \(\mathbf{n}\) on the surface: \( H= – \frac{1}{2} div ((\mathbf{n})) \).
  • Trace of the Shape Operator: In addition, \(H\) can be written as the trace of the shape operator (or Weingarten map) associated with the surface. This representation connects the mean curvature to the linear transformation that describes how the surface bends.

Why Study MCF?

Understanding MCF begins with a simple question: what happens when a surface evolves according to its own curvature? In practice, this flow helps smooth surfaces over time, driving them toward more regular shapes. This process is especially intriguing because it captures essential properties of geometric evolution without external forces, relying solely on the surface’s own shape.

A fundamental result by Huisken in 1984 in his paper titled “Flow by Mean Curvature of Convex Surfaces into Spheres” provides deep insight into this phenomenon. Huisken studied the evolution of a special type of hypersurfaces under MCF, and proved that any strictly convex smooth hypersurface, under MCF, evolves into a sphere before collapsing to a point in a finite, self-similar manner in a finite time. His work highlights how MCF transforms surfaces, regularizing their shape as they shrink and demonstrating the flow’s inherent tendency to round out irregularities. In mathematics, it is common practice to first test an idea or prove a theorem in simple and nice settings (such as convex compact spaces) before attempting to test/generalize to more abstract spaces.

While the general theory of MCF applies to a wide range of surfaces, it is particularly insightful to consider how a sphere—a familiar, symmetric object—shrinks under the flow. The sphere offers a clean, intuitive example where the underlying mathematics remains manageable, yet it showcases the essential features of MCF, such as shrinking to a point while maintaining symmetry. This serves as an ideal starting point for understanding more complex behaviors in less regular surfaces, grounded in the theoretical framework established by Huisken’s work.

Example: Shrinking Sphere Under MCF

Consider a sphere of radius \(R(t)\) at time \(t\). Under MCF, the sphere’s surface shrinks over time. We want to find how the volume of the sphere changes as it evolves under this flow.

For a sphere of radius \(R(t)\), the mean curvature \(H\) is given by: \( H=\frac{2}{R(t)}\), and the normal vector \(\mathbf{n}\) points radially inward, so: \(\mathbf{n}=\frac{−X}{R(t)}\)

Substituting these in the MCF equation: \( \frac{\partial X}{\partial t}=\frac{-2}{R(t)} \frac{-X}{R(t)}=\frac{2X}{R^2(t)}\)

To determine how the radius \(R(t)\) changes over time, observe that each point on the surface \(X\) can be written as: \(X=R(t) \mathbf{u}\) where \(\mathbf{u}\) is a unit vector in the direction of \(X\). Thus, \( \frac{\partial X}{\partial t} = \frac{\partial (R(t) \mathbf{u})}{\partial t} = \frac{2\mathbf{u}}{R(t)} \), which simplifies to \( \mathbf{u}\frac{dR(t)}{dt}= \frac{2 \mathbf{u}}{R(t)}\) so \(\frac{dR(t)}{dt}= \frac{2 }{R(t)} \).

To solve for \( R(t) \), we separate the variables and integrate \(\int R(t) \, dR(t) = \int 2 \, dt \) which gives \(\frac{R^2(t)}{2} = 2t + C \), where \( C \) is the constant of integration. Solving for \( R(t) \) gives \( R(t) = \sqrt{4t + C’} \) where \( C’ = 2C \). If we consider the initial condition \( R(0) = R_0 \), then \( R_0^2 = C’ \), and hence: \(R(t) = \sqrt{R_0^2 – 4t} \)

Volume Shrinkage:

The volume \( V(t) \) of the sphere at time \( t \) is \( V(t) = \frac{4}{3} \pi R^3(t) \). As the sphere shrinks over time, the volume decreases according to the following equation: \( V(t) = \frac{4}{3} \pi \left(R_0^2 – 4t \right)^{3/2} \), and it continues to shrink until \( t = \frac{R_0^2}{4} \), at which point the radius reaches zero and the sphere vanishes.

The accompanying animation illustrates the shrinking process, showing how both the radius and volume decrease over time under the influence of MCF until the sphere vanishes. The MCF drives the surface points of the sphere inward in the direction of the inward normals, leading to a uniform shrinking process. For a perfectly symmetric object like a sphere, the shrinking occurs uniformly, preserving the spherical shape until the radius reduces to zero.

The shrinking radius is computed via the derived formula above, \( R(t)= \sqrt{R_0^2 – 4t} \), and the shrinking becomes faster as time progress. The plots that are displayed at the end show the relation between \( R(t) \) and \( V(t) \) over time.

Shrinking sphere under MCF + evolution over time plots

Since machine computers operate in a discrete, digital environment, and GP is inherently computational, continuous geometric surfaces must be represented in a form that machine computers can process. This form is typically a mesh, a discrete approximation of a surface composed of vertices, edges, and faces, and given that MCF is intrinsically a continuous model described by a partial differential equation in space and time, directly applying it to meshes necessitates the discretization of its continuous operators (the Laplacian and time derivatives) as well.

Discretizing MCF

As implied by the above paragraph, discretization is the process of transforming a continuous mathematical model \(\mathcal{M}\), which is defined over a continuous domain \(\Omega \subset \mathbb{R}^n\) and governed by continuous variables \(u(x)\) and operators \(\mathcal{L}\), into a discrete model \(\mathcal{M}_h\) suitable for computational analysis. This involves:

  • Replacing the continuous domain \(\Omega\) with a finite set of discrete points or elements \(\Omega_h = {x_i}_{i=1}^N\).
  • Approximating continuous functions \(u(x)\) by discrete counterparts \(u_h(x_i)\).
  • Substituting continuous operators \(\mathcal{L}\) with discrete analogs \(\mathcal{L}_h\).

The goal is to ensure that the discrete model \(\mathcal{M}_h\) preserves key properties such as consistency, stability, and convergence, so that \(\mathcal{M}_h\) faithfully reflects the behavior of the continuous model \(\mathcal{M}\) as the discretization is refined, i.e., as \(h \rightarrow 0\).

While many continuous models can be discretized, the accuracy and efficiency of the approximation depend on the chosen discretization technique and its implementation. Common discretization techniques include: Finite difference methods (FDM), finite element methods (FEM), finite volume methods (FVM), and spectral methods.

Discretizing the MCF model can be systematically divided into two main components: spatial discretization and temporal discretization. In this article, we employ the FEM for spatial discretization and the FDM for temporal discretization.

Spatial Discretization via FEM

Spatial discretization involves representing the continuous spatial domain of the model as a mesh and approximating differential operators—primarily the Laplacian, and also gradient and divergence as quantities at the discrete mesh in question.

Assume the surface in question is discretized into a mesh which consists of a collection of triangular elements \(\{T_i\}\) with vertices \(v_j\). The first step to approximate the Laplace-Beltrami operator is to define basis functions \(\{\phi_j\}\) associated with the vertices and defined over \(T_i\). These are often linear or piecewise linear functions, \(\phi_j\) takes the value 1 at vertex \(v_j\)​ and 0 at all other vertices. These basis functions form a set that allows any function \(X\) on the surface to be approximated as a linear combination of them \( X≈ \sum_j X_j \phi_j\), where \(\mathbf{X_j}=\mathbf{X(v_j)}\) represents the value of the position function at vertex \(v_j\).

Next, we compute the gradient of each basis function within each triangle, which will be constant over the triangle because the functions are linear. The Laplacian operator, which in continuous terms is the divergence of the gradient, is discretized by integrating the product of the gradients of the basis functions over the surface. This leads to the construction of the stiffness matrix \(\mathbf{L}\), where each entry \( \mathbf{L_{ij}}\) is derived from the inner products of the gradients of the basis functions \(\phi_i\) and \(\phi_j\) ( \( \nabla \phi_i\) and \(\nabla \phi_j\)), weighted by the cotangent of the angles opposite the edge connecting vertices \(v_i\), and \(v_j\). The diagonal entries \(\mathbf{L_{ii}}\)​ sum the contributions from all adjacent vertices.

More precisely, for two vertices \(v_i\) and \(v_j\) connected by an edge, the stiffness matrix \(\mathbf{L}\) is computed using the following integral:

\[ \mathbf{L_{ij}}=\int_{\sigma} \nabla \phi_i \cdot \nabla \phi_j dA\]
where \(\sigma\) is the surface, and \(dA\) is the area element.

For a piecewise linear basis function on a triangular mesh, those gradients in question are constant within each triangle, so the integral simplifies to the sum over the triangles \(T_k\) that include the edge connecting \(v_i\) and \(v_j\):

\[ \mathbf{L_{ij}}= \sum_{T_k \in N(v_i,v_j)} (\nabla \phi_i \cdot \nabla \phi_j ) \text{Area} (T_k) \]

Where \( N(v_i,v_j)\) denotes the set of triangles sharing the edge between \(v_i\) and \(v_j\). Now the inner product \(\nabla \phi_i \cdot \nabla \phi_j\) can be computed using the cotangent formula
\[\nabla \phi_i \cdot \nabla \phi_j = \frac{-1}{2 \text{Area}(T_k)}(\cot{\alpha_{ij}}+\cot{\beta_{ij}}) \]

Thus, the entries \( \mathbf{L_{ij}} \) are:

\[\mathbf{L_{ij}}= – \frac{1}{2} (\cot{\alpha_{ij}}+\cot{\beta_{ij}}) \]

Finally the diangonal enteries \(\mathbf{L_{ii}}\) sum the contributions from all the adjacent verticies:
\[\mathbf{L_{ii}}=\sum_{i \neq j} \mathbf{L_{ij}} = \sum_{i \neq j} \frac{1}{2} (\cot{\alpha_{ij}}+\cot{\beta_{ij}})\]

To balance the discretization, we also need the mass matrix \(\mathbf{M}\), which arises from integrating the basis functions themselves. Formally, the entries of the mass matrix \(\mathbf{M}\) are given by:
\[ \mathbf{M_{ij}} = \int_{\sigma} \phi_i \phi_j dA\]

In simple cases1, the mass matrix is typically diagonal, with each entry written as:
\[ \mathbf{M}_{ii} = \frac{1}{3} \sum_{T \ni v_i} \text{Area}(T),\]

where the sum is over all triangles \(T\) that share the vertex \((v_i).\)
The final discretized Laplace-Beltrami operator is represented by the generalized eigenvalue problem \(\mathbf{L}X = \lambda \mathbf{M} X\), where \(\mathbf{L}\) encodes the differential operator and \( \mathbf{M}\) ensures that the discretization respects the geometry of the mesh.

Both matrices (the mass matrix \( \mathbf{M}\) and the stiffness matrix \( \mathbf{L}\)) are positive semi-definite, and sparse.

Intuitively,

  • The stiffness matrix \(\mathbf{L}\) captures the geometric and differential properties of the surface. It represents how the curvature is distributed over the mesh by measuring how the gradients of the basis functions (associated with each vertex) interact with each other. The entries of this matrix are weighted by the angles in the triangles, which essentially encode how “stiff” or “resistant” the mesh is to deformation. In other words, it determines how much the shape of the surface resists change when forces (like curvature flow) are applied.
  • The mass matrix \(\mathbf{M}\) accounts for the distribution of area or “mass” over the surface. It ensures that the discretization respects the surface’s geometry by distributing weight across the vertices according to the areas of the surrounding triangles. This matrix is often diagonal, with each entry corresponding to the area associated with a vertex, making sure that the mesh’s physical properties, like area and volume, are properly balanced in computations.

Solving the general eignvalue problem, we reach the following approximation of the Laplacian: \(\Delta \approx \mathbf{M}^{-1} \mathbf{L}\).

Remark: In discretizing the Laplacian operator, several other approaches can be employed too such as divided difference, higher-order elements, discrete exterior calculus,..etc. In any case, no discretization approach of the Laplacian could keep every natural property of its ideal continuous form.

Temporal Discretization via FDM

Temporal discretization refers to the approximation of the time-dependent aspects of the MCF model (the time derivative \(\frac{\partial X}{\partial t}\)). This step is critical for evolving the surface over time according to the mean curvature dynamics. Three common approaches are used for this type of discretization: Explicit methods (e.g. Explicit Euler method), Implicit methods (e.g. Implicit Euler method), and Semi-implicit methods (e.g Desbrun et al. (1999)).

  • Explicit Euler Method: This method (also known as forward Euler method) is where the new future positions of the mesh vertices are computed based on the current positions and the mean curvature at those positions. While this method is simple to implement, it may impose stability constraints on the time step size. The discretized vertex update rule is given by:
    \[\frac{X^{(k+1)} – X^{(k)}}{\Delta t} = \mathbf{M}^{-1} \mathbf{L} X^{(k)} \]
    Rearranging this, we get: \[ X^{(k+1)} = X^{(k)} + \Delta t \mathbf{M}^{-1}\mathbf{L} X^{(k)} \] where \(X^{(k)}\) and \(X^{(k+1)}\) are surface position matrix of the mesh vertices at time steps \(k\) and \(k+1\), respectively.
  • Implicit Euler Method: Also known as backward Euler, this method involves solving a system of equations at each time step, it might provide greater stability compared to the Explicit Euler method and allow for larger time steps. The discretized vertex update rule is given by: \[ X^{(k+1)} \approx X^{(k)} + \Delta t \mathbf{M}^{-1}\mathbf{L} X^{(k+1)} \] This equation is said to be fully implicit because the Laplacian depends on the vertex positions at the time step one is trying to solve for, making it nonlinear and difficult to solve.
  • Desbrun et al.’s Semi-Implicit Method: (Desbrun et al. 1999) proposed a semi-implicit method, which is a compromise between the simplicity of explicit methods and the stability of implicit methods. The idea is to treat the Laplacian and the vertex positions in the following manner: Instead of calculating the Laplacian at the next time step (which makes the equation nonlinear), they compute it at the current time step which simplifies the problem, and the vertex positions are still updated at the next time step. Their discretized vertex update rule is given by: \[ X^{(k+1)} \approx X^{(k)} + \Delta t\mathbf{M^{-1}}(\mathbf{L}X^{(k)} )X^{(k+1)} \] This equation is still implicit in \(X^{(k+1)}\), but the Laplacian is evaluated at the known positions \(X^{(k)}\), making the system linear and easier to solve. The update rule can be re-arranged into: \[X^{(k+1)} \approx (I- \Delta t\mathbf{M^{-1}L}X^{(k)} )^{-1}X^{(k)} \]. This method offers nice stability. However, it is not as accurate as fully implicit methods, because it only approximates the Laplacian based on the current positions. It might smooth the mesh, but not as precisely as solving the full nonlinear system such as the fully-implicit. In addition, it does not generalize to other geometric flows that require more complex handling of nonlinearities.
  • Adaptive Time Stepping: This approach adjusts the time step size dynamically based on the evolution of the surface, allowing for finer resolution during critical changes and coarser resolution during smoother phases. Progyan Das from my team worked on this approach.

Remark. The derivations for the above update rules are tacitly included in the part we derive their second-order forms.

Comparative Analysis and Empirical Validation: Euler Explicit vs. Desbrun et al.’s

  1. Accuracy
    • Euler Explicit Method:
      • Order of Accuracy: This method is first-order accurate in time, meaning that the global error in the solution decreases linearly with the time step \(\Delta t\). If the exact solution at time \(t\) is denoted by \(X_{t}\) and the numerical solution by at time \(t_n\) by \(X_{t_n}\)​, the error \(E_n = X(t_n) \) satisfies\(E_n \approx C \Delta t\), where \(C\) is a constant dependent on the problem. This linear relationship implies that halving the time step approximately halves the error.
      • Error Propagation: Errors tend to accumulate more rapidly, especially for larger time steps. Because this method updates the solution based only on information from the current time step so, if the time step \(\Delta t\) is too large, the method may not accurately capture the evolution of the curvature, leading to significant errors that compound over time. The error propagation can be expressed as \(X^{(k+1)}=X^{(k)}+ \Delta t \dot F(X^{(k)})\), where \(F(X^{(k)})\) is the update function. If\(\Delta t\) is too large, the local truncation error, which is \(O(\Delta t^2)\), becomes significant, causing larger cumulative errors.
      • Handling of Complex Geometries: Will probably struggle with highly irregular meshes, which is a direct consequence of the above bullet. leading to larger errors in curvature computation.
    • Desbrun et al. Semi-Implicit Method:
      • Order of Accuracy: This method is also first-order accurate in time because it is essentially a modified backward Euler scheme, where the implicit part is handled for spatial discretization, but the time discretization remains first-order.
      • Error Propagation Reduction: The method implicitly handles the curvature of the mesh by solving a linear system at every update, which incorporates more information about the solution at the next time step. This implicit approach effectively reduces errors, and stabilizes the solution especially when larger time steps are used compared to Euler’s explicit method.
      • Numerical Diffusion: Moreover, it has a better control over numerical diffusion —a phenomenon where fine details of the mesh are smoothed out excessively—compared to the explicit method, leading to more accurate smoothing. Numerical diffusion can be mathematically described by how the curvature smoothing term affects the higher-order modes of the solution and here is where the implicit nature of the method helps preserve these modes more effectively than Euler’s explicit method.
  2. Stability
    • Euler Explicit Method:
      • Conditionally Stable: The stability here depends on the time step size; it requires small time steps to maintain stability.
      • CFL Condition: The time step must satisfy the Courant-Friedrichs-Lewy (CFL) condition, which can severely restrict the time step size, especially for fine meshes. The CFL condition constrains the time step to be inversely proportional to the square of the mesh resolution. This means that as the mesh becomes finer, the time step must decrease quadratically, which significantly increases the number of iterations required for convergence.
    • Desbrun et al. Semi-Implicit Method:
      • Unconditionally Stable: Allows larger time steps without sacrificing stability. This is a key advantage for computational efficiency.
      • Robustness: More stable under large deformations or irregular meshes, making it suitable for a broader range of applications than the explicit method.
  3. Computational Efficiency and Memory Usage
    • Euler Explicit Method:
      • Efficiency: Simpler to implement and faster per iteration due to direct updates, but requires more iterations for convergence due to the small time steps needed.
      • Memory Usage: Lower memory requirements since it does not require solving linear systems.
      • Parallelization: Easier to parallelize due to the independence of the update steps.
    • Desbrun et al. Semi-Implicit Method:
      • Efficiency: More computationally intensive per iteration due to the need to solve linear systems, but fewer iterations may be needed due to larger permissible time steps.
      • Memory Usage: Higher memory consumption due to the storage of matrices for linear system solving.
      • Parallelization: More challenging to parallelize because of the dependencies introduced by solving the linear system.
  4. Implementation Complexity
    • Euler Explicit Method:
      • Complexity: Conceptually simpler and easier to implement. It involves straightforward updates without the need for solving linear systems.
      • Dependencies: Minimal dependencies between updates, making it a more accessible method for quick implementations.
    • Desbrun et al. Semi-Implicit Method:
      • Complexity: More complex to implement due to the need to solve large, sparse linear systems at each time step.
      • Dependencies: Involves matrix assembly and inversion, which can introduce additional challenges in implementation.
  5. Parameter Sensitivity
  • Euler Explicit Method:
    • Sensitivity: Highly sensitive to time step size. Small changes can significantly affect stability and accuracy.
  • Desbrun et al. Semi-Implicit Method:
    • Sensitivity: Less sensitive to time step size, allowing for greater flexibility in choosing time steps.

Overall Assessment:

  • Euler Explicit Method is advantageous for its simplicity, ease of implementation, and parallelization potential. However, it is limited by stability constraints, accuracy issues, and higher sensitivity to parameter choices.
  • Desbrun et al. Semi-Implicit Method offers superior stability, accuracy when compared to the explicit, and reduced numerical diffusion at the cost of increased computational complexity and memory usage. It is better suited for applications requiring robust and accurate smoothing, particularly in the context of complex or irregular meshes.

Oh… this felt like eating five horrible McDonald’s cheeseburgers. 🍔🍔🍔🍔🍔 Right? So, let’s compress this previous analysis into a nice compact table for quick reference.

AspectEuler Explicit MethodDesbrun et al. Semi-Implicit Method
AccuracyFirst-order accurate in time.
Higher error accumulation, especially for large time steps.
Struggles with complex geometries.
First-order accurate in time.
Better error reduction, especially for large time steps.
Better control over numerical diffusion.
StabilityConditionally stable.
Requires small time steps, dictated by the CFL condition.
Unconditionally stable.
Allows larger time steps without sacrificing stability.
Computational EfficiencySimple and fast per iteration.
Inefficient for fine meshes due to small time step requirement.
Computationally more expensive due to solving linear systems.
Efficient for larger time steps.
Memory UsageLower memory usage.Higher memory usage due to storing and solving linear systems.
Implementation ComplexityRelatively simple to implement.More complex due to the need to solve linear systems.
ParallelizationEasier to parallelize due to independent updates.More challenging to parallelize due to solving linear systems.
Parameter SensitivitySensitive to time step size (CFL condition).Less sensitive to time step size, allowing more flexibility.
Numerical DiffusionHigher numerical diffusion, especially for large time steps.Better control over numerical diffusion, preserving more detail.
Euler explicit vs Desbrun et al. semi-implicit methods compact comparative assessment

The next experiment aims to empirically illustrate the above tradeoffs in the stability and accuracy aspects presented above for the two methods in question for MCF since the main aim of the project is to come up with an approach to improving the accuracy and stability of MCF applied to triangle meshes by incorporating higher-order derivatives in the time integration process.

Experimental Design:

  1. Mesh Preparation:
    • Target Mesh: Load a 3D mesh model, and store it as the target mesh.
    • Noisy Mesh: Add a controlled amount of noise to simulate imperfections to the target mesh and store the output as the noisy mesh.
  2. Application of MCF Methods:
    • Apply each MCF method to the noisy mesh across a series of pre-determined range of time steps. The number of iterations for both methods is fixed (10 iterations) per each time step. Choose the range of the time steps such that it includes very small (those satisfying the CFL condition) to relatively large time steps to allow for a comprehensive analysis of the methods’ behavior under different conditions.
  3. Data Logging: For each time step, record the following data for each method:
    • Error Metric: Frobenius norm of the difference between the smoothed mesh and the original target mesh.
    • Stability Metric: the number of maximum \(\Delta t\) where the error stays less than a pre-determined error threshold by the user.
    • Computational Time: Time taken for execution.
  4. Output Plots
    • Error vs. Time Steps: The error is plotted against the time step size on a log-log scale. The slope of this curve will indicate the convergence rate:
      • Steeper Slope: Indicates faster convergence and higher accuracy.
      • Flatter Slope: Suggests slower convergence and potential inaccuracies.
    • Visual and Quantitative Evaluation: The final smoothed meshes are presented for visual comparison against the target mesh.
    • Time Step Limitations: The maximum time step for which each method remains stable is identified.
Frobenius norm vs time steps size for Euler’s explicit and Desbrun et al.’s semi-implicit methods
Runtime (seconds) vs time steps size for Euler’s explicit and Desbrun et al.’s semi-implicit methods

Interpretation of Results:

  • The accuracy plot shows that the Explicit method tends to produce much higher errors compared to the Desbrun’s especially at larger time steps, confirming that it requires small time steps to maintain stability. In addition, we can see that with Desbrun’s, the error remains constant across a wide range of varying time steps, this confirms that the method’s accuracy and stability are not affected by the choice of time step within the range tested, it also verifies its robustness, and tendency to reduce error propagation as time step sizes increase. In regards to the rate of convergence, the slope of the error curve vs time-steps in Euler’s explicit is \(m=1.29\), and Desbrun’s method rate of convergence is \(m=1\).
  • In this experiment, we set the error threshold to 3. Turns out the maximum step size \(\Delta t\) where the error is below this threshold for Euler’s explicit is 2.2204e-16 while for Desbrun et al.’s is (the final in our range).

Inspiration: This experimental design was inspired by a task I did during my second SGI project on 2D Differentiable Representation of Curve Networks, under Mikhail Bessmeltsev. 🙂

The Case for Higher-Order Time Discretization

As demonstrated by both theory and practice, even robust methods like Desbrun et al.’s semi-implicit method for MCF face limitations with first-order time discretization. While this category of methods offers a compromise between explicit and fully implicit methods, first-order discretization still imposes constraints on accuracy in numerical simulations. These limitations stem from the truncation errors inherent in first-order approximations.

First-order methods approximate the time derivative using only vertex velocities (the first derivative of position with respect to time) and disregard higher-order terms, such as acceleration (the second derivative). This omission means they fail to account for how the geometry of the surface might be changing or accelerating locally since higher derivatives encode information about local curvatures. If we conceptualize the next iterate as\(X^{(k+1)}_i = \mu (X_i^{(k)}, I^{(k)})\) where \(I^{(k)}\) is an information vector, local geometric properties of the surface at the current iterate \(X_i^{(k)}\), it becomes clear that the more detailed information we incorporate into \(I^{(k)}\), the more accurate the next state \(X^{(k+1)}_i \) will be.

In other words, higher-order discretizations, such as second-order ones, lead to a significant reduction in truncation errors, better convergence, and a more accurate representation of the geometry over time for larger time steps which contributes to a more economical utilization of computational resources (e.g reduced number of iterations).

Deriving Higher-order Discretizations for MCF

In this section, we derive the second-order accurate, in time, vertex update rules for the explicit forward Euler and the semi-implicit due to Desbrun et al. Starting from the continuous form of the MCF equation, we use a Taylor expansion to approximate the position of a point on the surface up to the second-order term, and for Desbrun et al.’s, we make use of Neumann series in our derivation.

Recall that from the first section, the MCF in its continuous form, is described following equation:

\[\frac{\partial X(u,t)}{\partial t} = \Delta X(u,t)\]

where: \( X(u,t) \in \mathbb{R}^3 \) is the position of a point on the surface at parameter \( u \) and time \( t \), and \(\Delta\) is the Laplace Beltrami operator.

Now, we apply a Taylor expansion of \( X(u,t) \) around time \( t=t_k\):

\[ X(u, t_k + \Delta t) = X(u, t_k) + \Delta t \frac{\partial X(u,t)}{\partial t}| _{t=t_k}+ \frac{\Delta t^2}{2} \frac{\partial^2 X(u,t)}{\partial t^2}| _{t=t_k} + \mathcal{O}(\Delta t^3) \]

Substituting the MCF equation \( \frac{\partial X}{\partial t} = \Delta X \): \[ X(u, t_k + \Delta t) = X(u, t_k) + \Delta t \Delta X(u,t) | _{t=t_k}+ \frac{\Delta t^2}{2} \frac{\partial}{\partial t} \left( \Delta X(u,t) \right)| _{t=t_k} + \mathcal{O}(\Delta t^3) \]

Where, \[\frac{\partial^2 X(u,t)}{\partial t^2}\Bigg |_{t=t_k} = \frac{\partial}{\partial t} \left( \Delta X(u,t) \right) \Bigg |_{t=t_k}\] This follows from the definition of Laplace Beltrami operator and Schwarz’s Theorem (also known as Clairaut’s Theorem on Equality of Mixed Partials)

Now, \[ \frac{\partial}{\partial t} \left( \nabla X(u,t) \right) \Bigg |_{t=t_k} = \Delta \left( \frac{\partial X(u,t)}{\partial t} \right) \Bigg |_{t=t_k}= \Delta \left( \Delta X(u,t) \right) \Bigg |_{t=t_k} = \Delta^2 X(u,t) \Bigg |_{t=t_k} \]

By substituting this in Taylor’s expansion, we get the continuous second-order expansion for \( X(u,t)\) at \(t_k\):

\[ X(u, t_k + \Delta t) = X(u, t_k) + \Delta t \Delta X(u,t) \Bigg |_{t=t_k} + \frac{\Delta t^2}{2} \Delta^2 X(u,t) \Bigg |_{t=t_k} + \mathcal{O}(\Delta t^3) (*) \]

Now let \( X^{(k)}_i\), \( X^{(k+1)}_i \) denote the position vector of vertex \(i\) at times \( t_k \), and\( t_{k+1} \) respectively, and \( \Delta t = t_{k+1}-t_k\). Using the Taylor expansion in \((*)\), the second-order vertex update rule becomes:

\[ X^{(k+1)}_i\approx X^{(k)}_i + \Delta t \Delta X^{(k)}_i + \frac{\Delta t^2}{2} \Delta^2 X^{(k)}_i \]

Now the spatial discrete approximation we use in this article is \(\Delta \approx \mathbf{ML}^{-1}\). Now we can write the vertex update rule for Forward Euler as follows:
\[ X^{(k+1)}_i\approx X^{(k)}_i + \Delta t \mathbf{ML}^{-1} X^{(k)}_i + \frac{\Delta t^2}{2} (\mathbf{ML}^{-1})^2 X^{(k)}_i \]

After doing some algebra, we reach the following matrix form:

\[ X^{(k+1)} \approx [ I + \Delta t \mathbf{ML}^{-1} + \frac{\Delta t^2}{2} (\mathbf{ML}^{-1})^2] X^{(k)} \]

For the semi-implicit form due to Desbrun et al.’s, things are a bit tricky. It should be easy by now to derive the first-order vertex update presented earlier in (ref) using Taylor expansion which takes the matrix form:

\[X^{k+1}_i \approx (I- \Delta t\mathbf{M^{-1}L} X^{k} )^{-1} X^{k} \]

To derive the second order term for this scheme, we expand the inverse matrix \( \left(I – \Delta t \mathbf{M}^{-1} L X^{(k)} \right)^{-1} \) using a Neumann series, but for this to work, we have to ensure that \( \mathbf{M}^{-1}L X^{(k)}\) satisfies the condition for convergence, i.e., its spectral radius is strictly less than 1: \[ \rho (\Delta t \mathbf{M}^{-1}L X^{(k)} )<1 \]

This means \(\Delta t\) must be chosen small enough, or the structure of \(\mathbf{M}^{-1}L X^{(k)}\) must ensure that its eigenvalues are small.

Assuming this is true, the Neumann series expansion for \( \left(I – \Delta t \mathbf{M}^{-1}L X^{(k)} \right)^{-1} \) can be written as: \[ \left(I – \Delta t \mathbf{M}^{-1} L X^{(k)}\right)^{-1} = I + \Delta t \mathbf{M}^{-1} L X^{(k)}+ \Delta t^2 \left(\mathbf{M}^{-1} L X^{(k)}\right)^2 + \mathcal{O}(\Delta t^3)\]

Substituting this approximation into the semi-implicit update rule, we get:

\[ X^{k+1} \approx \left( I + \Delta t \mathbf{M}^{-1} L X^{(k)} + \Delta t^2 \left( \mathbf{M}^{-1} L X^{(k)} \right)^2 \right) X^{k} \qquad (***)\]

Discussion. The advantage of using Neumann series in deriving the second-order time discretization is that it allows us to approximate \( \left(I – \Delta t \mathbf{M}^{-1} L X^{(k)} \right)^{-1} \) without having to directly compute the matrix inverse, which can be computationally expensive for large meshes. Instead, the expansion provides a series of manageable terms so with that we can economically exploit the accuracy benefits attained from adding the higher-order terms. With that said, the major disadvantage here is that if \(\Delta t\) becomes too large, the Neumann series may fail to converge or lead to unstable behavior, limiting a bit its effectiveness for semi-implicit schemes over larger intervals. However, it would be not correct to say that it is impossible to circumvent the stability issue. We talk about this in a subsequent article.

An Alternative Discretization based on (Huisken’s, 1984) MCF Evolution Equations

Earlier in the article, we mentioned the landmark result of the influence of MCF on strictly convex smooth hypersurface in Euclidean spaces due to (Huisken, 1984). To establish this result, Huisken derived several key equations that rigorously describe how various geometric quantities change over time as general surfaces evolve under MCF.

  1. \( \frac{\partial n}{\partial t}=\nabla H \)
  2. \( \frac{\partial H}{\partial t}=\Delta H+|A|^2 H =\Delta H + (H^2-2K)H. \) where \(A\) is the second fundamental form, and \(K\) is the Gaussian curvature.

These are called the surface evolution equations, The first equation describes how the unit normal vector \(n\) evolves over time, linking its rate of change to the gradient of the mean curvature \(H\). The second equation tracks the evolution of \(H\) itself as the surface changes.

From these, we can derive the following equation for the second derivative of the surface position \(X\) with respect to time:

\[\frac{\partial^2 X}{\partial t^2} = \frac{\partial}{\partial t}(H\mathbf{n}) = \frac{\partial H}{\partial t}\mathbf{n} + H\frac{\partial \mathbf{n}}{\partial t}= (\Delta H +(H^2-2K)H)\mathbf{n}+H\nabla H\]

This expression consists of geometric quantities that can be approximated on a mesh—though they tend to be noisy. We can also write the equation in an alternative form:

\[\frac{\partial^2 X}{\partial t^2}=(\Delta H)\mathbf{n} + (H^2-2K)\Delta x + H\nabla H\]

Why is this important here? Since we are discussing higher-order discretizations of MCF, we are interested in discovering new equivalent (and hopefully economical) ways to describe the temporal derivatives in question. (Huisken, 1984) provides some, and thus a natural question arises: can we discretize the components, \(\mathbf{n}, H, \Delta H, \nabla H, K\) of \( \frac{\partial^2 X}{\partial t^2}\) and \(\frac{\partial X}{\partial t}\) ​ to derive a second-order discretization for MCF using Taylor series? The answer is yes.

For example, the Gaussian curvature \(K\) can be discretized using the angle deficit method. The normal vector \(\mathbf{n}\) at a vertex can be estimated as the area-weighted average of the normals of the adjacent triangles. The mean curvature \(H\) as \(\mathbf{M^{-1}L}\) applied to the verticies of the mesh. The gradient \(\nabla H\) can be approximated using finite differences or based on the stiffness matrix \(mathbf{L}\) and adjacent vertex data, while \(\Delta H\) can be discretized using \(\mathbf{L}\) applied to the discrete mean curvature \(H\).

Using these discretized quantities, we arrive at the following vertex update formula: \[ X_i^{(k+1)} \approx X_i^{(k)}+\Delta t H_i \mathbf{n}_i + \frac{\Delta t^2}{2} ( (\mathbf{L} H_i) \mathbf{n}_i + (H_i^2-2K_i)\Delta X_i^{(k)} + H_i\nabla H_i)) \] However, as mentioned this discretization approach is often not preferred due to the significant noise in the quantities \(\mathbf{n}, H, \Delta H, \nabla H\), and \(K\) on meshes.

Visualizing \(n, H, \Delta H, \nabla H,\) and \(K\) on a Mesh

Visualizing mean curvature, and normals per vertices on a sphere mesh (r=0.5), with color mapping based on mean curvature.
Gradient of mean curvature H

The next experiment aims to visualize higher-order effects in MCF by plotting small arcs at each mesh vertex. These arcs are defined as

\(f(h) = X_i + (Hn) \big |_{X=X_i} \cdot h + \frac{1}{2} (\frac{\partial}{\partial t}(Hn))\big |_{X=X_i}= \frac{\partial H}{\partial t}n \big |_{X=X_i}+ H\frac{\partial n}{\partial t}\big |_{X=X_i}= (\Delta H +(H^2-2K)H)n \big |_{X=X_i}+H\nabla H) \big |_{X=X_i} \cdot h^2 \)

where \(\frac{\partial X}{\partial t}\) is the mean curvature normal and \(\frac{\partial^2X}{\partial t^2}\) is a second-order term from Huisken’s calculations.

The idea of this experiment was suggested by Prof. Justin Solomon on day two of the project!

The first-order term moves the surface in the direction of the normal, scaled by the mean curvature. This means regions with higher curvature see greater movement compared to flatter regions. The second-order term refines this by adding curvature-dependent corrections. It can enhance or counteract the displacement done via the first-order term, affecting the arc’s bending and potentially leading to different geometric changes. In addition, the second-order term can indeed add accuracy to the displacement, providing a more precise description of surface evolution. However, higher-order terms are also more sensitive to mesh noise and discretization errors, which can introduce potential instabilities or oscillations, particularly in regions with poor mesh quality and such instabilities can be amplified in regions with high curvature, where numerical errors from the second-order term might dominate. In our case, these instabilities are reflected in exaggerated displacements, resulting in disproportionately large polylines at certain vertices. With more trivial meshes, this instability problem will not be as amplified as it is the case with complex meshes.

Implementation of Second-Order Semi-implicit (Desbrun et al.’s, 1999):

The following are the output results of our second-order Desbrun et al’s semi-implicit method (\(***\)).

An animation showing the evolution of the surface under second-order Desbrun et. al’s MCF method, 50 iterations.

Yay or Nay: Circular Arc-Based Discretizations for Curvature-Driven Flows

In the Taylor expansion used to derive the vertex-update rule for MCF, the position \(X(u,t)\) of each vertex is typically approximated by a quadratic polynomial in time:

\[ X(u, t_k + \Delta t) = X(u, t_k) + \Delta t \frac{\partial X(u,t)}{\partial t}| _{t=t_k}+ \frac{\Delta t^2}{2} \frac{\partial^2 X(u,t)}{\partial t^2}| _{t=t_k} + \mathcal{O}(\Delta t^3) \]

where \( \frac{\partial X}{\partial t} = \Delta X \), the Laplacian of the position, is the driving term in MCF, and \( \frac{\partial^2 X}{\partial t^2} \) is obtained from differentiating this expression again. While this quadratic approximation is computationally straightforward and effective for small time steps, it does not inherently capture the geometric structure of the flow. 

For curvature-driven flows, such as the evolution of a sphere under MCF, where the curvature \( H \) is constant at each point, circular arcs may provide a more natural approximation. Circular arcs reflect the constant curvature evolution by following a trajectory where the velocity of each vertex aligns with the normal direction, and the path of the vertex forms part of a circle. This would involve approximating the update as:

\[X(u, t_k + \Delta t) \approx X(u, t_k) + r(\cos(\theta) – 1) \mathbf{n},\]

where \( r \) is the radius of curvature and \( \theta \) is the angle swept by the vertex in time \( \Delta t \), with \( \mathbf{n} \) being the surface normal. While circular arcs introduce more computational complexity, they better approximate the geometric behavior of curvature-dominated flows and may lead to improved accuracy and stability in such cases.

Key Takeaway:

  1. MCF is important in GP!
  2. Lots of discretization approaches for the Laplacian exist, none of them could keep every natural property of its ideal continuous form. You choose what is suitable for your problem, and application.
  3. Coming up with new equivalent formulations for MC \(H\), and the rate of change of the position vector-valued function \(X\) of points on the surface would open more doors for finding new economical discretizations

Future work: Will venture more into the math of MCF, focusing specifically on points 2 and 3 from the Key Takeaways. Additionally, I explore some tangentials in regards to higher-order integrators for MCF, and other geometric flows.


A Humorous Fail:

On the second day while coding the first-order discretizations, I forgot to include the mass matrix \(\mathbf{M}\) which resulted in a smoothed horribly deformed bear. This demonstrates the critical role of \(\mathbf{M}\) in ensuring that the discretization respects the surface’s geometry by appropriately distributing weight across the vertices according to the areas of the surrounding triangles.

  1. When the basis functions used are piecewise linear and the mesh structure is uniform ↩︎

Bibiliography:

  • Justin Solomon (Director). (2013, May 8). Lecture 12: Finite Elements and the Laplacian [Video recording]. https://www.youtube.com/watch?v=7_xDIg-pOC4
  • Huisken, G. (1984). Flow by mean curvature of convex surfaces into spheres. Journal of Differential Geometry20(1), 237-266.
  • Desbrun, M., Meyer, M., Schröder, P., & Barr, A. H. (1999, July). Implicit fairing of irregular meshes using diffusion and curvature flow. In Proceedings of the 26th annual conference on Computer graphics and interactive techniques (pp. 317-324).
  • Patanè, G. (2017). An introduction to Laplacian spectral distances and kernels: Theory, computation, and applications. In ACM SIGGRAPH 2017 Courses (pp. 1-54).
  • Hughes, T. J. R. (2000). The finite element method: Linear static and dynamic finite element analysis. Dover Publications.
  • Evans, L. C. (2010). Partial differential equations. American Mathematical Society.

Categories
Research

Global Intersection Analysis

Introduction

Simulations often become unstable as a result of self-intersection or intersection between two meshes. This instability can lead to wrong simulation results and incorrect outputs. At the same time, self and inter collisions of meshes are a necessity for artistic purposes in cases such as 3D animation. To resolve this issue, Baraff et al. (2003) reveals a method known as global intersection analysis, which specifies how to resolve these mesh intersections while allowing the simulation to run as it normally would. This ensures that we obtain simulation stability while also allowing for mesh intersections when required.

The goal

For the sake of our one week project time period, our project focussed on creating a simple implementation of the global intersection analysis algorithm along with a simple resolution method between two different meshes.

The global intersection analysis algorithm required a few steps:

  1. Identify a contour along which both meshes intersect each other
  2. Identify the points inside and outside this contour for both mesh
  3. Use these points to resolve the intersections by adding a constraint that “pulls” them out of each other (subsequently resolving the intersection)
  4. Allow the physics solver to run as it normally would after resolution of intersection.

Our entire implementation can be found here.

Development

We implemented this project in 3 parts.

  1. The global intersection analysis algorithm, which would perform the steps above
  2. Our physics solver, which would run an XPBD algorithm for the cloth to follow
  3. Integrating both these systems together in terms of the constraints needed for the cloth XPBD to take intersection resolution into account.

For our demo scene, we used a simple rigid ball and a tesselated plane to act as a cloth mesh.

Both test meshes with their intersection contour

Contour Identification and Flood Filling

For our GIA implementation, we created a scene which identifies a contour, identifies the points inside and outside and provides this data as output for further processing.

Points identified as inside and outside on the cloth mesh
Points identified as inside and outside the contour along the ball mesh

Our identification of points is done through a flood fill algorithm.

The original method mentioned by Baraff et al. (2003) discusses selecting the region with lesser area, determining that to be the area inside the contour and then flood filling from inside this smaller area and larger area with two separate colors to identify both regions on the mesh. It isn’t entirely clear how to choose a point in this area or how to determine this area solely from the contour. As a result, we created and implemented our own method for the flood fill which works well in our current test cases.

Our method follows standard flood filling techniques, where we identify a point on the surface with one color, then recursively repeat the coloring process to all points connected to the last one. To identify if a point is across the contour, we change the color if the edge that connects two points is identified as one that intersects with the edges on the contour.

To make this more efficient, we identify contour intersecting edges beforehand. We also colored any points along the contour beforehand to ensure the flood fill would not overwrite or mistake those points. As a heuristic, we also begin the flood fill along any identified points along the contour since pre-coloring a point means the flood fill algorithm will not continue if all points adjacent to a point are colored.

We utilised the IPC Toolkit [3] to perform intersection operations between the surface and contour edges.

XPBD Implementation

Using Matthias Müller’s implementation[2] as reference, we also implemented an XPBD cloth simulator.

Our cloth dangling by a few points on it’s top edges

We provide various parameters to change as well, such as the simulation time step and cloth properties.

Integration

To integrate both systems, the most important part was integrating the GIA collision resolution into the cloth mesh’s XPBD constraint formulation. For simplicity, we push away the cloth in a force whose direction is the vector from center of the cloth contour points and the ball mesh contour points. This way, we would have a close-enough approximation to push the cloth away. We then integrate this constraint into the system, where we identify and assign this force direction for the cloth after performing a GIA call. One GIA call would perform a contouring, vertex identification and force direction identification before assigning it to the cloth’s constraint. The cloth’s XPBD solver would then take this into account when determining the next position of points on the cloth.

Result

The results of our integration can be seen here.

The UI also allows user to play around with different simulation settings and options.

To get a look into the source code, our GitHub repository can be found here.

Scope for improvement

Given our 1 week period, there are some areas for improvement.

  1. Our current GIA implementation is too slow to be used in realtime and can be made faster.
  2. Our intersection resolution method is not what is exactly outlined in the original “Untangling cloth” paper but a rough approximation. Baraff et al. (2003) provide a much more precise method to determine per vertex resolution along the meshes.
  3. XPBD may not be the most accurate simulator of cloth (though for an application such as animation, it works well enough visually).
  4. Our GIA flood fill can sometimes identify vertices incorrectly, leading to inaccuracies in the mesh-intersection resolution.

References

[1] Baraff, D., Witkin, A., & Kass, M. (2003). Untangling cloth. ACM Transactions on Graphics, 22(3), 862–870. https://doi.org/10.1145/882262.882357

[2] Matthias Müller. pages/tenMinutePhysics/15-selfCollision.html at master · matthias-research/pages. GitHub. https://github.com/matthias-research/pages/blob/master/tenMinutePhysics/15-selfCollision.html

[3] IPC Toolkit. https://ipctk.xyz/

This blog was written by Sachin Kishan as one of the outcomes of a project during the SGI 2024 Fellowship under the mentorship of Zachary Ferguson.

Categories
Research

Arc-Length Splines

Introduction

A spline is a function that usually represents a piecewise polynomial. Splines have a variety of applications, from being used to visualize a curve to representing motion in animation, to model a surface, or for artistic visualizations. 

Our project on Arc Length Splines aimed to satisfy a new necessity in spline formulation- an analytically computable arc length. We aim to allow users to output the arc length and further constrain the spline using the arc length. Further, we want to ensure some amount of continuity on the spline.

Background

Continuity

In the context of splines, continuity can be defined as having no sudden change in value across the spline function. There are two kinds of continuity discussed here, namely parametric (C) and geometric (G) continuity. 

There are different levels of continuity dubbed C1 continuity, C2 continuity, and so on until Cn continuity. (The same applies to G as well)

A Ci continuous spline is defined as one where for a function f(t) that defines the spline, f i (t) is continuous for all points t, where i represents the order of the derivative. This is usually true for each curve that is in the spline. The points of interest where continuity may not occur are at the points where two curves meet, called the joints. To verify if there is continuity at the joints, we can use the following method: 

Assume we have two curves A and B where the endpoint of A is connected to the starting point of B. These curves are parameterized on the domain [0,1].

A and B form a Ci continuous spline if 

Ai (1) = Bi (0)

This can be used to verify any spline made of an arbitrary number of curves, given that for every two consecutive curves connected at a joint, they satisfy the above condition. 

Any Ci spline is also C0, C1 … Ci-2, Ci-1continuous.

Gi continuity is based on the “geometric” smoothness of the curve. Where G1 is the continuity at the tangents along the spline. G2 continuity refers to the continuty of the curvature along the spline. 

Continuity plays an important role in how a spline’s application may behave. For example, if a spline is acting as a path for an object moving along it, if the spline is not C1 continuous (meaning each consecutive curve has an endpoint that meets at the joint of the spline but with discontinuities for the first differential), the spline may not be useful. This is because objects would likely move at different speeds along different curves(or increase and decrease in speed at the joints.) This is what makes continuity a desirable property in spline applications. 

Approximating vs interpolating control points

Splines can be changed based on their control points. Splines that go through the entire set of control points are called interpolating splines (since they interpolate between the entire given set of control points). In an approximating spline, the spline approximates the polyline made by the control points into a smooth curve. 

Local and Global control

Based on the constraints on the spline formulation, the control points have influence over certain parts of the spline. In cases where a control point only effects the curve at which it acts as a joint (or is the point that directly decides how two curve segments are to be formed), the spline is said to have control points with local control. (Local being the curves before and after the control point at a joint). 

In some cases, constraints may influence larger parts of the spline beyond local curves, this is referred to as global control. 

It is usually desirable to have local control for applications like spline drawing, where we want to create a curve of a certain visual form without largely impacting the rest of the spline it is a part of.

Arc Length

Arc length is the distance along the curve given two points on the curve. For any arbitrary curves, arc length can be integrated between two points. For special curves such as circles or parabolas, there is an analytically computed arc length. Our choice of curve and interpolation schemes to explore were driven by where we could analytically calculate arc length for a curve. This also influenced other considerations of our spline.

Past Work

Most of our inspiration came from the paper “A Class of C^2 Interpolating Splines” by Cem Yuksel(2020). Yuksel was able to classify four different types of splines that maintained C^2 continuity and interpolation. These were splines using the Bézier interpolation function, circular interpolation function, elliptical interpolation function, and hybrid (circular-elliptical) interpolation function. Yuksel(2020) used a combination of the base curve and a blending function to maintain continuity. 

The use of a blending function is excellent for ensuring continuity between different kinds of curves. However, our goal is to have an exact arc length. The function Yuksel(2020) uses was a second order trigonometric function. Thus, there was no closed formula for the arc-length of for the blending function.

Our work

After initial review of past literature- we compiled a set of curves in which an analytical solution for arc length exists. We dived deeper into curves which have closed formula arc length. The main curves of this family are circles, parabolas, catenoids, and the logarithmic spiral. We then focussed on creating splines out of this initial set of curves. Through this process, we recognised that circular curves may be the best option, not only because their arc length is easy to calculate, but presented properties to create a locally controllable spline with C1 continuity.

A spline is formed when we connect an endpoint of two curves to each other. To formulate a circular spline, there are different ways to connect them. One option was to just draw line segments between two arcs. The problem is that this may not always be continuous if further constraints are not provided.

Another way to connect two circular arcs is to use the Dubins Path. Dubins Paths are usually used in robotics or car path planning. Since they are used to find the shortest path between any two curves which have constraints on curvature, they work perfectly to find a connecting path between two circular arcs.

Dubin’s path between two circular arcs, which turns out to be a line segment. Source: Salix Alba [2]

If a line segment between two curves is not a line segment, Dubin’s path instead form it’s own curve to draw between the two. These are known as CCC paths (where C stands for curve).

Dubin’s path between two circular arcs, which turns out to be another circular arc. Source: Salix Alba [3]

This ensures we no longer constrain user outcome and provide additional output for a user to adjust visually. The arc length of a Dubins path can also be computed analytically. This is because in all possible cases of a Dubins path, the path will either be a line segment or another circular arc. This allows for the analytical calculation of arc length across the entire spline while also ensuring C1 continuity. This fits with our requirements for an arc-length computable spline retaining at least a C1 continuity property.

We then worked on implementing the idea above: a set of circular arcs whose spline is formed with the dubins path connecting them. 

To implement circular arcs, we used two points and a tangent from one of them to determine a circular arc. This would equate to 3 points used per arc. Using a set of 6 points (which form 2 circular arcs), we then create a path between them. 

Our current implementation allows for the creation of circular arcs with straight paths formed between them.

Our implementation in python. Each arc is made of three points and is connected by a line segment

Future work 

Our current implementation is limited, we intend to improve it in the following ways:

  1. Ensure a Dubins path for all cases of circular arcs
  2. Create arc length constraints as well as outputs for users to view and edit
  3. Better interface to clearly demarcate between different curves, tangent points and circular arc points.
  4. Implementing higher order continuity while ensuring lesser constraints on spline editing.

References

[1] Cem Yuksel. 2020. A Class of C2 Interpolating Splines. ACM Trans. Graph. 39, 5, Article 160 (October 2020), 14 pages. https://doi.org/10.1145/3400301

[2] Salix Alba (2016, 6 February) Dubins3.svg https://upload.wikimedia.org/wikipedia/commons/e/e7/Dubins3.svg

[3] Salix Alba (2016, 6 February) Dubins2.svg https://upload.wikimedia.org/wikipedia/commons/e/e7/Dubins2.svg

[4] Freya Holmér. (2022, December 7). The continuity of Splines [Video]. YouTube. https://www.youtube.com/watch?v=jvPPXbo87ds

This blog was written by Sachin Kishan and Brittney Fahnestock during the SGI 2024 Fellowship as one of the outcomes of a project under the mentorship of Sofia Wyetzner and the support of Shanthika Naik as teaching assistant.

Categories
Research

Adaptive Isotropic Remeshing for TetSphere Splatting

Primary mentor: Minghao Guo

Volunteer assistant: Shanthika Naik

Student: Dianlun (Jennifer) Luo

1. Introduction

TetSphere Splatting, introduced by Guo et al. (2024), is a cutting-edge method for reconstructing 3D shapes with high-quality geometry. It employs an explicit, Lagrangian geometry representation to efficiently create high-quality meshes. Unlike conventional object reconstruction techniques, which typically use Eulerian representations and face challenges with high computational demands and suboptimal mesh quality, TetSphere Splatting leverages tetrahedral meshes, providing superior results without the need for neural networks or post-processing. This method is robust and versatile, making it suitable for a wide range of applications, including single-view 3D reconstruction and image-/text-to-3D content generation.

The primary objective of this project was to explore volumetric geometry reconstruction through the implementation of two key enhancements:

  1. Geometry Optimization via Subdivision/Remeshing: This involved integrating subdivision and remeshing techniques to refine the geometry optimization process in TetSphere Splatting, allowing for the capture of finer details and improving overall tessellation quality.
  2. Adaptive TetSphere Modification: Mechanisms were developed to adaptively split and merge tetrahedral spheres during the optimization phase, enhancing the flexibility of the method and improving geometric detail capture.

In this project, I implemented Adaptive Isotropic Remeshing from scratch. The algorithm was tested on two models:

  • final_surface_mesh.obj generated by TetSphere Splatting in the TetSphere Splatting codebase.
  • The Stanford Bunny mesh (bun_zipper_res2.ply).

2. Adaptive Isotropic Remeshing Algorithm

The remeshing procedure can be broken down into four major steps:

  1. Split all edges at their midpoint that are longer than \( \frac{4}{3} l \), where \( l \) is the local sizing field.
  2. Collapse all edges shorter than \( \frac{4}{5} l \) into their midpoint.
  3. Flip edges in order to minimize the deviation from the ideal vertex valence of 6 (or 4 on boundaries).
  4. Relocate vertices on the surface by applying tangential smoothing.
2.1 Adaptive Sizing Field

Instead of using a uniform target edge length, an adaptive sizing field was computed using the trimesh package. The constant target edge length \( L \) was replaced by an adaptive sizing field \( L(x) \), which is intuitive to control, simple to implement, and efficient to compute.

Splitting edges based on their length and the angle of their endpoint normals can lead to anisotropically stretched triangles, requiring an additional threshold parameter to control the allowed deviation of endpoint normals. The remeshing approach is based on a single, highly intuitive parameter: the approximation tolerance \( \epsilon \). This parameter controls the maximum allowed geometric deviation between the triangle mesh and the underlying smooth surface geometry.

The method first computes the curvature field of the input mesh and then derives the optimal local edge lengths (the sizing field \( L(x) \)) from the maximum curvature and the approximation tolerance. The process is as follows:

Algorithm: Compute Adaptive Sizing Field

Input: Mesh M, Points P, Parameter ε
Output: Sizing Field Values S

Step 1: Compute Mean Curvature
    H ← discrete_mean_curvature_measure(M, P, radius = 1.0)

Step 2: Compute Gaussian Curvature
    K ← discrete_gaussian_curvature_measure(M, P, radius = 1.0)

Step 3: Compute Maximum Curvature
    C_max ← |H| + sqrt(|H² - K|)

Step 4: Ensure Minimum Curvature Threshold
    C_max ← max(C_max, 1e-6)

Step 5: Compute Sizing Field Values
    S ← sqrt(max((6ε / C_max) - 3ε², 0))

Step 6: Return Sizing Field Values
    Return S
2.2 Edge Splitting

Edge splitting is a process where edges that exceed a certain length threshold are split at their midpoint. This step helps maintain consistent element sizes across the mesh and prevent the occurrence of overly stretched or irregular elements. By splitting longer edges, we can achieve better mesh quality, improve computational efficiency, and optimize the mesh for various applications.

Algorithm: Split Edges of a Mesh

Input: Mesh M, Parameter ε
Output: Modified Mesh with Split Edges

Step 1: Initialize new vertex and face lists
    new_vertices ← M.vertices
    new_faces ← M.faces

Step 2: Iterate over each edge in the mesh
    For each edge e = (v₀, v₁) in M.edges:
        - Compute edge length: l_edge ← ||M.vertices[v₀] - M.vertices[v₁]||
        - Compute sizing field for edge vertices: sizing_values ← sizing_field(M, [v₀, v₁], ε)
        - Calculate target length: l_target ← (sizing_values[0] + sizing_values[1]) / 2

        If l_edge > (4/3) * l_target:
            Step 3: Split the edge
                - Compute midpoint: midpoint ← (M.vertices[v₀] + M.vertices[v₁]) / 2
                - Add the new vertex to new_vertices
                - Find adjacent faces containing edge (v₀, v₁)
                - For each adjacent face f = (v₀, v₁, v₂):
                    - Remove face f from new_faces
                    - Add new faces [v₀, new_vertex, v₂] and [v₁, new_vertex, v₂] to new_faces

Step 4: Create and return new mesh
    Return M ← trimesh.Trimesh(vertices=new_vertices, faces=new_faces)
2.3 Edge Collapse

Edge collapse simplifies the mesh by merging the vertices of short edges. This step is particularly useful for eliminating unnecessary small elements that may introduce inefficiencies in the mesh. By collapsing short edges, the mesh is coarsened in areas where high detail is not needed.

Algorithm: Collapse Edges of a Mesh

Input: Mesh M, Parameter ε
Output: Modified Mesh with Collapsed Edges

Step 1: Initialize face list and sets
    new_faces ← M.faces
    vertices_to_remove ← []
    collapsed_vertices ← {}

Step 2: Iterate over each edge in the mesh
    For each edge e = (v₀, v₁) in M.edges:
        - Compute edge length: l_edge ← ||M.vertices[v₀] - M.vertices[v₁]||
        - Compute sizing field for edge vertices: sizing_values ← sizing_field(M, [v₀, v₁], ε)
        - Calculate target length: l_target ← (sizing_values[0] + sizing_values[1]) / 2

        If l_edge < (4/5) * l_target and v₀, v₁ not in collapsed_vertices:
            Step 3: Collapse the edge
                - Compute midpoint: midpoint ← (M.vertices[v₀] + M.vertices[v₁]) / 2
                - Set M.vertices[v₀] ← midpoint
                - Mark v₁ for removal: vertices_to_remove ← vertices_to_remove ∪ {v₁}
                - Mark v₀, v₁ as collapsed: collapsed_vertices ← collapsed_vertices ∪ {v₀, v₁}

                Find adjacent faces containing v₁ but not v₀:
                For each adjacent face f = (v₁, v₂, v₃):
                    - Remove face f from new_faces
                    - Add new face [v₀, v₂, v₃] to new_faces
                    - Mark v₂, v₃ as collapsed

Step 4: Update the mesh
    - Create index_map to remap vertex indices for remaining vertices
    - Remove vertices in vertices_to_remove from M.vertices
    - Remap the faces in new_faces based on index_map

Step 5: Create and return new mesh
    Return M ← trimesh.Trimesh(vertices=new_vertices, faces=new_faces)
    
2.4 Edge Flipping

Edge flipping is performed when it reduces the squared difference between the actual valence of the four vertices in the two adjacent triangles and their optimal values. For interior vertices, the ideal valence is 6, while for boundary vertices, it is 4.

To preserve key features of the mesh, sharp edges and material boundaries should not be flipped.

Algorithm: Check Edge Flip Condition

Input: Vertices V, Edges E, Faces F, Edge e = (v₁, v₂)
Output: Flip Condition, New Vertices v₃, v₄, Adjacent Faces f₁, f₂

Step 1: Find adjacent faces of edge e
    adjacent_faces ← [f ∈ F : v₁ ∈ f and v₂ ∈ f]

    If len(adjacent_faces) ≠ 2:
        Return False, None, None, None, None

Step 2: Identify third vertices of adjacent faces
    f₁, f₂ ← adjacent_faces
    v₃ ← [v ∈ f₁ : v ≠ v₁ and v ≠ v₂]
    v₄ ← [v ∈ f₂ : v ≠ v₁ and v ≠ v₂]

    If len(v₃) ≠ 1 or len(v₄) ≠ 1 or v₃[0] = v₄[0]:
        Return False, None, None, None, None

Step 3: Check edge conditions
    v₃, v₄ ← v₃[0], v₄[0]
    
    If [v₃, v₄] ∈ E:
        Return False, None, None, None, None

Step 4: Check normal vector alignment
    Compute normal₁ ← normal(f₁), normal₂ ← normal(f₂)
    Normalize both normal₁ and normal₂
    Compute dot product dot_product ← dot(normal₁, normal₂)

    If dot_product < 0.9:
        Return False, None, None, None, None

Step 5: Compute valence difference before and after edge flip
    valence ← compute_valence(V, E)
    Compute valence difference before and after the flip

Return difference_after < difference_before, v₃, v₄, f₁, f₂
Algorithm: Flip Edge of a Mesh

Input: Mesh M, Edge e = (v₁, v₂)
Output: Modified Mesh M

Step 1: Check if edge flip condition is satisfied
    condition, v₃, v₄, f₁, f₂ ← flip_edge_condition(M.vertices, M.edges_unique, M.faces, e)

    If condition = False:
        Return M

Step 2: Update the faces
    Remove f₁ and f₂ from the mesh faces
    Add new faces [v₁, v₃, v₄] and [v₂, v₃, v₄] to the mesh

Step 3: Create and return the updated mesh
    Return new_mesh ← trimesh.Trimesh(vertices=M.vertices, faces=new_faces)
2.5 Tangential Relaxation

Tangential smoothing relocates vertices to create a smoother surface without significantly altering the mesh geometry. The process averages the positions of the neighboring vertices, constrained to the surface, which maintains the mesh’s adherence to the original geometry while improving its smoothness.

Algorithm: Tangential Relaxation

Input: Mesh M, Sizing Field L
Output: Relaxed Mesh

Step 1: Initialize variables
    vertices ← M.vertices
    faces ← M.faces
    relaxed_vertices ← zeros_like(vertices)

Step 2: For each vertex i in vertices:
    Find incident faces for vertex i
        incident_faces ← where(faces == i)
        weighted_barycenter_sum ← zeros(3)
        weight_sum ← 0

    Step 3: For each face f in incident_faces:
        Compute barycenter and area
            face_vertices ← faces[f]
            triangle_vertices ← vertices[face_vertices]
            barycenter ← mean(triangle_vertices)
            
            v₀, v₁, v₂ ← triangle_vertices
            area ← ||cross(v₁ - v₀, v₂ - v₀)|| / 2

        Step 4: Compute weight based on sizing field
            sizing_at_barycenter ← mean([L[v] for v in face_vertices])
            weight ← area × sizing_at_barycenter

        Step 5: Accumulate weighted barycenter sum
            weighted_barycenter_sum ← weighted_barycenter_sum + weight × barycenter
            weight_sum ← weight_sum + weight

    Step 6: Compute relaxed vertex position
        relaxed_vertices[i] ← weighted_barycenter_sum / weight_sum if weight_sum ≠ 0 else vertices[i]

Step 7: Create and return relaxed mesh
    Return relaxed_mesh ← trimesh.Trimesh(vertices=relaxed_vertices, faces=faces)
2.6 Adaptive Remeshing

The adaptive remeshing algorithm iterates through the above steps to refine the mesh over multiple iterations.

Algorithm: Adaptive Remeshing

Input: Mesh M, Parameter ε, Number of Iterations: iteration
Output: Remeshed and Smoothed Mesh

For i = 1 to iteration:

    Step 1: Split edges
        M ← split_edges(M, ε)  

    Step 2: Collapse edges
        M ← collapse_edges(M, ε)

    Step 3: Flip edges
        For each edge e in M.edges_unique:
            M ← flip_edge(M, e)

    Step 4: Tangential relaxation
        Recalculate sizing values based on the updated mesh
        sizing_values ← sizing_field(M, M.vertices, ε)
        M ← tangential_relaxation(M, sizing_values)

Return M

3. Results

3.1 final_surface_mesh.obj

The TetSphere Splatting output (final_surface_mesh.obj) and its corresponding remeshed version (dog_remeshed_iter1.obj) were evaluated to assess the effectiveness of the remeshing algorithm.

  • final_surface_mesh.obj:
    • Vertices: 16,675
    • Faces: 33,266
  • dog_remeshed_iter1.obj:
    • Vertices: 26,603
    • Faces: 53,122

Comparison:
The remeshing process results in a significant increase in both vertices and faces, which produces a more refined and detailed mesh. This is clear in the comparison image below. The dog_remeshed_iter1.obj (right) shows a denser mesh structure compared to the final_surface_mesh.obj (left). The Adaptive Isotropic Remeshing algorithm enhances resolution and captures finer geometric details.

final_surface_mesh.obj
dog_remeshed_iter1.obj
3.2 bun_zipper_res2.ply

The Stanford Bunny mesh (bun_zipper_res2.ply) was similarly processed through multiple iterations of remeshing to evaluate the progressive refinement of the mesh.

  • bun_zipper_res2.ply:
    • Vertices: 8,171
    • Faces: 16,301
  • bunny_remeshed_iter1.obj:
    • Vertices: 6,476
    • Faces: 12,962
  • bunny_remeshed_iter2.obj:
    • Vertices: 5,196
    • Faces: 10,381
  • bunny_remeshed_iter3.obj:
    • Vertices: 4,159
    • Faces: 8,300

Comparison:
The progressive reduction in vertices and faces over the iterations demonstrates the remeshing algorithm’s ability to simplify the mesh while retaining the overall geometry. The remeshed iterations display a higher number of equilateral triangles, which creates a more uniform and well-proportioned mesh. This improvement is obvious when comparing the final iteration (bunny_remeshed_iter3.obj) with the original model (bun_zipper_res2.ply). The triangles become more equilateral and evenly distributed, resulting in a smoother and more consistent surface.

bun_zipper_res2.ply
bunny_remeshed_iter1.obj
bunny_remeshed_iter2.obj
bunny_remeshed_iter3.obj

4. Performance Analysis

The current implementation of the remeshing algorithm was developed in Python, utilizing the trimesh library. While Python is easy to use, performance limitations may arise when handling large-scale meshes or real-time rendering.

A potential solution to improve performance is to transition the core remeshing algorithms to C++.

Benefits of Using C++:
  • Speed: C++ allows for significantly faster execution times due to its lower-level memory management and optimization capabilities.
  • Parallelization: Advanced threading and parallel computing techniques in C++ can accelerate the remeshing process.
  • Memory Efficiency: C++ provides better control over memory allocation, which is crucial when working with large datasets.

5. References

M. Botsch and L. Kobbelt, “A remeshing approach to multiresolution modeling,” Proceedings of the 2004 Eurographics/ACM SIGGRAPH Symposium on Geometry Processing, Nice, France, 2004, pp. 185–192. doi: 10.1145/1057432.1057457.

M. Dunyach, D. Vanderhaeghe, L. Barthe, and M. Botsch, “Adaptive remeshing for real-time mesh deformation,” in Eurographics 2013 – Short Papers, M.-A. Otaduy and O. Sorkine, Eds. The Eurographics Association, 2013, pp. 29–32. doi: 10.2312/conf/EG2013/short/029-032.

“Remeshing,” Lecture Slides for CS468, Stanford University, 2012. [Online]. Available: https://graphics.stanford.edu/courses/cs468-12-spring/LectureSlides/13_Remeshing1.pdf. [Accessed: Sep. 12, 2024].

M. Guo, B. Wang, K. He, and W. Matusik, “TetSphere Splatting: Representing high-quality geometry with Lagrangian volumetric meshes,” arXiv, 2024. [Online]. Available: https://arxiv.org/abs/2405.20283. [Accessed: Sep. 12, 2024].

B. Kerbl, G. Kopanas, T. Leimkuehler, and G. Drettakis, “3D Gaussian splatting for real-time radiance field rendering,” ACM Trans. Graph., vol. 42, no. 4, Art. no. 139, Jul. 2023. doi: 10.1145/3592433.

Categories
Math Research

3D Chladni Patterns

Introduction

Chladni patterns are usually created by putting some light, scattered object like sand onto a metal plate. The metal plate is then made to vibrate, which forms different patterns on the plate depending on the frequency of the wave.

Skrodzki et al. (2016) introduce a method to bring Chladni patterns into the third dimension.

Chladni patterns represent the points along which multiple waves meet to form nodes. These nodes are points along the standing wave formed by the combination of waves where a particle has 0 displacement from its mean.

Depending on the boundary condition, the final solution for the Chladni formulation varies. We can choose between Dirichlet or Neumann conditions.

\(A · sin(u · π · x) · sin(v · π · y) · sin(w · π · z) + B · sin(u · π · x) · sin(v · π · z) · sin(w · π · y) \)
\(+ C · sin(u · π · y) · sin(v · π · x) · sin(w · π · z) + D · sin(u · π · y) · sin(v · π · z) · sin(w · π · x) \)
\(+E · sin(u · π · z) · sin(v · π · x) · sin(w · π · y) + F · sin(u · π · z) · sin(v · π · y) · sin(w · π · x) \)

The above is the solution given a Dirichlet boundary condition. To get the solution with a Nuemann boundary condition, it is the same as the above solution, where all the sin functions are cos instead. More details as to the use of amplitudes and wave number are discussed by Skrodzki et al. (2016).

Using the solution, we can use it as an implicit surface for rendering. This can be done using a standard cube marching algorithm.

Outcome

Through the above formulation and cube marching techniques, our group created two open source web versions. A shadertoy implementation as well as a 3js implementation.

Source code and weblinks to both implementations can be seen here.

Image 1: A screenshot of the rendering along with some of the options a user can modify

References

[1] Skrodzki, M., Reitebuch, U., & Polthier, K. (2016). Chladni Figures Revisited: A peek into the third dimension. Proceedings of Bridges 2016: Mathematics, Music, Art, Architecture, Education, Culture, 481–484. http://www.archive.bridgesmathart.org/2016/bridges2016-481.html

This blog was written by Sachin Kishan, Nicolas Pigadas and Bethlehem Tassew during the SGI 2024 Fellowship as one of the outcomes of a one week project under the mentorship of Martin Skrodzki and support of Alberto Tono as teaching assistant.

Categories
Research

Using Bicubic Curves for a Differentiable Representation of 2D Curve Networks

Project Mentor: Mikhail Bessmeltsev

SGI Volunteer: Daniel Perazzo

Raster vs Vector Images

Raster and vector graphics are two fundamental types of image representations used in digital design, each with distinct characteristics and applications. A raster image is a collection of colored pixels in a bitmap, making it well-suited for capturing detailed, complex images like photographs. However, raster graphics provide no additional information about the structure of the object they represent. This limitation can be problematic in scenarios where precision and scalability are essential. For example, representing an artist’s sketch as a raster image may result in jagged curves and the loss of fine details from the original sketch. In such cases, a vector graphics representation is preferred. Unlike raster images, vector graphics use mathematical equations to define shapes, lines, and curves, allowing unlimited scalability without losing quality. This makes vector graphics ideal for logos, illustrations, and other designs where clarity and precision are crucial, regardless of size.

Bitmap and vector images on magnifying
Image from https://en.wikipedia.org/wiki/Image_tracing

Bicubic Curves for Vectorization

In this project, we aimed to obtain a vector graphics representation of a given sketch in raster form using bicubic curves. A bicubic curve is a polynomial of the form

{\displaystyle p(x,y)=\sum \limits _{i=0}^{3}\sum _{j=0}^{3}a_{ij}x^{i}y^{j}.}

The value a bicubic function takes on a given grid cell is assumed to be uniquely determined by the value the function and its derivatives take at the vertices of the cell. Using bicubic curves allows us to more easily handle self-intersections and closed curves. Additionally, we can create discontinuities at the edges of grid cells by adding some noise to the given bicubic function. This will enable us to vectorize not only smooth but also discontinuous sketches.

A bicubic curve with a self-intersection at (0, 0)
Adding discontinuities at the grid cell edges

Preliminary Work

During the first week of the project, we implemented the following using Python:

  1. For any given grid cell, given its size and the function value p(v) as well as the derivatives px(v), py(v), and pxy(v) at each vertex v, we reconstructed and plotted the unique bicubic function satisfying these values. This reduces to solving a system of sixteen linear equations, following the steps mentioned in the “Bicubic Interpolation” Wikipedia page {4}.
  2. We added noise to the bicubic functions in each cell to create discontinuities. We also found the endpoints of curves in the grid due to the discontinuities formed above.
  3. At each point in the grid, we found the gradient of the bicubic functions using fsolve {5} and used them to plot the tangents at some points on the curves.
  4. We looked for self-intersections within a grid cell by finding saddle points of the bicubic function (such that the function value is 0 there). If the function is factorizable, it suffices to find the points of intersections between the 0-levels of two or more factors. Here too, we used fsolve.
Bicubic curves with tangents produced by picking four random values between -1 and 1 at each vertex

Optimization and Rendering

Vectorization

During the second week of the project, I attempted to vectorize the bitmap image by minimizing a function known as the L2 loss.

The L2 loss between the predicted and target bitmap images is computed as the sum of the squared differences for each pixel. To account for this, I gave the bicubic curves some thickness by replacing the curve Z = 0 (here Z = p(x, y) is a bicubic function) in each grid cell with the image of the function e-Z². This allowed me to determine the color at each point in the grid, defining the color of each pixel as the color at its center. Additionally, I calculated the gradient of the color at each point in a grid cell with respect to the coordinates of the corresponding bicubic curve. Using this approach, I defined the L2 loss function and its gradient as follows:

# Defining the loss function
def L2_loss(coeffs, d, width, height, x, y):
    c = np.zeros((height, width))
    for i in range(height):
        for j in range(width):
            c[i, j] = np.exp(-bicubic_function(x[j], y[i], coeffs) ** 2)
    return np.sum((c - d) ** 2)

# Derivative of c(x, y) wrt to the coeff of the x ** m * y ** n term
def colour_gradient(x, y, coeffs, m, n):
    return np.exp(-bicubic_function(x, y, coeffs) ** 2) * (-2 * bicubic_function(x, y, coeffs) * x ** m * y ** n)

# Gradient of the loss function
def gradient(coeffs, d, width, height, x, y):
    grad = np.zeros(16,)
    c = np.zeros((height, width))
    for i in range(height):
        for j in range(width):
            c[i, j] = np.exp(-bicubic_function(x[j], y[i], coeffs) ** 2)
            for n in range(4):
                for m in range(4):
                    k = m + 4 * n
                    grad[k] += 2 * (c[i, j] - d[i, j]) * colour_gradient(x[j], y[i], coeffs, m, n)
    return grad

In this code, d is the normalized numpy array of the target image and x and y are the lists of the x- and y-coordinates of the corresponding pixel centers.

Minimizing this L2 loss in each grid cell gives the vectorization that best approximates the target raster image.

Results

Below are some test bitmap images and their vectorized outputs.

Though the curves mimic the structure of the original sketch, they are not always in the same direction as in the actual sketch. It seems like the curves assume this orientation to depict the thickness of the sketch lines.

Conclusion

Our initial plan was to utilize the code developed in the first week to create a differentiable vectorization. This approach would involve differentiating the functions we wrote to find the endpoints and intersections with respect to the coefficients of the bicubic curve. We would also use the tangents we calculated to determine the curve that best fits the original sketch.

Throughout the project, we worked with bicubic curves to ensure intersections and continuity of isolines between adjacent grid cells. However, we suspect that biquadratic curves might be sufficient for our needs and could significantly improve the implementation’s effectiveness. We are yet to confirm this and also explore whether an even lower degree might be feasible.

This project has opened up exciting avenues for future work and the potential to refine our methods further promises valuable contributions to the field. I am deeply grateful to Mikhail and Daniel for their guidance and support throughout this project.

References:

{1} Guillard, Benoit, et al. “MeshUDF: Fast and Differentiable Meshing of Unsigned Distance Field Networks”, ECCV 2022, https://doi.org/10.48550/arXiv.2111.14549

{2} Yan, Chuan, et al. “Deep Sketch Vectorization via Implicit Surface Extraction” SIGGRAPH 2024, https://doi.org/10.1145/3658197.

{3} Liu, Shichen, et al. “Soft Rasterizer: A Differentiable Renderer for Image-based 3D Reasoning” ICCV 2019, https://doi.org/10.48550/arXiv.1904.01786.

{4} “Bicubic Interpolation”. Wikipedia, https://en.wikipedia.org/wiki/Bicubic_interpolation.

{5} SciPy – fsolve (scipy.optimize.fsolve), https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fsolve.html

Categories
Math Research

Redundant SDFs

Fellows: Aniket Rajnish, Ashlyn Lee, Megan Grosse, Taylor
Mentors: Ishit Mehta, Sina Nabizadeh

Repository: https://github.com/aniketrajnish/Redundant-SDFs

1.0 Problem Statement

Signed Distance Functions (SDFs) are widely used in computer graphics for rendering and animation. They provide an implicit representation of shapes, which can be converted to explicit representations using methods such as Marching Cubes (MC) and, more recently, Reach for the Spheres (RfS) [Sellan et al. 2023]. However, as we will demonstrate later in this blog post, SDFs have certain limitations in representing various types of surfaces, which can affect the quality of the resulting explicit representations.

Vector Distance Functions (VDFs) [Gomes and Faugeras, 2003] offer a promising alternative, with the potential to overcome some of these limitations. VDFs extend the concept of SDFs by incorporating directional information, allowing for more accurate representation of complex geometries.

This study explores both SDF and VDF-based reconstruction methods, comparing their effectiveness and identifying scenarios where each approach excels. By examining these techniques in detail, we aim to provide insights into their respective strengths and weaknesses, guiding future applications in computer graphics and geometric modeling.

2.0 Reconstruction from SDF

SDFs are one method of representing a surface implicitly. They involve assigning a scalar value to a set of points in a volume. The sign  of the SDF value indicates whether the point is inside (-ve) or outside (+ve) the shape, while the magnitude represents the distance to the nearest point on the surface

Mathematically, an SDF \(f(x)\) for a surface \(S\) can be defined as:

$$
f(x) = \begin{cases} -\min_{y \in S} |x – y| & \text{if } x \text{ is inside } S\\
0 & \text{if } x \text{ is on } S\\
\min_{y \in S} |x – y| & \text{if } x \text{ is outside } S \end{cases}
$$

Where, \(|x-y|\) denotes the Euclidean distance between points \(x\) and \(y\).

2.1 Extracting SDF

To reconstruct a surface from an SDF, we first need to obtain SDF values for a set of points in space. There are two main approaches we’ve implemented for sampling these points:

  • Random Sampling – Generating uniformly distributed random points within a bounding box that encompasses the mesh. This ensures a broad coverage of the space but may miss fine details of the surface.
  • Proximity-based Sampling – Generating points near the surface of the mesh. This method provides a higher density of samples near the surface, potentially capturing more detail. We implement this by first sampling the points and then adding Gaussian noise to these points.
def extract_sdf(sample_method: SampleMethod, mesh, num_samples=1250000, batch_size=50000, sigma=0.1):
    v, f = gpy.read_mesh(mesh)
    v = gpy.normalize_points(v)

    bbox_min = np.min(v, axis=0)
    bbox_max = np.max(v, axis=0)
    bbox_range = bbox_max - bbox_min    
    pad = 0.1 * bbox_range
    q_min = bbox_min - pad
    q_max = bbox_max + pad

    if sample_method == SampleMethod.RANDOM:
        P = np.random.uniform(q_min, q_max, (num_samples, 3))
    elif sample_method == SampleMethod.PROXIMITY:
        P = gpy.random_points_on_mesh(v, f, num_samples)
        noise = np.random.normal(0, sigma * np.mean(bbox_range), P.shape)
        P += noise   
    
    sdf, _, _ = gpy.signed_distance(P, v, f)

    return P, sdf

gpy.random_points_on_mesh() generates points on the mesh surface and gpy.signed_distance() computes the SDF values for the sampled points.

Here’s an example reconstruction using both random and proximity-based sampling for the bunny mesh using the Reach for the Arcs Algorithm:

2.2 Reconstruction using Marching Cubes

Marching Cubes is a classic algorithm for extracting a polygonal mesh from an implicit function, such as an SDF. The algorithm works by dividing the space into a regular grid of cubes and then “marching” through these cubes to construct triangles that approximate the surface.

The key steps of the Marching Cubes algorithm are:

  • Divide the 3D space into a regular grid of cubes.
  • For each cube:
    • Evaluate the SDF at each of the 8 vertices.
    • Determine which edges of the cube intersect the surface (where the SDF changes sign).
    • Calculate the exact intersection points on the edges using linear interpolation.
    • Connect these intersection points to form triangles based on a predefined table of cases.

How we implement the marching cubes algorithm in our project:

def marching_cubes(self, grid_res=128):
    bbox_min = np.min(self.pts, axis=0)
    bbox_max = np.max(self.pts, axis=0)
    bbox_range = bbox_max - bbox_min
    bbox_min -= 0.1 * bbox_range
    bbox_max += 0.1 * bbox_range

    dims = (grid_res, grid_res, grid_res)
    x, y, z = np.meshgrid(
        np.linspace(bbox_min[0], bbox_max[0], dims[0]),
        np.linspace(bbox_min[1], bbox_max[1], dims[1]),
        np.linspace(bbox_min[2], bbox_max[2], dims[2]),
        indexing='ij'
    )

    grid_pts = np.column_stack((x.ravel(), y.ravel(), z.ravel()))
    self.grid_sdf = griddata(self.pts, self.sdf, grid_pts, method='linear', fill_value=np.max(self.sdf))
    self.grid_sdf = self.grid_sdf.reshape(dims)

    verts, faces, _, _ = measure.marching_cubes(self.grid_sdf, level=0)
    verts = verts * (bbox_max - bbox_min) / (np.array(dims) - 1) + bbox_min

    return verts, faces

In this implementation:

  • We first define a bounding box that encompasses all sample points.
  • We create a regular 3D grid within this bounding box.
  • We interpolate SDF values onto this regular grid using scipy.interpolate.griddata.
  • We use the skimage.measure.marching_cubes function to perform the actual Marching Cubes algorithm.
  • Finally, we scale and translate the resulting vertices back to the original coordinate system.

The Marching Cubes algorithm uses a look-up table to determine how to triangulate each cube based on the sign of the SDF at its vertices. There are 256 possible configurations, which can be reduced to 15 unique cases due to symmetry.

One limitation of the Marching Cubes algorithm is that it can miss thin features or sharp edges due to the fixed grid resolution. Higher resolutions can capture more detail but at the cost of increased computation time and memory usage.

In the next sections, we’ll explore more advanced reconstruction techniques that aim to overcome some of these limitations.

2.3 Reconstruction using Reach for the Spheres

Reach for the Spheres (RFS) is a surface reconstruction method introduced by Sellan et al. in 2023. Unlike Marching Cubes, which operates on a fixed grid, RFS iteratively refines an initial mesh to better approximate the surface represented by the SDF. The key idea behind RFS is to use local sphere approximations to guide the mesh refinement process.

The key steps of the RFS  algorithm are:

  • Start with an initial mesh (typically a low-resolution sphere).
  • For each vertex in the mesh:
    • Compute the local reach (maximum sphere radius that doesn’t intersect the surface elsewhere).
    • Move the vertex to the center of this sphere.
  • Refine the mesh by splitting long edges and collapsing short ones.
  • Repeat steps 2-3 for a fixed number of iterations or until convergence.

For a vertex \(v\) and its normal \(n\), the local reach \(r\) is computed as:

$$
r = \arg\max_{r > 0} { \forall p \in \mathbb{R}^3 : |p – v| < r \Rightarrow f(p) > -r }
$$

Where \(f(p)\) is the SDF value at point \(p\). The new position of the vertex is then:

$$v_{new} = v + r \cdot n$$

How we implement the RFS algorithm in our project:

def reach_for_spheres(self, num_iters=5):
    v_init, f_init = gpy.icosphere(2)
    v_init = gpy.normalize_points(v_init)

    sdf_func = lambda x: griddata(self.pts, self.sdf, x, method='linear', fill_value=np.max(self.sdf))

    for _ in range(num_iters):
        v, f = gpy.reach_for_the_spheres(self.pts, sdf_func, v_init, f_init)
        v_init, f_init = v, f

    return v, f

In this implementation:

  • We start with an initial icosphere mesh with the subdivision level set to 2 which gives us  162 vertices and 320 faces.
  • We define an SDF function using interpolation from our sampled points.
  • We iteratively apply the reach_for_the_spheres function from the gpytoolbox library.
  • The resulting vertices and faces are used as input for the next iteration.

Based on empirical observations, we’ve found that the best results are achieved when the number of faces and vertices in the initial icosphere is close to, but slightly less than, the expected output mesh complexity. Keeping it greater just simply results in the icosphere as the output.

Subdivision Level (n)VerticesFaces
01220
14280
2162320
36421280
425625120
51024220480
64096281920
Subdivision levels of an icosphere and corresponding vertex and face counts generated using the gpy.icosphere function.

Advantages of RFS:

  • Adaptive resolution: RFS naturally adapts to surface features, using smaller spheres in areas of high curvature and larger spheres in flatter regions.
  • Feature preservation: The method is better at preserving sharp features compared to Marching Cubes
  • Topology handling: RFS can handle complex topologies and thin features more effectively than grid-based method.

2.4 Reconstruction using Reach for the Arcs

When reconstructing a mesh with Reach for the Arcs (RFA), this algorithm allows us to generate a mesh from a signed distance function. This method is a sort of extension of Reach for the Spheres in the way that it is focused on the curvature of our initial mesh as the focus for the reconstruction.  As mentioned in the paper, the run time is slower than reconstructing using Reach for the Sphere but it produces better outputs for closed surfaces with some curvature like spot and strawberry models in Section 2.5.

The key difference between RFA & RFS:

  • Instead of spheres, we use circular arcs to approximate the local surface.
  • The algorithm considers the principal curvatures of the surface when determining the arc parameters.
  • RFA includes additional steps for handling thin features and self-intersections.

For a vertex \(v\) with normal \(n\) and principal curvatures \(\kappa_1\) and \(\kappa_2\), the local arc is defined by:

$$ v(t) = v + r(\cos(t) – 1)n + r\sin(t)d $$

Where \(r = 1/\max(|\kappa_1|, |\kappa_2|)\) is the arc radius, and \(d\) is the direction of maximum curvature.

How we implement the RFA algorithm in our project:

def reach_for_arcs(self, **kwargs):
    return gpy.reach_for_the_arcs(
        self.pts, self.sdf,
        rng_seed=kwargs.get('rng_seed', 3452),
        fine_tune_iters=kwargs.get('fine_tune_iters', 3),
        batch_size=kwargs.get('batch_size', 10000),
        screening_weight=kwargs.get('screening_weight', 10.0),
        local_search_iters=kwargs.get('local_search_iters', 20),
        local_search_t=kwargs.get('local_search_t', 0.01),
        tol=kwargs.get('tol', 0.0001)
    )

This implementation uses the reach_for_the_arcs function from the gpytoolbox library. The key parameters are:

  • fine_tune_iters: Number of iterations for fine-tuning the mesh
  • batch_size: Number of points to process in each batch
  • screening_weight: Weight for the screening term in the optimization
  • local_search_iters and local_search_t: Parameters for the local search step.

Advantages of RFA:

  • Improved handling of high curvature: The use of arcs allows for better approximation of highly curved surfaces.
  • Better preservation of thin features: RFA includes specific steps to handle thin features that might be missed by RFS.
  • Reduced self-intersections: The algorithm includes checks to prevent self-intersections in the reconstructed mesh.

2.5 Validation of SDF reconstruction methods

For the validation we used the following params while performing the SDF based reconstruction:

  • SDF Computation
    • Method: SDFMethod.SAMPLE (sampling points)
    • Sampling Method: SampleMethod.PROXIMITY (sample points around the mesh surface)
    • Batch Size: 50000 (earlier we were extracting SDF in batches)
    • Sigma: 0.1 (for gaussian noise)
  • Marching Cubes
    • Grid resolution: 128 x 128 x 128
  • Reach for The Spheres
    • Number of Iterations: 10
    • Initial mesh: Icosphere with subdivision level 4 (2562 verts & 5120 faces)
  • Reach for The Arcs
    • Fine-tune iterations: 3
    • Batch size: 10,000
    • Screening weight: 10
    • Local search iterations: 20
    • Local search step size: 0.01
    • Tolerance: 0.0001

3.0 Reconstruction using VDF (Vector Distance Functions)

Reconstructions using Vector Distance Functions (VDFs) are similar to those using SDFs, except we have added a directional component. After creating a grid of points in a space, we can construct a vector stemming from each point that is both pointing in the direction that minimizes the distance to the mesh, as well as has the magnitude of the SDF at that point. In this way, the tip of each vector will lie on the real surface. From here, we can extract a point cloud of points on the surface and implement a reconstruction method.

Mathematically, VDF can be used for arbitrary smooth maniform \(M\) with \(k\) dimension.

The following proof is from Gomes and Faugeras’ 2003 paper:

Given a smooth manifold \(M\), which is a subset of \(\mathbb{R}^n\), for every point \(x\), \(\delta(x)\) is that distance of \(\textbf{x}\) to \(M\). Let \(\textbf{u(x)}\) be the vector length of \(\delta(x)\), since \(u(x) = \delta{(x)}D\delta{(x)}\) for \(||D\delta{x}||\)=1. At a point \(\textbf{x}\) where \(\delta \) is differeientable, let \(\textbf{y} = P_{M}(x)\) be an unique projection of \(\textbf{x}\) onto \(M\). This point is such that \(\delta(x) = || x-y||\). Assuming \(M\) is smooth at \(y\), \(\textbf{x-y}\) is normal to \(M\) and parallel to \(D\delta{(x)}\).

Thus:

$$
\textbf{u(x)}= \textbf{x-y}= x-P_{M}(x)
$$

Advantages of VDFs:

Due to their directional nature, VDFs have the potential to have more versatility and accuracy compared to SDFs.

There are four main benefits of using VDFs instead of SDFs:

  • VDFs can be used with shapes with arbitrary codimension.
  • VDFs perform better mesh reconstruction for closed surfaces.
  • VDFs can represent open surfaces that SDFs cannot.
  • VDFs can represent non-manifold surfaces like a Möbius strip.

Disadvantages of VDFs:

Despite these advantages, VDFs do have some drawbacks:

  • VDFs require more information to be known compared to SDFs
    • SDFs only need one number (distance) associated with each point in the grid
    • VDFs, on the other hand, require three numbers (assuming a 3D grid) associated with each point, representing the x, y, and z directions associated with each vector
  • VDFs can be harder to visualize/debug compared to SDFs
  • Depending on the method used to obtain them, VDFs can be inaccurate at some grid points
    • For example, singularities in the gradient can create issues when attempting to use the gradient for a VDF
    • These inaccurate VDFs can create inaccurate point clouds, which can result in poor reconstructions

Now, we will discuss the methods we used for creating VDFs and using them for reconstructions.

3.1 Reconstruction using Gradient

The gradient-based method leverages the fact that the gradient of a Signed Distance Function (SDF) points in the direction of maximum increase in distance. By reversing this direction and scaling by the SDF value, we obtain vectors pointing towards the surface.

Given an SDF \(\phi(x)\), the VDF \(\mathbf{u}(x)\) is defined as:

$$ \mathbf{u}(x) = -\phi(x) \frac{\nabla \phi(x)}{|\nabla \phi(x)|} $$

Where \(\nabla \phi(x)\) is the gradient of the SDF at point \(x\).

Implementation:

def compute_gradient(self):
    def finite_difference(f, axis):
        f_pos = np.roll(f, shift=-1, axis=axis)
        f_neg = np.roll(f, shift=1, axis=axis)
        f_pos[0] = f[0]  # neumann boundary condition
        f_neg[-1] = f[-1]
        return (f_pos - f_neg) / 2.0
    
    distance_grid = self.signed_distance.reshape(self.grid_shape)
    gradients = np.zeros((np.prod(self.grid_shape), 3))
    for dim in range(3):
        grad = finite_difference(distance_grid, axis=dim).flatten()
        gradients[:, dim] = grad
    
    magnitudes = np.linalg.norm(gradients, axis=1, keepdims=True)
    magnitudes = np.clip(magnitudes, a_min=1e-10, a_max=None)
    normalized_gradients = gradients / magnitudes

    signed_distance_reshaped = self.signed_distance.reshape(-1, 1)
    self.vector_distance = -1 * normalized_gradients * signed_distance_reshaped
    self.vector_distance[:, [0, 1, 2]] = self.vector_distance[:, [1, 0, 2]]

def compute_surface_points(self):
    self.vdf_pts = self.grid_vertices + self.vector_distance * self.signed_distance.reshape(-1, 1)

This method computes gradients using finite differences, normalizes them, and then scales by the SDF value to obtain the VDF. The surface points are then computed by adding the VDF to the grid vertices.

However, the gradient VDF reconstruction method, while theoretically sound, presented significant practical challenges in our experiments. As evident in the image above, this method produced numerous artifacts and distortions:.These issues likely arise from singularities and numerical instabilities in the gradient field, especially in complex geometries.

Given these results, we opted to focus on the barycentric coordinate method for our VDF reconstructions. This approach proved more robust and reliable, especially for complex 3D models with intricate surface details.

3.2 Reconstruction using Barycentric Coordinates

The barycentric coordinate method utilizes the signed_distance function from gpytoolbox to find the closest points on the mesh for each grid point. These closest points, represented in barycentric coordinates, are then converted to Cartesian coordinates to form a point cloud.

For a point \(p\) with barycentric coordinates \((u, v, w)\) relative to a triangle with vertices \((v_1, v_2, v_3)\), the Cartesian coordinates are:

$$ p = u v_1 + v v_2 + w v_3 $$

Implementation:

def compute_barycentric(self):
    self.signed_distance, self.face_indices, self.barycentric_coords = gpy.signed_distance(self.grid_vertices, self.V, self.F)
    self.vdf_pts = self.barycentric_to_cartesian()
    self.vector_distance = self.grid_vertices - self.vdf_pts
    magnitudes = np.linalg.norm(self.vector_distance, axis=1, keepdims=True)
    self.vector_distance = self.vector_distance / magnitudes

def barycentric_to_cartesian(self):
    cartesian_coords = np.zeros((self.barycentric_coords.shape[0], 3))
    face_vertex_coordinates = self.V[self.F[self.face_indices]]
    for i, (bary_coords, face_vertices) in enumerate(zip(self.barycentric_coords, face_vertex_coordinates)):
        cartesian_coords[i] = np.dot(bary_coords, face_vertices)
    return cartesian_coords

This method computes the closest points on the mesh for each grid point and converts them to Cartesian coordinates. The VDF is then computed as the normalized vector from the grid point to its closest point on the mesh.

3.3 Reconstruction using one point per face

The one-point-per-face method aims to create a more uniform distribution of points on the surface by using the centroid of each face as a sample point. This approach is particularly effective for meshes with similarly-sized faces.

For a triangle face with vertices \((v_1, v_2, v_3)\), the centroid \(c\) is:

$$ c = \frac{v_1 + v_2 + v_3}{3} $$

The normal vector \(\mathbf{n}\) for the face is:

$$ \mathbf{n} = \frac{(v_2 – v_1) \times (v_3 – v_1)}{|(v_2 – v_1) \times (v_3 – v_1)|} $$

Implementation:

def compute_centroid_normal(self):
    print('Computing centroid-normal reconstruction...')
    
    def compute_triangle_centroids(vertices, faces):
        return np.mean(vertices[faces], axis=1)

    def compute_face_normals(vertices, faces):
        v0, v1, v2 = vertices[faces[:, 0]], vertices[faces[:, 1]], vertices[faces[:, 2]]
        normals = np.cross(v1 - v0, v2 - v0)
        return normals / np.linalg.norm(normals, axis=1, keepdims=True)

    centroids = compute_triangle_centroids(self.V, self.F)
    normals = compute_face_normals(self.V, self.F)

    # Double-sided representation
    self.vdf_pts = np.concatenate((centroids, centroids), axis=0)
    self.vector_distance = np.concatenate((normals, -1 * normals), axis=0)

This method computes the centroids and normals for each face of the mesh. It creates a double-sided representation by duplicating the centroids and flipping the normals, which can help in creating a watertight reconstruction.

4.0 Analyzing the Results

After performing various reconstruction methods using Signed Distance Functions (SDFs) and Vector Distance Functions (VDFs), it’s crucial to analyze and compare the results. This section covers the visualization, and quantitative evaluation of the reconstructed meshes.

4.1 Visualizing the Results

The project uses the Polyscope library for visualizing the original and reconstructed meshes. The visualize method in all SDFReconstructor, VDFReconstructor and RayReconstructor classes handles the visualization.

This method creates a visual comparison between the original mesh and the reconstructed meshes or point clouds from different methods. The meshes are positioned side by side for easy comparison.

Additionally, the project includes a rendering functionality to export the visualization results as images. This is implemented in the render.py file

4.2 Fitness Score using Chamfer Distance

The project uses the Chamfer distance to compute a fitness score for the reconstructed meshes. The implementation is in the fitness.py file.

The Chamfer distance measures the average nearest-neighbor distance between two point sets. The fitness score is computed as:

$$\text{Fitness} = \frac{1}{1 + \text{Chamfer distance}}$$

This results in a value between 0 and 1, where higher values indicate better reconstruction.

5.0 Some Naive experimentation to get the results

Our initial experiments with the gradient VDF method produced unexpected results due to singularities, leading to inaccurate point clouds. This prompted us to explore alternative VDF methods that directly sample points on the surface, resulting in more reliable reconstructions.

We implemented a ray-based reconstruction technique, the RayReconstructor class, as an additional method. This method generates a point cloud by shooting random rays and finding their intersections with the surface, then uses Poisson Surface Reconstruction to create the final mesh.

6.0 Takeaways & Further Exploration

Overall, VDFs are a promising alternative to SDFs and with the correct reconstruction method, can provide almost identical resultant meshes. They are able to provide more detail at the same resolution when compared to Marching Cubes and Reach for the Spheres. While Reach for the Arcs gives an impressive amount of detail, it often adds too much detail to where it results in details that are not present in the original mesh. VDFs, when constructed/used correctly, do not have this problem. If time permitted, we would have experimented with more non-manifold surfaces and compared the results of SDF and VDF reconstructions of these surfaces. VDFs do not require a clear outside and inside like SDFs do, so we could also explore more such surfaces manifold and non-manifold alike.

References

William E. Lorensen and Harvey E. Cline. 1987. Marching cubes: A high resolution 3D surface construction algorithm. SIGGRAPH Comput. Graph. 21, 4 (July 1987), 163–169. https://doi.org/10.11enecccevitchejdkevcjvrgvjhjhhkgcbggtkhtl45/37402.37422

Silvia Sellán, Christopher Batty, and Oded Stein. 2023. Reach For the Spheres: Tangency-aware surface reconstruction of SDFs. In SIGGRAPH Asia 2023 Conference Papers (SA ’23). Association for Computing Machinery, New York, NY, USA, Article 73, 1–11.
https://doi.org/10.1145/3610548.3618196

Silvia Sellán, Yingying Ren, Christopher Batty, and Oded Stein. 2024. Reach for the Arcs: Reconstructing Surfaces from SDFs via Tangent Points. In ACM SIGGRAPH 2024 Conference Papers (SIGGRAPH ’24). Association for Computing Machinery, New York, NY, USA, Article 23, 1–12.
https://doi.org/10.1145/3641519.3657419

José Gomes and Olivier Faugeras. The Vector Distance Functions. International Journal of Computer Vision, vol. 52, no. 2–3, 2003.
https://doi.org/10.1023/A:1022956108418

Categories
Research

Implementing Spectral Conformal Parameterization using RXMesh

Fellows: Sachin Kishan and Diany Pressato

Teaching Assistant: Supriya Gadi Patil

Mentor: Ahmed Mahmoud

Introduction

RXMesh is a framework for processing triangle mesh on the GPU. It allows developers to easily use the GPU’s massive parallelism to perform common geometry processing tasks on triangle mesh. For a preliminary introduction to RXMesh, check out this tutorial blog here.

Our task was to use RXMesh to implement some classic algorithms in geometry processing. This blog details one of the algorithms we implemented: Spectral Conformal Parameterization (SCP).

This blog will discuss the algorithm in terms of its basic concepts, the steps of the algorithm, the implementation using RXMesh, and our results.

Our entire implementation of SCP application in RXMesh can be found here.

Spectral Conformal Parameterization

Parameterization refers to the mapping of a triangle surface mesh in \(R^{3}\) space to a planar \(R^{2}\) space. This mapping can then be subsequently used in different application such as texture mapping, where we map the values from some other \(R^{2}\) space (e.g., image) to the surface of the mesh. In this case, we map color values to the mesh surface.

Often when we perform mesh parameterization, there is a possibility to change some properties of the geometric elements of the mesh, such as the size, shape, or angle between edges. It is desired to preserve these properties since we usually want to create a direct mapping of our intended values (e.g., texture color) to be in the “correct” positions on the mesh. What is correct is usually up to artistic opinion, and thus a mapping that preserves these properties is desirable. A mapping that conserves angles (and thus shape) is called a conformal mapping.

Mullen et al. (2008) provides a method to perform conformal parameterization of an open surface mesh by modelling the problem as a generalized Eigenvalue problem. Using an initial guess, we can perform multiple iterations of a general eigen vector solving method (i.e., the power method [2] ) to determine the largest eigen vector. This eigen vector, in our case, is the parameterized coordinates of the given mesh, giving us our desired output.

The algorithm for the power method is as follows:

Given \((A-\lambda B)x=0 \)

We want to find \(x\).

Given an initial guess of eigen-vector \(x_{0}\)

For some number of iterations, on iteration k

  • \(B y_{k+1} = A x_{k}\)
  • \(x_{k+1}= \frac{y_{k+1}}{\left| y_{k+1} \right|} \)

where \(x\) is the desired eigenvector solution.

Implementation

We applied the power method to our application in the following way, based on the formulation given by Mullen et al. (2008) in equation 7:

\([B-\frac{1}{V_{b}}e_{b}e_{b}^{t}]u=\mu L_{c}u \)

To find the eigenvector \(u\), we can the power method since the problem is in the form \(B y_{k+1} = A x_{k}\). This can be done in the following way:

\(T2=dot(eb,u) \)
\(T1=Bu – T2 * eb \)
\(L_{c}u=T1 \)
\(u=u/normalise (u) \)

Where \(u\) is the parameterized coordinates. \(T1\) and \(T2\) are temporary variables to store certain scalar or vector values. Other variables are explained by Mullen et al. (2008). This loop can be implemented like this using RXMesh:

for (int i = 0; i < iterations; i++) 
{
    cuComplex T2 = eb.dot(uv_mat);
    rx.for_each_vertex(rxmesh::DEVICE,
                       [eb, T2, T1, uv_mat, B] __device__(
                       const rxmesh::VertexHandle vh) mutable
                       {
                           T1(vh, 0) =
                       cuCsubf(cuCmulf(B(vh, vh), uv_mat(vh, 0)),
                       cuCmulf(eb(vh, 0), T2));
                       });

    Lc.solve(T1, uv_mat);
    float norm = uv_mat.norm2();
    uv_mat.multiply(1.0f / norm);
    if (std::abs(prv_norm - norm) < 0.0001)//convergence check
    {
        break; 
    }
    prv_norm = norm;
}

\(T1\) is a vector whose values are per-vertex, which is why we perform a per vertex computation. \(T2\) is a scalar in which we store the value of the dot product dot(eb,uv).

Since our system matrix \(L_{c}\) has dimensions 2V x 2V, we use complex representations of value entries in the vectors and matrices we use. So, for example, our \(L_{c}\) matrix now has dimensions V x V where each entry is a complex number (with a real and complex part).

This means subsequent computations must use CUDA functions that operate on complex numbers (e.g., cuCsubf which subtracts complex numbers).

All initial pre-computations only allocate to the real part of the entry, except for the calculation of the area term \(A\) where we calculate \(L_{c}=L_{D}-A\) . For the area term, we assign values to the complex part of the area term. \(L_{D}\) represents the standard Laplace-Beltrami operator, which we calculate using the cotangent edge weights and assign to the real part.

After we obtain the eigenvector \(u\) in complex form, we can convert it to a standard V x 2 matrix form and use that for further processing (such as rendering).

Results

Using the SCP implementation, we can now apply texturing to our open meshes as shown in the examples below.

The bunny with a checker texture.
The beetle mesh with a checker texture.

References

[1] Mullen, P., Tong, Y., Alliez, P., & Desbrun, M. (2008). Spectral Conformal Parameterization. Computer Graphics Forum, 27(5), 1487–1494. https://doi.org/10.1111/j.1467-8659.2008.01289.x

[2] Kozłowski, W. (1992). The power method for the generalized eigenvalue problem. Roczniki Polskiego Towarzystwa Matematycznego. Seria 3, Matematyka Stosowana/Matematyka Stosowana/Mathematica Applicanda, 21(35). https://doi.org/10.14708/ma.v21i35.1790

This blog was written by Sachin Kishan and Diany Pressato during the SGI 2024 Fellowship as one of the outcomes of a two-week project under the mentorship of Ahmed Mahmoud and the support of Supriya Gadi Patil as a teaching assistant.

Categories
Research

Implementing ARAP Deformations using RXMesh

Fellows: Sachin Kishan

Teaching Assistant: Supriya Gadi Patil

Mentor: Ahmed Mahmoud

Introduction

RXMesh is a framework for processing triangle mesh on the GPU. It allows developers to easily use the GPU’s massive parallelism to perform common geometry processing tasks. For a preliminary introduction to RXMesh, check out this tutorial blog here.

Our task was to use RXMesh to implement some classic algorithms in geometry processing. This blog details one of the algorithms we implemented: As-Rigid-As-Possible(ARAP) Surface Modeling [1].

This blog will discuss the algorithm in terms of its basic concepts, the steps of the algorithm, the implementation using RXMesh, and the final results.

Our entire implementation of ARAP deformation in RXMesh can be found here.

As-Rigid-As-Possible Surface Modeling

Sorkine and Alexa (2007) formulated a method where we can retain the shape of an object as much as possible despite deforming one (or a set of) its vertices. This retention of shape is achieved by ensuring that subsequent transformations to the mesh after deforming its vertices are rigid. This property of rigidity means that a mesh can only undergo rigid transformations such as translation or rotation, but not stretching or shearing (since these are not rigid transformations).

The algorithm’s focus is to minimize the energy which measures the extent of deformation (dubbed rigidity energy). This energy is minimized when given a specific set of positions or a specific set of rotations. Notice that we can not minimize both at once—since we need the minimized rotational transformation for minimizing the energy from deformed positions. Thus, we instead solve both alternatively, i.e., solve for deformed positions assuming rotation is fixed, then solve for rotation assuming the position is fixed.

Using a Laplacian, we can form a sparse system of linear equations which we can then solve to obtain a new set of positions for the mesh. The algorithm works as follows for a single iteration

  1. Identify the rotation matrix for the given deformed positions which minimizes deformation energy.
  2. Solve for the new positions to minimize the energy, given a set of linear equations.

In this case, the set of linear equations to solve is in the form :

\(Lx=b\)

Where \(L\) is the Laplace-Beltrami operator, which can be formed by finding the cotangent edge weights of the mesh, \(x\) represents the set of new vertex positions we solve for, and \(b\) represents the right-hand side of Equation 8 from Sorkine and Alexa (2007). The derivation is detailed further in the original paper.

\(b = \sum_{j\in N(i)}^{} \frac{w_{ij}}{2} (R_{i}+R_{j})(p_{i}-p_{j})\)

To further minimize the energy, we can perform multiple iterations of the same steps, using updated vertex positions and consequent rotations per iteration.

Implementation

The following code snippets show how the algorithm is implemented in RXMesh.

We first deform the input mesh and set constraints on each vertex depending on whether it is fixed, deformed, or free to move. This will be taken into account when calculating the system matrix \(L\).

 // set constraints
    const vec3<float> sphere_center(0.1818329, -0.99023, 0.325066);
    rx.for_each_vertex(DEVICE, [=] __device__(const VertexHandle& vh) {
        const vec3<float> p(deformed_vertex_pos(vh, 0),
                            deformed_vertex_pos(vh, 1),
                            deformed_vertex_pos(vh, 2));

        // fix the bottom
        if (p[2] < -0.63) {
            constraints(vh) = 2;
        }

        // move the jaw
        if (glm::distance(p, sphere_center) < 0.1) {
            constraints(vh) = 1;
        }
    });

In the code above, we set the mesh to deform the dragon’s jaw and fix its bottom vertices (legs). This can vary depending on the type of input we want to have. Another possibility is to add a feature for users to click and select vertices to deform or fix them in space interactively.

Before we perform the iteration step, we calculate the Laplacian. It can be calculated as the cotangent edge weight matrix. So, we create a kernel using RXMesh.

    auto calc_mat = [&](VertexHandle v_id, VertexIterator& vv) {
        if (constraints(v_id, 0) == 0) {
            for (int i = 0; i < vv.size(); i++) {
          laplace_mat(v_id, v_id) += weight_matrix(v_id, vv[i]);
          laplace_mat(v_id, vv[i]) -= weight_matrix(v_id, vv[i]);
            }
        } else {
            for (int i = 0; i < vv.size(); i++) {
                laplace_mat(v_id, vv[i]) = 0;
            }
            laplace_mat(v_id, v_id) = 1;
        }
    };

The weight matrix is calculated using a special query called “EVDiamond”, where we obtain the two end vertices and the two opposite vertices of an edge. This gives us exactly the data we need to find the cotangent edge weight.

    auto calc_weights = [&](EdgeHandle edge_id, VertexIterator& vv) {
        // the edge goes from p-r while the q and s are the opposite vertices
        const rxmesh::VertexHandle p_id = vv[0];
        const rxmesh::VertexHandle r_id = vv[2];
        const rxmesh::VertexHandle q_id = vv[1];
        const rxmesh::VertexHandle s_id = vv[3];

        if (!p_id.is_valid() || !r_id.is_valid() || !q_id.is_valid() ||
            !s_id.is_valid()) {
            return;
        }

        const vec3<T> p(coords(p_id, 0), coords(p_id, 1), coords(p_id, 2));
        const vec3<T> r(coords(r_id, 0), coords(r_id, 1), coords(r_id, 2));
        const vec3<T> q(coords(q_id, 0), coords(q_id, 1), coords(q_id, 2));
        const vec3<T> s(coords(s_id, 0), coords(s_id, 1), coords(s_id, 2));


        T weight = 0;
   weight += dot((p - q), (r - q)) / length(cross(p - q, r - q));
   weight += dot((p - s), (r - s)) / length(cross(p - s, r - s));

        weight /= 2;
        weight = std::max(0.f, weight);

        edge_weights(p_id, r_id) = weight;
        edge_weights(r_id, p_id) = weight;
    };

We now have all the building blocks necessary to perform an iteration step. This is what an iteration loop looks like:

    // process step
    for (int i = 0; i < iterations; i++) {
        // solver for rotation
        calculate_rotation_matrix<float, CUDABlockSize>
            <<<lb_rot.blocks, lb_rot.num_threads, lb_rot.smem_bytes_dyn>>>(
                rx.get_context(),
                ref_vertex_pos,
                deformed_vertex_pos,
                rotations,
                weight_matrix);

        // solve for position
        calculate_b<float, CUDABlockSize>
            <<<lb_b_mat.blocks,
               lb_b_mat.num_threads,
               lb_b_mat.smem_bytes_dyn>>>(rx.get_context(),
                                          ref_vertex_pos,
                                          deformed_vertex_pos,
                                          rotations,
                                          weight_matrix,
                                          b_mat,
                                          constraints);

        laplace_mat.solve(b_mat, *deformed_vertex_pos_mat);
    }

We perform multiple iterations of the algorithm to effectively minimize the energy. We can see in the code above that we alternatively solve for rotation and then for positions. We first calculate the rotation matrix given a set of deformed positions, then calculate the \(b\) matrix (as defined earlier) and then solve the system of equations to identify the new vertex positions.

The following kernel calculates the rotation matrix per vertex to minimize the rigidity energy given their positions.

auto cal_rot = [&](VertexHandle v_id, VertexIterator& vv) {
        Eigen::Matrix3f S = Eigen::Matrix3f::Zero();

        for (int j = 0; j < vv.size(); j++) {

            float w = weight_matrix(v_id, vv[j]);

            Eigen::Vector<float, 3> pi_vector = {
              ref_vertex_pos(v_id, 0) - ref_vertex_pos(vv[j], 0),
              ref_vertex_pos(v_id, 1) - ref_vertex_pos(vv[j], 1),
              ref_vertex_pos(v_id, 2) - ref_vertex_pos(vv[j], 2)
};


            Eigen::Vector<float, 3> pi_dash_vector = {
                deformed_vertex_pos(v_id, 0) - deformed_vertex_pos(vv[j], 0),
                deformed_vertex_pos(v_id, 1) - deformed_vertex_pos(vv[j], 1),
                deformed_vertex_pos(v_id, 2) - deformed_vertex_pos(vv[j], 2)};

            S = S + w * pi_vector * pi_dash_vector.transpose();
        }

        Eigen::Matrix3f U;         // left singular vectors
        Eigen::Matrix3f V;         // right singular vectors
        Eigen::Vector3f sing_val;  // singular values

        svd(S, U, sing_val, V);

        const float smallest_singular_value = sing_val.minCoeff();


        Eigen::Matrix3f R = V * U.transpose();

        if (R.determinant() < 0) {
            U.col(smallest_singular_value) =
                U.col(smallest_singular_value) * -1;
            R = V * U.transpose();
        }

        // Matrix R to vector attribute R

        for (int i = 0; i < 3; i++)
            for (int j = 0; j < 3; j++)
                rotations(v_id, i * 3 + j) = R(i, j);
    };

We first find the Single Value Decomposition (SVD) of the sum (\(S_{i} \)) where

\(S_{i}=\sum_{j\in N(i)} w_{ij}e_{ij}e^{‘T}_{ij} \) (Equation 5 in Sorkine and Alexa(2007))

We then use that to calculate rotation and assign that as the vertice’s rotation. After we calculate rotation, we calculate the right-hand side of Equation 8 as stated earlier.

auto init_lambda = [&](VertexHandle v_id, VertexIterator& vv) {
        // variable to store ith entry of b_mat
        Eigen::Vector3f bi(0.0f, 0.0f, 0.0f);

        // get rotation matrix for ith vertex
        Eigen::Matrix3f Ri = Eigen::Matrix3f::Zero(3, 3);

        for (int i = 0; i < 3; i++) {
            for (int j = 0; j < 3; j++)
                Ri(i, j) = rotations(v_id, i * 3 + j);
        }

        for (int nei_index = 0; nei_index < vv.size(); nei_index++) {
            // get rotation matrix for neightbor j
            Eigen::Matrix3f Rj = Eigen::Matrix3f::Zero(3, 3);
            for (int i = 0; i < 3; i++)
                for (int j = 0; j < 3; j++)
                    Rj(i, j) = rotations(vv[nei_index], i * 3 + j);

            // find rotation addition
            Eigen::Matrix3f rot_add = Ri + Rj;
            // find coord difference
            Eigen::Vector3f vert_diff = {
                ref_vertex_pos(v_id, 0) - ref_vertex_pos(vv[nei_index], 0),
                ref_vertex_pos(v_id, 1) - ref_vertex_pos(vv[nei_index], 1),
                ref_vertex_pos(v_id, 2) - ref_vertex_pos(vv[nei_index], 2)};

            // update bi
            bi = bi +
                 0.5 * weight_mat(v_id, vv[nei_index]) * rot_add * vert_diff;
        }

        if (constraints(v_id, 0) == 0) {
            b_mat(v_id, 0) = bi[0];
            b_mat(v_id, 1) = bi[1];
            b_mat(v_id, 2) = bi[2];
        } else {
            b_mat(v_id, 0) = deformed_vertex_pos(v_id, 0);
            b_mat(v_id, 1) = deformed_vertex_pos(v_id, 1);
            b_mat(v_id, 2) = deformed_vertex_pos(v_id, 2);
        }
    };

After calling both these kernels, we use the solve function to solve the equations. We then update the mesh with this new set of vertex positions.

Results

Using the implementation above, we can create a per-frame ARAP implementation to animate movement for a mesh. In the example below, we fix the upper half vertices of Spot the Cow and allow the bottom half to freely deform. We then move each of the bottom “feet” vertices in every frame to animate a swimming or walking-like movement in Spot.

The bottom feet vertices are user-deformed , the upper vertices are fixed, and the bottom half (except the feet) are free to move—and are determined by the ARAP algorithm.

Here’s another example with the dragon mesh.

Deform a vertex of the dragon’s jaw while keeping its feet fixed.

Further work

Polyscope shows us that the framerate of these demos is very poor (below 30 fps). We would expect it to do better, especially when using the GPU. A performance analysis as well as optimizing of the solver will be required to make the application reach desired levels of real-time processing.

References

[1] Olga Sorkine and Marc Alexa. 2007. As-rigid-as-possible surface modeling. In Proceedings of the fifth Eurographics symposium on Geometry processing (SGP ’07). Eurographics Association, Goslar, DEU, 109–116. https://doi.org/10.2312/SGP/SGP07/109-116

This blog was written by Sachin Kishan during the SGI 2024 Fellowship as one of the outcomes of a two-week project under the mentorship of Ahmed Mahmoud and the support of Supriya Gadi Patil as a teaching assistant.