Triangle flipping Python snippet

triangle-flipping.py

def find_flips(mesh, start_el_nr=None, flip_flags=None, remaining_els=None, thresh=1e-5):
    if start_el_nr is None:
        start_el_nr = 0

    if flip_flags is None:
        flip_flags = np.zeros(len(mesh), dtype=np.bool)

    if remaining_els is None:
        remaining_els = set(xrange(len(mesh)))
        remaining_els.remove(start_el_nr)

    t2t = mesh.triangle_to_neighbor_triangles

    queue = [ # pairs of (from, to)
            (start_el_nr, neigh_el) for neigh_el in t2t[start_el_nr]
            if neigh_el in remaining_els
            ]

    pts = mesh.vertex_coordinates

    while queue:
        from_el, to_el = queue.pop(-1)

        if to_el not in remaining_els:
            continue

        #       v_from
        #        o._
        #       /   ~._
        #      /       ~._
        #     /           ~._
        # va o---------------o vb
        #     \          _.~
        #      \      _.~
        #       \  _.~
        #        o~
        #       v_to

        from_vertices = set(mesh.elements[from_el])
        to_vertices = set(mesh.elements[to_el])
        shared_edge = from_vertices & to_vertices
        va, vb = shared_edge
        v_to, = to_vertices - shared_edge
        v_from, = from_vertices - shared_edge

        from_normal = np.cross(pts[v_from]-pts[va], pts[vb]-pts[va])
        from_normal /= la.norm(from_normal)

        to_normal = np.cross(pts[vb]-pts[va], pts[v_to]-pts[va])
        to_normal /= la.norm(to_normal)

        from_matches_normal = (
                (la.norm(from_normal - mesh.normals[from_el]) < thresh)
                ^ flip_flags[from_el])

        to_matches_normal = la.norm(to_normal - mesh.normals[to_el]) < thresh

        flip_flags[to_el] = to_matches_normal ^ from_matches_normal
        remaining_els.remove(to_el)

        queue.extend(
                (to_el, neigh_el) for neigh_el in t2t[to_el]
                if neigh_el in remaining_els
                )

def perform_flips(mesh, flip_flags):
    new_elements = np.empty_like(mesh.elements)

    for el in xrange(len(mesh.elements)):
        if flip_flags[el]:
            a,b,c = mesh.elements[el]
            new_elements[el] = a,c,b
        else:
            new_elements[el] = mesh.elements[el]

    return Mesh(mesh.vertex_coordinates, new_elements)