Is there a way to calculate mean curvature of a triangular mesh?

Question

So I've written a script which does this:

Before:

After:

The code is:

#fname - filename of imported .stl
#thic - required thickness of output model
#txt - ascii or binary stl
#cyc - cycles of smoothing
#cut - times triangle divided
#fac - smoothing factor
#per - smoothing reps
#bbX - bounding box X
#bbY - bounding box Y
#bbZ - bounding box Z
#trX, trY, trZ - count of elements in array
def s_surface(fname, thic, txt, cyc, cut, fac, rep, bbX, bbY, bbZ, trX, trY, trZ ):
    bb = [[(bbX/2, 0, 0), (1,0,0), True, False], 
          [(-1*bbX/2, 0, 0), (1,0,0), False, True],
          [(0, bbY/2, 0), (0,1,0), True, False], 
          [(0,-1*bbY/2, 0), (0,1,0), False, True],
          [(0, 0, bbZ/2), (0,0,1), True, False], 
          [(0, 0, -1*bbZ/2), (0,0,1), False, True]
         ]
    tr = [trX, trY, trZ]
    bpy.ops.import_mesh.stl(filepath=fname)
    ob_new = bpy.context.selected_objects[0]
    bpy.context.scene.objects.active = ob_new
    bpy.ops.object.origin_set(type='GEOMETRY_ORIGIN')
    bpy.ops.object.mode_set(mode='EDIT')
    bpy.ops.mesh.select_all(action='SELECT')
    bpy.ops.mesh.normals_make_consistent(inside=False)
    bpy.ops.mesh.remove_doubles(threshold=0.02)
    for i in range(cyc):
        bpy.ops.mesh.subdivide(number_cuts=cut)
        bpy.ops.mesh.vertices_smooth(factor=fac, repeat=rep, xaxis=True, yaxis=True, zaxis=True)
    bpy.ops.mesh.normals_make_consistent(inside=False)
    for i in range(len(bb)):
        bpy.ops.mesh.bisect(plane_co=bb[i][0],plane_no=bb[i][1], clear_outer=bb[i][2],clear_inner=bb[i][3]) 
        bpy.ops.mesh.select_all(action='SELECT')
    bpy.ops.mesh.normals_make_consistent(inside=False)
    bpy.ops.object.mode_set(mode='OBJECT')
    for i in range(3):
        bpy.ops.object.origin_set(type='GEOMETRY_ORIGIN')
        bpy.ops.object.modifier_add(type='ARRAY')
        bpy.context.object.modifiers["Array"].count = tr[i]
        bpy.context.object.modifiers["Array"].relative_offset_displace[0] = 0
        bpy.context.object.modifiers["Array"].relative_offset_displace[i] = 1
        bpy.context.object.modifiers["Array"].use_merge_vertices = True
        bpy.context.object.modifiers["Array"].merge_threshold = 0.01
        bpy.ops.object.modifier_apply(apply_as='DATA', modifier="Array")
    bpy.ops.object.origin_set(type='GEOMETRY_ORIGIN')
    bpy.ops.object.modifier_add(type='SOLIDIFY')
    bpy.context.object.modifiers["Solidify"].thickness = thic
    bpy.context.object.modifiers["Solidify"].use_quality_normals = True
    bpy.context.object.modifiers["Solidify"].use_even_offset = True
    bpy.ops.object.modifier_apply(apply_as='DATA', modifier="Solidify")
    nfname = ""
    nfname = fname[0:-4] + ".stl"
    bpy.ops.export_mesh.stl(filepath=nfname, ascii = txt)

Surfaces I am trying to smooth are supposed to be triply periodic minimal surfaces, so it means that after smoothing is applied mean curvature is meant to be 0 in all vertexes of a mesh. How can I calculate mean curvature to check if smoothing is done right?

Answer 1

This is an answer based on this Computer Graphics SE indicated in the comments (the one by Nathan Reed).

The math is described in this answer, but in short the curvature is calculated by vertex so:

• Get all edges from this vertex
• And for each edge compare the projection along the edge of the normals at it extremities
• Take the mean of all that

Now, as didn't want to just code a copy of it, I propose a little enhancement:

The calculation is also vertex centered, but we go along each of its surrounding face (need to be triangulated), to calculate a curvature weighted by the angle of this face at the vertex we are calculating around.

Once the sum done over all faces, we mean this sum by the total angle around the vertex.

Making so, the curvature of a large angle will be more important than for a small angle. This is like taking the integral of the curvatures around the vertex.

I think this may be close to the theoretical calculus as it consist of taking the curvature for each cutting plane turning around the vertex normal.

In both algorithms, the base calculus is the same and this corresponds to the CGSE answer:

def curvature_along_edge( vert, other ):
    normal_diff = other.normal - vert.normal
    vert_diff = other.co - vert.co
    return normal_diff.dot( vert_diff ) / vert_diff.length_squared

But in the case we want to use angles, we need to organize a surrounding loop (ring) over the vertex for which we want the curvature, to have the several angles with the good signs.

For instance here:

We want the ring/loop:

• For vertex 0: 5, 8, 4
• For vertex 8: 1, 6, 3, 7, 2, 4, 0, 5, (1)

etc.

These loops are counterclockwise, as vertices of Blender's polygons are.

This is done by:

# Get vertices in the face order but starting from a given vert
def following_verts_of_vert( vert, face ):
    i0 = index_of( vert, face.verts )
    i1 = (i0 + 1) % 3
    i2 = (i0 + 2) % 3
    return face.verts[i0], face.verts[i1], face.verts[i2]

# Create the oriented ring around vert
def ring_from_vert( vert ):
    vertices = []
    for face in vert.link_faces:
        i0, i1, i2 = following_verts_of_vert( vert, face )
        vertices.append( [i1, i2] )
    result = vertices[0]    
    added = True
    while added and len(vertices):
        added = False
        prev = search_link( result[0], vertices, 1 )
        if prev:
            result = [prev[0]] + result
            vertices.remove( prev )
            added = True
        next = search_link( result[-1], vertices, 0 )
        if next and next[1] not in result:
            result.append( next[1] )
            vertices.remove( next )
            added = True
    return result

(this code above is not really needed, or could be optimized, but was my first though to have a constant orientation reference for the curvature)

So that finally, the mean curvature around a vertex is calculated by:

def angle_between_edges( vert, other1, other2 ):
    edge1 = other1.co - vert.co
    edge2 = other2.co - vert.co
    product = edge1.cross( edge2 )
    sinus = product.length / (edge1.length * edge2.length)
    return asin( min(1.0, sinus) )

def mean_curvature_vert( vert ):
    ring = ring_from_vert( vert )
    ring_curvatures = [curvature_along_edge( vert, other ) for other in ring]
    total_angle = 0.0
    curvature = 0.0
    for i in range(len(ring)-1):
        angle = angle_between_edges( vert, ring[i], ring[i+1] )
        total_angle += angle
        curvature += angle * (ring_curvatures[i] + ring_curvatures[i+1])

    return curvature / (2.0 * total_angle)

Here is the compared results:

They are close, but the comparison I've used needs some explanation: the resulting curvatures are normalized so that they fit to the interval for vertex groups [0, 1]:

In short, this vertex group assignment compares the contrasts but not the values.

This is done so:

def assign_to_vertex_group( obj, group_name, curvatures ):
    vertex_group = ensure_vertex_group( obj, group_name )

    curvatures = [abs(c) for c in curvatures]

    min_curvature = min( curvatures )
    max_curvature = max( curvatures )
    vg_fac = 1.0 / (max_curvature - min_curvature) if max_curvature != min_curvature else 1.0

    for i, vert in enumerate( obj.data.vertices ):
        vertex_group.add( [vert.index], (curvatures[i] - min_curvature) * vg_fac, 'REPLACE' )

Here is the blend file with the 2 scripts named : cgseSimple and cgseSimple2.

Note:

Also included Mike's implementation (mike text) + the same rewritten a bit (mikesPaper), for understanding and comparison purpose. But from all that, it is hard to know what is the more accurate (how to determinate which is true?). I was not able to understand how this algo is able to determinate curvature orientation (convex vs. concave inclinations).

Also tested all that on big meshes (not included due to the blend file size limit). But that accentuate the diff between the different approaches.

Answer 2

So I found an article on computing curvatures of triangular meshes.

And wrote this piece:

import bpy
import math
import mathutils 

def create_pairs(k):
    tr = []
    bpy.ops.object.mode_set(mode='OBJECT')
    obj = bpy.context.active_object
    bpy.ops.object.mode_set(mode='EDIT')
    bpy.ops.mesh.select_mode(type="VERT")
    bpy.ops.mesh.select_all(action='DESELECT')
    bpy.ops.object.mode_set(mode='OBJECT')
    obj.data.vertices[k].select = True
    bpy.ops.object.mode_set(mode='EDIT')
    bpy.ops.mesh.select_more()
    bpy.ops.object.mode_set(mode='OBJECT')
    polys = [i.index for i in bpy.context.active_object.data.polygons if i.select]

    polys2 =[]
    for i in polys:
        tr = [obj.data.polygons[i].vertices[0],
              obj.data.polygons[i].vertices[1],
              obj.data.polygons[i].vertices[2]]
        tr.insert(0,tr.pop(tr.index(k)))      
        polys2.append(tr)
    triang = []
    tr = []
    for i in range(len(polys2)):
        for j in range(len(polys2)):
            if len(set(polys2[i])&set(polys2[j])) == 2 and i != j:
               if  i not in triang: 
                   triang.append(i)
                   triang.append(j)
                   tr.append([polys2[i], polys2[j]])
    return tr

def cot(pair):
    obj = bpy.context.active_object
    p = list(set(pair[0]) & set(pair[1]))
    ab = list(set(pair[0]) - set(pair[1])) + list(set(pair[1]) - set(pair[0]))
    vec_a1 = obj.data.vertices[p[0]].co - obj.data.vertices[ab[0]].co
    vec_a2 = obj.data.vertices[ab[0]].co - obj.data.vertices[p[1]].co
    vec_b1 = obj.data.vertices[p[0]].co - obj.data.vertices[ab[1]].co
    vec_b2 = obj.data.vertices[ab[1]].co - obj.data.vertices[p[1]].co
    cos_a = (vec_a1.x * vec_a2.x + vec_a1.y * vec_a2.y + vec_a1.z * vec_a2.z)/(math.sqrt(vec_a1.x**2 + vec_a1.y**2 + vec_a1.z**2)* math.sqrt(vec_a2.x**2 + vec_a2.y**2 + vec_a2.z**2))
    cos_b = (vec_b1.x * vec_b2.x + vec_b1.y * vec_b2.y + vec_b1.z * vec_b2.z)/(math.sqrt(vec_b1.x**2 + vec_b1.y**2 + vec_b1.z**2)* math.sqrt(vec_b2.x**2 + vec_b2.y**2 + vec_b2.z**2))
    alpha = cos_a/(math.sqrt(1-cos_a**2))
    beta = cos_b/(math.sqrt(1-cos_b**2))
    return alpha + beta

def sq_norm(pair):
    obj = bpy.context.active_object
    p = list(set(pair[0]) & set(pair[1]))
    return (obj.matrix_world * (obj.data.vertices[p[0]].co - obj.data.vertices[p[1]].co)).length**2

def com_edge(pair):
    obj = bpy.context.active_object
    p = list(set(pair[0]) & set(pair[1]))
    return (obj.matrix_world * (obj.data.vertices[p[0]].co - obj.data.vertices[p[1]].co))

def v_area(ring):
    v_a = 0
    for i in range(len(ring)):
        v_a = v_a + cot(ring[i]) * sq_norm(ring[i])
    v_a = 0.125*v_a
    return v_a

def mean_curvature(k): 
    ring = create_pairs(k)
    v_area(ring)
    mean = mathutils.Vector((0,0,0))
    for i in range(len(ring)):
        mean = mean + cot(ring[i]) * com_edge(ring[i])
    mean = 0.5 * (0.5 * v_area(ring) * mean).length 
    return mean

obj = bpy.context.active_object
for k in range(len(obj.data.vertices)):
    print(mean_curvature(k))

Computes mean curvature at each vertex of a mesh, needs tweaking - not working with obtuse triangles, other than that seems to get the job done.

Edit: fixed bugs, beatified code