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