Computer Graphics Stack Exchange is a question and answer site for computer graphics researchers and programmers. It only takes a minute to sign up.
Sign up to join this community
I have a mesh and in the region around each triangle, I want to compute an estimate of the principal curvature directions. I have never done this sort of thing before and Wikipedia does not help a lot. Can you describe or point me to a simple algorithm that can help me compute this estimate?
Assume that I know positions and normals of all vertices.
When I needed an estimate of mesh curvature for a skin shader, the algorithm I ended up settling on was this:
First, I computed a scalar curvature for each edge in the mesh. If the edge has positions $p_1, p_2$ and normals $n_1, n_2$, then I estimated its curvature as:
$$\text{curvature} = \frac{(n_2 - n_1) \cdot (p_2 - p_1)}{|p_2 - p_1|^2}$$
This calculates the difference in normals, projected along the edge, as a fraction of the length of the edge. (See below for how I came up with this formula.)
Then, for each vertex I looked at the curvatures of all the edges touching it. In my case, I just wanted a scalar estimate of "average curvature", so I ended up taking the geometric mean of the absolute values of all the edge curvatures at each vertex. For your case, you might find the minimum and maximum curvatures, and take those edges to be the principal curvature directions (maybe orthonormalizing them with the vertex normal). That's a bit rough, but it might give you a good enough result for what you want to do.
The motivation for this formula is looking at what happens in 2D when applied to a circle:
Suppose you have a circle of radius $r$ (so its curvature is $1/r$), and you have two points on the circle, with their normals $n_1, n_2$. The positions of the points, relative to the circle's center, are going to be $p_1 = rn_1$ and $p_2 = rn_2$, due to the property that a circle or sphere's normals always point directly out from its center.
Therefore you can recover the radius as $r = |p_1| / |n_1|$ or $|p_2| / |n_2|$. But in general, the vertex positions won't be relative to the circle's center. We can work around this by subtracting the two: $$\begin{aligned} p_2 - p_1 &= rn_2 - rn_1 \\ &= r(n_2 - n_1) \\ r &= \frac{|p_2 - p_1|}{|n_2 - n_1|} \\ \text{curvature} = \frac{1}{r} &= \frac{|n_2 - n_1|}{|p_2 - p_1|} \end{aligned}$$
The result is exact only for circles and spheres. However, we can extend it to make it a bit more "tolerant", and use it on arbitrary 3D meshes, and it seems to work reasonably well. We can make the formula more "tolerant" by first projecting the vector $n_2 - n_1$ onto the direction of the edge, $p_2 - p_1$. This allows for these two vectors not being exactly parallel (as they are in the circle case); we'll just project away any component that's not parallel. We can do this by dotting with the normalized edge vector: $$\begin{aligned} \text{curvature} &= \frac{(n_2 - n_1) \cdot \text{normalize}(p_2 - p_1)}{|p_2 - p_1|} \\ &= \frac{(n_2 - n_1) \cdot (p_2 - p_1)/|p_2 - p_1|}{|p_2 - p_1|} \\ &= \frac{(n_2 - n_1) \cdot (p_2 - p_1)}{|p_2 - p_1|^2} \end{aligned}$$
Et voilà, there's the formula that appeared at the top of this answer. By the way, a nice side benefit of using the signed projection (the dot product) is that the formula then gives a signed curvature: positive for convex, and negative for concave surfaces.
Another approach I can imagine using, but haven't tried, would be to estimate the second fundamental form of the surface at each vertex. This could be done by setting up a tangent basis at the vertex, then converting all neighboring vertices into that tangent space, and using least-squares to find the best-fit 2FF matrix. Then the principal curvature directions would be the eigenvectors of that matrix. This seems interesting as it could let you find curvature directions "implied" by the neighboring vertices without any edges explicitly pointing in those directions, but on the other hand is a lot more code, more computation, and perhaps less numerically robust.
A paper that takes this approach is Rusinkiewicz, "Estimating Curvatures and Their Derivatives on Triangle Meshes". It works by estimating the best-fit 2FF matrix per triangle, then averaging the matrices per-vertex (similar to how smooth normals are calculated).
Just to an add another way to the excellent @NathanReed answer, you can use mean and gaussian curvature that can be obtained with a discrete Laplace-Beltrami.
So suppose that the 1-ring neighbourhood of $v_i$ in your mesh looks like this
$A(v_i)$ can be simply a $\frac{1}{3}$ of the areas of the triangles that form this ring and the indicated $v_j$ is one of the neighbouring vertices.
Now let's call $f(v_i)$ the function defined by your mesh (must be a differentiable manifold) at a certain point. The most popular discretization of the Laplace-Beltrami operator that I know is the cotangent discretization and is given by:
$$\Delta_S f(v_i) = \frac{1}{2A(v_i)} \sum_{v_j \in N_1(v_i)} (cot \alpha_{ij} + cot \beta_{ij}) (f(v_j) - f(v_i)) $$
Where $v_j \in N_1(v_i)$ means every vertex in the one ring neighbourhood of $v_i$.
With this is pretty simple to compute the mean curvature (now for simplicity let's call the function of your mesh at the vertex of interest simply $v$ ) is
$$H = \frac{1}{2} || \Delta_S v || $$
Now let's introduce the angle $\theta_j$ as
The Gaussian curvature is:
$$K = (2\pi - \sum_j \theta_j) / A$$
After all of this pain, the principal discrete curvatures are given by:
$$k_1 = H + \sqrt{H^2 - K} \ \ \text{and} \ \ k_2 = H - \sqrt{H^2 - K}$$
If you are interested in the subject (and to add some reference to this post) an excellent read is: Discrete Differential-Geometry Operators for Triangulated 2-Manifolds [Meyer et al. 2003].
For the images I thank my ex-professor Niloy Mitra as I found them in some notes I took for his lectures.
