diff options
Diffstat (limited to 'mesh_inset/model.py')
-rw-r--r-- | mesh_inset/model.py | 507 |
1 files changed, 507 insertions, 0 deletions
diff --git a/mesh_inset/model.py b/mesh_inset/model.py new file mode 100644 index 00000000..4e3b0f66 --- /dev/null +++ b/mesh_inset/model.py @@ -0,0 +1,507 @@ +# ##### BEGIN GPL LICENSE BLOCK ##### +# +# This program is free software; you can redistribute it and/or +# modify it under the terms of the GNU General Public License +# as published by the Free Software Foundation; either version 2 +# of the License, or (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program; if not, write to the Free Software Foundation, +# Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +# +# ##### END GPL LICENSE BLOCK ##### + +"""Manipulations of Models. +""" + +__author__ = "howard.trickey@gmail.com" + +from . import geom +from . import triquad +from . import offset +import math + + +def PolyAreasToModel(polyareas, bevel_amount, bevel_pitch, quadrangulate): + """Convert a PolyAreas into a Model object. + + Assumes polyareas are in xy plane. + + Args: + polyareas: geom.PolyAreas + bevel_amount: float - if > 0, amount of bevel + bevel_pitch: float - if > 0, angle in radians of bevel + quadrangulate: bool - should n-gons be quadrangulated? + Returns: + geom.Model + """ + + m = geom.Model() + if not polyareas: + return m + polyareas.points.AddZCoord(0.0) + m.points = polyareas.points + for pa in polyareas.polyareas: + PolyAreaToModel(m, pa, bevel_amount, bevel_pitch, quadrangulate) + return m + + +def PolyAreaToModel(m, pa, bevel_amount, bevel_pitch, quadrangulate): + if bevel_amount > 0.0: + BevelPolyAreaInModel(m, pa, bevel_amount, bevel_pitch, quadrangulate) + elif quadrangulate: + if len(pa.poly) == 0: + return + qpa = triquad.QuadrangulateFaceWithHoles(pa.poly, pa.holes, pa.points) + m.faces.extend(qpa) + m.colors.extend([ pa.color ] * len(qpa)) + else: + m.faces.append(pa.poly) + # TODO: just the first part of QuadrangulateFaceWithHoles, to join + # holes to outer poly + m.colors.append(pa.color) + + +def ExtrudePolyAreasInModel(mdl, polyareas, depth, cap_back): + """Extrude the boundaries given by polyareas by -depth in z. + + Assumes polyareas are in xy plane. + + Arguments: + mdl: geom.Model - where to do extrusion + polyareas: geom.Polyareas + depth: float + cap_back: bool - if True, cap off the back + Side Effects: + For all edges in polys in polyareas, make quads in Model + extending those edges by depth in the negative z direction. + The color will be the color of the face that the edge is part of. + """ + + for pa in polyareas.polyareas: + back_poly = _ExtrudePoly(mdl, pa.poly, depth, pa.color, True) + back_holes = [] + for p in pa.holes: + back_holes.append(_ExtrudePoly(mdl, p, depth, pa.color, False)) + if cap_back: + qpa = triquad.QuadrangulateFaceWithHoles(back_poly, back_holes, + polyareas.points) + # need to reverse each poly to get normals pointing down + for i, p in enumerate(qpa): + t = list(p) + t.reverse() + qpa[i] = tuple(t) + model.faces.extend(qpa) + model.colors.extend([pa.color] * len(qpa)) + + +def _ExtrudePoly(mdl, poly, depth, color, isccw): + """Extrude the poly by -depth in z + + Arguments: + mdl: geom.Model - where to do extrusion + poly: list of vertex indices + depth: float + color: tuple of three floats + isccw: True if counter-clockwise + Side Effects + For all edges in poly, make quads in Model + extending those edges by depth in the negative z direction. + The color will be the color of the face that the edge is part of. + Returns: + list of int - vertices for extruded poly + """ + + if len(poly) < 2: + return + extruded_poly = [] + points = mdl.points + if isccw: + incr = 1 + else: + incr = -1 + for i, v in enumerate(poly): + vnext = poly[(i+incr) % len(poly)] + (x0,y0,z0) = points.pos[v] + (x1,y1,z1) = points.pos[vnext] + vextrude = points.AddPoint((x0,y0,z0-depth)) + vnextextrude = points.AddPoint((x1,y1,z1-depth)) + if isccw: + sideface = [v, vextrude, vnextextrude, vnext] + else: + sideface = [v, vnext, vnextextrude, vextrude] + mdl.faces.append(sideface) + mdl.colors.append(color) + extruded_poly.append(vextrude) + return extruded_poly + + +def BevelPolyAreaInModel(mdl, polyarea, + bevel_amount, bevel_pitch, quadrangulate): + """Bevel the interior of polyarea in model. + + This does smart beveling: advancing edges are merged + rather than doing an 'overlap'. Advancing edges that + hit an opposite edge result in a split into two beveled areas. + + If the polyarea is not in the xy plane, do the work in a + transformed model, and then transfer the changes back. + + Arguments: + mdl: geom.Model - where to do bevel + polyarea geom.PolyArea - area to bevel into + bevel_amount: float - if > 0, amount of bevel + bevel_pitch: float - if > 0, angle in radians of bevel + quadrangulate: bool - should n-gons be quadrangulated? + Side Effects: + Faces and points are added to model to model the + bevel and the interior of the polyareas. + """ + + pa_norm = polyarea.Normal() + if pa_norm == (0.0, 0.0, 1.0): + m = mdl + pa_rot = polyarea + else: + (pa_rot, inv_rot, inv_map) = _RotatedPolyAreaToXY(polyarea, pa_norm) + # don't actually have to add the original faces into model, just their points. + m = geom.Model() + m.points = pa_rot.points + vspeed = math.tan(bevel_pitch) + off = offset.Offset(pa_rot, 0.0, vspeed) + off.Build(bevel_amount) + inner_pas = AddOffsetFacesToModel(m, off, polyarea.color) + for pa in inner_pas.polyareas: + if quadrangulate: + if len(pa.poly) == 0: + continue + qpa = triquad.QuadrangulateFaceWithHoles(pa.poly, pa.holes, pa.points) + m.faces.extend(qpa) + m.colors.extend([ pa.color ] * len(qpa)) + else: + m.faces.append(pa.poly) + m.colors.append(pa.color) + if m != mdl: + _AddTransformedPolysToModel(mdl, m.faces, m.points, inv_rot, inv_map) + + +def AddOffsetFacesToModel(mdl, off, color = (0.0, 0.0, 0.0)): + """Add the faces due to an offset into model. + + Returns the remaining interiors of the offset as a PolyAreas. + + Args: + mdl: geom.Model - model to add offset faces into + off: offset.Offset + color: (float, float, float) - color to make the faces + Returns: + geom.PolyAreas + """ + + mdl.points = off.polyarea.points + assert(len(mdl.points.pos) == 0 or len(mdl.points.pos[0]) == 3) + o = off + ostack = [ ] + while o: + if o.endtime != 0.0: + for face in o.facespokes: + n = len(face) + for i, spoke in enumerate(face): + nextspoke = face[(i+1) % n] + v0 = spoke.origin + v1 = nextspoke.origin + v2 = nextspoke.dest + v3 = spoke.dest + if v2 == v3: + mface = [v0, v1, v2] + else: + mface = [v0, v1, v2, v3] + mdl.faces.append(mface) + mdl.colors.append(color) + ostack.extend(o.inneroffsets) + if ostack: + o = ostack.pop() + else: + o = None + return off.InnerPolyAreas() + + +def BevelSelectionInModel(mdl, selected_faces, + bevel_amount, bevel_pitch, quadrangulate, as_region): + """Bevel the selected faces in the model. + + If as_region is False, each face is beveled individually, + otherwise regions of contiguous faces are merged into + PolyAreas and beveled as a whole. + + TODO: something if extracted PolyAreas are not approximately + planar. + + Args: + mdl: geom.Model + selected_faces: list of list of int + bevel_amount: float - amount to inset + bevel_pitch: float - angle of bevel side + quadrangulate: bool - should insides be quadrangulated? + as_region: bool - should faces be merged into regions? + Side effect: + Beveling faces will be added to the model + """ + + pas = [] + if as_region: + pas = RegionToPolyAreas(selected_faces, mdl.points) + else: + for face in selected_faces: + pas.append(geom.PolyArea(mdl.points, face)) + for pa in pas: + BevelPolyAreaInModel(mdl, pa, + bevel_amount, bevel_pitch, quadrangulate) + + +def RegionToPolyAreas(faces, points): + """Find polygonal outlines induced by union of faces. + + Finds the polygons formed by boundary edges (those not + sharing an edge with another face in region_faces), and + turns those into PolyAreas. + In the general case, there will be holes inside. + + Args: + faces: list of list of int - each sublist is a face (indices into points) + points: geom.Points - gives coordinates for vertices + Returns: + list of geom.PolyArea + """ + + ans = [] + (edges, vtoe) = _GetEdgeData(faces) + (face_adj, is_interior_edge) = _GetFaceGraph(faces, edges, vtoe, points) + (components, ftoc) = _FindFaceGraphComponents(faces, face_adj) + for c in range(len(components)): + boundary_edges = set() + vstobe = dict() + for e, ((vs, ve), f) in enumerate(edges): + if ftoc[f] != c or is_interior_edge[e]: + continue + boundary_edges.add(e) + vstobe[vs] = e + polys = [] + while boundary_edges: + e = boundary_edges.pop() + ((vstart, ve), _) = edges[e] + poly = [ vstart, ve ] + while ve != vstart: + if ve not in vstobe: + print("whoops, couldn't close boundary") + break + nexte = vstobe[ve] + ((_, ve), _) = edges[nexte] + boundary_edges.remove(nexte) + if ve != vstart: + poly.append(ve) + polys.append(poly) + if len(polys) == 0: + # can happen if an entire closed polytope is given + # at least until we do an edge check + return [] + elif len(polys) == 1: + ans.append(geom.PolyArea(points, polys[0])) + else: + outerf = _FindOuterPoly(polys, points) + pa = geom.PolyArea(points, polys[outerf]) + pa.holes = [ polys[i] for i in range(len(polys)) if i != outerf ] + ans.append(pa) + return ans + + +def _GetEdgeData(faces): + """Find edges from faces, and some lookup dictionaries. + + Args: + faces: list of list of int - each a closed CCW polygon of vertex indices + Returns: + (list of ((int, int), int), dict{ int->list of int}) - + list elements are ((startv, endv), face index) + dict maps vertices to edge indices + """ + + edges = [] + vtoe = dict() + for findex, f in enumerate(faces): + nf = len(f) + for i, v in enumerate(f): + endv = f[(i+1) % nf] + edges.append(((v, endv), findex)) + eindex = len(edges)-1 + if v in vtoe: + vtoe[v].append(eindex) + else: + vtoe[v] = [ eindex ] + return (edges, vtoe) + + +def _GetFaceGraph(faces, edges, vtoe, points): + """Find the face adjacency graph. + + Faces are adjacent if they share an edge, + and the shared edge goes in the reverse direction, + and if the angle between them isn't too large. + + Args: + faces: list of list of int + edges: list of ((int, int), int) - see _GetEdgeData + vtoe: dict{ int->list of int } - see _GetEdgeData + points: geom.Points + Returns: + (list of list of int, list of bool) - + first list: each sublist is adjacent face indices for each face + second list: maps edge index to True if it separates adjacent faces + """ + + face_adj = [ [] for i in range(len(faces)) ] + is_interior_edge = [ False ] * len(edges) + for e, ((vs, ve), f) in enumerate(edges): + for othere in vtoe[ve]: + ((_, we), g) = edges[othere] + if we == vs: + # face g is adjacent to face f + # TODO: angle check + if g not in face_adj[f]: + face_adj[f].append(g) + is_interior_edge[e] = True + # Don't bother with mirror relations, will catch later + return (face_adj, is_interior_edge) + + +def _FindFaceGraphComponents(faces, face_adj): + """Partition faces into connected components. + + Args: + faces: list of list of int + face_adj: list of list of int - see _GetFaceGraph + Returns: + (list of list of int, list of int) - + first list partitions face indices into separate lists, each a component + second list maps face indices into their component index + """ + + if not faces: + return ([], []) + components = [] + ftoc = [ -1 ] * len(faces) + for i in range(len(faces)): + if ftoc[i] == -1: + compi = len(components) + comp = [] + _FFGCSearch(i, faces, face_adj, ftoc, compi, comp) + components.append(comp) + return (components, ftoc) + + +def _FFGCSearch(findex, faces, face_adj, ftoc, compi, comp): + """Depth first search helper function for _FindFaceGraphComponents + + Searches recursively through all faces connected to findex, adding + each face found to comp and setting ftoc for that face to compi. + """ + + comp.append(findex) + ftoc[findex] = compi + for otherf in face_adj[findex]: + if ftoc[otherf] == -1: + _FFGCSearch(otherf, faces, face_adj, ftoc, compi, comp) + +def _FindOuterPoly(polys, points): + """Assuming polys has one that contains the rest, find that one. + + Args: + polys: list of list of int - list of polys given by vertex indices + points: geom.Points + Returns: + int - the index in polys of the outermost one + """ + + if len(polys) < 2: + return 0 + for i, poly in enumerate(polys): + otherpoly = polys[(i+1) % len(polys)] + if geom.PointInside(points.pos[otherpoly[0]], poly, points) == 1: + return i + print("whoops, couldn't find an outermost poly") + return 0 + + +def _RotatedPolyAreaToXY(polyarea, norm): + """Return a PolyArea rotated to xy plane. + + Only the points in polyarea will be transferred. + + Args: + polyarea: geom.PolyArea + norm: the normal for polyarea + Returns: + (geom.PolyArea, (float, ..., float), dict{ int -> int }) - new PolyArea, + 4x3 inverse transform, dict mapping new verts to old ones + """ + + # find rotation matrix that takes norm to (0,0,1) + (nx, ny, nz) = norm + if abs(nx) < abs(ny) and abs(nx) < abs(nz): + v = (vx, vy, vz) = geom.Norm3(0.0, nz, - ny) + elif abs(ny) < abs(nz): + v = (vx, vy, vz) = geom.Norm3(nz, 0.0, - nx) + else: + v = (vx, vy, vz) = geom.Norm3(ny, - nx, 0.0) + (ux, uy, uz) = geom.Cross3(v, norm) + rotmat = [ux, vx, nx, uy, vy, ny, uz, vz, nz, 0.0, 0.0, 0.0] + # rotation matrices are orthogonal, so inverse is transpose + invrotmat = [ux, uy, uz, vx, vy, vz, nx, ny, nz, 0.0, 0.0, 0.0] + pointmap = dict() + invpointmap = dict() + newpoints = geom.Points() + for poly in [polyarea.poly] + polyarea.holes: + for v in poly: + vcoords = polyarea.points.pos[v] + newvcoords = geom.MulPoint3(vcoords, rotmat) + newv = newpoints.AddPoint(newvcoords) + pointmap[v] = newv + invpointmap[newv] = v + pa = geom.PolyArea(newpoints) + pa.poly = [ pointmap[v] for v in polyarea.poly ] + pa.holes = [ [ pointmap[v] for v in hole ] for hole in polyarea.holes ] + return (pa, invrotmat, invpointmap) + + +def _AddTransformedPolysToModel(mdl, polys, points, transform, pointmap): + """Add (transformed) the points and faces to a model. + + Add polys to mdl. The polys have coordinates given by indices into points.pos; + those need to be transformed by multiplying by the transform matrix. + The vertices may already exist in mdl. Rather than relying on AddPoint to detect + the duplicate (transform rounding error makes that dicey), the pointmap dictionary + is used to map vertex indices in polys into those in mdl - if they exist already. + + Args: + mdl: geom.Model - where to put new vertices, faces + polys: list of list of int - each sublist a poly + points: geom.Points - coords for vertices in polys + transform: (float, ..., float) - 12-tuple, a 4x3 transform matrix + pointmap: dict { int -> int } - maps new vertex indices to old ones + Side Effects: + The model gets new faces and vertices, based on those in polys. + We are allowed to modify pointmap, as it will be discarded after call. + """ + + for i, coords in enumerate(points.pos): + if i not in pointmap: + p = geom.MulPoint3(coords, transform) + pointmap[i] = mdl.points.AddPoint(p) + for poly in polys: + mpoly = [ pointmap[v] for v in poly ] + mdl.faces.append(mpoly) |