mesh_inset/model.py

   1 # ##### BEGIN GPL LICENSE BLOCK #####
   2 #
   3 #  This program is free software; you can redistribute it and/or
   4 #  modify it under the terms of the GNU General Public License
   5 #  as published by the Free Software Foundation; either version 2
   6 #  of the License, or (at your option) any later version.
   7 #
   8 #  This program is distributed in the hope that it will be useful,
   9 #  but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 #  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 #  GNU General Public License for more details.
  12 #
  13 #  You should have received a copy of the GNU General Public License
  14 #  along with this program; if not, write to the Free Software Foundation,
  15 #  Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  16 #
  17 # ##### END GPL LICENSE BLOCK #####
  18
  19 # <pep8 compliant>
  20
  21 """Manipulations of Models.
  22 """
  23
  24 __author__ = "howard.trickey@gmail.com"
  25
  26 from . import geom
  27 from . import triquad
  28 from . import offset
  29 import math
  30
  31
  32 def PolyAreasToModel(polyareas, bevel_amount, bevel_pitch, quadrangulate):
  33     """Convert a PolyAreas into a Model object.
  34
  35     Assumes polyareas are in xy plane.
  36
  37     Args:
  38       polyareas: geom.PolyAreas
  39       bevel_amount: float - if > 0, amount of bevel
  40       bevel_pitch: float - if > 0, angle in radians of bevel
  41       quadrangulate: bool - should n-gons be quadrangulated?
  42     Returns:
  43       geom.Model
  44     """
  45
  46     m = geom.Model()
  47     if not polyareas:
  48         return m
  49     polyareas.points.AddZCoord(0.0)
  50     m.points = polyareas.points
  51     for pa in polyareas.polyareas:
  52         PolyAreaToModel(m, pa, bevel_amount, bevel_pitch, quadrangulate)
  53     return m
  54
  55
  56 def PolyAreaToModel(m, pa, bevel_amount, bevel_pitch, quadrangulate):
  57     if bevel_amount > 0.0:
  58         BevelPolyAreaInModel(m, pa, bevel_amount, bevel_pitch, quadrangulate,
  59             False)
  60     elif quadrangulate:
  61         if len(pa.poly) == 0:
  62             return
  63         qpa = triquad.QuadrangulateFaceWithHoles(pa.poly, pa.holes, pa.points)
  64         m.faces.extend(qpa)
  65         m.face_data.extend([pa.data] * len(qpa))
  66     else:
  67         m.faces.append(pa.poly)
  68         # TODO: just the first part of QuadrangulateFaceWithHoles, to join
  69         # holes to outer poly
  70         m.face_data.append(pa.data)
  71
  72
  73 def ExtrudePolyAreasInModel(mdl, polyareas, depth, cap_back):
  74     """Extrude the boundaries given by polyareas by -depth in z.
  75
  76     Assumes polyareas are in xy plane.
  77
  78     Arguments:
  79       mdl: geom.Model - where to do extrusion
  80       polyareas: geom.Polyareas
  81       depth: float
  82       cap_back: bool - if True, cap off the back
  83     Side Effects:
  84       For all edges in polys in polyareas, make quads in Model
  85       extending those edges by depth in the negative z direction.
  86       The application data will be the data of the face that the edge
  87       is part of.
  88     """
  89
  90     for pa in polyareas.polyareas:
  91         back_poly = _ExtrudePoly(mdl, pa.poly, depth, pa.data, True)
  92         back_holes = []
  93         for p in pa.holes:
  94             back_holes.append(_ExtrudePoly(mdl, p, depth, pa.data, False))
  95         if cap_back:
  96             qpa = triquad.QuadrangulateFaceWithHoles(back_poly, back_holes,
  97               polyareas.points)
  98             # need to reverse each poly to get normals pointing down
  99             for i, p in enumerate(qpa):
 100                 t = list(p)
 101                 t.reverse()
 102                 qpa[i] = tuple(t)
 103             mdl.faces.extend(qpa)
 104             mdl.face_data.extend([pa.data] * len(qpa))
 105
 106
 107 def _ExtrudePoly(mdl, poly, depth, data, isccw):
 108     """Extrude the poly by -depth in z
 109
 110     Arguments:
 111       mdl: geom.Model - where to do extrusion
 112       poly: list of vertex indices
 113       depth: float
 114       data: application data
 115       isccw: True if counter-clockwise
 116     Side Effects
 117       For all edges in poly, make quads in Model
 118       extending those edges by depth in the negative z direction.
 119       The application data will be the data of the face that the edge
 120       is part of.
 121     Returns:
 122       list of int - vertices for extruded poly
 123     """
 124
 125     if len(poly) < 2:
 126         return
 127     extruded_poly = []
 128     points = mdl.points
 129     if isccw:
 130         incr = 1
 131     else:
 132         incr = -1
 133     for i, v in enumerate(poly):
 134         vnext = poly[(i + incr) % len(poly)]
 135         (x0, y0, z0) = points.pos[v]
 136         (x1, y1, z1) = points.pos[vnext]
 137         vextrude = points.AddPoint((x0, y0, z0 - depth))
 138         vnextextrude = points.AddPoint((x1, y1, z1 - depth))
 139         if isccw:
 140             sideface = [v, vextrude, vnextextrude, vnext]
 141         else:
 142             sideface = [v, vnext, vnextextrude, vextrude]
 143         mdl.faces.append(sideface)
 144         mdl.face_data.append(data)
 145         extruded_poly.append(vextrude)
 146     return extruded_poly
 147
 148
 149 def BevelPolyAreaInModel(mdl, polyarea,
 150     bevel_amount, bevel_pitch, quadrangulate, as_percent):
 151     """Bevel the interior of polyarea in model.
 152
 153     This does smart beveling: advancing edges are merged
 154     rather than doing an 'overlap'.  Advancing edges that
 155     hit an opposite edge result in a split into two beveled areas.
 156
 157     If the polyarea is not in the xy plane, do the work in a
 158     transformed model, and then transfer the changes back.
 159
 160     Arguments:
 161       mdl: geom.Model - where to do bevel
 162       polyarea geom.PolyArea - area to bevel into
 163       bevel_amount: float - if > 0, amount of bevel
 164       bevel_pitch: float - if > 0, angle in radians of bevel
 165       quadrangulate: bool - should n-gons be quadrangulated?
 166       as_percent: bool - if True, interpret amount as percent of max
 167     Side Effects:
 168       Faces and points are added to model to model the
 169       bevel and the interior of the polyareas.
 170     """
 171
 172     pa_norm = polyarea.Normal()
 173     if pa_norm == (0.0, 0.0, 1.0):
 174         m = mdl
 175         pa_rot = polyarea
 176     else:
 177         (pa_rot, inv_rot, inv_map) = _RotatedPolyAreaToXY(polyarea, pa_norm)
 178         # don't have to add the original faces into model, just their points.
 179         m = geom.Model()
 180         m.points = pa_rot.points
 181     vspeed = math.tan(bevel_pitch)
 182     off = offset.Offset(pa_rot, 0.0, vspeed)
 183     if as_percent:
 184         bevel_amount = bevel_amount * off.MaxAmount() / 100.0
 185     off.Build(bevel_amount)
 186     inner_pas = AddOffsetFacesToModel(m, off, polyarea.data)
 187     for pa in inner_pas.polyareas:
 188         if quadrangulate:
 189             if len(pa.poly) == 0:
 190                 continue
 191             qpa = triquad.QuadrangulateFaceWithHoles(pa.poly, pa.holes,
 192                 pa.points)
 193             m.faces.extend(qpa)
 194             m.face_data.extend([pa.data] * len(qpa))
 195         else:
 196             m.faces.append(pa.poly)
 197             m.face_data.append(pa.data)
 198     if m != mdl:
 199         _AddTransformedPolysToModel(mdl, m.faces, m.points, m.face_data,
 200              inv_rot, inv_map)
 201
 202
 203 def AddOffsetFacesToModel(mdl, off, data=None):
 204     """Add the faces due to an offset into model.
 205
 206     Returns the remaining interiors of the offset as a PolyAreas.
 207
 208     Args:
 209       mdl: geom.Model - model to add offset faces into
 210       off: offset.Offset
 211       data: any - application data to be copied to the faces
 212     Returns:
 213       geom.PolyAreas
 214     """
 215
 216     mdl.points = off.polyarea.points
 217     assert(len(mdl.points.pos) == 0 or len(mdl.points.pos[0]) == 3)
 218     o = off
 219     ostack = []
 220     while o:
 221         if o.endtime != 0.0:
 222             for face in o.facespokes:
 223                 n = len(face)
 224                 for i, spoke in enumerate(face):
 225                     nextspoke = face[(i + 1) % n]
 226                     v0 = spoke.origin
 227                     v1 = nextspoke.origin
 228                     v2 = nextspoke.dest
 229                     v3 = spoke.dest
 230                     if v2 == v3:
 231                         mface = [v0, v1, v2]
 232                     else:
 233                         mface = [v0, v1, v2, v3]
 234                     mdl.faces.append(mface)
 235                     mdl.face_data.append(data)
 236         ostack.extend(o.inneroffsets)
 237         if ostack:
 238             o = ostack.pop()
 239         else:
 240             o = None
 241     return off.InnerPolyAreas()
 242
 243
 244 def BevelSelectionInModel(mdl, bevel_amount, bevel_pitch, quadrangulate,
 245         as_region, as_percent):
 246     """Bevel all the faces in the model, perhaps as one region.
 247
 248     If as_region is False, each face is beveled individually,
 249     otherwise regions of contiguous faces are merged into
 250     PolyAreas and beveled as a whole.
 251
 252     TODO: something if extracted PolyAreas are not approximately
 253     planar.
 254
 255     Args:
 256       mdl: geom.Model
 257       bevel_amount: float - amount to inset
 258       bevel_pitch: float - angle of bevel side
 259       quadrangulate: bool - should insides be quadrangulated?
 260       as_region: bool - should faces be merged into regions?
 261       as_percent: bool - should amount be interpreted as a percent
 262           of the maximum amount (if True) or an absolute amount?
 263     Side effect:
 264       Beveling faces will be added to the model
 265     """
 266
 267     pas = []
 268     if as_region:
 269         pas = RegionToPolyAreas(mdl.faces, mdl.points, mdl.face_data)
 270     else:
 271         for f, face in enumerate(mdl.faces):
 272             pas.append(geom.PolyArea(mdl.points, face, [],
 273                 mdl.face_data[f]))
 274     for pa in pas:
 275         BevelPolyAreaInModel(mdl, pa,
 276             bevel_amount, bevel_pitch, quadrangulate, as_percent)
 277
 278
 279 def RegionToPolyAreas(faces, points, data):
 280     """Find polygonal outlines induced by union of faces.
 281
 282     Finds the polygons formed by boundary edges (those not
 283     sharing an edge with another face in region_faces), and
 284     turns those into PolyAreas.
 285     In the general case, there will be holes inside.
 286     We want to associate data with the region PolyAreas.
 287     Just choose a representative element of data[] when
 288     more than one face is combined into a PolyArea.
 289
 290     Args:
 291       faces: list of list of int - each sublist is a face (indices into points)
 292       points: geom.Points - gives coordinates for vertices
 293       data: list of any - parallel to faces, app data to put in PolyAreas
 294     Returns:
 295       list of geom.PolyArea
 296     """
 297
 298     ans = []
 299     (edges, vtoe) = _GetEdgeData(faces)
 300     (face_adj, is_interior_edge) = _GetFaceGraph(faces, edges, vtoe, points)
 301     (components, ftoc) = _FindFaceGraphComponents(faces, face_adj)
 302     for c in range(len(components)):
 303         boundary_edges = set()
 304         betodata = dict()
 305         vstobe = dict()
 306         for e, ((vs, ve), f) in enumerate(edges):
 307             if ftoc[f] != c or is_interior_edge[e]:
 308                 continue
 309             boundary_edges.add(e)
 310             # vstobe[v] is set of edges leaving v
 311             # (could be more than one if boundary touches itself at a vertex)
 312             if vs in vstobe:
 313                 vstobe[vs].append(e)
 314             else:
 315                 vstobe[vs] = [e]
 316             betodata[(vs, ve)] = data[f]
 317         polys = []
 318         poly_data = []
 319         while boundary_edges:
 320             e = boundary_edges.pop()
 321             ((vstart, ve), face_i) = edges[e]
 322             poly = [vstart, ve]
 323             datum = betodata[(vstart, ve)]
 324             while ve != vstart:
 325                 if ve not in vstobe:
 326                     print("whoops, couldn't close boundary")
 327                     break
 328                 nextes = vstobe[ve]
 329                 if len(nextes) == 1:
 330                     nexte = nextes[0]
 331                 else:
 332                     # find a next edge with face index face_i
 333                     # TODO: this is not guaranteed to work,
 334                     # as continuation edge may have been for a different
 335                     # face that is now combined with face_i by erasing
 336                     # interior edges. Find a better algorithm here.
 337                     nexte = -1
 338                     for ne_cand in nextes:
 339                         if edges[ne_cand][1] == face_i:
 340                             nexte = ne_cand
 341                             break
 342                     if nexte == -1:
 343                         # case mentioned in TODO may have happened;
 344                         # just choose any nexte - may mess things up
 345                         nexte = nextes[0]
 346                 ((_, ve), face_i) = edges[nexte]
 347                 if nexte not in boundary_edges:
 348                     print("whoops, nexte not a boundary edge", nexte)
 349                     break
 350                 boundary_edges.remove(nexte)
 351                 if ve != vstart:
 352                     poly.append(ve)
 353             polys.append(poly)
 354             poly_data.append(datum)
 355         if len(polys) == 0:
 356             # can happen if an entire closed polytope is given
 357             # at least until we do an edge check
 358             return []
 359         elif len(polys) == 1:
 360             ans.append(geom.PolyArea(points, polys[0], [], poly_data[0]))
 361         else:
 362             outerf = _FindOuterPoly(polys, points, faces)
 363             pa = geom.PolyArea(points, polys[outerf], [], poly_data[outerf])
 364             pa.holes = [polys[i] for i in range(len(polys)) if i != outerf]
 365             ans.append(pa)
 366     return ans
 367
 368
 369 def _GetEdgeData(faces):
 370     """Find edges from faces, and some lookup dictionaries.
 371
 372     Args:
 373       faces: list of list of int - each a closed CCW polygon of vertex indices
 374     Returns:
 375       (list of ((int, int), int), dict{ int->list of int}) -
 376         list elements are ((startv, endv), face index)
 377         dict maps vertices to edge indices
 378     """
 379
 380     edges = []
 381     vtoe = dict()
 382     for findex, f in enumerate(faces):
 383         nf = len(f)
 384         for i, v in enumerate(f):
 385             endv = f[(i + 1) % nf]
 386             edges.append(((v, endv), findex))
 387             eindex = len(edges) - 1
 388             if v in vtoe:
 389                 vtoe[v].append(eindex)
 390             else:
 391                 vtoe[v] = [eindex]
 392     return (edges, vtoe)
 393
 394
 395 def _GetFaceGraph(faces, edges, vtoe, points):
 396     """Find the face adjacency graph.
 397
 398     Faces are adjacent if they share an edge,
 399     and the shared edge goes in the reverse direction,
 400     and if the angle between them isn't too large.
 401
 402     Args:
 403       faces: list of list of int
 404       edges: list of ((int, int), int) - see _GetEdgeData
 405       vtoe: dict{ int->list of int } - see _GetEdgeData
 406       points: geom.Points
 407     Returns:
 408       (list of  list of int, list of bool) -
 409         first list: each sublist is adjacent face indices for each face
 410         second list: maps edge index to True if it separates adjacent faces
 411     """
 412
 413     face_adj = [[] for i in range(len(faces))]
 414     is_interior_edge = [False] * len(edges)
 415     for e, ((vs, ve), f) in enumerate(edges):
 416         for othere in vtoe[ve]:
 417             ((_, we), g) = edges[othere]
 418             if we == vs:
 419                 # face g is adjacent to face f
 420                 # TODO: angle check
 421                 if g not in face_adj[f]:
 422                     face_adj[f].append(g)
 423                     is_interior_edge[e] = True
 424                 # Don't bother with mirror relations, will catch later
 425     return (face_adj, is_interior_edge)
 426
 427
 428 def _FindFaceGraphComponents(faces, face_adj):
 429     """Partition faces into connected components.
 430
 431     Args:
 432       faces: list of list of int
 433       face_adj: list of list of int - see _GetFaceGraph
 434     Returns:
 435       (list of list of int, list of int) -
 436         first list partitions face indices into separate lists,
 437             each a component
 438         second list maps face indices into their component index
 439     """
 440
 441     if not faces:
 442         return ([], [])
 443     components = []
 444     ftoc = [-1] * len(faces)
 445     for i in range(len(faces)):
 446         if ftoc[i] == -1:
 447             compi = len(components)
 448             comp = []
 449             _FFGCSearch(i, faces, face_adj, ftoc, compi, comp)
 450             components.append(comp)
 451     return (components, ftoc)
 452
 453
 454 def _FFGCSearch(findex, faces, face_adj, ftoc, compi, comp):
 455     """Depth first search helper function for _FindFaceGraphComponents
 456
 457     Searches recursively through all faces connected to findex, adding
 458     each face found to comp and setting ftoc for that face to compi.
 459     """
 460
 461     comp.append(findex)
 462     ftoc[findex] = compi
 463     for otherf in face_adj[findex]:
 464         if ftoc[otherf] == -1:
 465             _FFGCSearch(otherf, faces, face_adj, ftoc, compi, comp)
 466
 467
 468 def _FindOuterPoly(polys, points, faces):
 469     """Assuming polys has one CCW-oriented face when looking
 470     down average normal of faces, return that one.
 471
 472     Only one of the faces should have a normal whose dot product
 473     with the average normal of faces is positive.
 474
 475     Args:
 476       polys: list of list of int - list of polys given by vertex indices
 477       points: geom.Points
 478       faces: list of list of int - original selected region, used to find
 479           average normal
 480     Returns:
 481       int - the index in polys of the outermost one
 482     """
 483
 484     if len(polys) < 2:
 485         return 0
 486     fnorm = (0.0, 0.0, 0.0)
 487     for face in faces:
 488         if len(face) > 2:
 489             fnorm = geom.VecAdd(fnorm, geom.Newell(face, points))
 490     if fnorm == (0.0, 0.0, 0.0):
 491         return 0
 492     # fnorm is really a multiple of the normal, but fine for test below
 493     for i, poly in enumerate(polys):
 494         if len(poly) > 2:
 495             pnorm = geom.Newell(poly, points)
 496             if geom.VecDot(fnorm, pnorm) > 0:
 497                 return i
 498     print("whoops, couldn't find an outermost poly")
 499     return 0
 500
 501
 502 def _RotatedPolyAreaToXY(polyarea, norm):
 503     """Return a  PolyArea rotated to xy plane.
 504
 505     Only the points in polyarea will be transferred.
 506
 507     Args:
 508       polyarea: geom.PolyArea
 509       norm: the normal for polyarea
 510     Returns:
 511       (geom.PolyArea, (float, ..., float), dict{ int -> int }) - new PolyArea,
 512           4x3 inverse transform, dict mapping new verts to old ones
 513     """
 514
 515     # find rotation matrix that takes norm to (0,0,1)
 516     (nx, ny, nz) = norm
 517     if abs(nx) < abs(ny) and abs(nx) < abs(nz):
 518         v = (vx, vy, vz) = geom.Norm3(0.0, nz, - ny)
 519     elif abs(ny) < abs(nz):
 520         v = (vx, vy, vz) = geom.Norm3(nz, 0.0, - nx)
 521     else:
 522         v = (vx, vy, vz) = geom.Norm3(ny, - nx, 0.0)
 523     (ux, uy, uz) = geom.Cross3(v, norm)
 524     rotmat = [ux, vx, nx, uy, vy, ny, uz, vz, nz, 0.0, 0.0, 0.0]
 525     # rotation matrices are orthogonal, so inverse is transpose
 526     invrotmat = [ux, uy, uz, vx, vy, vz, nx, ny, nz, 0.0, 0.0, 0.0]
 527     pointmap = dict()
 528     invpointmap = dict()
 529     newpoints = geom.Points()
 530     for poly in [polyarea.poly] + polyarea.holes:
 531         for v in poly:
 532             vcoords = polyarea.points.pos[v]
 533             newvcoords = geom.MulPoint3(vcoords, rotmat)
 534             newv = newpoints.AddPoint(newvcoords)
 535             pointmap[v] = newv
 536             invpointmap[newv] = v
 537     pa = geom.PolyArea(newpoints)
 538     pa.poly = [pointmap[v] for v in polyarea.poly]
 539     pa.holes = [[pointmap[v] for v in hole] for hole in polyarea.holes]
 540     pa.data = polyarea.data
 541     return (pa, invrotmat, invpointmap)
 542
 543
 544 def _AddTransformedPolysToModel(mdl, polys, points, poly_data,
 545         transform, pointmap):
 546     """Add (transformed) the points and faces to a model.
 547
 548     Add polys to mdl.  The polys have coordinates given by indices
 549     into points.pos; those need to be transformed by multiplying by
 550     the transform matrix.
 551     The vertices may already exist in mdl.  Rather than relying on
 552     AddPoint to detect the duplicate (transform rounding error makes
 553     that dicey), the pointmap dictionar is used to map vertex indices
 554     in polys into those in mdl - if they exist already.
 555
 556     Args:
 557       mdl: geom.Model - where to put new vertices, faces
 558       polys: list of list of int - each sublist a poly
 559       points: geom.Points - coords for vertices in polys
 560       poly_data: list of any - parallel to polys
 561       transform: (float, ..., float) - 12-tuple, a 4x3 transform matrix
 562       pointmap: dict { int -> int } - maps new vertex indices to old ones
 563     Side Effects:
 564       The model gets new faces and vertices, based on those in polys.
 565       We are allowed to modify pointmap, as it will be discarded after call.
 566     """
 567
 568     for i, coords in enumerate(points.pos):
 569         if i not in pointmap:
 570             p = geom.MulPoint3(coords, transform)
 571             pointmap[i] = mdl.points.AddPoint(p)
 572     for i, poly in enumerate(polys):
 573         mpoly = [pointmap[v] for v in poly]
 574         mdl.faces.append(mpoly)
 575         mdl.face_data.append(poly_data[i])