<p>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.  </p>

<p>So suppose that the 1-ring neighbourhood of $v_i$ in your mesh looks like this</p>

<p>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; <a href="https://i.sstatic.net/aWaSm.png" rel="noreferrer"><img src="https://i.sstatic.net/aWaSm.png" alt="enter image description here"></a> </p>

<p>$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. </p>

<p>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: </p>

<p>$$\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)) $$</p>

<p>Where $v_j \in N_1(v_i)$ means every vertex in the one ring neighbourhood of $v_i$. </p>

<p>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</p>

<p>$$H = \frac{1}{2} || \Delta_S v || $$</p>

<p>Now let's introduce the angle $\theta_j$ as </p>

<p>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<a href="https://i.sstatic.net/ysAHy.png" rel="noreferrer"><img src="https://i.sstatic.net/ysAHy.png" alt="enter image description here"></a></p>

<p>The Gaussian curvature is:</p>

<p>$$K = (2\pi - \sum_j \theta_j) / A$$</p>

<p>After all of this pain, the principal discrete curvatures are given by: </p>

<p>$$k_1 = H + \sqrt{H^2 - K} \ \  \text{and} \ \  k_2 = H - \sqrt{H^2 - K}$$</p>

<hr>

<p>If you are interested in the subject (and to add some reference to this post) an excellent read is:   <a href="http://http.cs.berkeley.edu/~jrs/meshpapers/MeyerDesbrunSchroderBarr.pdf" rel="noreferrer">Discrete Differential-Geometry Operators
for Triangulated 2-Manifolds</a> [Meyer et al. 2003]. </p>

<p>For the images I thank my ex-professor Niloy Mitra as I found them in some notes I took for his lectures.</p>

CC BY-SA 3.0
2015-11-17T17:28:29.103Z
Both answers are really good, it was hard for me to pick. Since I had asked about the simplest way, I think Nathan takes the cake.
ap_
CC BY-SA 3.0
2016-05-28T09:33:15.520Z
What Meyer et al. 2003 did not (perhaps explicitly) mention was how to compute the curvature for border vertices. Since they have taken an angle deficit approach, the Gaussian curvature for border vertices should read $K = (\pi - \sum_j{\theta_j})/A_{mixed}$.
teodron
CC BY-SA 3.0
2016-11-18T12:00:45.223Z
@teodron Might you have any insights on mean curvature for border vertices? Can such a thing be defined?
Museful
CC BY-SA 3.0
2016-11-18T12:50:42.583Z
@Museful I am a bit worried about the mean curvature not being negative, no matter the surface type. If the Laplacian-like operator is defined on a border vertex, then it's a matter of evaluating the same expression including only the triangles that make up the faces of the surface incident at $v_i$. There are more recent papers on discrete curvatures, however..
teodron
CC BY-SA 3.0
2015-11-14T09:36:56.117Z
2015-11-14T09:54:46.410Z
2015-11-14T09:54:46.410Z
cifz
<p>When I needed an estimate of mesh curvature for a skin shader, the algorithm I ended up settling on was this:</p>
<p>First, I computed a scalar curvature for each edge in the mesh. If the edge has positions <span class="math-container">$p_1, p_2$</span> and normals <span class="math-container">$n_1, n_2$</span>, then I estimated its curvature as:</p>
<p><span class="math-container">$$\text{curvature} = \frac{(n_2 - n_1) \cdot (p_2 - p_1)}{|p_2 - p_1|^2}$$</span></p>
<p>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.)</p>
<p>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 &quot;average curvature&quot;, 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.</p>
<hr />
<p>The motivation for this formula is looking at what happens in 2D when applied to a circle:</p>
<p><img src="https://i.sstatic.net/PL8jY.png" alt="curvature formula applied to two points on a circle" /></p>
<p>Suppose you have a circle of radius <span class="math-container">$r$</span> (so its curvature is <span class="math-container">$1/r$</span>), and you have two points on the circle, with their normals <span class="math-container">$n_1, n_2$</span>. The positions of the points, relative to the circle's center, are going to be <span class="math-container">$p_1 = rn_1$</span> and <span class="math-container">$p_2 = rn_2$</span>, due to the property that a circle or sphere's normals always point directly out from its center.</p>
<p>Therefore you can recover the radius as <span class="math-container">$r = |p_1| / |n_1|$</span> or <span class="math-container">$|p_2| / |n_2|$</span>. But in general, the vertex positions won't be relative to the circle's center. We can work around this by subtracting the two:
<span class="math-container">$$\begin{aligned}
p_2 - p_1 &amp;= rn_2 - rn_1 \\
&amp;= r(n_2 - n_1) \\
r &amp;= \frac{|p_2 - p_1|}{|n_2 - n_1|} \\
\text{curvature} = \frac{1}{r} &amp;= \frac{|n_2 - n_1|}{|p_2 - p_1|}
\end{aligned}$$</span></p>
<p>The result is exact only for circles and spheres. However, we can extend it to make it a bit more &quot;tolerant&quot;, and use it on arbitrary 3D meshes, and it seems to work reasonably well. We can make the formula more &quot;tolerant&quot; by first projecting the vector <span class="math-container">$n_2 - n_1$</span> onto the direction of the edge, <span class="math-container">$p_2 - p_1$</span>. 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:
<span class="math-container">$$\begin{aligned}
\text{curvature} &amp;= \frac{(n_2 - n_1) \cdot \text{normalize}(p_2 - p_1)}{|p_2 - p_1|} \\
&amp;= \frac{(n_2 - n_1) \cdot (p_2 - p_1)/|p_2 - p_1|}{|p_2 - p_1|} \\
&amp;= \frac{(n_2 - n_1) \cdot (p_2 - p_1)}{|p_2 - p_1|^2}
\end{aligned}$$</span></p>
<p>Et voilà, there's the formula that appeared at the top of this answer. By the way, a nice side benefit of using the <em>signed</em> projection (the dot product) is that the formula then gives a signed curvature: positive for convex, and negative for concave surfaces.</p>
<hr />
<p>Another approach I can imagine using, but haven't tried, would be to estimate the <a href="https://en.wikipedia.org/wiki/Second_fundamental_form" rel="noreferrer">second fundamental form</a> 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 &quot;implied&quot; 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.</p>
<p>A paper that takes this approach is <a href="https://pixl.cs.princeton.edu/pubs/Rusinkiewicz_2004_ECA/index.php" rel="noreferrer">Rusinkiewicz, &quot;Estimating Curvatures and Their Derivatives on Triangle Meshes&quot;</a>. It works by estimating the best-fit 2FF matrix per triangle, then averaging the matrices per-vertex (similar to how smooth normals are calculated).</p>

CC BY-SA 4.0
2019-08-05T17:26:04.553Z
FYI if that matters, I've used your answer here https://blender.stackexchange.com/questions/146819/is-there-a-way-to-calculate-mean-curvature-of-a-triangular-mesh/147371#147371 but adding a weighting using the angle around p1. Don't know if you find that valuable? Anyway feel free to comment. Thanks.
lemon
CC BY-SA 4.0
2015-11-14T00:43:02.683Z
2022-09-14T05:59:45.937Z
2022-09-14T05:59:45.937Z
Nathan Reed
<p>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 <a href="https://en.wikipedia.org/wiki/Principal_curvature">Wikipedia</a> does not help a lot. Can you describe or point me to a simple algorithm that can help me compute this estimate?</p>

<p>Assume that I know positions and normals of all vertices.</p>

<mesh>
What is the simplest way to compute principal curvature for a mesh triangle?
ap_
# What is the simplest way to compute principal curvature for a mesh triangle?

(Asked by: username_0 on 2015-11-13)

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](https://en.wikipedia.org/wiki/Principal_curvature) 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.

**Tags:** mesh


---

## 2 answers

### Answer by username_4 on 2015-11-14. Score: 36 | Accepted answer

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:

![curvature formula applied to two points on a circle](https://i.sstatic.net/PL8jY.png)

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](https://en.wikipedia.org/wiki/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"](https://pixl.cs.princeton.edu/pubs/Rusinkiewicz_2004_ECA/index.php). It works by estimating the best-fit 2FF matrix per triangle, then averaging the matrices per-vertex (similar to how smooth normals are calculated).

> **Comments:**
> **username_5**: FYI if that matters, I've used your answer here https://blender.stackexchange.com/questions/146819/is-there-a-way-to-calculate-mean-curvature-of-a-triangular-mesh/147371#147371 but adding a weighting using the angle around p1. Don't know if you find that valuable? Anyway feel free to comment. Thanks.  

---

### Answer by username_1 on 2015-11-14. Score: 25

Just to an add another way to the excellent @username_4 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

[![enter image description here](https://i.sstatic.net/aWaSm.png)](https://i.sstatic.net/aWaSm.png)

$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

[![enter image description here](https://i.sstatic.net/ysAHy.png)](https://i.sstatic.net/ysAHy.png)

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](http://http.cs.berkeley.edu/~jrs/meshpapers/MeyerDesbrunSchroderBarr.pdf) [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.

> **Comments:**
> **username_0**: Both answers are really good, it was hard for me to pick. Since I had asked about the simplest way, I think username_4 takes the cake.  
> **username_2**: What Meyer et al. 2003 did not (perhaps explicitly) mention was how to compute the curvature for border vertices. Since they have taken an angle deficit approach, the Gaussian curvature for border vertices should read $K = (\pi - \sum_j{\theta_j})/A_{mixed}$.  
> **username_3**: @username_2 Might you have any insights on mean curvature for border vertices? Can such a thing be defined?  
> **username_2**: @username_3 I am a bit worried about the mean curvature not being negative, no matter the surface type. If the Laplacian-like operator is defined on a border vertex, then it's a matter of evaluating the same expression including only the triangles that make up the faces of the surface incident at $v_i$. There are more recent papers on discrete curvatures, however..  

---

