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?

If the code is not relevant, he could simply strip it from the question. — Leander, Aug 1 '19 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 '19 at 22:09
@gandalf3 Thank you! This sounds like a reasonable method: computergraphics.stackexchange.com/questions/1718/… — Robert Roth, Aug 2 '19 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 '19 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 '19 at 10:32

lemon · Accepted Answer · 2019-08-04 18:00:54Z

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.

@RobertRoth, you're welcome. Do you have meshes with known curvatures for testing? — lemon, Aug 4 '19 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 '19 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 '19 at 15:43
Excellent answer and extremely useful code, thank you @lemon — nantille, Oct 12 '19 at 15:26

Community · Accepted Answer · 2020-06-10 12:57:44Z

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 '19 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 '19 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 '19 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 '19 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 '19 at 17:13

Stack Exchange Network

current community

your communities

more stack exchange communities

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

2 Answers 2

Your Answer

Not the answer you're looking for? Browse other questions tagged python mesh scripting geometry mathematics or ask your own question.

Linked

Hot Network Questions

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

2 Answers 2

Your Answer

Sign up or log in

Post as a guest

Not the answer you're looking for? Browse other questions tagged python mesh scripting geometry mathematics or ask your own question.

Linked

Related

Hot Network Questions