mesh_inset/model.py

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