# Computation of the tensor of curvature on surface mesh

Is there a formula which enables the computation the tensor of curvature knowing the following at each vertex and cell of a triangulated mesh:

• Normal vector
• Two arbitrary vectors in the tangent space
• Mean and Gauss curvature

I am trying to script it using vtk/python. The final goal is to compute the principal directions, eigenvectors of this tensor.

• Can you clarify the following question: Of what would you like to know the curvature? I assume that you mean "the curvature of the manifold which the triangulation (surface mesh) describes", but I can think of other interpretations of your question as well. Sep 27 at 19:18
• So if I understand correctly, you assume that mean and gauss curvatures are known? If this is the case, you know the sum and product of the principal curvature, so you can find them. I guess you want to find these starting from the mesh. Nov 1 at 11:56

As suggested in Kindlmann the curvature of a surface is defined by the relationship between positional changes in the neighborhood of a point placed on the surface and the change in the surface normal. Given a level-set $$\Phi(\mathbf{x})$$, we consider that the value of the level-set is positive inside the object, negative outside. Hence, we define the normal unit vector: $$$$\mathbf{n}=-\frac{\nabla \Phi}{|\nabla \Phi|}$$$$ The curvature information is contained in the 3x3 matrix $$\nabla{\bf{n}^T}$$. Considering the Hessian matrix: $$$$\bf{H} = \begin{bmatrix} \frac{\partial^2\Phi}{\partial x^2}& \frac{\partial \Phi}{\partial xy} & \frac{\partial \Phi}{\partial xz} \\[1ex] % <-- 1ex more space between rows of matrix \frac{\partial \Phi}{\partial xy} & \frac{\partial^2\Phi}{\partial y^2} & \frac{\partial \Phi}{\partial yz} \\[1ex] \frac{\partial \Phi}{\partial xz} & \frac{\partial \Phi}{\partial yz} & \frac{\partial^2\Phi}{\partial z^2} \end{bmatrix}$$$$ The projection matrix is defined as $$\bf{P}=\bf{I}-\bf{n}\bf{n}^T$$ and allows us to project the hessian on the tangent plane (we are interested only on the changement in the direction of the normal vector, not of the magnitude). We multiply $$\nabla{\bf{n}^T}$$ by the projection matrix $$\bf{P}$$, obtaining the Geometric tensor: $$$$\bf{G}=\nabla{\bf{n}}^T\bf{P}$$$$ The projection matrix, defined as $$\bf{P}=\bf{I}-\bf{n}\bf{n}^T$$, projects the matrix on the tangent plane to the surface described by the function $$\Phi(\bf{x})=0$$. As described in Kindlmann it is possible to write the relationship: $$$$\nabla {\mathbf{n}^T}=-\frac{1}{|{\nabla \Phi}|}(\bf{P}\bf{H})$$$$ The Hessian matrix describes how the gradient changes around the neighborhood of the points placed on an iso-surface of the function $$\Phi(x)$$. In order to describe the curvature we are interested only in changes of direction of the gradient. Hence we project $$\bf{H}$$ on the tangent plane. The restriction of the Hessian to the tangent plane is a symmetric matrix and it is possible to find an orthonormal basis $$\{\bf{p}_1,\bf{p}_2,\bf{n}\}$$ able to diagonalize the matrix. In this basis we will obtain: $$$$\nabla {\bf{n}^T}= \begin{bmatrix} k_1& 0 & \sigma_1 \\[1ex] % <-- 1ex more space between rows of matrix 0 & k_2 & \sigma_2 \\[1ex] 0 & 0 & 0 \end{bmatrix}$$$$ $$\bf{p}_1$$ and $$\bf{p}_2$$ are the two eigenvectors associated to the principal curvatures, with eigenvalues $$k_1$$ and $$k_2$$. The other two values $$\sigma_1$$ and $$\sigma_2$$ describe how the normal tilts. This aspect is called flowline curvature.

All of this to say that I would search numerically the eigenvectors of the projected Hessian Matrix through an eigenvalues-eigenvector identification algorithm. For example in python you have the function eig of the library numpy. This would give you the eigenvectors of the principal direction.

• Can you explain the important points on those references? Sep 25 at 21:35
• edited with details Sep 26 at 9:04
• But the OP already knows the Gaussian and mean curvatures. The principal curvatures can be computed solving a 2 by 2 system. Sep 26 at 13:53
• The problem in the FEM context is that the vector field of normal vectors is not differentiable. Sep 27 at 19:16
• Ok, I did not know. I supposed we are in a differentiable normal field case, hence not in a FEM framework. It was not clear from the question (or maybe it is me that another time did not understand it?) Sep 27 at 19:58

You could also try an approach I have outlined below based on mimicking the discrete analogues of the Gauss map, the shape operator and its link to the principle directions and principle curvatures.

Step 1. At a given point or cell $$i_0$$ from the grid, take the unit normal $$\vec{n}_{i_0}$$. Then determine the set $$N(i_0)$$ of all vertices and cells immediately adjacent to $$i_0$$.

Step 2. For each $$j \in N(i_0)$$ there is the unique unit normal vector $$\vec{n}_{j}$$ of $$j$$. Thus, you obtain the set of all unit normals $$\vec{n}_j$$ where $$j \in N(i_0)$$.

Step 3. For each each unit normal $$\vec{n}_j$$, where $$j \in N(i_0)$$, define the orthogonal projection $$P_{i_0} \vec{n}_j \, =\, \vec{n}_j \, -\, (\vec{n}_j \cdot \vec{n}_{i_0})\, \vec{n}_{i_0}$$ onto the plane perpendicular to $$\vec{n}_{i_0}$$, which can be interpreted as the tangent plane of the surface at $${i_0}$$.

Step 4. Interpreting all projected 2D vector $$P_{i_0} \vec{n}_j \, : \, j \in N(i_{0})$$ as column vectors of size $$2 \times 1$$, define the $$2\times 2$$ symmetric matrix $$S_{i_{0}} \, =\, \sum_{j \, \in \, N(i_0)} \, \Big(\,P_{i_0} \vec{n}_j \Big)\, \Big(\,P_{i_0} \vec{n}_j \Big)^T$$

Step 5. Calculate the eigenvectors of $$S_{i_0}$$. These are your dsicrete principle directions at $$i_0$$. The eigenvalues would be the discrete principal curvatures at $$i_0$$.