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