Computation of the tensor of curvature on surface mesh

Question

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.

Answer 1

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: \begin{equation} \mathbf{n}=-\frac{\nabla \Phi}{|\nabla \Phi|} \end{equation} The curvature information is contained in the 3x3 matrix $\nabla{\bf{n}^T}$. Considering the Hessian matrix: \begin{equation} \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} \end{equation} 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: \begin{equation} \bf{G}=\nabla{\bf{n}}^T\bf{P} \end{equation} 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: \begin{equation} \nabla {\mathbf{n}^T}=-\frac{1}{|{\nabla \Phi}|}(\bf{P}\bf{H}) \end{equation} 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: \begin{equation} \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} \end{equation} $\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.

Answer 2

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