mesh_inset/geom.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 """Geometry classes and operations.
  22 Also, vector file representation (Art).
  23 """
  24
  25 __author__ = "howard.trickey@gmail.com"
  26
  27 import math
  28
  29 # distances less than about DISTTOL will be considered
  30 # essentially zero
  31 DISTTOL = 1e-3
  32 INVDISTTOL = 1e3
  33
  34
  35 class Points(object):
  36     """Container of points without duplication, each mapped to an int.
  37
  38     Points are either have dimension at least 2, maybe more.
  39
  40     Implementation:
  41     In order to efficiently find duplicates, we quantize the points
  42     to triples of ints and map from quantized triples to vertex
  43     index.
  44
  45     Attributes:
  46       pos: list of tuple of float - coordinates indexed by
  47           vertex number
  48       invmap: dict of (int, int, int) to int - quantized coordinates
  49           to vertex number map
  50     """
  51
  52     def __init__(self, initlist=[]):
  53         self.pos = []
  54         self.invmap = dict()
  55         for p in initlist:
  56             self.AddPoint(p)
  57
  58     @staticmethod
  59     def Quantize(p):
  60         """Quantize the float tuple into an int tuple.
  61
  62         Args:
  63           p: tuple of float
  64         Returns:
  65           tuple of int - scaled by INVDISTTOL and rounded p
  66         """
  67
  68         return tuple([int(round(v * INVDISTTOL)) for v in p])
  69
  70     def AddPoint(self, p, allowdups = False):
  71         """Add point p to the Points set and return vertex number.
  72
  73         If there is an existing point which quantizes the same,,
  74         don't add a new one but instead return existing index.
  75         Except if allowdups is True, don't do that deduping.
  76
  77         Args:
  78           p: tuple of float - coordinates (2-tuple or 3-tuple)
  79         Returns:
  80           int - the vertex number of added (or existing) point
  81         """
  82
  83         qp = Points.Quantize(p)
  84         if qp in self.invmap and not allowdups:
  85             return self.invmap[qp]
  86         else:
  87             self.invmap[qp] = len(self.pos)
  88             self.pos.append(p)
  89             return len(self.pos) - 1
  90
  91     def AddPoints(self, points, allowdups = False):
  92         """Add another set of points to this set.
  93
  94         We need to return a mapping from indices
  95         in the argument points space into indices
  96         in this point space.
  97
  98         Args:
  99           points: Points - to union into this set
 100         Returns:
 101           list of int: maps added indices to new ones
 102         """
 103
 104         vmap = [0] * len(points.pos)
 105         for i in range(len(points.pos)):
 106             vmap[i] = self.AddPoint(points.pos[i], allowdups)
 107         return vmap
 108
 109     def AddZCoord(self, z):
 110         """Change this in place to have a z coordinate, with value z.
 111
 112         Assumes the coordinates are currently 2d.
 113
 114         Args:
 115           z: the value of the z coordinate to add
 116         Side Effect:
 117           self now has a z-coordinate added
 118         """
 119
 120         assert(len(self.pos) == 0 or len(self.pos[0]) == 2)
 121         newinvmap = dict()
 122         for i, (x, y) in enumerate(self.pos):
 123             newp = (x, y, z)
 124             self.pos[i] = newp
 125             newinvmap[self.Quantize(newp)] = i
 126         self.invmap = newinvmap
 127
 128     def AddToZCoord(self, i, delta):
 129         """Change the z-coordinate of point with index i to add delta.
 130
 131         Assumes the coordinates are currently 3d.
 132
 133         Args:
 134           i: int - index of a point
 135           delta: float - value to add to z-coord
 136         """
 137
 138         (x, y, z) = self.pos[i]
 139         self.pos[i] = (x, y, z + delta)
 140
 141
 142 class PolyArea(object):
 143     """Contains a Polygonal Area (polygon with possible holes).
 144
 145     A polygon is a list of vertex ids, each an index given by
 146     a Points object. The list represents a CCW-oriented
 147     outer boundary (implicitly closed).
 148     If there are holes, they are lists of CW-oriented vertices
 149     that should be contained in the outer boundary.
 150     (So the left face of both the poly and the holes is
 151     the filled part.)
 152
 153     Attributes:
 154       points: Points
 155       poly: list of vertex ids
 156       holes: list of lists of vertex ids (each a hole in poly)
 157       data: any - application data (can hold color, e.g.)
 158     """
 159
 160     def __init__(self, points=None, poly=None, holes=None, data=None):
 161         self.points = points if points else Points()
 162         self.poly = poly if poly else []
 163         self.holes = holes if holes else []
 164         self.data = data
 165
 166     def AddHole(self, holepa):
 167         """Add a PolyArea's poly as a hole of self.
 168
 169         Need to reverse the contour and
 170         adjust the the point indexes and self.points.
 171
 172         Args:
 173           holepa: PolyArea
 174         """
 175
 176         vmap = self.points.AddPoints(holepa.points)
 177         holepoly = [vmap[i] for i in holepa.poly]
 178         holepoly.reverse()
 179         self.holes.append(holepoly)
 180
 181     def ContainsPoly(self, poly, points):
 182         """Tests if poly is contained within self.poly.
 183
 184         Args:
 185           poly: list of int - indices into points
 186           points: Points - maps to coords
 187         Returns:
 188           bool - True if poly is fully contained within self.poly
 189         """
 190
 191         for v in poly:
 192             if PointInside(points.pos[v], self.poly, self.points) == -1:
 193                 return False
 194         return True
 195
 196     def Normal(self):
 197         """Returns the normal of the polyarea's main poly."""
 198
 199         pos = self.points.pos
 200         poly = self.poly
 201         if len(pos) == 0 or len(pos[0]) == 2 or len(poly) == 0:
 202             print("whoops, not enough info to calculate normal")
 203             return (0.0, 0.0, 1.0)
 204         return Newell(poly, self.points)
 205
 206
 207 class PolyAreas(object):
 208     """Contains a list of PolyAreas and a shared Points.
 209
 210     Attributes:
 211       polyareas: list of PolyArea
 212       points: Points
 213     """
 214
 215     def __init__(self):
 216         self.polyareas = []
 217         self.points = Points()
 218
 219     def scale_and_center(self, scaled_side_target):
 220         """Adjust the coordinates of the polyareas so that
 221         it is centered at the origin and has its longest
 222         dimension scaled to be scaled_side_target."""
 223
 224         if len(self.points.pos) == 0:
 225             return
 226         (minv, maxv) = self.bounds()
 227         maxside = max([maxv[i] - minv[i] for i in range(2)])
 228         if maxside > 0.0:
 229             scale = scaled_side_target / maxside
 230         else:
 231             scale = 1.0
 232         translate = [-0.5 * (maxv[i] + minv[i]) for i in range(2)]
 233         dim = len(self.points.pos[0])
 234         if dim == 3:
 235             translate.append([0.0])
 236         for v in range(len(self.points.pos)):
 237             self.points.pos[v] = tuple([scale * (self.points.pos[v][i] + \
 238                 translate[i]) for i in range(dim)])
 239
 240     def bounds(self):
 241         """Find bounding box of polyareas in xy.
 242
 243         Returns:
 244           ([minx,miny],[maxx,maxy]) - all floats
 245         """
 246
 247         huge = 1e100
 248         minv = [huge, huge]
 249         maxv = [-huge, -huge]
 250         for pa in self.polyareas:
 251             for face in [pa.poly] + pa.holes:
 252                 for v in face:
 253                     vcoords = self.points.pos[v]
 254                     for i in range(2):
 255                         if vcoords[i] < minv[i]:
 256                             minv[i] = vcoords[i]
 257                         if vcoords[i] > maxv[i]:
 258                             maxv[i] = vcoords[i]
 259         if minv[0] == huge:
 260             minv = [0.0, 0.0]
 261         if maxv[0] == huge:
 262             maxv = [0.0, 0.0]
 263         return (minv, maxv)
 264
 265
 266 class Model(object):
 267     """Contains a generic 3d model.
 268
 269     A generic 3d model has vertices with 3d coordinates.
 270     Each vertex gets a 'vertex id', which is an index that
 271     can be used to refer to the vertex and can be used
 272     to retrieve the 3d coordinates of the point.
 273
 274     The actual visible part of the geometry are the faces,
 275     which are n-gons (n>2), specified by a vector of the
 276     n corner vertices.
 277     Faces may also have data associated with them,
 278     and the data will be copied into newly created faces
 279     from the most likely neighbor faces..
 280
 281     Attributes:
 282       points: geom.Points - the 3d vertices
 283       faces: list of list of indices (each a CCW traversal of a face)
 284       face_data: list of any - if present, is parallel to
 285           faces list and holds arbitrary data
 286     """
 287
 288     def __init__(self):
 289         self.points = Points()
 290         self.faces = []
 291         self.face_data = []
 292
 293
 294 class Art(object):
 295     """Contains a vector art diagram.
 296
 297     Attributes:
 298       paths: list of Path objects
 299     """
 300
 301     def __init__(self):
 302         self.paths = []
 303
 304
 305 class Paint(object):
 306     """A color or pattern to fill or stroke with.
 307
 308     For now, just do colors, but could later do
 309     patterns or images too.
 310
 311     Attributes:
 312       color: (r,g,b) triple of floats, 0.0=no color, 1.0=max color
 313     """
 314
 315     def __init__(self, r=0.0, g=0.0, b=0.0):
 316         self.color = (r, g, b)
 317
 318     @staticmethod
 319     def CMYK(c, m, y, k):
 320         """Return Paint specified in CMYK model.
 321
 322         Uses formula from 6.2.4 of PDF Reference.
 323
 324         Args:
 325           c, m, y, k: float - in range [0, 1]
 326         Returns:
 327           Paint - with components in rgb form now
 328         """
 329
 330         return Paint(1.0 - min(1.0, c + k),
 331             1.0 - min(1.0, m + k), 1.0 - min(1.0, y + k))
 332
 333 black_paint = Paint()
 334 white_paint = Paint(1.0, 1.0, 1.0)
 335
 336 ColorDict = {
 337     'aqua': Paint(0.0, 1.0, 1.0),
 338     'black': Paint(0.0, 0.0, 0.0),
 339     'blue': Paint(0.0, 0.0, 1.0),
 340     'fuchsia': Paint(1.0, 0.0, 1.0),
 341     'gray': Paint(0.5, 0.5, 0.5),
 342     'green': Paint(0.0, 0.5, 0.0),
 343     'lime': Paint(0.0, 1.0, 0.0),
 344     'maroon': Paint(0.5, 0.0, 0.0),
 345     'navy': Paint(0.0, 0.0, 0.5),
 346     'olive': Paint(0.5, 0.5, 0.0),
 347     'purple': Paint(0.5, 0.0, 0.5),
 348     'red': Paint(1.0, 0.0, 0.0),
 349     'silver': Paint(0.75, 0.75, 0.75),
 350     'teal': Paint(0.0, 0.5, 0.5),
 351     'white': Paint(1.0, 1.0, 1.0),
 352     'yellow': Paint(1.0, 1.0, 0.0)
 353 }
 354
 355
 356 class Path(object):
 357     """Represents a path in the PDF sense, with painting instructions.
 358
 359     Attributes:
 360       subpaths: list of Subpath objects
 361       filled: True if path is to be filled
 362       fillevenodd: True if use even-odd rule to fill (else non-zero winding)
 363       stroked: True if path is to be stroked
 364       fillpaint: Paint to fill with
 365       strokepaint: Paint to stroke with
 366     """
 367
 368     def __init__(self):
 369         self.subpaths = []
 370         self.filled = False
 371         self.fillevenodd = False
 372         self.stroked = False
 373         self.fillpaint = black_paint
 374         self.strokepaint = black_paint
 375
 376     def AddSubpath(self, subpath):
 377         """"Add a subpath."""
 378
 379         self.subpaths.append(subpath)
 380
 381     def Empty(self):
 382         """Returns True if this Path as no subpaths."""
 383
 384         return not self.subpaths
 385
 386
 387 class Subpath(object):
 388     """Represents a subpath in PDF sense, either open or closed.
 389
 390     We'll represent lines, bezier pieces, circular arc pieces
 391     as tuples with letters giving segment type in first position
 392     and coordinates (2-tuples of floats) in the other positions.
 393
 394     Segment types:
 395      ('L', a, b)       - line from a to b
 396      ('B', a, b, c, d) - cubic bezier from a to b, with control points c,d
 397      ('Q', a, b, c)    - quadratic bezier from a to b, with 1 control point c
 398      ('A', a, b, rad, xrot, large-arc, ccw) - elliptical arc from a to b,
 399        with rad=(rx, ry) as radii, xrot is x-axis rotation in degrees,
 400        large-arc is True if arc should be >= 180 degrees,
 401        ccw is True if start->end follows counter-clockwise direction
 402        (see SVG spec); note that after rad,
 403        the rest are floats or bools, not coordinate pairs
 404     Note that s[1] and s[2] are the start and end points for any segment s.
 405
 406     Attributes:
 407       segments: list of segment tuples (see above)
 408       closed: True if closed
 409     """
 410
 411     def __init__(self):
 412         self.segments = []
 413         self.closed = False
 414
 415     def Empty(self):
 416         """Returns True if this subpath as no segments."""
 417
 418         return not self.segments
 419
 420     def AddSegment(self, seg):
 421         """Add a segment."""
 422
 423         self.segments.append(seg)
 424
 425     @staticmethod
 426     def SegStart(s):
 427         """Return start point for segment.
 428
 429         Args:
 430           s: a segment tuple
 431         Returns:
 432           (float, float): the coordinates of the segment's start point
 433         """
 434
 435         return s[1]
 436
 437     @staticmethod
 438     def SegEnd(s):
 439         """Return end point for segment.
 440
 441         Args:
 442           s: a segment tuple
 443         Returns:
 444           (float, float): the coordinates of the segment's end point
 445         """
 446
 447         return s[2]
 448
 449
 450 class TransformMatrix(object):
 451     """Transformation matrix for 2d coordinates.
 452
 453     The transform matrix is:
 454       [ a b 0 ]
 455       [ c d 0 ]
 456       [ e f 1 ]
 457     and coordinate transformation is defined by:
 458       [x' y' 1] = [x y 1] x TransformMatrix
 459
 460     Attributes:
 461       a, b, c, d, e, f: floats
 462     """
 463
 464     def __init__(self, a=1.0, b=0.0, c=0.0, d=1.0, e=0.0, f=0.0):
 465         self.a = a
 466         self.b = b
 467         self.c = c
 468         self.d = d
 469         self.e = e
 470         self.f = f
 471
 472     def __str__(self):
 473         return str([self.a, self.b, self.c, self.d, self.e, self.f])
 474
 475     def Copy(self):
 476         """Return a copy of this matrix."""
 477
 478         return TransformMatrix(self.a, self.b, self.c, self.d, self.e, self.f)
 479
 480     def ComposeTransform(self, a, b, c, d, e, f):
 481         """Apply the transform given the the arguments on top of this one.
 482
 483         This is accomplished by returning t x sel
 484         where t is the transform matrix that would be formed from the args.
 485
 486         Arguments:
 487           a, b, c, d, e, f: float - defines a composing TransformMatrix
 488         """
 489
 490         newa = a * self.a + b * self.c
 491         newb = a * self.b + b * self.d
 492         newc = c * self.a + d * self.c
 493         newd = c * self.b + d * self.d
 494         newe = e * self.a + f * self.c + self.e
 495         newf = e * self.b + f * self.d + self.f
 496         self.a = newa
 497         self.b = newb
 498         self.c = newc
 499         self.d = newd
 500         self.e = newe
 501         self.f = newf
 502
 503     def Apply(self, pt):
 504         """Return the result of applying this transform to pt = (x,y).
 505
 506         Arguments:
 507           (x, y) : (float, float)
 508         Returns:
 509           (x', y'): 2-tuple of floats, the result of [x y 1] x self
 510         """
 511
 512         (x, y) = pt
 513         return (self.a * x + self.c * y + self.e, \
 514             self.b * x + self.d * y + self.f)
 515
 516
 517 def ApproxEqualPoints(p, q):
 518     """Return True if p and q are approximately the same points.
 519
 520     Args:
 521       p: n-tuple of float
 522       q: n-tuple of float
 523     Returns:
 524       bool - True if the 1-norm <= DISTTOL
 525     """
 526
 527     for i in range(len(p)):
 528         if abs(p[i] - q[i]) > DISTTOL:
 529             return False
 530         return True
 531
 532
 533 def PointInside(v, a, points):
 534     """Return 1, 0, or -1 as v is inside, on, or outside polygon.
 535
 536     Cf. Eric Haines ptinpoly in Graphics Gems IV.
 537
 538     Args:
 539       v : (float, float) or (float, float, float) - coordinates of a point
 540       a : list of vertex indices defining polygon (assumed CCW)
 541       points: Points - to get coordinates for polygon
 542     Returns:
 543       1, 0, -1: as v is inside, on, or outside polygon a
 544     """
 545
 546     (xv, yv) = (v[0], v[1])
 547     vlast = points.pos[a[-1]]
 548     (x0, y0) = (vlast[0], vlast[1])
 549     if x0 == xv and y0 == yv:
 550         return 0
 551     yflag0 = y0 > yv
 552     inside = False
 553     n = len(a)
 554     for i in range(0, n):
 555         vi = points.pos[a[i]]
 556         (x1, y1) = (vi[0], vi[1])
 557         if x1 == xv and y1 == yv:
 558             return 0
 559         yflag1 = y1 > yv
 560         if yflag0 != yflag1:
 561             xflag0 = x0 > xv
 562             xflag1 = x1 > xv
 563             if xflag0 == xflag1:
 564                 if xflag0:
 565                     inside = not inside
 566             else:
 567                 z = x1 - (y1 - yv) * (x0 - x1) / (y0 - y1)
 568                 if z >= xv:
 569                     inside = not inside
 570         x0 = x1
 571         y0 = y1
 572         yflag0 = yflag1
 573     if inside:
 574         return 1
 575     else:
 576         return -1
 577
 578
 579 def SignedArea(polygon, points):
 580     """Return the area of the polgon, positive if CCW, negative if CW.
 581
 582     Args:
 583       polygon: list of vertex indices
 584       points: Points
 585     Returns:
 586       float - area of polygon, positive if it was CCW, else negative
 587     """
 588
 589     a = 0.0
 590     n = len(polygon)
 591     for i in range(0, n):
 592         u = points.pos[polygon[i]]
 593         v = points.pos[polygon[(i + 1) % n]]
 594         a += u[0] * v[1] - u[1] * v[0]
 595     return 0.5 * a
 596
 597
 598 def VecAdd(a, b):
 599     """Return vector a-b.
 600
 601     Args:
 602       a: n-tuple of floats
 603       b: n-tuple of floats
 604     Returns:
 605       n-tuple of floats - pairwise addition a+b
 606     """
 607
 608     n = len(a)
 609     assert(n == len(b))
 610     return tuple([a[i] + b[i] for i in range(n)])
 611
 612
 613 def VecSub(a, b):
 614     """Return vector a-b.
 615
 616     Args:
 617       a: n-tuple of floats
 618       b: n-tuple of floats
 619     Returns:
 620       n-tuple of floats - pairwise subtraction a-b
 621     """
 622
 623     n = len(a)
 624     assert(n == len(b))
 625     return tuple([a[i] - b[i] for i in range(n)])
 626
 627
 628 def VecDot(a, b):
 629     """Return the dot product of two vectors.
 630
 631     Args:
 632       a: n-tuple of floats
 633       b: n-tuple of floats
 634     Returns:
 635       n-tuple of floats - dot product of a and b
 636     """
 637
 638     n = len(a)
 639     assert(n == len(b))
 640     sum = 0.0
 641     for i in range(n):
 642         sum += a[i] * b[i]
 643     return sum
 644
 645
 646 def VecLen(a):
 647     """Return the Euclidean length of the argument vector.
 648
 649     Args:
 650       a: n-tuple of floats
 651     Returns:
 652       float: the 2-norm of a
 653     """
 654
 655     s = 0.0
 656     for v in a:
 657         s += v * v
 658     return math.sqrt(s)
 659
 660
 661 def Newell(poly, points):
 662     """Use Newell method to find polygon normal.
 663
 664     Assume poly has length at least 3 and points are 3d.
 665
 666     Args:
 667       poly: list of int - indices into points.pos
 668       points: Points - assumed 3d
 669     Returns:
 670       (float, float, float) - the average normal
 671     """
 672
 673     sumx = 0.0
 674     sumy = 0.0
 675     sumz = 0.0
 676     n = len(poly)
 677     pos = points.pos
 678     for i, ai in enumerate(poly):
 679         bi = poly[(i + 1) % n]
 680         a = pos[ai]
 681         b = pos[bi]
 682         sumx += (a[1] - b[1]) * (a[2] + b[2])
 683         sumy += (a[2] - b[2]) * (a[0] + b[0])
 684         sumz += (a[0] - b[0]) * (a[1] + b[1])
 685     return Norm3(sumx, sumy, sumz)
 686
 687
 688 def Norm3(x, y, z):
 689     """Return vector (x,y,z) normalized by dividing by squared length.
 690     Return (0.0, 0.0, 1.0) if the result is undefined."""
 691     sqrlen = x * x + y * y + z * z
 692     if sqrlen < 1e-100:
 693         return (0.0, 0.0, 1.0)
 694     else:
 695         try:
 696             d = math.sqrt(sqrlen)
 697             return (x / d, y / d, z / d)
 698         except:
 699             return (0.0, 0.0, 1.0)
 700
 701
 702 # We're using right-hand coord system, where
 703 # forefinger=x, middle=y, thumb=z on right hand.
 704 # Then, e.g., (1,0,0) x (0,1,0) = (0,0,1)
 705 def Cross3(a, b):
 706     """Return the cross product of two vectors, a x b."""
 707
 708     (ax, ay, az) = a
 709     (bx, by, bz) = b
 710     return (ay * bz - az * by, az * bx - ax * bz, ax * by - ay * bx)
 711
 712
 713 def MulPoint3(p, m):
 714     """Return matrix multiplication of p times m
 715     where m is a 4x3 matrix and p is a 3d point, extended with 1."""
 716
 717     (x, y, z) = p
 718     return (x * m[0] + y * m[3] + z * m[6] + m[9],
 719         x * m[1] + y * m[4] + z * m[7] + m[10],
 720         x * m[2] + y * m[5] + z * m[8] + m[11])