diff options
Diffstat (limited to 'mesh_inset/geom.py')
| -rw-r--r-- | mesh_inset/geom.py | 719 |
1 files changed, 719 insertions, 0 deletions
diff --git a/mesh_inset/geom.py b/mesh_inset/geom.py new file mode 100644 index 00000000..98059dcf --- /dev/null +++ b/mesh_inset/geom.py @@ -0,0 +1,719 @@ +# ##### BEGIN GPL LICENSE BLOCK ##### +# +# This program is free software; you can redistribute it and/or +# modify it under the terms of the GNU General Public License +# as published by the Free Software Foundation; either version 2 +# of the License, or (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program; if not, write to the Free Software Foundation, +# Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +# +# ##### END GPL LICENSE BLOCK ##### + +# <pep8 compliant> + +"""Geometry classes and operations. +Also, vector file representation (Art). +""" + +__author__ = "howard.trickey@gmail.com" + +import math + +# distances less than about DISTTOL will be considered +# essentially zero +DISTTOL = 1e-3 +INVDISTTOL = 1e3 + + +class Points(object): + """Container of points without duplication, each mapped to an int. + + Points are either have dimension at least 2, maybe more. + + Implementation: + In order to efficiently find duplicates, we quantize the points + to triples of ints and map from quantized triples to vertex + index. + + Attributes: + pos: list of tuple of float - coordinates indexed by + vertex number + invmap: dict of (int, int, int) to int - quantized coordinates + to vertex number map + """ + + def __init__(self, initlist=[]): + self.pos = [] + self.invmap = dict() + for p in initlist: + self.AddPoint(p) + + @staticmethod + def Quantize(p): + """Quantize the float tuple into an int tuple. + + Args: + p: tuple of float + Returns: + tuple of int - scaled by INVDISTTOL and rounded p + """ + + return tuple([int(round(v * INVDISTTOL)) for v in p]) + + def AddPoint(self, p): + """Add point p to the Points set and return vertex number. + + If there is an existing point which quantizes the same,, + don't add a new one but instead return existing index. + + Args: + p: tuple of float - coordinates (2-tuple or 3-tuple) + Returns: + int - the vertex number of added (or existing) point + """ + + qp = Points.Quantize(p) + if qp in self.invmap: + return self.invmap[qp] + else: + self.invmap[qp] = len(self.pos) + self.pos.append(p) + return len(self.pos) - 1 + + def AddPoints(self, points): + """Add another set of points to this set. + + We need to return a mapping from indices + in the argument points space into indices + in this point space. + + Args: + points: Points - to union into this set + Returns: + list of int: maps added indices to new ones + """ + + vmap = [0] * len(points.pos) + for i in range(len(points.pos)): + vmap[i] = self.AddPoint(points.pos[i]) + return vmap + + def AddZCoord(self, z): + """Change this in place to have a z coordinate, with value z. + + Assumes the coordinates are currently 2d. + + Args: + z: the value of the z coordinate to add + Side Effect: + self now has a z-coordinate added + """ + + assert(len(self.pos) == 0 or len(self.pos[0]) == 2) + newinvmap = dict() + for i, (x, y) in enumerate(self.pos): + newp = (x, y, z) + self.pos[i] = newp + newinvmap[self.Quantize(newp)] = i + self.invmap = newinvmap + + def AddToZCoord(self, i, delta): + """Change the z-coordinate of point with index i to add delta. + + Assumes the coordinates are currently 3d. + + Args: + i: int - index of a point + delta: float - value to add to z-coord + """ + + (x, y, z) = self.pos[i] + self.pos[i] = (x, y, z + delta) + + +class PolyArea(object): + """Contains a Polygonal Area (polygon with possible holes). + + A polygon is a list of vertex ids, each an index given by + a Points object. The list represents a CCW-oriented + outer boundary (implicitly closed). + If there are holes, they are lists of CW-oriented vertices + that should be contained in the outer boundary. + (So the left face of both the poly and the holes is + the filled part.) + + Attributes: + points: Points + poly: list of vertex ids + holes: list of lists of vertex ids (each a hole in poly) + data: any - application data (can hold color, e.g.) + """ + + def __init__(self, points=None, poly=None, holes=None, data=None): + self.points = points if points else Points() + self.poly = poly if poly else [] + self.holes = holes if holes else [] + self.data = data + + def AddHole(self, holepa): + """Add a PolyArea's poly as a hole of self. + + Need to reverse the contour and + adjust the the point indexes and self.points. + + Args: + holepa: PolyArea + """ + + vmap = self.points.AddPoints(holepa.points) + holepoly = [vmap[i] for i in holepa.poly] + holepoly.reverse() + self.holes.append(holepoly) + + def ContainsPoly(self, poly, points): + """Tests if poly is contained within self.poly. + + Args: + poly: list of int - indices into points + points: Points - maps to coords + Returns: + bool - True if poly is fully contained within self.poly + """ + + for v in poly: + if PointInside(points.pos[v], self.poly, self.points) == -1: + return False + return True + + def Normal(self): + """Returns the normal of the polyarea's main poly.""" + + pos = self.points.pos + poly = self.poly + if len(pos) == 0 or len(pos[0]) == 2 or len(poly) == 0: + print("whoops, not enough info to calculate normal") + return (0.0, 0.0, 1.0) + return Newell(poly, self.points) + + +class PolyAreas(object): + """Contains a list of PolyAreas and a shared Points. + + Attributes: + polyareas: list of PolyArea + points: Points + """ + + def __init__(self): + self.polyareas = [] + self.points = Points() + + def scale_and_center(self, scaled_side_target): + """Adjust the coordinates of the polyareas so that + it is centered at the origin and has its longest + dimension scaled to be scaled_side_target.""" + + if len(self.points.pos) == 0: + return + (minv, maxv) = self.bounds() + maxside = max([maxv[i] - minv[i] for i in range(2)]) + if maxside > 0.0: + scale = scaled_side_target / maxside + else: + scale = 1.0 + translate = [-0.5 * (maxv[i] + minv[i]) for i in range(2)] + dim = len(self.points.pos[0]) + if dim == 3: + translate.append([0.0]) + for v in range(len(self.points.pos)): + self.points.pos[v] = tuple([scale * (self.points.pos[v][i] + \ + translate[i]) for i in range(dim)]) + + def bounds(self): + """Find bounding box of polyareas in xy. + + Returns: + ([minx,miny],[maxx,maxy]) - all floats + """ + + huge = 1e100 + minv = [huge, huge] + maxv = [-huge, -huge] + for pa in self.polyareas: + for face in [pa.poly] + pa.holes: + for v in face: + vcoords = self.points.pos[v] + for i in range(2): + if vcoords[i] < minv[i]: + minv[i] = vcoords[i] + if vcoords[i] > maxv[i]: + maxv[i] = vcoords[i] + if minv[0] == huge: + minv = [0.0, 0.0] + if maxv[0] == huge: + maxv = [0.0, 0.0] + return (minv, maxv) + + +class Model(object): + """Contains a generic 3d model. + + A generic 3d model has vertices with 3d coordinates. + Each vertex gets a 'vertex id', which is an index that + can be used to refer to the vertex and can be used + to retrieve the 3d coordinates of the point. + + The actual visible part of the geometry are the faces, + which are n-gons (n>2), specified by a vector of the + n corner vertices. + Faces may also have data associated with them, + and the data will be copied into newly created faces + from the most likely neighbor faces.. + + Attributes: + points: geom.Points - the 3d vertices + faces: list of list of indices (each a CCW traversal of a face) + face_data: list of any - if present, is parallel to + faces list and holds arbitrary data + """ + + def __init__(self): + self.points = Points() + self.faces = [] + self.face_data = [] + + +class Art(object): + """Contains a vector art diagram. + + Attributes: + paths: list of Path objects + """ + + def __init__(self): + self.paths = [] + + +class Paint(object): + """A color or pattern to fill or stroke with. + + For now, just do colors, but could later do + patterns or images too. + + Attributes: + color: (r,g,b) triple of floats, 0.0=no color, 1.0=max color + """ + + def __init__(self, r=0.0, g=0.0, b=0.0): + self.color = (r, g, b) + + @staticmethod + def CMYK(c, m, y, k): + """Return Paint specified in CMYK model. + + Uses formula from 6.2.4 of PDF Reference. + + Args: + c, m, y, k: float - in range [0, 1] + Returns: + Paint - with components in rgb form now + """ + + return Paint(1.0 - min(1.0, c + k), + 1.0 - min(1.0, m + k), 1.0 - min(1.0, y + k)) + +black_paint = Paint() +white_paint = Paint(1.0, 1.0, 1.0) + +ColorDict = { + 'aqua': Paint(0.0, 1.0, 1.0), + 'black': Paint(0.0, 0.0, 0.0), + 'blue': Paint(0.0, 0.0, 1.0), + 'fuchsia': Paint(1.0, 0.0, 1.0), + 'gray': Paint(0.5, 0.5, 0.5), + 'green': Paint(0.0, 0.5, 0.0), + 'lime': Paint(0.0, 1.0, 0.0), + 'maroon': Paint(0.5, 0.0, 0.0), + 'navy': Paint(0.0, 0.0, 0.5), + 'olive': Paint(0.5, 0.5, 0.0), + 'purple': Paint(0.5, 0.0, 0.5), + 'red': Paint(1.0, 0.0, 0.0), + 'silver': Paint(0.75, 0.75, 0.75), + 'teal': Paint(0.0, 0.5, 0.5), + 'white': Paint(1.0, 1.0, 1.0), + 'yellow': Paint(1.0, 1.0, 0.0) +} + + +class Path(object): + """Represents a path in the PDF sense, with painting instructions. + + Attributes: + subpaths: list of Subpath objects + filled: True if path is to be filled + fillevenodd: True if use even-odd rule to fill (else non-zero winding) + stroked: True if path is to be stroked + fillpaint: Paint to fill with + strokepaint: Paint to stroke with + """ + + def __init__(self): + self.subpaths = [] + self.filled = False + self.fillevenodd = False + self.stroked = False + self.fillpaint = black_paint + self.strokepaint = black_paint + + def AddSubpath(self, subpath): + """"Add a subpath.""" + + self.subpaths.append(subpath) + + def Empty(self): + """Returns True if this Path as no subpaths.""" + + return not self.subpaths + + +class Subpath(object): + """Represents a subpath in PDF sense, either open or closed. + + We'll represent lines, bezier pieces, circular arc pieces + as tuples with letters giving segment type in first position + and coordinates (2-tuples of floats) in the other positions. + + Segment types: + ('L', a, b) - line from a to b + ('B', a, b, c, d) - cubic bezier from a to b, with control points c,d + ('Q', a, b, c) - quadratic bezier from a to b, with 1 control point c + ('A', a, b, rad, xrot, large-arc, ccw) - elliptical arc from a to b, + with rad=(rx, ry) as radii, xrot is x-axis rotation in degrees, + large-arc is True if arc should be >= 180 degrees, + ccw is True if start->end follows counter-clockwise direction + (see SVG spec); note that after rad, + the rest are floats or bools, not coordinate pairs + Note that s[1] and s[2] are the start and end points for any segment s. + + Attributes: + segments: list of segment tuples (see above) + closed: True if closed + """ + + def __init__(self): + self.segments = [] + self.closed = False + + def Empty(self): + """Returns True if this subpath as no segments.""" + + return not self.segments + + def AddSegment(self, seg): + """Add a segment.""" + + self.segments.append(seg) + + @staticmethod + def SegStart(s): + """Return start point for segment. + + Args: + s: a segment tuple + Returns: + (float, float): the coordinates of the segment's start point + """ + + return s[1] + + @staticmethod + def SegEnd(s): + """Return end point for segment. + + Args: + s: a segment tuple + Returns: + (float, float): the coordinates of the segment's end point + """ + + return s[2] + + +class TransformMatrix(object): + """Transformation matrix for 2d coordinates. + + The transform matrix is: + [ a b 0 ] + [ c d 0 ] + [ e f 1 ] + and coordinate tranformation is defined by: + [x' y' 1] = [x y 1] x TransformMatrix + + Attributes: + a, b, c, d, e, f: floats + """ + + def __init__(self, a=1.0, b=0.0, c=0.0, d=1.0, e=0.0, f=0.0): + self.a = a + self.b = b + self.c = c + self.d = d + self.e = e + self.f = f + + def __str__(self): + return str([self.a, self.b, self.c, self.d, self.e, self.f]) + + def Copy(self): + """Return a copy of this matrix.""" + + return TransformMatrix(self.a, self.b, self.c, self.d, self.e, self.f) + + def ComposeTransform(self, a, b, c, d, e, f): + """Apply the transform given the the arguments on top of this one. + + This is accomplished by returning t x sel + where t is the transform matrix that would be formed from the args. + + Arguments: + a, b, c, d, e, f: float - defines a composing TransformMatrix + """ + + newa = a * self.a + b * self.c + newb = a * self.b + b * self.d + newc = c * self.a + d * self.c + newd = c * self.b + d * self.d + newe = e * self.a + f * self.c + self.e + newf = e * self.b + f * self.d + self.f + self.a = newa + self.b = newb + self.c = newc + self.d = newd + self.e = newe + self.f = newf + + def Apply(self, pt): + """Return the result of applying this tranform to pt = (x,y). + + Arguments: + (x, y) : (float, float) + Returns: + (x', y'): 2-tuple of floats, the result of [x y 1] x self + """ + + (x, y) = pt + return (self.a * x + self.c * y + self.e, \ + self.b * x + self.d * y + self.f) + + +def ApproxEqualPoints(p, q): + """Return True if p and q are approximately the same points. + + Args: + p: n-tuple of float + q: n-tuple of float + Returns: + bool - True if the 1-norm <= DISTTOL + """ + + for i in range(len(p)): + if abs(p[i] - q[i]) > DISTTOL: + return False + return True + + +def PointInside(v, a, points): + """Return 1, 0, or -1 as v is inside, on, or outside polygon. + + Cf. Eric Haines ptinpoly in Graphics Gems IV. + + Args: + v : (float, float) or (float, float, float) - coordinates of a point + a : list of vertex indices defining polygon (assumed CCW) + points: Points - to get coordinates for polygon + Returns: + 1, 0, -1: as v is inside, on, or outside polygon a + """ + + (xv, yv) = (v[0], v[1]) + vlast = points.pos[a[-1]] + (x0, y0) = (vlast[0], vlast[1]) + if x0 == xv and y0 == yv: + return 0 + yflag0 = y0 > yv + inside = False + n = len(a) + for i in range(0, n): + vi = points.pos[a[i]] + (x1, y1) = (vi[0], vi[1]) + if x1 == xv and y1 == yv: + return 0 + yflag1 = y1 > yv + if yflag0 != yflag1: + xflag0 = x0 > xv + xflag1 = x1 > xv + if xflag0 == xflag1: + if xflag0: + inside = not inside + else: + z = x1 - (y1 - yv) * (x0 - x1) / (y0 - y1) + if z >= xv: + inside = not inside + x0 = x1 + y0 = y1 + yflag0 = yflag1 + if inside: + return 1 + else: + return -1 + + +def SignedArea(polygon, points): + """Return the area of the polgon, positive if CCW, negative if CW. + + Args: + polygon: list of vertex indices + points: Points + Returns: + float - area of polygon, positive if it was CCW, else negative + """ + + a = 0.0 + n = len(polygon) + for i in range(0, n): + u = points.pos[polygon[i]] + v = points.pos[polygon[(i + 1) % n]] + a += u[0] * v[1] - u[1] * v[0] + return 0.5 * a + + +def VecAdd(a, b): + """Return vector a-b. + + Args: + a: n-tuple of floats + b: n-tuple of floats + Returns: + n-tuple of floats - pairwise addition a+b + """ + + n = len(a) + assert(n == len(b)) + return tuple([a[i] + b[i] for i in range(n)]) + + +def VecSub(a, b): + """Return vector a-b. + + Args: + a: n-tuple of floats + b: n-tuple of floats + Returns: + n-tuple of floats - pairwise subtraction a-b + """ + + n = len(a) + assert(n == len(b)) + return tuple([a[i] - b[i] for i in range(n)]) + + +def VecDot(a, b): + """Return the dot product of two vectors. + + Args: + a: n-tuple of floats + b: n-tuple of floats + Returns: + n-tuple of floats - dot product of a and b + """ + + n = len(a) + assert(n == len(b)) + sum = 0.0 + for i in range(n): + sum += a[i] * b[i] + return sum + + +def VecLen(a): + """Return the Euclidean lenght of the argument vector. + + Args: + a: n-tuple of floats + Returns: + float: the 2-norm of a + """ + + s = 0.0 + for v in a: + s += v * v + return math.sqrt(s) + + +def Newell(poly, points): + """Use Newell method to find polygon normal. + + Assume poly has length at least 3 and points are 3d. + + Args: + poly: list of int - indices into points.pos + points: Points - assumed 3d + Returns: + (float, float, float) - the average normal + """ + + sumx = 0.0 + sumy = 0.0 + sumz = 0.0 + n = len(poly) + pos = points.pos + for i, ai in enumerate(poly): + bi = poly[(i + 1) % n] + a = pos[ai] + b = pos[bi] + sumx += (a[1] - b[1]) * (a[2] + b[2]) + sumy += (a[2] - b[2]) * (a[0] + b[0]) + sumz += (a[0] - b[0]) * (a[1] + b[1]) + return Norm3(sumx, sumy, sumz) + + +def Norm3(x, y, z): + """Return vector (x,y,z) normalized by dividing by squared length. + Return (0.0, 0.0, 1.0) if the result is undefined.""" + sqrlen = x * x + y * y + z * z + if sqrlen < 1e-100: + return (0.0, 0.0, 1.0) + else: + try: + d = math.sqrt(sqrlen) + return (x / d, y / d, z / d) + except: + return (0.0, 0.0, 1.0) + + +# We're using right-hand coord system, where +# forefinger=x, middle=y, thumb=z on right hand. +# Then, e.g., (1,0,0) x (0,1,0) = (0,0,1) +def Cross3(a, b): + """Return the cross product of two vectors, a x b.""" + + (ax, ay, az) = a + (bx, by, bz) = b + return (ay * bz - az * by, az * bx - ax * bz, ax * by - ay * bx) + + +def MulPoint3(p, m): + """Return matrix multiplication of p times m + where m is a 4x3 matrix and p is a 3d point, extended with 1.""" + + (x, y, z) = p + return (x * m[0] + y * m[3] + z * m[6] + m[9], + x * m[1] + y * m[4] + z * m[7] + m[10], + x * m[2] + y * m[5] + z * m[8] + m[11]) |