mesh_inset/triquad.py

   1 # SPDX-FileCopyrightText: 2011-2022 Blender Foundation
   2 #
   3 # SPDX-License-Identifier: GPL-2.0-or-later
   4
   5
   6 from . import geom
   7 import math
   8 import random
   9 from math import sqrt, hypot
  10
  11 # Points are 3-tuples or 2-tuples of reals: (x,y,z) or (x,y)
  12 # Faces are lists of integers (vertex indices into coord lists)
  13 # After triangulation/quadrangulation, the tris and quads will
  14 # be tuples instead of lists.
  15 # Vmaps are lists taking vertex index -> Point
  16
  17 TOL = 1e-7     # a tolerance for fuzzy equality
  18 GTHRESH = 75   # threshold above which use greedy to _Quandrangulate
  19 ANGFAC = 1.0   # weighting for angles in quad goodness measure
  20 DEGFAC = 10.0  # weighting for degree in quad goodness measure
  21
  22 # Angle kind constants
  23 Ang0 = 1
  24 Angconvex = 2
  25 Angreflex = 3
  26 Angtangential = 4
  27 Ang360 = 5
  28
  29
  30 def TriangulateFace(face, points):
  31     """Triangulate the given face.
  32
  33     Uses an easy triangulation first, followed by a constrained delauney
  34     triangulation to get better shaped triangles.
  35
  36     Args:
  37       face: list of int - indices in points, assumed CCW-oriented
  38       points: geom.Points - holds coordinates for vertices
  39     Returns:
  40       list of (int, int, int) - 3-tuples are CCW-oriented vertices of
  41           triangles making up the triangulation
  42     """
  43
  44     if len(face) <= 3:
  45         return [tuple(face)]
  46     tris = EarChopTriFace(face, points)
  47     bord = _BorderEdges([face])
  48     triscdt = _CDT(tris, bord, points)
  49     return triscdt
  50
  51
  52 def TriangulateFaceWithHoles(face, holes, points):
  53     """Like TriangulateFace, but with holes inside the face.
  54
  55     Works by making one complex polygon that has segments to
  56     and from the holes ("islands"), and then using the same method
  57     as TriangulateFace.
  58
  59     Args:
  60       face: list of int - indices in points, assumed CCW-oriented
  61       holes: list of list of int - each sublist is like face
  62           but CW-oriented and assumed to be inside face
  63       points: geom.Points - holds coordinates for vertices
  64     Returns:
  65       list of (int, int, int) - 3-tuples are CCW-oriented vertices of
  66           triangles making up the triangulation
  67     """
  68
  69     if len(holes) == 0:
  70         return TriangulateFace(face, points)
  71     allfaces = [face] + holes
  72     sholes = [_SortFace(h, points) for h in holes]
  73     joinedface = _JoinIslands(face, sholes, points)
  74     tris = EarChopTriFace(joinedface, points)
  75     bord = _BorderEdges(allfaces)
  76     triscdt = _CDT(tris, bord, points)
  77     return triscdt
  78
  79
  80 def QuadrangulateFace(face, points):
  81     """Quadrangulate the face (subdivide into convex quads and tris).
  82
  83     Like TriangulateFace, but after triangulating, join as many pairs
  84     of triangles as possible into convex quadrilaterals.
  85
  86     Args:
  87       face: list of int - indices in points, assumed CCW-oriented
  88       points: geom.Points - holds coordinates for vertices
  89     Returns:
  90       list of 3-tuples or 4-tuples of ints - CCW-oriented vertices of
  91           quadrilaterals and triangles making up the quadrangulation.
  92     """
  93
  94     if len(face) <= 3:
  95         return [tuple(face)]
  96     tris = EarChopTriFace(face, points)
  97     bord = _BorderEdges([face])
  98     triscdt = _CDT(tris, bord, points)
  99     qs = _Quandrangulate(triscdt, bord, points)
 100     return qs
 101
 102
 103 def QuadrangulateFaceWithHoles(face, holes, points):
 104     """Like QuadrangulateFace, but with holes inside the faces.
 105
 106     Args:
 107       face: list of int - indices in points, assumed CCW-oriented
 108       holes: list of list of int - each sublist is like face
 109           but CW-oriented and assumed to be inside face
 110       points: geom.Points - holds coordinates for vertices
 111     Returns:
 112       list of 3-tuples or 4-tuples of ints - CCW-oriented vertices of
 113           quadrilaterals and triangles making up the quadrangulation.
 114     """
 115
 116     if len(holes) == 0:
 117         return QuadrangulateFace(face, points)
 118     allfaces = [face] + holes
 119     sholes = [_SortFace(h, points) for h in holes]
 120     joinedface = _JoinIslands(face, sholes, points)
 121     tris = EarChopTriFace(joinedface, points)
 122     bord = _BorderEdges(allfaces)
 123     triscdt = _CDT(tris, bord, points)
 124     qs = _Quandrangulate(triscdt, bord, points)
 125     return qs
 126
 127
 128 def _SortFace(face, points):
 129     """Rotate face so leftmost vertex is first, where face is
 130     list of indices in points."""
 131
 132     n = len(face)
 133     if n <= 1:
 134         return face
 135     lefti = 0
 136     leftv = face[0]
 137     for i in range(1, n):
 138         # following comparison is lexicographic on n-tuple
 139         # so sorts on x first, using lower y as tie breaker.
 140         if points.pos[face[i]] < points.pos[leftv]:
 141             lefti = i
 142             leftv = face[i]
 143     return face[lefti:] + face[0:lefti]
 144
 145
 146 def EarChopTriFace(face, points):
 147     """Triangulate given face, with coords given by indexing into points.
 148     Return list of faces, each of which will be a triangle.
 149     Use the ear-chopping method."""
 150
 151     # start with lowest coord in 2d space to try
 152     # to get a pleasing uniform triangulation if starting with
 153     # a regular structure (like a grid)
 154     start = _GetLeastIndex(face, points)
 155     ans = []
 156     incr = 1
 157     n = len(face)
 158     while n > 3:
 159         i = _FindEar(face, n, start, incr, points)
 160         vm1 = face[(i - 1) % n]
 161         v0 = face[i]
 162         v1 = face[(i + 1) % n]
 163         face = _ChopEar(face, i)
 164         n = len(face)
 165         incr = - incr
 166         if incr == 1:
 167             start = i % n
 168         else:
 169             start = (i - 1) % n
 170         ans.append((vm1, v0, v1))
 171     ans.append(tuple(face))
 172     return ans
 173
 174
 175 def _GetLeastIndex(face, points):
 176     """Return index of coordinate that is leftmost, lowest in face."""
 177
 178     bestindex = 0
 179     bestpos = points.pos[face[0]]
 180     for i in range(1, len(face)):
 181         pos = points.pos[face[i]]
 182         if pos[0] < bestpos[0] or \
 183                 (pos[0] == bestpos[0] and pos[1] < bestpos[1]):
 184             bestindex = i
 185             bestpos = pos
 186     return bestindex
 187
 188
 189 def _FindEar(face, n, start, incr, points):
 190     """An ear of a polygon consists of three consecutive vertices
 191     v(-1), v0, v1 such that v(-1) can connect to v1 without intersecting
 192     the polygon.
 193     Finds an ear, starting at index 'start' and moving
 194     in direction incr. (We attempt to alternate directions, to find
 195     'nice' triangulations for simple convex polygons.)
 196     Returns index into faces of v0 (will always find one, because
 197     uses a desperation mode if fails to find one with above rule)."""
 198
 199     angk = _ClassifyAngles(face, n, points)
 200     for mode in range(0, 5):
 201         i = start
 202         while True:
 203             if _IsEar(face, i, n, angk, points, mode):
 204                 return i
 205             i = (i + incr) % n
 206             if i == start:
 207                 break  # try next higher desperation mode
 208
 209
 210 def _IsEar(face, i, n, angk, points, mode):
 211     """Return true, false depending on ear status of vertices
 212     with indices i-1, i, i+1.
 213     mode is amount of desperation: 0 is Normal mode,
 214     mode 1 allows degenerate triangles (with repeated vertices)
 215     mode 2 allows local self crossing (folded) ears
 216     mode 3 allows any convex vertex (should always be one)
 217     mode 4 allows anything (just to be sure loop terminates!)"""
 218
 219     k = angk[i]
 220     vm2 = face[(i - 2) % n]
 221     vm1 = face[(i - 1) % n]
 222     v0 = face[i]
 223     v1 = face[(i + 1) % n]
 224     v2 = face[(i + 2) % n]
 225     if vm1 == v0 or v0 == v1:
 226         return (mode > 0)
 227     b = (k == Angconvex or k == Angtangential or k == Ang0)
 228     c = _InCone(vm1, v0, v1, v2, angk[(i + 1) % n], points) and \
 229         _InCone(v1, vm2, vm1, v0, angk[(i - 1) % n], points)
 230     if b and c:
 231         return _EarCheck(face, n, angk, vm1, v0, v1, points)
 232     if mode < 2:
 233         return False
 234     if mode == 3:
 235         return SegsIntersect(vm2, vm1, v0, v1, points)
 236     if mode == 4:
 237         return b
 238     return True
 239
 240
 241 def _EarCheck(face, n, angk, vm1, v0, v1, points):
 242     """Return True if the successive vertices vm1, v0, v1
 243     forms an ear.  We already know that it is not a reflex
 244     Angle, and that the local cone containment is ok.
 245     What remains to check is that the edge vm1-v1 doesn't
 246     intersect any other edge of the face (besides vm1-v0
 247     and v0-v1).  Equivalently, there can't be a reflex Angle
 248     inside the triangle vm1-v0-v1.  (Well, there are
 249     messy cases when other points of the face coincide with
 250     v0 or touch various lines involved in the ear.)"""
 251     for j in range(0, n):
 252         fv = face[j]
 253         k = angk[j]
 254         b = (k == Angreflex or k == Ang360) \
 255             and not(fv == vm1 or fv == v0 or fv == v1)
 256         if b:
 257             # Is fv inside closure of triangle (vm1,v0,v1)?
 258             c = not(Ccw(v0, vm1, fv, points) \
 259                           or Ccw(vm1, v1, fv, points) \
 260                           or Ccw(v1, v0, fv, points))
 261             fvm1 = face[(j - 1) % n]
 262             fv1 = face[(j + 1) % n]
 263             # To try to deal with some degenerate cases,
 264             # also check to see if either segment attached to fv
 265             # intersects either segment of potential ear.
 266             d = SegsIntersect(fvm1, fv, vm1, v0, points) or \
 267                       SegsIntersect(fvm1, fv, v0, v1, points) or \
 268                       SegsIntersect(fv, fv1, vm1, v0, points) or \
 269                       SegsIntersect(fv, fv1, v0, v1, points)
 270             if c or d:
 271                 return False
 272     return True
 273
 274
 275 def _ChopEar(face, i):
 276     """Return a copy of face (of length n), omitting element i."""
 277
 278     return face[0:i] + face[i + 1:]
 279
 280
 281 def _InCone(vtest, a, b, c, bkind, points):
 282     """Return true if point with index vtest is in Cone of points with
 283     indices a, b, c, where Angle ABC has AngleKind Bkind.
 284     The Cone is the set of points inside the left face defined by
 285     segments ab and bc, disregarding all other segments of polygon for
 286     purposes of inside test."""
 287
 288     if bkind == Angreflex or bkind == Ang360:
 289         if _InCone(vtest, c, b, a, Angconvex, points):
 290             return False
 291         return not((not(Ccw(b, a, vtest, points)) \
 292                              and not(Ccw(b, vtest, a, points)) \
 293                              and Ccw(b, a, vtest, points))
 294                             or
 295                             (not(Ccw(b, c, vtest, points)) \
 296                              and not(Ccw(b, vtest, c, points)) \
 297                              and Ccw(b, a, vtest, points)))
 298     else:
 299         return Ccw(a, b, vtest, points) and Ccw(b, c, vtest, points)
 300
 301
 302 def _JoinIslands(face, holes, points):
 303     """face is a CCW face containing the CW faces in the holes list,
 304     where each hole is sorted so the leftmost-lowest vertex is first.
 305     faces and holes are given as lists of indices into points.
 306     The holes should be sorted by softface.
 307     Add edges to make a new face that includes the holes (a Ccw traversal
 308     of the new face will have the inside always on the left),
 309     and return the new face."""
 310
 311     while len(holes) > 0:
 312         (hole, holeindex) = _LeftMostFace(holes, points)
 313         holes = holes[0:holeindex] + holes[holeindex + 1:]
 314         face = _JoinIsland(face, hole, points)
 315     return face
 316
 317
 318 def _JoinIsland(face, hole, points):
 319     """Return a modified version of face that splices in the
 320     vertices of hole (which should be sorted)."""
 321
 322     if len(hole) == 0:
 323         return face
 324     hv0 = hole[0]
 325     d = _FindDiag(face, hv0, points)
 326     newface = face[0:d + 1] + hole + [hv0] + face[d:]
 327     return newface
 328
 329
 330 def _LeftMostFace(holes, points):
 331     """Return (hole,index of hole in holes) where hole has
 332     the leftmost first vertex.  To be able to handle empty
 333     holes gracefully, call an empty hole 'leftmost'.
 334     Assumes holes are sorted by softface."""
 335
 336     assert(len(holes) > 0)
 337     lefti = 0
 338     lefthole = holes[0]
 339     if len(lefthole) == 0:
 340         return (lefthole, lefti)
 341     leftv = lefthole[0]
 342     for i in range(1, len(holes)):
 343         ihole = holes[i]
 344         if len(ihole) == 0:
 345             return (ihole, i)
 346         iv = ihole[0]
 347         if points.pos[iv] < points.pos[leftv]:
 348             (lefti, lefthole, leftv) = (i, ihole, iv)
 349     return (lefthole, lefti)
 350
 351
 352 def _FindDiag(face, hv, points):
 353     """Find a vertex in face that can see vertex hv, if possible,
 354     and return the index into face of that vertex.
 355     Should be able to find a diagonal that connects a vertex of face
 356     left of v to hv without crossing face, but try two
 357     more desperation passes after that to get SOME diagonal, even if
 358     it might cross some edge somewhere.
 359     First desperation pass (mode == 1): allow points right of hv.
 360     Second desperation pass (mode == 2): allow crossing boundary poly"""
 361
 362     besti = - 1
 363     bestdist = 1e30
 364     for mode in range(0, 3):
 365         for i in range(0, len(face)):
 366             v = face[i]
 367             if mode == 0 and points.pos[v] > points.pos[hv]:
 368                 continue  # in mode 0, only want points left of hv
 369             dist = _DistSq(v, hv, points)
 370             if dist < bestdist:
 371                 if _IsDiag(i, v, hv, face, points) or mode == 2:
 372                     (besti, bestdist) = (i, dist)
 373         if besti >= 0:
 374             break  # found one, so don't need other modes
 375     assert(besti >= 0)
 376     return besti
 377
 378
 379 def _IsDiag(i, v, hv, face, points):
 380     """Return True if vertex v (at index i in face) can see vertex hv.
 381     v and hv are indices into points.
 382     (v, hv) is a diagonal if hv is in the cone of the Angle at index i on face
 383     and no segment in face intersects (h, hv).
 384     """
 385
 386     n = len(face)
 387     vm1 = face[(i - 1) % n]
 388     v1 = face[(i + 1) % n]
 389     k = _AngleKind(vm1, v, v1, points)
 390     if not _InCone(hv, vm1, v, v1, k, points):
 391         return False
 392     for j in range(0, n):
 393         vj = face[j]
 394         vj1 = face[(j + 1) % n]
 395         if SegsIntersect(v, hv, vj, vj1, points):
 396             return False
 397     return True
 398
 399
 400 def _DistSq(a, b, points):
 401     """Return distance squared between coords with indices a and b in points.
 402     """
 403
 404     diff = Sub2(points.pos[a], points.pos[b])
 405     return Dot2(diff, diff)
 406
 407
 408 def _BorderEdges(facelist):
 409     """Return a set of (u,v) where u and v are successive vertex indices
 410     in some face in the list in facelist."""
 411
 412     ans = set()
 413     for i in range(0, len(facelist)):
 414         f = facelist[i]
 415         for j in range(1, len(f)):
 416             ans.add((f[j - 1], f[j]))
 417         ans.add((f[-1], f[0]))
 418     return ans
 419
 420
 421 def _CDT(tris, bord, points):
 422     """Tris is a list of triangles ((a,b,c), CCW-oriented indices into points)
 423     Bord is a set of border edges (u,v), oriented so that tris
 424     is a triangulation of the left face of the border(s).
 425     Make the triangulation "Constrained Delaunay" by flipping "reversed"
 426     quadrangulaterals until can flip no more.
 427     Return list of triangles in new triangulation."""
 428
 429     td = _TriDict(tris)
 430     re = _ReveresedEdges(tris, td, bord, points)
 431     ts = set(tris)
 432     # reverse the reversed edges until done.
 433     # reversing and edge adds new edges, which may or
 434     # may not be reversed or border edges, to re for
 435     # consideration, but the process will stop eventually.
 436     while len(re) > 0:
 437         (a, b) = e = re.pop()
 438         if e in bord or not _IsReversed(e, td, points):
 439             continue
 440         # rotate e in quad adbc to get other diagonal
 441         erev = (b, a)
 442         tl = td.get(e)
 443         tr = td.get(erev)
 444         if not tl or not tr:
 445             continue  # shouldn't happen
 446         c = _OtherVert(tl, a, b)
 447         d = _OtherVert(tr, a, b)
 448         if c is None or d is None:
 449             continue  # shouldn't happen
 450         newt1 = (c, d, b)
 451         newt2 = (c, a, d)
 452         del td[e]
 453         del td[erev]
 454         td[(c, d)] = newt1
 455         td[(d, b)] = newt1
 456         td[(b, c)] = newt1
 457         td[(d, c)] = newt2
 458         td[(c, a)] = newt2
 459         td[(a, d)] = newt2
 460         if tl in ts:
 461             ts.remove(tl)
 462         if tr in ts:
 463             ts.remove(tr)
 464         ts.add(newt1)
 465         ts.add(newt2)
 466         re.extend([(d, b), (b, c), (c, a), (a, d)])
 467     return list(ts)
 468
 469
 470 def _TriDict(tris):
 471     """tris is a list of triangles (a,b,c), CCW-oriented indices.
 472     Return dict mapping all edges in the triangles to the containing
 473     triangle list."""
 474
 475     ans = dict()
 476     for i in range(0, len(tris)):
 477         (a, b, c) = t = tris[i]
 478         ans[(a, b)] = t
 479         ans[(b, c)] = t
 480         ans[(c, a)] = t
 481     return ans
 482
 483
 484 def _ReveresedEdges(tris, td, bord, points):
 485     """Return list of reversed edges in tris.
 486     Only want edges not in bord, and only need one representative
 487     of (u,v)/(v,u), so choose the one with u < v.
 488     td is dictionary from _TriDict, and is used to find left and right
 489     triangles of edges."""
 490
 491     ans = []
 492     for i in range(0, len(tris)):
 493         (a, b, c) = tris[i]
 494         for e in [(a, b), (b, c), (c, a)]:
 495             if e in bord:
 496                 continue
 497             (u, v) = e
 498             if u < v:
 499                 if _IsReversed(e, td, points):
 500                     ans.append(e)
 501     return ans
 502
 503
 504 def _IsReversed(e, td, points):
 505     """If e=(a,b) is a non-border edge, with left-face triangle tl and
 506     right-face triangle tr, then it is 'reversed' if the circle through
 507     a, b, and (say) the other vertex of tl contains the other vertex of tr.
 508     td is a _TriDict, for finding triangles containing edges, and points
 509     gives the coordinates for vertex indices used in edges."""
 510
 511     tl = td.get(e)
 512     if not tl:
 513         return False
 514     (a, b) = e
 515     tr = td.get((b, a))
 516     if not tr:
 517         return False
 518     c = _OtherVert(tl, a, b)
 519     d = _OtherVert(tr, a, b)
 520     if c is None or d is None:
 521         return False
 522     return InCircle(a, b, c, d, points)
 523
 524
 525 def _OtherVert(tri, a, b):
 526     """tri should be a tuple of 3 vertex indices, two of which are a and b.
 527     Return the third index, or None if all vertices are a or b"""
 528
 529     for v in tri:
 530         if v != a and v != b:
 531             return v
 532     return None
 533
 534
 535 def _ClassifyAngles(face, n, points):
 536     """Return vector of anglekinds of the Angle around each point in face."""
 537
 538     return [_AngleKind(face[(i - 1) % n], face[i], face[(i + 1) % n], points) \
 539         for i in list(range(0, n))]
 540
 541
 542 def _AngleKind(a, b, c, points):
 543     """Return one of the Ang... constants to classify Angle formed by ABC,
 544     in a counterclockwise traversal of a face,
 545     where a, b, c are indices into points."""
 546
 547     if Ccw(a, b, c, points):
 548         return Angconvex
 549     elif Ccw(a, c, b, points):
 550         return Angreflex
 551     else:
 552         vb = points.pos[b]
 553         udotv = Dot2(Sub2(vb, points.pos[a]), Sub2(points.pos[c], vb))
 554         if udotv > 0.0:
 555             return Angtangential
 556         else:
 557             return Ang0   # to fix: return Ang360 if "inside" spur
 558
 559
 560 def _Quandrangulate(tris, bord, points):
 561     """Tris is list of triangles, forming a triangulation of region whose
 562     border edges are in set bord.
 563     Combine adjacent triangles to make quads, trying for "good" quads where
 564     possible. Some triangles will probably remain uncombined"""
 565
 566     (er, td) = _ERGraph(tris, bord, points)
 567     if len(er) == 0:
 568         return tris
 569     if len(er) > GTHRESH:
 570         match = _GreedyMatch(er)
 571     else:
 572         match = _MaxMatch(er)
 573     return _RemoveEdges(tris, match)
 574
 575
 576 def _RemoveEdges(tris, match):
 577     """tris is list of triangles.
 578     er is as returned from _MaxMatch or _GreedyMatch.
 579
 580     Return list of (A,D,B,C) resulting from deleting edge (A,B) causing a merge
 581     of two triangles; append to that list the remaining unmatched triangles."""
 582
 583     ans = []
 584     triset = set(tris)
 585     while len(match) > 0:
 586         (_, e, tl, tr) = match.pop()
 587         (a, b) = e
 588         if tl in triset:
 589             triset.remove(tl)
 590         if tr in triset:
 591             triset.remove(tr)
 592         c = _OtherVert(tl, a, b)
 593         d = _OtherVert(tr, a, b)
 594         if c is None or d is None:
 595             continue
 596         ans.append((a, d, b, c))
 597     return ans + list(triset)
 598
 599
 600 def _ERGraph(tris, bord, points):
 601     """Make an 'Edge Removal Graph'.
 602
 603     Given a list of triangles, the 'Edge Removal Graph' is a graph whose
 604     nodes are the triangles (think of a point in the center of them),
 605     and whose edges go between adjacent triangles (they share a non-border
 606     edge), such that it would be possible to remove the shared edge
 607     and form a convex quadrilateral.  Forming a quadrilateralization
 608     is then a matter of finding a matching (set of edges that don't
 609     share a vertex - remember, these are the 'face' vertices).
 610     For better quadrilaterlization, we'll make the Edge Removal Graph
 611     edges have weights, with higher weights going to the edges that
 612     are more desirable to remove.  Then we want a maximum weight matching
 613     in this graph.
 614
 615     We'll return the graph in a kind of implicit form, using edges of
 616     the original triangles as a proxy for the edges between the faces
 617     (i.e., the edge of the triangle is the shared edge). We'll arbitrarily
 618     pick the triangle graph edge with lower-index start vertex.
 619     Also, to aid in traversing the implicit graph, we'll keep the left
 620     and right triangle triples with edge 'ER edge'.
 621     Finally, since we calculate it anyway, we'll return a dictionary
 622     mapping edges of the triangles to the triangle triples they're in.
 623
 624     Args:
 625       tris: list of (int, int, int) giving a triple of vertex indices for
 626           triangles, assumed CCW oriented
 627       bord: set of (int, int) giving vertex indices for border edges
 628       points: geom.Points - for mapping vertex indices to coords
 629     Returns:
 630       (list of (weight,e,tl,tr), dict)
 631         where edge e=(a,b) is non-border edge
 632         with left face tl and right face tr (each a triple (i,j,k)),
 633         where removing the edge would form an "OK" quad (no concave angles),
 634         with weight representing the desirability of removing the edge
 635         The dict maps int pairs (a,b) to int triples (i,j,k), that is,
 636         mapping edges to their containing triangles.
 637     """
 638
 639     td = _TriDict(tris)
 640     dd = _DegreeDict(tris)
 641     ans = []
 642     ctris = tris[:]  # copy, so argument not affected
 643     while len(ctris) > 0:
 644         (i, j, k) = tl = ctris.pop()
 645         for e in [(i, j), (j, k), (k, i)]:
 646             if e in bord:
 647                 continue
 648             (a, b) = e
 649             # just consider one of (a,b) and (b,a), to avoid dups
 650             if a > b:
 651                 continue
 652             erev = (b, a)
 653             tr = td.get(erev)
 654             if not tr:
 655                 continue
 656             c = _OtherVert(tl, a, b)
 657             d = _OtherVert(tr, a, b)
 658             if c is None or d is None:
 659                 continue
 660             # calculate amax, the max of the new angles that would
 661             # be formed at a and b if tl and tr were combined
 662             amax = max(Angle(c, a, b, points) + Angle(d, a, b, points),
 663                        Angle(c, b, a, points) + Angle(d, b, a, points))
 664             if amax > 180.0:
 665                 continue
 666             weight = ANGFAC * (180.0 - amax) + DEGFAC * (dd[a] + dd[b])
 667             ans.append((weight, e, tl, tr))
 668     return (ans, td)
 669
 670
 671 def _GreedyMatch(er):
 672     """er is list of (weight,e,tl,tr).
 673
 674     Find maximal set so that each triangle appears in at most
 675     one member of set"""
 676
 677     # sort in order of decreasing weight
 678     er.sort(key=lambda v: v[0], reverse=True)
 679     match = set()
 680     ans = []
 681     while len(er) > 0:
 682         (_, _, tl, tr) = q = er.pop()
 683         if tl not in match and tr not in match:
 684             match.add(tl)
 685             match.add(tr)
 686             ans.append(q)
 687     return ans
 688
 689
 690 def _MaxMatch(er):
 691     """Like _GreedyMatch, but use divide and conquer to find best possible set.
 692
 693     Args:
 694       er: list of (weight,e,tl,tr)  - see _ERGraph
 695     Returns:
 696       list that is a subset of er giving a maximum weight match
 697     """
 698
 699     (ans, _) = _DCMatch(er)
 700     return ans
 701
 702
 703 def _DCMatch(er):
 704     """Recursive helper for _MaxMatch.
 705
 706     Divide and Conquer approach to finding max weight matching.
 707     If we're lucky, there's an edge in er that separates the edge removal
 708     graph into (at least) two separate components.  Then the max weight
 709     is either one that includes that edge or excludes it - and we can
 710     use a recursive call to _DCMatch to handle each component separately
 711     on what remains of the graph after including/excluding the separating edge.
 712     If we're not lucky, we fall back on _EMatch (see below).
 713
 714     Args:
 715       er: list of (weight, e, tl, tr) (see _ERGraph)
 716     Returns:
 717       (list of (weight, e, tl, tr), float) - the subset forming a maximum
 718           matching, and the total weight of the match.
 719     """
 720
 721     if not er:
 722         return ([], 0.0)
 723     if len(er) == 1:
 724         return (er, er[0][0])
 725     match = []
 726     matchw = 0.0
 727     for i in range(0, len(er)):
 728         (nc, comp) = _FindComponents(er, i)
 729         if nc == 1:
 730             # er[i] doesn't separate er
 731             continue
 732         (wi, _, tl, tr) = er[i]
 733         if comp[tl] != comp[tr]:
 734             # case 1: er separates graph
 735             # compare the matches that include er[i] versus
 736             # those that exclude it
 737             (a, b) = _PartitionComps(er, comp, i, comp[tl], comp[tr])
 738             ax = _CopyExcluding(a, tl, tr)
 739             bx = _CopyExcluding(b, tl, tr)
 740             (axmatch, wax) = _DCMatch(ax)
 741             (bxmatch, wbx) = _DCMatch(bx)
 742             if len(ax) == len(a):
 743                 wa = wax
 744                 amatch = axmatch
 745             else:
 746                 (amatch, wa) = _DCMatch(a)
 747             if len(bx) == len(b):
 748                 wb = wbx
 749                 bmatch = bxmatch
 750             else:
 751                 (bmatch, wb) = _DCMatch(b)
 752             w = wa + wb
 753             wx = wax + wbx + wi
 754             if w > wx:
 755                 match = amatch + bmatch
 756                 matchw = w
 757             else:
 758                 match = [er[i]] + axmatch + bxmatch
 759                 matchw = wx
 760         else:
 761             # case 2: er not needed to separate graph
 762             (a, b) = _PartitionComps(er, comp, -1, 0, 0)
 763             (amatch, wa) = _DCMatch(a)
 764             (bmatch, wb) = _DCMatch(b)
 765             match = amatch + bmatch
 766             matchw = wa + wb
 767         if match:
 768             break
 769     if not match:
 770         return _EMatch(er)
 771     return (match, matchw)
 772
 773
 774 def _EMatch(er):
 775     """Exhaustive match helper for _MaxMatch.
 776
 777     This is the case when we were unable to find a single edge
 778     separating the edge removal graph into two components.
 779     So pick a single edge and try _DCMatch on the two cases of
 780     including or excluding that edge.  We may be lucky in these
 781     subcases (say, if the graph is currently a simple cycle, so
 782     only needs one more edge after the one we pick here to separate
 783     it into components).  Otherwise, we'll end up back in _EMatch
 784     again, and the worse case will be exponential.
 785
 786     Pick a random edge rather than say, the first, to hopefully
 787     avoid some pathological cases.
 788
 789     Args:
 790       er: list of (weight, el, tl, tr) (see _ERGraph)
 791     Returns:
 792        (list of (weight, e, tl, tr), float) - the subset forming a maximum
 793           matching, and the total weight of the match.
 794     """
 795
 796     if not er:
 797         return ([], 0.0)
 798     if len(er) == 1:
 799         return (er, er[1][1])
 800     i = random.randint(0, len(er) - 1)
 801     eri = (wi, _, tl, tr) = er[i]
 802     # case a: include eri. exclude other edges that touch tl or tr
 803     a = _CopyExcluding(er, tl, tr)
 804     a.append(eri)
 805     (amatch, wa) = _DCMatch(a)
 806     wa += wi
 807     if len(a) == len(er) - 1:
 808         # if a excludes only eri, then er didn't touch anything else
 809         # in the graph, and the best match will always include er
 810         # and we can skip the call for case b
 811         wb = -1.0
 812         bmatch = []
 813     else:
 814         b = er[:i] + er[i + 1:]
 815         (bmatch, wb) = _DCMatch(b)
 816     if wa > wb:
 817         match = amatch
 818         match.append(eri)
 819         matchw = wa
 820     else:
 821         match = bmatch
 822         matchw = wb
 823     return (match, matchw)
 824
 825
 826 def _FindComponents(er, excepti):
 827     """Find connected components induced by edges, excluding one edge.
 828
 829     Args:
 830       er: list of (weight, el, tl, tr) (see _ERGraph)
 831       excepti: index in er of edge to be excluded
 832     Returns:
 833       (int, dict): int is number of connected components found,
 834           dict maps triangle triple ->
 835               connected component index (starting at 1)
 836      """
 837
 838     ncomps = 0
 839     comp = dict()
 840     for i in range(0, len(er)):
 841         (_, _, tl, tr) = er[i]
 842         for t in [tl, tr]:
 843             if t not in comp:
 844                 ncomps += 1
 845                 _FCVisit(er, excepti, comp, t, ncomps)
 846     return (ncomps, comp)
 847
 848
 849 def _FCVisit(er, excepti, comp, t, compnum):
 850     """Helper for _FindComponents depth-first-search."""
 851
 852     comp[t] = compnum
 853     for i in range(0, len(er)):
 854         if i == excepti:
 855             continue
 856         (_, _, tl, tr) = er[i]
 857         if tl == t or tr == t:
 858             s = tl
 859             if s == t:
 860                 s = tr
 861             if s not in comp:
 862                 _FCVisit(er, excepti, comp, s, compnum)
 863
 864
 865 def _PartitionComps(er, comp, excepti, compa, compb):
 866     """Partition the edges of er by component number, into two lists.
 867
 868     Generally, put odd components into first list and even into second,
 869     except that component compa goes in the first and compb goes in the second,
 870     and we ignore edge er[excepti].
 871
 872     Args:
 873       er: list of (weight, el, tl, tr) (see _ERGraph)
 874       comp: dict - mapping triangle triple -> connected component index
 875       excepti: int - index in er to ignore (unless excepti==-1)
 876       compa: int - component to go in first list of answer (unless 0)
 877       compb: int - component to go in second list of answer (unless 0)
 878     Returns:
 879       (list, list) - a partition of er according to above rules
 880     """
 881
 882     parta = []
 883     partb = []
 884     for i in range(0, len(er)):
 885
 886         if i == excepti:
 887             continue
 888         tl = er[i][2]
 889         c = comp[tl]
 890         if c == compa or (c != compb and (c & 1) == 1):
 891             parta.append(er[i])
 892         else:
 893             partb.append(er[i])
 894     return (parta, partb)
 895
 896
 897 def _CopyExcluding(er, s, t):
 898     """Return a copy of er, excluding all those involving triangles s and t.
 899
 900     Args:
 901       er: list of (weight, e, tl, tr) - see _ERGraph
 902       s: 3-tuple of int - a triangle
 903       t: 3-tuple of int - a triangle
 904     Returns:
 905       Copy of er excluding those with tl or tr == s or t
 906     """
 907
 908     ans = []
 909     for e in er:
 910         (_, _, tl, tr) = e
 911         if tl == s or tr == s or tl == t or tr == t:
 912             continue
 913         ans.append(e)
 914     return ans
 915
 916
 917 def _DegreeDict(tris):
 918     """Return a dictionary mapping vertices in tris to the number of triangles
 919     that they are touch."""
 920
 921     ans = dict()
 922     for t in tris:
 923         for v in t:
 924             if v in ans:
 925                 ans[v] = ans[v] + 1
 926             else:
 927                 ans[v] = 1
 928     return ans
 929
 930
 931 def PolygonPlane(face, points):
 932     """Return a Normal vector for the face with 3d coords given by indexing
 933     into points."""
 934
 935     if len(face) < 3:
 936         return (0.0, 0.0, 1.0)    # arbitrary, we really have no idea
 937     else:
 938         coords = [points.pos[i] for i in face]
 939         return Normal(coords)
 940
 941
 942 # This Normal appears to be on the CCW-traversing side of a polygon
 943 def Normal(coords):
 944     """Return an average Normal vector for the point list, 3d coords."""
 945
 946     if len(coords) < 3:
 947         return (0.0, 0.0, 1.0)    # arbitrary
 948
 949     (ax, ay, az) = coords[0]
 950     (bx, by, bz) = coords[1]
 951     (cx, cy, cz) = coords[2]
 952
 953     if len(coords) == 3:
 954         sx = (ay - by) * (az + bz) + \
 955             (by - cy) * (bz + cz) + \
 956             (cy - ay) * (cz + az)
 957         sy = (az - bz) * (ax + bx) + \
 958             (bz - cz) * (bx + cx) + \
 959             (cz - az) * (cx + ax)
 960         sz = (ax - bx) * (by + by) + \
 961             (bx - cx) * (by + cy) + \
 962             (cx - ax) * (cy + ay)
 963         return Norm3(sx, sy, sz)
 964     else:
 965         sx = (ay - by) * (az + bz) + (by - cy) * (bz + cz)
 966         sy = (az - bz) * (ax + bx) + (bz - cz) * (bx + cx)
 967         sz = (ax - bx) * (ay + by) + (bx - cx) * (by + cy)
 968         return _NormalAux(coords[3:], coords[0], sx, sy, sz)
 969
 970
 971 def _NormalAux(rest, first, sx, sy, sz):
 972     (ax, ay, az) = rest[0]
 973     if len(rest) == 1:
 974         (bx, by, bz) = first
 975     else:
 976         (bx, by, bz) = rest[1]
 977     nx = sx + (ay - by) * (az + bz)
 978     ny = sy + (az - bz) * (ax + bx)
 979     nz = sz + (ax - bx) * (ay + by)
 980     if len(rest) == 1:
 981         return Norm3(nx, ny, nz)
 982     else:
 983         return _NormalAux(rest[1:], first, nx, ny, nz)
 984
 985
 986 def Norm3(x, y, z):
 987     """Return vector (x,y,z) normalized by dividing by squared length.
 988     Return (0.0, 0.0, 1.0) if the result is undefined."""
 989     sqrlen = x * x + y * y + z * z
 990     if sqrlen < 1e-100:
 991         return (0.0, 0.0, 1.0)
 992     else:
 993         try:
 994             d = sqrt(sqrlen)
 995             return (x / d, y / d, z / d)
 996         except:
 997             return (0.0, 0.0, 1.0)
 998
 999
1000 # We're using right-hand coord system, where
1001 # forefinger=x, middle=y, thumb=z on right hand.
1002 # Then, e.g., (1,0,0) x (0,1,0) = (0,0,1)
1003 def Cross3(a, b):
1004     """Return the cross product of two vectors, a x b."""
1005
1006     (ax, ay, az) = a
1007     (bx, by, bz) = b
1008     return (ay * bz - az * by, az * bx - ax * bz, ax * by - ay * bx)
1009
1010
1011 def Dot2(a, b):
1012     """Return the dot product of two 2d vectors, a . b."""
1013
1014     return a[0] * b[0] + a[1] * b[1]
1015
1016
1017 def Perp2(a, b):
1018     """Return a sort of 2d cross product."""
1019
1020     return a[0] * b[1] - a[1] * b[0]
1021
1022
1023 def Sub2(a, b):
1024     """Return difference of 2d vectors, a-b."""
1025
1026     return (a[0] - b[0], a[1] - b[1])
1027
1028
1029 def Add2(a, b):
1030     """Return the sum of 2d vectors, a+b."""
1031
1032     return (a[0] + b[0], a[1] + b[1])
1033
1034
1035 def Length2(v):
1036     """Return length of vector v=(x,y)."""
1037
1038     return hypot(v[0], v[1])
1039
1040
1041 def LinInterp2(a, b, alpha):
1042     """Return the point alpha of the way from a to b."""
1043
1044     beta = 1 - alpha
1045     return (beta * a[0] + alpha * b[0], beta * a[1] + alpha * b[1])
1046
1047
1048 def Normalized2(p):
1049     """Return vector p normlized by dividing by its squared length.
1050     Return (0.0, 1.0) if the result is undefined."""
1051
1052     (x, y) = p
1053     sqrlen = x * x + y * y
1054     if sqrlen < 1e-100:
1055         return (0.0, 1.0)
1056     else:
1057         try:
1058             d = sqrt(sqrlen)
1059             return (x / d, y / d)
1060         except:
1061             return (0.0, 1.0)
1062
1063
1064 def Angle(a, b, c, points):
1065     """Return Angle abc in degrees, in range [0,180),
1066     where a,b,c are indices into points."""
1067
1068     u = Sub2(points.pos[c], points.pos[b])
1069     v = Sub2(points.pos[a], points.pos[b])
1070     n1 = Length2(u)
1071     n2 = Length2(v)
1072     if n1 == 0.0 or n2 == 0.0:
1073         return 0.0
1074     else:
1075         costheta = Dot2(u, v) / (n1 * n2)
1076         if costheta > 1.0:
1077             costheta = 1.0
1078         if costheta < - 1.0:
1079             costheta = - 1.0
1080         return math.acos(costheta) * 180.0 / math.pi
1081
1082
1083 def SegsIntersect(ixa, ixb, ixc, ixd, points):
1084     """Return true if segment AB intersects CD,
1085     false if they just touch.  ixa, ixb, ixc, ixd are indices
1086     into points."""
1087
1088     a = points.pos[ixa]
1089     b = points.pos[ixb]
1090     c = points.pos[ixc]
1091     d = points.pos[ixd]
1092     u = Sub2(b, a)
1093     v = Sub2(d, c)
1094     w = Sub2(a, c)
1095     pp = Perp2(u, v)
1096     if abs(pp) > TOL:
1097         si = Perp2(v, w) / pp
1098         ti = Perp2(u, w) / pp
1099         return 0.0 < si < 1.0 and 0.0 < ti < 1.0
1100     else:
1101         # parallel or overlapping
1102         if Dot2(u, u) == 0.0 or Dot2(v, v) == 0.0:
1103             return False
1104         else:
1105             pp2 = Perp2(w, v)
1106             if abs(pp2) > TOL:
1107                 return False  # parallel, not collinear
1108             z = Sub2(b, c)
1109             (vx, vy) = v
1110             (wx, wy) = w
1111             (zx, zy) = z
1112             if vx == 0.0:
1113                 (t0, t1) = (wy / vy, zy / vy)
1114             else:
1115                 (t0, t1) = (wx / vx, zx / vx)
1116             return 0.0 < t0 < 1.0 and 0.0 < t1 < 1.0
1117
1118
1119 def Ccw(a, b, c, points):
1120     """Return true if ABC is a counterclockwise-oriented triangle,
1121     where a, b, and c are indices into points.
1122     Returns false if not, or if colinear within TOL."""
1123
1124     (ax, ay) = (points.pos[a][0], points.pos[a][1])
1125     (bx, by) = (points.pos[b][0], points.pos[b][1])
1126     (cx, cy) = (points.pos[c][0], points.pos[c][1])
1127     d = ax * by - bx * ay - ax * cy + cx * ay + bx * cy - cx * by
1128     return d > TOL
1129
1130
1131 def InCircle(a, b, c, d, points):
1132     """Return true if circle through points with indices a, b, c
1133     contains point with index d (indices into points).
1134     Except: if ABC forms a counterclockwise oriented triangle
1135     then the test is reversed: return true if d is outside the circle.
1136     Will get false, no matter what orientation, if d is cocircular, with TOL^2.
1137       | xa ya xa^2+ya^2 1 |
1138       | xb yb xb^2+yb^2 1 | > 0
1139       | xc yc xc^2+yc^2 1 |
1140       | xd yd xd^2+yd^2 1 |
1141     """
1142
1143     (xa, ya, za) = _Icc(points.pos[a])
1144     (xb, yb, zb) = _Icc(points.pos[b])
1145     (xc, yc, zc) = _Icc(points.pos[c])
1146     (xd, yd, zd) = _Icc(points.pos[d])
1147     det = xa * (yb * zc - yc * zb - yb * zd + yd * zb + yc * zd - yd * zc) \
1148           - xb * (ya * zc - yc * za - ya * zd + yd * za + yc * zd - yd * zc) \
1149           + xc * (ya * zb - yb * za - ya * zd + yd * za + yb * zd - yd * zb) \
1150           - xd * (ya * zb - yb * za - ya * zc + yc * za + yb * zc - yc * zb)
1151     return det > TOL * TOL
1152
1153
1154 def _Icc(p):
1155     (x, y) = (p[0], p[1])
1156     return (x, y, x * x + y * y)