Estimating Principal Curvature Directions on Discrete Surfaces

Question

I am continuing with working on implementing curvature estimations on triangulated surfaces (See: here), using the algorithm of : Meyer, M., Desbrun, M., Schröder, P., & Barr, A. H. (2003). Discrete differential-geometry operators for triangulated 2-manifolds. In Visualization and mathematics III (pp. 35-57). Springer Berlin Heidelberg.

After getting to a stage where I can use Mathematica to calculate mean and Gaussian curvature, I am now trying to estimate principal curvatures and their directions on a triangulated surface.

The estimation of the values of principal curvatures is straightforward from the values of Gaussian and mean curvature, however the estimation of the principal directions seems to be a bit more involved. The goal is to find at each vertex i on the mesh, the local symmetric curvature tensor $B=\left(\begin{smallmatrix}B_1&B_2\\B_2&B_3\end{smallmatrix}\right)$ (defined with respect to a 2D coordinate system tangent to i), that best describes the normal curvatures calculated on all edges attached to vertex i.

As far as I understand the algorithm of Meyer et al., I need to do the following steps at each vertex on a mesh.

1. Calculate a local 2D orthonormal coordinate system that is tangent to the surface at the vertex of interest. This is done by using a) the local normal vector b) some arbitrary choice of edge to one of the nearest neighbours which is then projected into the tangent plane at the vertex of interest, and then c) a vector perpendicular to both.

2. Estimate the normal curvature along each of the edges surrounding the vertex. This can be done directly from the edge vector and the normal vector.

3. Project all edges into the local 2D coordinate system and use this to calculate normal curvatures as a function of the parameters $B_1,B_2,B_3$.

4. Calculate the weighted sum of the square of the differences between normal curvatures along each edge, and the curvatures calculated as a function of $B_1,B_2,B_3$. Weights are related to the area around in each vertex and the angles of the triangles.

5. Minimise this weighted sum to find the values of $B_1,B_2,B_3$ under the constraints that the sum of $B_1$ and $B_2$ gives the mean curvature previously calculated, and $B_1 B_3-B_2^2$ gives the Gaussian curvature also previously calculated.

I have tried implementing this in the following code:

For completeness I also add my code from the previous post (See: here, see especially J.M.s implementation which is quite a bit faster than mine):

Estimation of mean and gaussian curvatures on the mesh

(*Angles at each triangle*)

va = VectorAngle[#1 - #2, #3 - #2] & @@@ 
     Partition[#, 3, 1, {2, -2}] & /@ 
   MeshPrimitives[mesh, {2}][[All, 1]];
(*List of labels of triangles and positions in the list at which the \
vertices are obtuse*)
obttrilist = Position[va, n_ /; n > \[Pi]/2]; 
(*Coordinates of Vertices on Mesh*)
mc = MeshCoordinates[mesh];
(*Number of vertices*)
nvert = MeshCellCount[mesh, 0];
(*Number of edges*)
nedge = MeshCellCount[mesh, 1];
(*Number of faces*)
nfaces = MeshCellCount[mesh, 2];
(*List of Edges, consisting of a list of pairs of vertex numbers*)

edgeList = MeshCells[mesh, 1][[All, 1]];
(*List of Triangles consisting of a list of triples of vertex numbers*)

triangleList = MeshCells[mesh, 2][[All, 1]];
(*Triangle Areas*)

areaTriangles = PropertyValue[{mesh, 2}, MeshCellMeasure]; 
(*Length of Edges*)

edgeLengths = PropertyValue[{mesh, 1}, MeshCellMeasure];
(*Positions of vertex i in the edge list (*SLOW*), Note this gives \
the edge index and either 1 or 2 depending on the order inside that \
edge*)
posinEdgeList = Position[edgeList, #] & /@ Range[nvert];
(*Positions of vertex i in the triangle list (*SLOW*), Note this \
gives the triangle index and either 1, 2 or 3 depending on the order \
inside that triangle *)

posinTriangleList = Position[triangleList, #] & /@ Range[nvert];
(*Number of nearest neighbours at each vertex*)

nearestneighbourList = Length[posinEdgeList[[#]]] & /@ Range[nvert];
(*Function that calculates for a given pair of vertex indices from a \
line i,j, what triangles in the mesh indices also contain these \
indices, output is the triangle index*)

trilistfunc[line_] := 
  Intersection[posinTriangleList[[line[[1]]]][[All, 1]], 
   posinTriangleList[[line[[2]]]][[All, 1]]];
(*List of triangle indices that are attached to each line, This means \
that trianglesAtLines[[k]] will return the indices of the triangles \
containing the line k (If only one index is returned we are on the \
boundary!)*)
trianglesAtLines = Map[trilistfunc, edgeList];
(*List of indices of edges that are on the boundary*)

boundaryedges = 
  Flatten[Position[Length[trianglesAtLines[[#]]] & /@ Range[nedge], 
    1]];
(*List of indices of vertices that are on the boundary*)

boundaryvertices = 
  Flatten[edgeList[[#]] & /@ boundaryedges] // DeleteDuplicates;
(*Function that calculates which vertices in the attached triangles \
to a given edge are opposite to this edge, vertices are given as \
indices*)
oppcornerfunction[i_] := 
  If[MemberQ[boundaryedges, i], {0, 
    0}, {Cases[
      Cases[triangleList[[trianglesAtLines[[i, 1]]]], 
       Except[edgeList[[i, 1]]]], Except[edgeList[[i, 2]]]][[1]], 
    Cases[Cases[triangleList[[trianglesAtLines[[i, 2]]]], 
       Except[edgeList[[i, 1]]]], Except[edgeList[[i, 2]]]][[1]]}];
(*List of pairs of vertex indices m and n, that are opposite to edge \
k, pairs are ordered according to the edge number, if {0,0} then on \
boundary*)
oppcornerList = Map[oppcornerfunction, Range[nedge]];
(*Function that calculates the cotangents of the angles of the \
corners of "both" triangles opposite to a given edge ordered \
according to edge number (gives 0,0 for edge)*)

cotanglepairfunc[i_] := 
  If[MemberQ[boundaryedges, i], {0, 
    0}, {Cot[
     va[[trianglesAtLines[[i, 1]]]][[
      Position[triangleList[[trianglesAtLines[[i, 1]]]], 
        oppcornerList[[i, 1]]][[1, 1]]]]], 
    Cot[va[[trianglesAtLines[[i, 2]]]][[
      Position[triangleList[[trianglesAtLines[[i, 2]]]], 
        oppcornerList[[i, 2]]][[1, 1]]]]]}];
(*List of pairs of the cotangents of the opposite corner angles to \
each edge k*)
canglepairs = Map[cotanglepairfunc, Range[nedge]];
(*Function so we choose vertex j attached to vertex i*)

sw12func[a_] := If[a[[2]] == 1, 2, 1];
(*Function to calculate the list of oriented ij vectors attached to \
vertex i*)
ijvectfunc[i_] := 
  Map[(mc[[i]] - 
      mc[[edgeList[[posinEdgeList[[i, #, 1]], 
        sw12func[posinEdgeList[[i, #]]]]]]]) &, 
   Range[Length[posinEdgeList[[i]]]]];
(*List of oriented ij vectors attached to each vertex i *)

ijvectlist = Map[ijvectfunc, Range[nvert]];
(*Function to calculate the Mean curvature vector at each vertex*)

ijCVfunc[i_] := 
  Total[ijvectlist[[i]]*
    Map[Total[canglepairs[[posinEdgeList[[i, #, 1]]]]] &, 
     Range[Length[posinEdgeList[[i]]]]]];
(*List of Mean Curvature vectors at each vertex*)

ijCVlist = Map[ijCVfunc, Range[nvert]];
(*Now we need to calculate the Voronoi Area, modified such that \
obtuse triangles are taken into account*)
(*Create Function to \
calculate mixed Voronoi Area (see paper for explanation)*)

trianglecoords = Map[mc[[#]] &@triangleList[[#]] &, Range[nfaces]];

faceNormalfunc[tricoords_] := 
  Cross[tricoords[[1]] - tricoords[[2]], 
   tricoords[[3]] - tricoords[[2]]];

facenormals = Map[faceNormalfunc, trianglecoords];

mcnewcalc = 
  Map[Total[
     Map[(facenormals[[#]]*areaTriangles[[#]]) &, 
      posinTriangleList[[All, All, 1]][[#]]]] &, Range[nvert]];

meancvnew = -Sign[MapThread[Dot, {mcnewcalc, ijCVlist}]] (Norm /@ 
     ijCVlist);

areaMixedfunction[i_] := 
  If[MemberQ[boundaryvertices, i], 0, Total[Map[Do[

       edgenumber = posinEdgeList[[i, #, 1]];

       d1 = trianglesAtLines[[edgenumber]][[1]];
       d2 = trianglesAtLines[[edgenumber]][[2]];
       AMixed = 0;
       If[MemberQ[obttrilist[[All, 1]], d1],
        (*Now do test to see which triangle area we add correcting \
for whether the triangles are obtuse or not*)

        ObtVnum = Position[obttrilist[[All, 1]], d1][[1, 1]];

        (*Vertex index of the obtuse part of the triangle*)

        Vnum = triangleList[[obttrilist[[ObtVnum, 1]], 
          obttrilist[[ObtVnum, 2]]]];
        If[Vnum == i,
         (*Triangle Obtuse at i, part of area T/2*)

         AMixed += (1/4)*areaTriangles[[d1]];
         ,
         (*Triangle Obtuse but not at i, area T/4*)

         AMixed += (1/8)*areaTriangles[[d1]];
         ]
        ,
        AMixed += (1/8)*(canglepairs[[edgenumber]][[
            1]])*(edgeLengths[[edgenumber]])^2
        (*If False we add the normal voronoi*)
        ];
       (*Repeat the test for the other triangle*)

       If[MemberQ[obttrilist[[All, 1]], d2],
        (*Now do test to see which triangle area we add*)

        ObtVnum = Position[obttrilist[[All, 1]], d2][[1, 1]];
        Vnum = 
         triangleList[[obttrilist[[ObtVnum, 1]], 
          obttrilist[[ObtVnum, 2]]]];;
        If[Vnum == i,
         (*Triangle Obtuse at i, therefore  area T/2*)

         AMixed += (1/4)*areaTriangles[[d2]];
         ,
         (*Triangle Obtuse but not at i, therefore add half of area T/
         4*)
         AMixed += (1/8)*areaTriangles[[d2]];
         ]
        ,
        AMixed += (1/8)*(canglepairs[[edgenumber]][[
            2]])*(edgeLengths[[edgenumber]])^2
        (*If False we add the normal voronoi*)
        ];
       Return[AMixed]
       , 1] &, Range[Length[posinEdgeList[[i]]]]]]];
(*Create a list of the Mixed areas per vertex*)

AmixList = Map[areaMixedfunction, Range[nvert]];
(*Gaussian Curvature*)

gaussCurv = 
  Map[If[MemberQ[boundaryvertices, #], 
     0, (2*\[Pi] - 
        Total[Extract[va, 
          Position[MeshCells[mesh, 2][[All, 1]], #]]])/
      AmixList[[#]]] &, Range[nvert]];
(*Mean Curvature*)

meanCurv = 
  Map[If[MemberQ[boundaryvertices, #], 
     0, (meancvnew[[#]]/AmixList[[#]])/4] &, Range[nvert]];

Estimation of principal curvature directions on the mesh

(*List of normalised normal vectors at each vertex (Vectors set to \
{0,0,0} on boundaries)*)

Nvectors = 
  Map[If[MemberQ[boundaryvertices, #], {0, 0, 
      0}, (mcnewcalc[[#]]/Norm[mcnewcalc[[#]]])] &, Range[nvert]];
(*Function to calculate weighting factor(s) for vertex i*)

wijfunc[i_] := 
  Map[(1/8)*(canglepairs[[posinEdgeList[[i, #, 1]]]][[1]] + 
       canglepairs[[posinEdgeList[[i, #, 1]]]][[2]])*
     Norm[ijvectlist[[i, #]]]^2 &, 
   Range[Length[posinEdgeList[[i]]]]];
(*Weighting factors for each edge *)

wij = Map[If[MemberQ[boundaryvertices, #], 0, wijfunc[#]] &, 
   Range[nvert]];
(*Calculate first the local orthonormal coordinate system*)
(*Use the \
first vector in the nearest neighbour vector list as the first vector \
coordinate*)
(*Calculate projection of this vector in the normal \
plane and normalise (Boundaries set to {0,0,0})*)

tgtdirn = 
  Map[If[MemberQ[boundaryvertices, #], {0, 0, 
      0}, (ijvectlist[[#, 
         1]] - (ijvectlist[[#, 1]].Nvectors[[#]]) Nvectors[[#]])] &, 
   Range[nvert]];
(*Normalise tgt vector to give one of the basis vectors in the \
tangent plane*)

dij1 = Map[
   If[MemberQ[boundaryvertices, #], {0, 0, 
      0}, (tgtdirn[[#]]/Norm[tgtdirn[[#]]])] &, Range[nvert]];
(*Calculate perpendicular vector in tangent plane and normalise*)

perpdirn = 
  Map[If[MemberQ[boundaryvertices, #], {0, 0, 
      0}, (Cross[dij1[[#]], Nvectors[[#]]])] &, Range[nvert]];
(*Normalise perpendicular vector to give the otherbasis vectors in \
the tangent plane*)

dij2 = Map[
   If[MemberQ[boundaryvertices, #], {0, 0, 
      0}, (perpdirn[[#]]/Norm[perpdirn[[#]]])] &, Range[nvert]];
(*Now we have an orthonormal coordinate system at each vertex on the \
surface (except for boundaries)*)

(*Function to estimate  normal curvature along each edge ij attached \
to vertex i*)
\[Kappa]ijnfunc[i_] := 
  Map[2*ijvectlist[[i, #]].Nvectors[[i]]/(Norm[
         ijvectlist[[i, #]]])^2 &, 
   Range[Length[posinEdgeList[[i]]]]];
(*Estimates of normal curvatures *)
\[Kappa]ijn = 
  Map[If[MemberQ[boundaryvertices, #], 0, \[Kappa]ijnfunc[#]] &, 
   Range[nvert]];
(*Function to calculate projections of edges into the tangent plane*)

ijprojfunc[i_] := 
  Map[ijvectlist[[i, #]] - (-ijvectlist[[i, #]].Nvectors[[i]]) \
Nvectors[[i]] &, Range[Length[posinEdgeList[[i]]]]];
(*List of Projections of edge vectors in the tangent plane*)

ijprojlist = Map[ijprojfunc[#] &, Range[nvert]];
(*Calculation of the weighted square of the differences between \
estimated curvatures and those calculated using the curvature tensor \
(to be found)*)

Efunc[i_] := 
  Map[wij[[i, #]]*({(ijprojlist[[i, #]]).dij1[[
            i]], (ijprojlist[[i, #]]).dij2[[i]]}.{{B1, B2}, {B2, 
           B3}}.{ijprojlist[[i, #]].dij1[[i]], 
          ijprojlist[[i, #]].dij2[[i]]} - \[Kappa]ijn[[i, #]])^2 &, 
   Range[Length[posinEdgeList[[i]]]]];
(*Functions to be minimized*)

Evals = Map[
   If[MemberQ[boundaryvertices, #], 0, Simplify[Total[Efunc[#]]]] &, 
   Range[nvert]];
(*Calculation of eigenvectors of the local curvature tensor*)

eigs = Map[
   If[MemberQ[boundaryvertices, #], 
     0, (Eigenvectors[{{B1, B2}, {B2, B3}} /. 
        NMinimize[{Evals[[#]], {B1 + B2 == 2*meanCurv[[#]], 
            B1*B3 - B2^2 == gaussCurv[[#]]}}, {B1, B2, B3}][[2]]])] &,
    Range[nvert]];
(*Calculation of the principal curvature direction vectors*)

principaldirectionvectors = 
  Map[If[MemberQ[boundaryvertices, #], {0, 0, 
      0}, {eigs[[#, 1, 1]]*dij2[[#]] + eigs[[#, 1, 2]]*dij1[[#]], 
      eigs[[#, 2, 1]]*dij2[[#]] + eigs[[#, 2, 2]]*dij1[[#]]}] &, 
   Range[nvert]];

The code (at least) runs. For the surface defined by:

P[m_, d_] := (1 - 2 m) d^2/(2 m );
Q[m_, d_, a_] := (1 - m) d^4/16 a^2 m;
R[m_, d_] := (m - 1) d^4/16 m;
eryth[m_, d_, a_] := (x^2 + y^2)^2 + P[m, d] (x^2 + y^2) + 
   Q[m, d, a] z^2 + R[m, d];
(*Larkin,T.J.& Kuchel,P.W.Bull.Math.Biol.(2010) 72:1323. Mathematical \
Models of Naturally "Morphed" Human Erythrocytes:Stomatocytes and \
Echinocytes.*)
mesh = 
 BoundaryDiscretizeRegion[
  ImplicitRegion[eryth[0.7, 10, 0.6] <= 0, {x, y, z}], 
  MaxCellMeasure -> {"Length" -> 0.5}]

we get the output

The code I used for plotting is:

tubelengthf = 5;
Lineplot4 = Table[{}, {j, 1, nvert}];
Lineplot5 = Table[{}, {j, 1, nvert}];

For[j = 1, j < nvert + 1, j++, 
 Lineplot4[[j]] = 
  Tube[{mc[[j]] - principaldirectionvectors[[j, 1]]/tubelengthf, 
    mc[[j]] + principaldirectionvectors[[j, 1]]/tubelengthf}, 0.01]]
For[j = 1, j < nvert + 1, j++, 
 Lineplot5[[j]] = 
  Tube[{mc[[j]] - principaldirectionvectors[[j, 2]]/tubelengthf, 
    mc[[j]] + principaldirectionvectors[[j, 2]]/tubelengthf}, 0.01]]

v = MeshCoordinates[mesh];
w = MeshCells[mesh, 2];
v1 = Graphics3D[{EdgeForm[], GraphicsComplex[v, w]}];
v2 = Graphics3D[{{Black, Lineplot4}, {Black, Lineplot5}}, {EdgeForm[],
     GraphicsComplex[v, w]}];
Show[v1, v2]

The code is somewhat slow, and this seems to be due to the numerical minimisation at each vertex on the surface, but the code is still testable for reasonably sized meshes. What can be seen however in the plot is that although I am getting in some locations a reasonable estimate, this is not very well defined over the entire surface, where we know from symmetry that the principal curvature directions must be aligned along the "longitudes and latitudes".

According to "Garimella R.V., Swartz B.K.: Curvature Estimation for Unstructured Triangulations of Surfaces. Tech. Rep. LAUR-03-8240, Los Alamos National Laboratory , Nov. 2003. "...the Meyer, Desbrun, Schröder, Barr algorithm doesn't perform as well as approaches which use NURBS-style surfaces that are fitted to the triangulated surface and then used to calculate local curvature tensors." Software such as Rhino can be used for this, however I would like to do this in Mathematica.

After this rather long preamble I finally come to the question.

If this poor estimation of principal curvature directions is due to the limitations of the algorithm (and not mistakes in my coding) does anyone have any hints of other ways of calculating in Mathematica the full curvature tensor?

Any suggestions would be most welcome.

Answer 1

For this answer, I shall be doing something slightly more ambitious. In particular, I will be computing the so-called curvature tensor, which encodes information on the normal vector $\mathbf n$, the principal curvatures $\kappa_1,\kappa_2$, and the principal directions $\mathbf v_1,\mathbf v_2$ as a symmetric matrix. More concretely, a curvature tensor $\mathbf E$ has the eigenvalues $\kappa_1,\kappa_2,0$, with corresponding eigenvectors $\mathbf v_1,\mathbf v_2,\mathbf n$. (Of course, one can then compute the Gaussian and mean curvature from the principal curvatures.)

Okay, I lied a bit in the first paragraph. What I will actually be computing is the Taubin curvature tensor $\mathbf M$, which has the same set of eigenvectors, but instead has the corresponding eigenvalues $\frac{3\kappa_1+\kappa_2}8,\frac{\kappa_1+3\kappa_2}8,0$, using the estimation procedure presented in the linked paper.

As with my previous answer, the following computations are only applicable for closed meshes; they will need to be modified for meshes with a boundary.

Again, I'll be using this algebraic surface with the symmetries of a dodecahedron:

dodeq = z^6 - 5 (x^2 + y^2) z^4 + 5 (x^2 + y^2)^2 z^2 - 2 (x^4 - 10 x^2 y^2 + 5 y^4) x z +
        (x^2 + y^2 + z^2)^3 - (x^2 + y^2 + z^2)^2 + (x^2 + y^2 + z^2) - 1;

dod = BoundaryDiscretizeRegion[ImplicitRegion[dodeq < 0, {x, y, z}],
                               MaxCellMeasure -> {"Length" -> 0.1}]

Extract the vertices, triangles, and neighboring vertex indices:

pts = MeshCoordinates[dod];
tri = MeshCells[dod, 2] /. Polygon[p_] :> p;

nbrs = Table[DeleteDuplicates[Flatten[List @@@ First[FindCycle[
             Extract[tri, Drop[SparseArray[Unitize[tri - k],
                                           Automatic, 1]["NonzeroPositions"], None, -1],
                     # /. {k, a_, b_} | {b_, k, a_} | {a_, b_, k} :> (a -> b) &]]]]],
             {k, Length[pts]}];

We then use Max's method to estimate the vertex normals:

nrms = Table[Normalize[Total[With[{c = pts[[k]], vl = pts[[#]]}, 
                                  Cross[vl[[1]] - c, vl[[2]] - c]/
                                  ((#.# &[vl[[1]] - c]) (#.# &[vl[[2]] - c]))] & /@
                             Partition[nbrs[[k]], 2, 1, 1],
                             Method -> "CompensatedSummation"]], {k, Length[pts]}];

Here then is the Taubin curvature tensor computation:

ctl = Table[With[{v = pts[[k]], n = nrms[[k]], nl = nbrs[[k]]}, 
                 Normalize[With[{vl = pts[[#]]},
                                (Norm[Cross[vl[[1]] - v, vl[[2]] - v]] + 
                                 Norm[Cross[vl[[2]] - v, vl[[3]] - v]])/2] & /@ 
                           Partition[nl, 3, 1, -2], Total].
                 Table[With[{d = v - vj}, -2 d.n/d.d
                            (Outer[Times, #, #] &[Normalize[d - Projection[d, n]]])],
                       {vj, pts[[nl]]}]], {k, Length[pts]}];

We can then use Eigensystem[] to extract the principal directions:

pdl = (Pick[#2, Unitize[Chop[#1]], 1] & @@ Eigensystem[#]) & /@ ctl;

Visualize the principal directions on the surface:

With[{h = 0.02}, 
     Graphics3D[{GraphicsComplex[pts, {EdgeForm[], Polygon[tri]}, VertexNormals -> nrms], 
                 MapThread[{Line[{#1 - h #2[[1]], #1 + h #2[[1]]}], 
                            Line[{#1 - h #2[[2]], #1 + h #2[[2]]}]} &, {pts, pdl}]}, 
                Boxed -> False]]

Instead of showing the image generated by that last one, here is the result of a side-by-side comparison of the estimated and true (derived from the implicit Cartesian formula) principal directions:

As previously noted, one can also compute the Gaussian and mean curvatures from the estimated Taubin tensor as well:

gc = Composition[Apply[Times], AffineTransform[{{3, -1}, {-1, 3}}],
                 DeleteCases[0], Chop, Eigenvalues] /@ ctl;

mc = Composition[Mean, AffineTransform[{{3, -1}, {-1, 3}}], 
                 DeleteCases[0], Chop, Eigenvalues] /@ ctl;

---

Here are the corresponding results for the erythrocyte model featured in this paper. The mesh used for these pictures was generated using the following:

With[{m = 7/10, d = 10, a = 3/5},
     erythro = BoundaryDiscretizeRegion[ImplicitRegion[
                       (x^2 + y^2)^2 - (1 - 1/(2 m)) d^2 (x^2 + y^2) + 
                       z^2 (1 - m) d^4/16 a^2 m + (m - 1) d^4/16 m < 0, {x, y, z}], 
                       MaxCellMeasure -> {"Length" -> 0.4}]];

Compare the Taubin-estimated and true principal directions:

Compare the true and Taubin-estimated Gaussian (left) and mean (right) curvatures:

Answer 2

OK, at least here is an attempt to solve my problem. Hopefully these thoughts and code may be useful to others.

It seems like there is no single ideal algorithm to solve this problem. Some work better than others for certain qualities of meshes and others working better for other types of meshes (see Vasa et al. 2016 for a quantitative comparison).

There are a variety of methods available to do this (many are implemented in other codes and softwares (e.g. Rhino, MeshLab etc). It may be possible to use these other codes inside Mathematica, I have not explored this yet.

The original idea of this question was to see whether one could fit a smooth surface (NURBS for example) in the neighbourhood of a vertex and then use the analytical equations for the fitted surface to calculate the local curvature tensor, and hence the principal directions, mean curvature and gaussian curvature.

It seems that this is not so trivial. In Mathematica NURBS are implemented using BSplineFunction. It takes a rectangular array of control points as input, and does not fit through an arbitrary number of points as in Interpolation. Unfortunately using Interpolation for unstructured grid it does not seem to allow one to use the "Spline" option, which would be great in order to get out the parameters.

This has been covered somewhat already see e.g. (Produce a spline from a set of {{x, y}, z} points and get its parameters/expression) or (How to make BSplineFunction pass each data point and naturally smooth?), but I haven’t yet managed to work out a solution for non-regular meshes or meshes in 3D in which there is no rectangular $uv$ parametrisation of the surface.

In the process of trying to find a solution I cam across another option in which quadric surfaces are fitted to the local neighbourhood of the vertex of interest (See discussion in Petitjean, S., ACM Computing Surveys 2002, A Survey of Methods for Recovering Quadrics in Triangle Meshes). This is relatively straightforward to implement; however, it seems (in my implementation) to be a bit slow, especially for larger meshes and 2nd nearest neighbourhood implementations.

The algorithm in its simplest form has the following steps.

1. Pick a vertex of interest i
2. Calculate an estimate for the normal vector to the surface at this vertex
3. select a group of nearest neighbours to this vertex and rotate them into a local coordinate system such that the local z-axis is aligned with the normal vector at vertex i
4. Fit the quadric equation to the surface.
5. The parameters of this equation can then be used to directly calculate the local curvature tensor at each point

There are variations on this algorithm in which more and more nearest neighbours can be added to the fit, or to improve the fit by using the fit to create a new estimate for the normal vector which is then used to calculate a new fit. (see Petitjean 2002, for details).

In the following I have implemented this using just a simple quadric fit with two choices for the neighbourhood of vertices to be fitted (see below).

The code for the fit is as follows:

(*Simple quadric fitting model to estimate curvatures on surfaces*)

curvatureEstimatesimplequadric[mesh0_, nextnn_] :=

 Module[{mesh = mesh0},
  (*Coordinates of Vertices on Mesh*)
  mc = MeshCoordinates[mesh];
  (*Number of vertices*)
  nvert = MeshCellCount[mesh, 0];
  (*Number of edges*)
  nedge = MeshCellCount[mesh, 1];
  (*Number of faces*)
  nfaces = MeshCellCount[mesh, 2];
  (*List of Edges, consisting of a list of pairs of vertex numbers*)

   edgeList = MeshCells[mesh, 1][[All, 1]];
  (*List of Triangles consisting of a list of triples of vertex \
numbers*)
  triangleList = MeshCells[mesh, 2][[All, 1]];
  (*Triangle Areas*)

  areaTriangles = PropertyValue[{mesh, 2}, MeshCellMeasure]; 
  (*Length of Edges*)

  edgeLengths = PropertyValue[{mesh, 1}, MeshCellMeasure];
  (*Positions of vertex i in the edge list (*SLOW*), 
  Note this gives the edge index and either 1 or 2 depending on the \
order inside that edge*)

  posinEdgeList = Position[edgeList, #] & /@ Range[nvert];
  (*Positions of vertex i in the triangle list (*SLOW*), 
  Note this gives the triangle index and either 1, 
  2 or 3 depending on the order inside that triangle *)

  posinTriangleList = Position[triangleList, #] & /@ Range[nvert];
  (*Number of nearest neighbours at each vertex*)

  nearestneighbourList = 
   Length[posinEdgeList[[#]]] & /@ Range[nvert];

  (*Function that calculates for a given pair of vertex indices from \
a line i,j, 
  what triangles in the mesh indices also contain these indices, 
  output is the triangle index*)

  trilistfunc[line_] := 
   Intersection[posinTriangleList[[line[[1]]]][[All, 1]], 
    posinTriangleList[[line[[2]]]][[All, 1]]];
  (*List of triangle indices that are attached to each line, 
  This means that trianglesAtLines[[
  k]] will return the indices of the triangles containing the line k \
(If only one index is returned we are on the boundary!)*)

  trianglesAtLines = Map[trilistfunc, edgeList];
  (*List of indices of edges that are on the boundary*)

  boundaryedges = 
   Flatten[Position[Length[trianglesAtLines[[#]]] & /@ Range[nedge], 
     1]];
  (*List of indices of vertices that are on the boundary*)

  boundaryvertices = 
   Flatten[edgeList[[#]] & /@ boundaryedges] // DeleteDuplicates;
  (*Function that calculates which vertices in the attached triangles \
to a given edge are opposite to this edge, 
  vertices are given as indices*)

  (*Create Function to calculate mixed Voronoi Area (see paper for \
explanation)*)

  trianglecoords = Map[mc[[#]] &@triangleList[[#]] &, Range[nfaces]];

  faceNormalfunc[tricoords_] := 
   Cross[tricoords[[1]] - tricoords[[2]], 
    tricoords[[3]] - tricoords[[2]]];

  facenormals = Map[faceNormalfunc, trianglecoords];

  mcnewcalc = 
   Map[Total[
      Map[(facenormals[[#]]*areaTriangles[[#]]) &, 
       posinTriangleList[[All, All, 1]][[#]]]] &, Range[nvert]];

  (*List of normalised normal vectors at each vertex*)

  Nvectors = 
   Map[(mcnewcalc[[#]]/Norm[mcnewcalc[[#]]]) &, Range[nvert]];

  (*Function to give the vertex indices of all the nearest neighbours \
j attached to vertex i by edges ij*)

  nneighbindexes[i_] := 
   Cases[Flatten[Map[edgeList[[#]] &, posinEdgeList[[i]][[All, 1]]]], 
    Except[i]];
  nextnneighbourindexes[i_] := 
   DeleteDuplicates[
    Flatten[Map[nneighbindexes[#] &, nneighbindexes[i]]]];
  (*List of points to be fitted around vertex i*)

  ptstofit[i_] := 
   If[nextnn == 1, Join[{mc[[i]]}, Map[mc[[#]] &, nneighbindexes[i]]],
     Join[{mc[[i]]}, Map[mc[[#]] &, nextnneighbourindexes[i]]]];
  (*The following calculates on next nearest neighbours (need to \
introduce code though inside this module?)*)
  (*ptstofit[i_]:=
  Join[{mc[[i]]},Map[mc[[#]]&,nextnneighbourindexes[
  i]]];*)
  (*calculation of points to fit in a rotated coordinate \
system aligned with the estimated normal and translated such that \
vertex i is at the origin*)

  localcoordpointslist[i_] := 
   Map[RotationMatrix[{Nvectors[[i]], {0, 0, 1}}] .(# - mc[[i]]) &, 
    ptstofit[i]];
  lmmodelfit[i_] := 
   LinearModelFit[localcoordpointslist[i], {1, x^2, x y, y^2}, {x, y}];
  lmmodelfits = 
   Map[If[MemberQ[boundaryvertices, #], 0, lmmodelfit[#]] &, 
    Range[nvert]];
  gaussc = 
   Map[If[MemberQ[boundaryvertices, #], 
      0, (4*lmmodelfits[[#]]["BestFitParameters"][[2]]*
         lmmodelfits[[#]]["BestFitParameters"][[4]] - 
        lmmodelfits[[#]]["BestFitParameters"][[3]]^2)] &, 
    Range[nvert]];
  meanc = 
   Map[If[MemberQ[boundaryvertices, #], 
      0, (lmmodelfits[[#]]["BestFitParameters"][[2]] + 
        lmmodelfits[[#]]["BestFitParameters"][[4]])] &, 
    Range[nvert]];

  eigenveccalcs = 
   Map[If[MemberQ[boundaryvertices, #], 0, 
      Eigenvectors[{{2*lmmodelfits[[#]]["BestFitParameters"][[2]], 
         lmmodelfits[[#]]["BestFitParameters"][[
          3]]}, {lmmodelfits[[#]]["BestFitParameters"][[3]], 
         2*lmmodelfits[[#]]["BestFitParameters"][[4]]}}]] &, 
    Range[nvert]];
  ev1 = Map[
    If[MemberQ[boundaryvertices, #], {0, 0, 0}, 
      RotationMatrix[{{0, 0, 1}, 
         Nvectors[[#]]}].(eigenveccalcs[[#]][[1, 1]]*{1, 0, 0} + 
         eigenveccalcs[[#]][[1, 2]]*{0, 1, 0})] &, Range[nvert]];
  ev2 = Map[
    If[MemberQ[boundaryvertices, #], {0, 0, 0}, 
      RotationMatrix[{{0, 0, 1}, 
         Nvectors[[#]]}].(eigenveccalcs[[#]][[2, 1]]*{1, 0, 0} + 
         eigenveccalcs[[#]][[2, 2]]*{0, 1, 0})] &, Range[nvert]];
  (*Perhaps do this in the eigensystem to speed up*)

  evals = Map[
    If[MemberQ[boundaryvertices, #], {0, 0}, 
      Eigenvalues[{{2*lmmodelfits[[#]]["BestFitParameters"][[2]], 
         lmmodelfits[[#]]["BestFitParameters"][[
          3]]}, {lmmodelfits[[#]]["BestFitParameters"][[3]], 
         2*lmmodelfits[[#]]["BestFitParameters"][[4]]}}]] &, 
    Range[nvert]];


  {nvert, mc, triangleList, meanc, gaussc, ev1, ev2, evals}
  ]

The results for simple surfaces like those defined above are:

For the quadric fitting there are two examples one with fitting only the first nearest neighbours and the other with the second nearest neighbours

Interestingly for the mean and gaussian curvatures the MDSB algorithm seems better but for the principal curvature directions we get much better results from the quadric fitting:

I will be interested to see what other solutions are out there and at the very least I hope this is helpful to someone else.