Scripts: curvature.py

# Color selected surface pieces by mean curvature.
#
# Gray at the average curvature value over the surface, and blue and red
# at +/- 3 standard deviations of curvature values across the surface.
#
# The curvature is estimated from the vertices and normals of the triangulated
# surface in simple way which will show artifacts from non-isotropic meshes.
# For each triangle edge it computes the normal vector rotation from one vertex
# to the other divided by the edge length.  The vertex mean curvature is the
# mean of the curvatures computed for each edge.
#

def color_by_curvature(p):
    k = vertex_mean_curvatures(p)
    r = 1.0/k
    print ('radius min %.4g, max %.4g, mean %.4g, std %.4g'
           % (r.min(), r.max(), r.mean(), r.std()))

    km = k.mean()
    cval = (km - 3*k.std(), km, km + 3*k.std())
    # Colors are red, green, blue, opacity in range 0 to 1.
    low_color, mean_color, high_color = (0,0,1,1),(.7,.7,.7,1),(1,0,0,1) # blue, gray, red
    colors = (low_color, mean_color, high_color)
    above_color = below_color = (1,1,0,1)
    import SurfaceColor
    rgba = SurfaceColor.interpolate_colormap(k, cval, colors, low_color, high_color)
    p.vertexColors = rgba

def vertex_mean_curvatures(p):
    va, ta = p.geometry         # Vertices and triangles
    na = p.normals
    nv, nt = len(va), len(ta)
    from numpy import zeros, float64
    k = zeros((nv,), float64)   # Vertex curvatures
    c = zeros((nv,), float64)   # Count of edges per vertex
    from math import acos, sqrt
    for t in range(nt):
        j0, j1, j2 = ta[t]
        for i0,i1 in ((j0,j1),(j1,j2),(j2,j0)):
            v0,v1 = va[i0],va[i1]
            n0,n1 = na[i0],na[i1]
            e = v1-v0           # vector between two vertices
            elen = sqrt((e*e).sum())
            e /= elen           # Normalize edge length.
            da = acos((n0*e).sum()) - acos((n1*e).sum())
            k01 = da / elen
            k[i0] += k01
            k[i1] += k01
            c[i0] += 1
            c[i1] += 1
    k /= c
    return k

import Surface
plist = Surface.selected_surface_pieces()
for p in plist:
    color_by_curvature(p)