# 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. Commented Aug 1, 2019 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. Commented Aug 1, 2019 at 22:09
• @gandalf3 Thank you! This sounds like a reasonable method: computergraphics.stackexchange.com/questions/1718/… Commented Aug 2, 2019 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. Commented Aug 2, 2019 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. Commented Aug 2, 2019 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! Commented Aug 4, 2019 at 16:41
• @RobertRoth, you're welcome. Do you have meshes with known curvatures for testing? Commented Aug 4, 2019 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 Commented Aug 5, 2019 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). Commented Aug 5, 2019 at 15:43
• Excellent answer and extremely useful code, thank you @lemon Commented Oct 12, 2019 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? Commented Aug 2, 2019 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") Commented Aug 2, 2019 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. Commented Aug 3, 2019 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. Commented Aug 3, 2019 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
Commented Aug 3, 2019 at 17:13

OpenMesh version of Nathan Reed and lemon:

import numpy as np
import openmesh as om
from vedo import *

def array_angle(array1, array2):
"""
INPUT: N x 2/3/...
"""

res = np.sum(array1 * array2, axis=1)
res /= np.linalg.norm(array1, axis=1)
res /= np.linalg.norm(array2, axis=1)
res = np.clip(res, -1.0, 1.0)
res = np.arccos(res)
return res

def run(path_mesh):
v = mesh_om.points()
f = mesh_om.face_vertex_indices()
n = mesh_om.vertex_normals()

if True:
"""
https://computergraphics.stackexchange.com/a/1719/17639
"""
e = mesh_om.ev_indices()
ev = v[e[:, 1]] - v[e[:, 0]]
en = n[e[:, 1]] - n[e[:, 0]]
e_curv = np.sum(ev * en, axis=1) / np.linalg.norm(ev, axis=1)
v_curv = np.zeros(v.shape[0])
for idx_e in range(e.shape[0]):
v_curv[e[idx_e, 0]] += e_curv[idx_e]
v_curv[e[idx_e, 1]] += e_curv[idx_e]
else:
"""
https://blender.stackexchange.com/a/147371/82691
"""
vf_indices = mesh_om.vf_indices()
ev_indices = mesh_om.ev_indices()
fe_indices = mesh_om.fe_indices()

ev_01 = v[ev_indices[:, 1]] - v[ev_indices[:, 0]]
en_01 = n[ev_indices[:, 1]] - n[ev_indices[:, 0]]
e_curv = np.sum(ev_01 * en_01, axis=1) / np.linalg.norm(ev_01, axis=1)

e_01 = v[f[:, 1]] - v[f[:, 0]]
e_12 = v[f[:, 2]] - v[f[:, 1]]
e_20 = v[f[:, 0]] - v[f[:, 2]]

f_angles = np.zeros_like(f).astype(np.float32)
f_angles[:, 0] = array_angle(e_01, -e_20)
f_angles[:, 1] = array_angle(-e_01, e_12)
f_angles[:, 2] = array_angle(e_20, -e_12)

v_curv = np.zeros(v.shape[0])
for idx_v in range(v.shape[0]):
va = 0
for i in range(vf_indices.shape[1]):
idx_f = vf_indices[idx_v, i]
if idx_f == -1:
break
idx_v_in_f = np.where(f[idx_f] == idx_v)[0][0]
a = f_angles[idx_f, idx_v_in_f]
va += a
idx_e0 = fe_indices[idx_f, (0+idx_v_in_f) % 3]
idx_e2 = fe_indices[idx_f, (2+idx_v_in_f) % 3]
v_curv[idx_v] += a * (e_curv[idx_e0] + e_curv[idx_e2])
v_curv[idx_v] /= (2 * va)

if True:
#! vis
mean_v_curv = np.mean(v_curv)
std_v_curv = np.std(v_curv)
k_min = mean_v_curv - 3 * std_v_curv
k_max = mean_v_curv + 3 * std_v_curv
v_curv = np.clip(v_curv, k_min, k_max)

M = Mesh([v, f])
M.show(resetcam=True)

if __name__ == "__main__":
path_mesh = "/home/lab9/Mesh/bunny.obj"
run(path_mesh)



When reading with Openmesh, the following message appears, indicating that the mesh may be flawed. At this time, the values of vertices and faces read by Openmesh and Blender may be different, and preprocessing may be required.

PolyMeshT::add_face: complex vertex


PolyMeshT::add_face: complex vertex is usually caused by the "wrong" faces, so we just need to delete some of the faces.

One possible solution is to use igl to read the mesh(keep the content of the file) with texcoord and normal(excluded here), and use openmesh to form the mesh, which will automatically delete the bad faces, we can calculate the mask and select valid f/ft/fn from the original array.

import numpy as np
import igl # conda install -c conda-forge igl
import openmesh as om

def write_mesh(path_mesh, v, vc=None, vt=None, vn=None, f=None, ft=None, fn=None, str_mtllib="", str_usemtl=""):
assert v.ndim == 2 and v.shape[1] == 3, ("\033[1;31mv.shape=%s\n\033[0m") % str(v.shape)
if vc is not None and len(vc) == 0:
vc = None #! vc = []
if vt is not None and len(vt) == 0:
vt = None
if vn is not None and len(vn) == 0:
vn = None
if ft is not None and len(ft) == 0:
ft = None
if fn is not None and len(fn) == 0:
fn = None

with open(path_mesh, "w") as fp:
if str_mtllib != "":
fp.write("mtllib " + str_mtllib + "\n")

if vc is not None:
v = np.concatenate([v, vc / 255.0], 1).reshape(-1, 6)
fp.write(("v {:f} {:f} {:f} {:.3f} {:.3f} {:.3f}\n" * v.shape[0]).format(*v.reshape(-1)))
else:
fp.write(("v {:f} {:f} {:f}\n" * v.shape[0]).format(*v.reshape(-1)))

if vt is not None:
fp.write(("vt {:f} {:f}\n" * vt.shape[0]).format(*vt.reshape(-1)))

if vn is not None:
fp.write(("vn {:f} {:f} {:f}\n" * vn.shape[0]).format(*vn.reshape(-1)))

if str_usemtl != "":
fp.write("usemtl " + str_usemtl + "\n")

if (ft is not None) and (fn is not None) and (f is not None):
num_f = f.shape[0]
# f v1/vt1/vn1 v2/vt2/vn2 v3/vt3/vn3 ...
for idx_f in range(num_f):
fp.write(("f {:d}/{:d}/{:d} {:d}/{:d}/{:d} {:d}/{:d}/{:d}\n").format(
f[idx_f, 0] + 1, ft[idx_f, 0] + 1, fn[idx_f, 0] + 1,
f[idx_f, 1] + 1, ft[idx_f, 1] + 1, fn[idx_f, 1] + 1,
f[idx_f, 2] + 1, ft[idx_f, 2] + 1, fn[idx_f, 2] + 1,
))
elif (ft is None) and (fn is not None) and (f is not None):
num_f = f.shape[0]
# f v1//vn1 v2//vn2 v3//vn3 ...
for idx_f in range(num_f):
fp.write(("f {:d}//{:d} {:d}//{:d} {:d}//{:d}\n").format(
f[idx_f, 0] + 1, fn[idx_f, 0] + 1,
f[idx_f, 1] + 1, fn[idx_f, 1] + 1,
f[idx_f, 2] + 1, fn[idx_f, 2] + 1,
))
elif (ft is not None) and (fn is None) and (f is not None):
num_f = f.shape[0]
# f v1/vt1 v2/vt2 v3/vt3 ...
for idx_f in range(num_f):
fp.write(("f {:d}/{:d} {:d}/{:d} {:d}/{:d}\n").format(
f[idx_f, 0] + 1, ft[idx_f, 0] + 1,
f[idx_f, 1] + 1, ft[idx_f, 1] + 1,
f[idx_f, 2] + 1, ft[idx_f, 2] + 1,
))
elif (ft is None) and (fn is None) and (f is not None):
# f v1 v2 v3 ....
fp.write(("f {:d} {:d} {:d}\n" * f.shape[0]).format(*(f.reshape(-1) + 1)))

if __name__ == "__main__":
path_mesh = "XXX.obj"
path_mesh_iglom = path_mesh[:-4] + "_IGLOM.obj"

v_igl, vt_igl, _, f_igl, ft_igl, _ = igl.read_obj(path_mesh)

#! check
if True:
dbl_area = igl.doublearea(v_igl, f_igl)
f_mask0 = (dbl_area > 2e-9).astype(np.bool_)  # ! 2e-9
if len(vt_igl) != 0:

mesh_iglom = om.PolyMesh(v_igl, f_igl)
v_iglom = mesh_iglom.points()
f_iglom = mesh_iglom.face_vertex_indices()

offset_f_igl = 0
for idx_f_iglom in range(f_iglom.shape[0]):
is_true = True
while is_true:
is_equ = True
is_equ &= f_igl[offset_f_igl, 0] == f_iglom[idx_f_iglom, 0]
is_equ &= f_igl[offset_f_igl, 1] == f_iglom[idx_f_iglom, 1]
is_equ &= f_igl[offset_f_igl, 2] == f_iglom[idx_f_iglom, 2]
if is_equ: