surfaces/surf/subd/butterfly.py

from __future__ import division

from itertools import combinations

from surf.geometry import PolygonMesh, Vertex


def _get_third_point(mesh, eid, fid):
    evids = set(mesh.edges[eid])
    fvids = set(mesh.faces[fid])
    other_vid = fvids - evids
    assert len(other_vid) == 1
    other_vid = other_vid.pop()
    return other_vid


def _make_edge_point(mesh, eid):
    final_vert = Vertex()

    a = 1 / 2
    b = 1 / 8
    c = -1 / 16

    a_verts = [mesh.vertices[v] for v in mesh.edges[eid]]
    for v in a_verts:
        final_vert += a * v

    for fid in mesh.faces_for_edge[eid]:
        other_vid = _get_third_point(mesh, eid, fid)
        final_vert += b * mesh.vertices[other_vid]

        other_edges = [e for e in mesh.edges_for_face[fid] if e != eid]
        for other_edge in other_edges:
            wing_face = [f for f in mesh.faces_for_edge[other_edge]
                    if f != fid]
            assert len(wing_face) == 1
            vid = _get_third_point(mesh, other_edge, wing_face[0])
            final_vert += c * mesh.vertices[vid]

    return final_vert


def refine(mesh):
    new_verts = list(mesh.vertices)
    nv_offset = len(new_verts)
    new_edges = []
    new_faces = []

    for eid, verts_for_edge in enumerate(mesh.edges):
        new_vert = _make_edge_point(mesh, eid)
        new_verts.append(new_vert)

        # populate *some* of the new edges (missing all edges on a face ...)
        new_vert_id = eid + nv_offset
        for evid in mesh.edges[eid]:
            new_edge = [evid, new_vert_id]
            new_edges.append(new_edge)

    for fid in range(len(mesh.faces)):
        # join all new edge points on a face:
        new_center_face = [i + nv_offset for i in mesh.edges_for_face[fid]]
        inner_edges = list(combinations(new_center_face, 2))
        new_edges.extend(inner_edges)

        # make new faces: center face
        new_faces.append(new_center_face)

        # three new faces, one per old vertex:
        for fvid in mesh.faces[fid]:
            possible_edges = mesh.edges_for_vert[fvid]
            face_edges = mesh.edges_for_face[fid]
            edges_for_fvid = list(set(possible_edges) & set(face_edges))
            assert len(edges_for_fvid) == 2
            new_faces.append([
                edges_for_fvid[0] + nv_offset,
                fvid,
                edges_for_fvid[1] + nv_offset
            ])

    return PolygonMesh(new_verts, new_faces, new_edges)