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. 22Also, vector file representation (Art). 23""" 24 25__author__ = "howard.trickey@gmail.com" 26 27import math 28 29# distances less than about DISTTOL will be considered 30# essentially zero 31DISTTOL = 1e-3 32INVDISTTOL = 1e3 33 34 35class 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 142class 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 207class 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 266class 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 294class 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 305class 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 333black_paint = Paint() 334white_paint = Paint(1.0, 1.0, 1.0) 335 336ColorDict = { 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 356class 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 387class 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 450class 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 517def 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 533def 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 579def 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 598def 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 613def 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 628def 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 646def 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 661def 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 688def 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) 705def 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 713def 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]) 721