Mean Discrete Curvature + Connectivity

[ideas-man]

Hi,

I am trying to modify the mean discrete curvature measure method. I am trying to refine the candidates this method uses to have only those candidates that are generated using the current point (x) or its immediate neighbors, in another words, to use connectivity information of the mesh.

Right now I came up with this code, but its super slow (around 200 times slower than the original given rad of 0.05). I was wondering if someone can help me out since I am pretty sure I am missing something here.

def dmc_measure(mesh, points, radius):

    points = np.asanyarray(points, dtype=np.float64)
    if not util.is_shape(points, (-1, 3)):
        raise ValueError('points must be (n,3)!')

    # axis aligned bounds
    bounds = np.column_stack((points - radius,
                              points + radius))

    # line segments that intersect axis aligned bounding box
    candidates = [list(mesh.face_adjacency_tree.intersection(b))
                  for b in bounds]

    neighbors = [mesh.vertices[mesh.vertex_neighbors[i]]
                 for i in range(len(points))]

    mean_curv = np.empty(len(points))
    for i, (x, x_candidates, x_neighbors) in enumerate(zip(points, candidates, neighbors)):

        ref_x_candidates = [c for c in x_candidates if
                            np.any([mesh.vertices[mesh.face_adjacency_edges[c]][0] in x_neighbors,
                                    mesh.vertices[mesh.face_adjacency_edges[c]][1] in x_neighbors])]

        endpoints = mesh.vertices[mesh.face_adjacency_edges[ref_x_candidates]]
        lengths = line_ball_intersection(
            endpoints[:, 0],
            endpoints[:, 1],
            center=x,
            radius=radius)
        angles = mesh.face_adjacency_angles[ref_x_candidates]
        signs = np.where(mesh.face_adjacency_convex[ref_x_candidates], 1, -1)
        mean_curv[i] = (lengths * angles * signs).sum() / 2

    return mean_curv

It looks like ref_x_candidates list comprehension is slow. Maybe I could simply replace candidates with something that could directly give me line segments generated using neighbors of the current point, then the whole bounds could be dropped completely, but I do not know how to do that.

Cheers

[davidackerman]

if it is useful, here is my attempt at a discrete mean curvature in a one-ring neighborhood. it is based on this post: https://computergraphics.stackexchange.com/a/1721. note there is an area approximation using the (total area)/3:

import trimesh
import numpy as np
import networkx as nx

def my_discrete_mean_curvature_measure(mesh):
    """Calculate discrete mean curvature of mesh using one-ring neighborhood."""

    # one-rings (immediate neighbors of) each vertex
    g = nx.from_edgelist(mesh.edges_unique)
    one_rings = [list(g[i].keys()) for i in range(len(mesh.vertices))]

    # cotangents of angles and store in dictionary based on corresponding vertex and face
    face_angles = mesh.face_angles_sparse
    cotangents = { f"{vertex},{face}":1/np.tan(angle) for vertex,face,angle in zip(face_angles.row, face_angles.col, face_angles.data)}

    # discrete Laplace-Beltrami contribution of the shared edge of adjacent faces:
    #        /*\
    #       / * \
    #      /  *  \
    #    vi___*___vj
    #
    # store results in dictionary with vertex ids as keys     
    fa = mesh.face_adjacency
    fae = mesh.face_adjacency_edges
    edge_measure = {f"{fae[i][0]},{fae[i][1]}":(mesh.vertices[fae[i][1]] - mesh.vertices[fae[i][0]])*(cotangents[f"{v[0]},{fa[i][0]}"]+cotangents[f"{v[1]},{fa[i][1]}"]) for i,v in enumerate(mesh.face_adjacency_unshared) }

    # calculate mean curvature using one-ring
    mean_curv = [0]*len(mesh.vertices)
    for vertex_id,face_ids in enumerate(mesh.vertex_faces):
        face_ids = face_ids[face_ids!=-1] #faces associated with vertex_id
        one_ring = one_rings[vertex_id]
        delta_s = 0;

        for one_ring_vertex_id in one_ring:
            if f"{vertex_id},{one_ring_vertex_id}" in edge_measure:
                delta_s += edge_measure[f"{vertex_id},{one_ring_vertex_id}"]
            elif f"{one_ring_vertex_id},{vertex_id}"  in edge_measure:
                delta_s -= edge_measure[f"{one_ring_vertex_id},{vertex_id}"]

        delta_s *= 1/(2*sum(mesh.area_faces[face_ids])/3) #use 1/3 of the areas
        mean_curv[vertex_id] = 0.5*np.linalg.norm(delta_s)

    return np.array(mean_curv)

def main():
    # test on sphere of radius 5
    radius = 5
    m = trimesh.creation.icosphere(radius=radius)
    print(my_discrete_mean_curvature_measure(m)/(0.5*(1/radius + 1/radius)))

if __name__ == "__main__":
    main()`