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

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

• If the code is not relevant, he could simply strip it from the question. – Leander Aug 1 at 13:17
• Sounds like you're looking for an algorithm, in which case you may be better off asking on math.se or computergraphics.se. As far as I know Blender doesn't come with any tool which does what you're looking for. – gandalf3 Aug 1 at 22:09
• @gandalf3 Thank you! This sounds like a reasonable method: computergraphics.stackexchange.com/questions/1718/… – Robert Roth Aug 2 at 7:00
• As with all things "bad" is probably not the most correct description, there are cases where you want to use bpy.ops. Most coders write blender scripts without bpy.ops, since it is slower, context-sensitive and more difficult to read. This answer provides some generalized explanation. Every bpy.ops triggers a scene update. – Leander Aug 2 at 7:29
• It seems that you have found a solution already. Please add it as an answer instead of into the question. It is not uncommon to answer your own question once you have found your solution. Adding it as an answer makes it easier for others to find it and enables them to upvote it. – Leander Aug 2 at 10:32

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 = []
i0, i1, i2 = following_verts_of_vert( vert, face )
vertices.append( [i1, i2] )
result = vertices[0]
prev = search_link( result[0], vertices, 1 )
if prev:
result = [prev[0]] + result
vertices.remove( prev )
next = search_link( result[-1], vertices, 0 )
if next and next[1] not in result:
result.append( next[1] )
vertices.remove( next )
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.

• Thank you @lemon! :) I will try it out tomorrow! – Robert Roth Aug 4 at 16:41
• @RobertRoth, you're welcome. Do you have meshes with known curvatures for testing? – lemon Aug 4 at 17:32
• Your code works extremely well! It is also very fast! You might want to post a comment on CG-SE with a link to your answer, as i'm certain that they will also like it! Here are some surfaces with known curvatures: imgur.com/a/UYxGIoP – Robert Roth Aug 5 at 15:14
• @RobertRoth, thanks. Could be optimized a bit (a part is unnecessary, that's indicated in the answer). Thank you for the image. Will try to check that when possible. For CGSE... maybe will ask them if this 'angle' is of some mean in their opinion... (I think all that is really dependent of the accuracy we need to have). – lemon Aug 5 at 15:43
• Excellent answer and extremely useful code, thank you @lemon – nantille Oct 12 at 15:26

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

• Two questions... what is the mean using obj.matrix_world (all is fine if the object is centered to origin, but if not...)? second point but more important: in a mesh like in the question, borders should have a big mean curvature, I think (in the way this algorithm is done), but how are they (border vertices) considered here, I don't see the point? – lemon Aug 2 at 17:42
• The code didn't work for me just yet. I added a small part which sets vertex weights depending on the calculated value but all i get is noise as seen in the picture. imgur.com/a/878bwcH Am i using the code wrongly? This is the code i used with your functions: obj = bpy.context.active_object vg = obj.vertex_groups.new("VertexGroup") for k in range(len(obj.data.vertices)): meanCurvature = mean_curvature(k) print(meanCurvature) bpy.ops.object.mode_set(mode='OBJECT') vg.add([k], meanCurvature/10, "ADD") – Robert Roth Aug 2 at 18:39
• @RobertRoth, the randoms you have are due to 'obj.matrix_world *' which should be removed. No mean here to consider the object transformation matrix (except for scales but?) as we talk about vectors here, not locations. – lemon Aug 3 at 14:21
• @RobertRoth, ok, noticed matrix changes all is the object is not at the center. I've done few experiments and the DMSB_III paper is not so simple (for me). But CGSE answers are more simple and the first one gives coherent results. Mike, all this is nothing against your answer, just try to understand things. – lemon Aug 3 at 17:00
• @RobertRoth I don't understand what causing this noise, and I get the same result: here I am looking for what causing this issue – Mike Aug 3 at 17:13