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) # TODO: new_faces = [] # TODO: new_edges = [] 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_edges.extend( [[i[0] + nv_offset, i[1] + nv_offset] for i in combinations(mesh.edges_for_face[fid], 2)]) return PolygonMesh(new_verts, new_faces, new_edges)