How calculate the surface curvature for each vertex of a mesh?

Question

I want to calculate the curvature value for each vertex of my mesh. All the scripts suggested here Is there a way to calculate mean curvature of a triangular mesh? don't work for me. These are the script (written by lemon) and the mesh mentioned above after running the script.

import bpy
import bmesh
from mathutils import Vector
from math import sqrt, asin
Helper function to select vertices
def select_vertices( obj, indices, end_with_mode ):
    bpy.ops.object.mode_set(mode = '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')
    for index in indices:
        obj.data.vertices[index].select = True
    bpy.ops.object.mode_set(mode = end_with_mode)
def ensure_vertex_group( obj, group_name ):
    vertex_group = obj.vertex_groups.get( group_name )
    if not vertex_group:
        vertex_group = obj.vertex_groups.new( name = group_name )
    for i in range( len(obj.data.vertices) ):
        vertex_group.add( [i], 0, 'ADD' )
    return vertex_group
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' )


def index_of( element, sequence ):
    for i, e in enumerate( sequence ):
        if e == element: return i
    return -1
def search_link( value, links, position ):
    for l in links:
        if l[position] == value: return l
    return None
def rotate( l, n ):
    return l[n:] + l[:n]
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
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
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)


def mean_curvature( obj, threshold ):
# Get bmesh access on the mesh
bm = bmesh.new()
bm.from_mesh( obj.data )

bm.verts.ensure_lookup_table()
bm.edges.ensure_lookup_table()
bm.faces.ensure_lookup_table()

above_threshold = {}
curvatures = []
for vert in bm.verts:
    curvature = mean_curvature_vert( vert )
    curvatures.append( curvature )
    if abs(curvature) > threshold:
        above_threshold[vert.index] = True

return curvatures, above_threshold


print( '-------------' )
bpy.ops.object.mode_set(mode = 'OBJECT')
obj = bpy.context.active_object
threshold = 0.01
curvatures, above_threshold = mean_curvature( obj, threshold )
assign_to_vertex_group( obj, 'Curvatures', curvatures )
select_vertices( obj, above_threshold, 'OBJECT' )

Answer 1

The script here is correct but you have few high value so this is not really visible, or they are in hidden parts.

I've added a magnify factor to compensate if wanted:

Updated file: