from collections import defaultdict from surf.geometry import PolygonMesh, centroid def refine(mesh): new_vertices = list(mesh.vertices) # this can be used to both get to fvid in new_vertices and to offest # mesh.faces to be dealing with the fvid generated from that face f_vert_offset = len(new_vertices) # this doesn't need to be calculated, but it is conceptually convenient to # have along the way. Alternatively one could get this by doing a: # >>> [e + e_vert_offset for e in mesh.edges_for_face[fvid] # (fvid for all faces) evid_for_fvid = defaultdict(list) new_edges = [] new_faces = [] # For each face, add a face Vertex for old_fid, old_face in enumerate(mesh.faces): # Set each face Vertex to be the centroid of all original Vertices for # the respective face. # here we have already created a list that is a copy of the original # verts. we append the new face Vertices to that same list, and simply # keep track of offsets into the new vert list ... see f_vert_offset # above and e_vert_offset below new_face_vertex = mesh.centroid(old_face) new_vertices.append(new_face_vertex) # edge verts will be new_vertices[e_vert_offset:]. Also, this can be used # to both get to evid in new_vertices and to offest mesh.edges to be # dealing with the evid generated from that edge e_vert_offset = len(new_vertices) # For each edge, add an edge Vertex. for old_eid, old_edge in enumerate(mesh.edges): tmp_verts = [] # Set each edge Vertex to be the average of the two neighbouring, (very # recently calculated) face Vertices ... tmp_verts.extend([new_vertices[f + f_vert_offset] for f in mesh.faces_for_edge[old_eid]]) # ... and its two original endpoints. tmp_verts.extend([mesh.vertices[vid] for vid in old_edge]) # centroid == average new_vertices.append(centroid(tmp_verts)) # make mapping from evid -> fvid for later for old_fid in mesh.faces_for_edge[old_eid]: fvid = old_fid + f_vert_offset evid = old_eid + e_vert_offset evid_for_fvid[fvid].append(evid) # For each original Vertex P, move the original Vertex to a new location # Here we mimic the F, R, and P spelling in the wiki article for old_vid in xrange(f_vert_offset): # take the average F of all n face Vertices for faces touching P ... F_verts = [] for old_fid in mesh.faces_for_vert[old_vid]: F_verts.append(new_vertices[f_vert_offset + old_fid]) F = centroid(F_verts) # and take the average R of all n new edge Vertices for edges adjacent # to P edges = [mesh.edges[eid] for eid in mesh.edges_for_vert[old_vid]] e_verts = [] for edge in edges: e_verts.extend([mesh.vertices[vid] for vid in edge]) R = centroid(e_verts) # where each edge midpoint is the average of its two endpoint vertices. # *Move* each original point to the point (or replace it in new_verts) new_vertex_point = (F + 2 * R + (len(edges) - 3) * new_vertices[old_vid] ) / len(edges) new_vertices[old_vid] = new_vertex_point # Connect each new Vertex point (with original id) to the new edge # points of all original edges incident on the original vertex. for old_eid in mesh.edges_for_vert[old_vid]: new_edges.append([old_vid, old_eid + e_vert_offset]) # For each face point, add an edge connecting the new face point to each # new edge point for fvid, evids in evid_for_fvid.iteritems(): for evid in evids: new_edges.append([fvid, evid]) # to define new faces need to know the following for each new face: for fvid in xrange(f_vert_offset, e_vert_offset): # new face_vid (1 vid) # vid connected to fvid (iterate over, use 1 vid periteration) connected_vids = mesh.faces[fvid - f_vert_offset] evids_for_current_face = evid_for_fvid[fvid] for connected_vid in connected_vids: # evid connected to both the vid and the fvid (2 vids) evids_for_vid = [e + e_vert_offset for e in mesh.edges_for_vert[connected_vid]] common_ids = list(set(evids_for_vid) & set(evids_for_current_face)) assert len(common_ids) == 2 # then connect the (4) vids together ... new_faces.append([fvid, common_ids[0], connected_vid, common_ids[1]]) return PolygonMesh(vertices=new_vertices, edges=new_edges, faces=new_faces) if __name__ == '__main__': import sys import json from surf.subd.cc import refine input_file_name = sys.argv[1] cube = json.load(open(input_file_name, 'r')) p = PolygonMesh(**cube) q = refine(p) print q