A RhinoPython script which computes the principal curvature directions at any given vertex on a mesh.

meshCurvatureTensors.py

import Rhino
import scriptcontext
import System.Guid
#import System.Drawing
import math
from matrix import Matrix

def Torus(R, r, toroidal, poloidal):
    """Returns a parametrically constructs a torus mesh."""
    mesh = Rhino.Geometry.Mesh()
    for i in range(poloidal):
        for j in range(toroidal):
            theta = i * math.pi * 2 / poloidal
            phi = j * math.pi * 2 / toroidal
            x = (R + r * math.cos(phi)) * math.cos(theta)
            y = (R + r * math.cos(phi)) * math.sin(theta)
            z = r * math.sin(phi)
            mesh.Vertices.Add(x, y, z)
    n = toroidal * poloidal
    for i in range(poloidal):
        for j in range(toroidal):
            mesh.Faces.AddFace(
            ((i+1) * toroidal)%n + j,
            ((i+1) * toroidal)%n + (j+1)%toroidal,
            i*toroidal + (j+1)%toroidal,
            i*toroidal + j
            )
    mesh.Normals.ComputeNormals()
    mesh.Compact()
    return mesh
    
def CurvatureTensors(mesh):
    """Returns the curvature tensors for each vertex of a mesh."""
    tensors = [] # list of Vector3f
    for i in range(mesh.Vertices.Count):
        v = mesh.TopologyVertices.TopologyVertexIndex(i)
        position = mesh.Vertices[i]
        # get vertex normal as a matrix
        normal = mesh.Normals[i]
        Nvi = Matrix([[normal.X], [normal.Y], [normal.Z]])
        # get sorted 1-ring
        ring = mesh.TopologyVertices.ConnectedTopologyVertices(v, True) # CW-sort (not CCW!)
        # calculate face weightings, wij
        wij = []
        n = len(ring)
        for j in range(n):
            vec0 = mesh.TopologyVertices[ring[(j+(n-1))%n]] - position
            vec1 = mesh.TopologyVertices[ring[j]] - position
            vec2 = mesh.TopologyVertices[ring[(j+1)%n]] - position
            # Assumes closed manifold
            # TODO: handle boundaries
            wij.append(0.5 * (Rhino.Geometry.Vector3f.CrossProduct(vec0, vec1).Length + 
                     Rhino.Geometry.Vector3f.CrossProduct(vec1, vec2).Length))
        wijSum = sum(wij)
        # calculate matrix, Mvi
        Mvi = Matrix([[0, 0, 0], [0, 0, 0], [0, 0, 0]])
        I = Matrix.identity(3)
        for j in range(n):
            vec = mesh.TopologyVertices[ring[j]] - position
            edgeAsMatrix = Matrix([[vec.X], [vec.Y], [vec.Z]])
            Tij = (I - (Nvi * Nvi.transpose())) * edgeAsMatrix
            Tij *= (1 / Tij.getFrobeniusNorm())
            kij = (Nvi.transpose() * 2 * edgeAsMatrix).getEntry(0, 0) / math.pow(vec.Length, 2)
            Mvi += Tij.multiply(Tij.transpose()).scalarMultiply((wij[j]/wijSum)*kij)
        # get eigenvalues and eigenvectors for Mvi
        evals, evecs = Mvi.eig()
        # replace eigenvector matrix with list of Vector3f
        evecs = [Rhino.Geometry.Vector3f(evecs.getEntry(0,i), evecs.getEntry(1,i), evecs.getEntry(2,i)) for i in range(3)]
        # scale eigenvectors by corresponding eigenvalues
        [e.Unitize() for e in evecs]
        evecs = [Rhino.Geometry.Vector3f.Multiply(evecs[i], evals[i]) for i in range(3)]
        # sort by absolute value of eigenvalues (norm < min < max)
        # sortv: abs curvature, curvature, Vector3f dir
        sortv = zip([abs(i) for i in evals], evals, evecs)
        sortv.sort()
        # add min/max principal vectors to list
        tensors.append([sortv[1][-1], sortv[2][-1]])
    return tensors
    
def drawCurvatureTensors():
    """Draws curvature tensors for each vertex of a selected mesh."""
    # select mesh
    gp = Rhino.Input.Custom.GetObject()
    gp.SetCommandPrompt("GetPoint with options")
    gp.GeometryFilter = Rhino.DocObjects.ObjectType.Mesh
    # set up the options
    dblOption = Rhino.Input.Custom.OptionDouble(1.0, 0.01, 5)
    gp.AddOptionDouble("Scale", dblOption)
    scale = 1
    objref = None
    while True:
        # perform the get operation. This will prompt the user to
        # input a point, but also allow for command line options
        # defined above
        get_rc = gp.Get()
        if gp.CommandResult() != Rhino.Commands.Result.Success:
            return gp.CommandResult()
        if get_rc == Rhino.Input.GetResult.Object:
            objref = gp.Object(0)
            scale = dblOption.CurrentValue
        elif get_rc==Rhino.Input.GetResult.Option:
            continue
        break
    # get mesh
    mesh = objref.Mesh()
    # get curvature tensors
    tensors = CurvatureTensors(mesh)
    # set up layers (for coloring min/max lines)
    layers = []
    layers.append(scriptcontext.doc.Layers.Add("minimum", System.Drawing.Color.Blue))
    layers.append(scriptcontext.doc.Layers.Add("maximum", System.Drawing.Color.Red))
    # find layer if already created
    if layers[0] == -1: layers[0] = scriptcontext.doc.Layers.Find("minimum", True)
    if layers[1] == -1: layers[1] = scriptcontext.doc.Layers.Find("maximum", True)
    # draw curvature tensors
    for i in range(mesh.Vertices.Count):
        #t = tensors[i]
        # scale lines
        t = [Rhino.Geometry.Vector3f.Multiply(tensors[i][j], scale) for j in range(2)]
        position = mesh.Vertices[i]
        for j in range(2): # 0 - min, 1 - max
            # set layer
            scriptcontext.doc.Layers.SetCurrentLayerIndex(layers[j], False)
            # create lines
            line1 = Rhino.Geometry.Line(position, t[j])
            t[j].Reverse()
            line2 = Rhino.Geometry.Line(position, t[j])
            # add lines
            scriptcontext.doc.Objects.AddLine(line1)
            scriptcontext.doc.Objects.AddLine(line2)
    # redraw and return
    scriptcontext.doc.Views.Redraw()
    return Rhino.Commands.Result.Success
 
if __name__ == "__main__":
    #torus = Torus(3, 2, 30, 30)
    #scriptcontext.doc.Objects.AddMesh(torus)
    #CurvatureTensors(torus)
    drawCurvatureTensors()