diff options
Diffstat (limited to 'mesh_inset/offset.py')
| -rw-r--r-- | mesh_inset/offset.py | 755 |
1 files changed, 755 insertions, 0 deletions
diff --git a/mesh_inset/offset.py b/mesh_inset/offset.py new file mode 100644 index 00000000..3a44b95b --- /dev/null +++ b/mesh_inset/offset.py @@ -0,0 +1,755 @@ +# ##### 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> + +"""Creating offset polygons inside faces.""" + +__author__ = "howard.trickey@gmail.com" + +import math +from . import triquad +from . import geom +from .triquad import Sub2, Add2, Angle, Ccw, Normalized2, Perp2, Length2, \ + LinInterp2, TOL +from .geom import Points + +AREATOL = 1e-4 + + +class Spoke(object): + """A Spoke is a line growing from an outer vertex to an inner one. + + A Spoke is contained in an Offset (see below). + + Attributes: + origin: int - index of origin point in a Points + dest: int - index of dest point + is_reflex: bool - True if spoke grows from a reflex angle + dir: (float, float, float) - direction vector (normalized) + speed: float - at time t, other end of spoke is + origin + t*dir. Speed is such that the wavefront + from the face edges moves at speed 1. + face: int - index of face containing this Spoke, in Offset + index: int - index of this Spoke in its face + destindex: int - index of Spoke dest in its face + """ + + def __init__(self, v, prev, next, face, index, points): + """Set attribute of spoke from points making up initial angle. + + The spoke grows from an angle inside a face along the bisector + of that angle. Its speed is 1/sin(.5a), where a is the angle + formed by (prev, v, next). That speed means that the perpendicular + from the end of the spoke to either of the prev->v or v->prev + edges will grow at speed 1. + + Args: + v: int - index of point spoke grows from + prev: int - index of point before v on boundary (in CCW order) + next: int - index of point after v on boundary (in CCW order) + face: int - index of face containing this spoke, in containing offset + index: int - index of this spoke in its face + points: geom.Points - maps vertex indices to 3d coords + """ + + self.origin = v + self.dest = v + self.face = face + self.index = index + self.destindex = -1 + vmap = points.pos + vp = vmap[v] + prevp = vmap[prev] + nextp = vmap[next] + uin = Normalized2(Sub2(vp, prevp)) + uout = Normalized2(Sub2(nextp, vp)) + uavg = Normalized2((0.5 * (uin[0] + uout[0]), \ + 0.5 * (uin[1] + uout[1]))) + if abs(Length2(uavg)) < TOL: + # in and out vectors are reverse of each other + self.dir = (uout[0], uout[1], 0.0) + self.is_reflex = False + self.speed = 1e7 + else: + # bisector direction is 90 degree CCW rotation of + # average incoming/outgoing + self.dir = (-uavg[1], uavg[0], 0.0) + self.is_reflex = Ccw(next, v, prev, points) + ang = Angle(prev, v, next, points) # in range [0, 180) + sin_half_ang = math.sin(math.pi * ang / 360.0) + if abs(sin_half_ang) < TOL: + self.speed = 1e7 + else: + self.speed = 1.0 / sin_half_ang + + def __repr__(self): + """Printing representation of a Spoke.""" + + return "@%d+%gt%s <%d,%d>" % (self.origin, \ + self.speed, str(self.dir), \ + self.face, self.index) + + def EndPoint(self, t, points, vspeed): + """Return the coordinates of the non-origin point at time t. + + Args: + t: float - time to end of spoke + points: geom.Points - coordinate map + vspeed: float - speed in z direction + Returns: + (float, float, float) - coords of spoke's endpoint at time t + """ + + p = points.pos[self.origin] + d = self.dir + v = self.speed + return (p[0] + v * t * d[0], p[1] + v * t * d[1], p[2] + vspeed * t) + + def VertexEvent(self, other, points): + """Intersect self with other spoke, and return the OffsetEvent, if any. + + A vertex event is with one advancing spoke intersects an adjacent + adavancing spoke, forming a new vertex. + + Args: + other: Spoke - other spoke to intersect with + points: Geom.points + Returns: + None or OffsetEvent - if there's an intersection in the growing + directions of the spokes, will return the OffsetEvent for + the intersection; + if lines are collinear or parallel, return None + """ + + vmap = points.pos + a = vmap[self.origin] + b = Add2(a, self.dir) + c = vmap[other.origin] + d = Add2(c, other.dir) + # find intersection of line ab with line cd + u = Sub2(b, a) + v = Sub2(d, c) + w = Sub2(a, c) + pp = Perp2(u, v) + if abs(pp) > TOL: + # lines or neither parallel nor collinear + si = Perp2(v, w) / pp + ti = Perp2(u, w) / pp + if si >= 0 and ti >= 0: + p = LinInterp2(a, b, si) + dist_ab = si * Length2(u) + dist_cd = ti * Length2(v) + time_ab = dist_ab / self.speed + time_cd = dist_cd / other.speed + time = max(time_ab, time_cd) + return OffsetEvent(True, time, p, self, other) + return None + + def EdgeEvent(self, other, offset): + """Intersect self with advancing edge and return OffsetEvent, if any. + + An edge event is when one advancing spoke intersects an advancing + edge. Advancing edges start out as face edges and move perpendicular + to them, at a rate of 1. The endpoints of the edge are the advancing + spokes on either end of the edge (so the edge shrinks or grows as + it advances). At some time, the edge may shrink to nothing and there + will be no EdgeEvent after that time. + + We represent an advancing edge by the first spoke (in CCW order + of face) of the pair of defining spokes. + + At time t, end of this spoke is at + o + d*s*t + where o=self.origin, d=self.dir, s= self.speed. + The advancing edge line has this equation: + oo + od*os*t + p*a + where oo, od, os are o, d, s for other spoke, and p is direction + vector parallel to advancing edge, and a is a real parameter. + Equating x and y of intersection point: + + o.x + d.x*s*t = oo.x + od.x*os*t + p.x*w + o.y + d.y*s*t = oo.y + od.y*os*t + p.y*w + + which can be rearranged into the form + + a = bt + cw + d = et + fw + + and solved for t, w. + + Args: + other: Spoke - the edge out of this spoke's origin is the advancing + edge to be checked for intersection + offset: Offset - the containing Offset + Returns: + None or OffsetEvent - with data about the intersection, if any + """ + + vmap = offset.polyarea.points.pos + o = vmap[self.origin] + oo = vmap[other.origin] + otherface = offset.facespokes[other.face] + othernext = otherface[(other.index + 1) % len(otherface)] + oonext = vmap[othernext.origin] + p = Normalized2(Sub2(oonext, oo)) + a = o[0] - oo[0] + d = o[1] - oo[1] + b = other.dir[0] * other.speed - self.dir[0] * self.speed + e = other.dir[1] * other.speed - self.dir[1] * self.speed + c = p[0] + f = p[1] + if abs(c) > TOL: + dem = e - f * b / c + if abs(dem) > TOL: + t = (d - f * a / c) / dem + w = (a - b * t) / c + else: + return None + elif abs(f) > TOL: + dem = b - c * e / f + if abs(dem) > TOL: + t = (a - c * d / f) / dem + w = (d - e * t) / f + else: + return None + else: + return None + if t < 0.0: + # intersection is in backward direction along self spoke + return None + if w < 0.0: + # intersection on wrong side of first end of advancing line segment + return None + # calculate the equivalent of w for the other end + aa = o[0] - oonext[0] + dd = o[1] - oonext[1] + bb = othernext.dir[0] * othernext.speed - self.dir[0] * self.speed + ee = othernext.dir[1] * othernext.speed - self.dir[1] * self.speed + cc = -p[0] + ff = -p[1] + if abs(cc) > TOL: + ww = (aa - bb * t) / cc + elif abs(ff) > TOL: + ww = (dd - ee * t) / ff + else: + return None + if ww < 0.0: + return None + evertex = (o[0] + self.dir[0] * self.speed * t, \ + o[1] + self.dir[1] * self.speed * t) + return OffsetEvent(False, t, evertex, self, other) + + +class OffsetEvent(object): + """An event involving a spoke during offset computation. + + The events kinds are: + vertex event: the spoke intersects an adjacent spoke and makes a new + vertex + edge event: the spoke hits an advancing edge and splits it + + Attributes: + is_vertex_event: True if this is a vertex event (else it is edge event) + time: float - time at which it happens (edges advance at speed 1) + event_vertex: (float, float) - intersection point of event + spoke: Spoke - the spoke that this event is for + other: Spoke - other spoke involved in event; if vertex event, this will + be an adjacent spoke that intersects; if an edge event, this is the + spoke whose origin's outgoing edge grows to hit this event's spoke + """ + + def __init__(self, isv, time, evertex, spoke, other): + """Creates and initializes attributes of an OffsetEvent.""" + + self.is_vertex_event = isv + self.time = time + self.event_vertex = evertex + self.spoke = spoke + self.other = other + + def __repr__(self): + """Printing representation of an event.""" + + if self.is_vertex_event: + c = "V" + else: + c = "E" + return "%s t=%5f %s %s %s" % (c, self.time, str(self.event_vertex), \ + repr(self.spoke), repr(self.other)) + + +class Offset(object): + """Represents an offset polygonal area, and used to construct one. + + Currently, the polygonal area must lie approximately in the XY plane. + As well as growing inwards in that plane, the advancing lines also + move in the Z direction at the rate of vspeed. + + Attributes: + polyarea: geom.PolyArea - the area we are offsetting from. + We share the polyarea.points, and add to it as points in + the offset polygonal area are computed. + facespokes: list of list of Spoke - each sublist is a closed face + (oriented CCW); the faces may mutually interfere. + These lists are spokes for polyarea.poly + polyarea.holes. + endtime: float - time when this offset hits its first + event (relative to beginning of this offset), or the amount + that takes this offset to the end of the total Build time + timesofar: float - sum of times taken by all containing Offsets + vspeed: float - speed that edges move perpendicular to offset plane + inneroffsets: list of Offset - the offsets that take over after this + (inside it) + """ + + def __init__(self, polyarea, time, vspeed): + """Set up initial state of Offset from a polyarea. + + Args: + polyarea: geom.PolyArea + time: float - time so far + """ + + self.polyarea = polyarea + self.facespokes = [] + self.endtime = 0.0 + self.timesofar = time + self.vspeed = vspeed + self.inneroffsets = [] + self.InitFaceSpokes(polyarea.poly) + for f in polyarea.holes: + self.InitFaceSpokes(f) + + def __repr__(self): + ans = ["Offset: endtime=%g" % self.endtime] + for i, face in enumerate(self.facespokes): + ans.append(("<%d>" % i) + str([str(spoke) for spoke in face])) + return '\n'.join(ans) + + def PrintNest(self, indent_level=0): + indent = " " * indent_level * 4 + print(indent + "Offset timesofar=", self.timesofar, "endtime=", + self.endtime) + print(indent + " polyarea=", self.polyarea.poly, self.polyarea.holes) + for o in self.inneroffsets: + o.PrintNest(indent_level + 1) + + def InitFaceSpokes(self, face_vertices): + """Initialize the offset representation of a face from vertex list. + + If the face has no area or too small an area, don't bother making it. + + Args: + face_vertices: list of int - point indices for boundary of face + Side effect: + A new face (list of spokes) may be added to self.facespokes + """ + + n = len(face_vertices) + if n <= 2: + return + points = self.polyarea.points + area = abs(geom.SignedArea(face_vertices, points)) + if area < AREATOL: + return + findex = len(self.facespokes) + fspokes = [Spoke(v, face_vertices[(i - 1) % n], \ + face_vertices[(i + 1) % n], findex, i, points) \ + for i, v in enumerate(face_vertices)] + self.facespokes.append(fspokes) + + def NextSpokeEvents(self, spoke): + """Return the OffsetEvents that will next happen for a given spoke. + + It might happen that some events happen essentially simultaneously, + and also it is convenient to separate Edge and Vertex events, so + we return two lists. + But, for vertex events, only look at the event with the next Spoke, + as the event with the previous spoke will be accounted for when we + consider that previous spoke. + + Args: + spoke: Spoke - a spoke in one of the faces of this object + Returns: + (float, list of OffsetEvent, list of OffsetEvent) - + time of next event, + next Vertex event list and next Edge event list + """ + + facespokes = self.facespokes[spoke.face] + n = len(facespokes) + bestt = 1e100 + bestv = [] + beste = [] + # First find vertex event (only the one with next spoke) + next_spoke = facespokes[(spoke.index + 1) % n] + ev = spoke.VertexEvent(next_spoke, self.polyarea.points) + if ev: + bestv = [ev] + bestt = ev.time + # Now find edge events, if this is a reflex vertex + if spoke.is_reflex: + prev_spoke = facespokes[(spoke.index - 1) % n] + for f in self.facespokes: + for other in f: + if other == spoke or other == prev_spoke: + continue + ev = spoke.EdgeEvent(other, self) + if ev: + if ev.time < bestt - TOL: + beste = [] + bestv = [] + bestt = ev.time + if abs(ev.time - bestt) < TOL: + beste.append(ev) + return (bestt, bestv, beste) + + def Build(self, target=2e100): + """Build the complete Offset structure or up until target time. + + Find the next event(s), makes the appropriate inner Offsets + that are inside this one, and calls Build on those Offsets to continue + the process until only a single point is left or time reaches target. + """ + + bestt = 1e100 + bestevs = [[], []] + for f in self.facespokes: + for s in f: + (t, ve, ee) = self.NextSpokeEvents(s) + if t < bestt - TOL: + bestevs = [[], []] + bestt = t + if abs(t - bestt) < TOL: + bestevs[0].extend(ve) + bestevs[1].extend(ee) + if bestt == 1e100: + # could happen if polygon is oriented wrong + # or in other special cases + return + if abs(bestt) < TOL: + # seems to be in a loop, so quit + return + self.endtime = bestt + (ve, ee) = bestevs + newfaces = [] + splitjoin = None + if target < self.endtime: + self.endtime = target + newfaces = self.MakeNewFaces(self.endtime) + elif ve and not ee: + # Only vertex events. + # Merging of successive vertices in inset face will + # take care of the vertex events + newfaces = self.MakeNewFaces(self.endtime) + else: + # Edge events too + # First make the new faces (handles all vertex events) + newfaces = self.MakeNewFaces(self.endtime) + # Only do one edge event (handle other simultaneous edge + # events in subsequent recursive Build calls) + splitjoin = self.SplitJoinFaces(newfaces, ee[0]) + nexttarget = target - self.endtime + if len(newfaces) > 0: + pa = geom.PolyArea(points=self.polyarea.points) + pa.data = self.polyarea.data + newt = self.timesofar + self.endtime + pa2 = None # may make another + if not splitjoin: + pa.poly = newfaces[0] + pa.holes = newfaces[1:] + elif splitjoin[0] == 'split': + (_, findex, newface0, newface1) = splitjoin + if findex == 0: + # Outer poly of polyarea was split. + # Now there will be two polyareas. + # If there were holes, need to allocate according to + # which one contains the holes. + pa.poly = newface0 + pa2 = geom.PolyArea(points=self.polyarea.points) + pa2.data = self.polyarea.data + pa2.poly = newface1 + if len(newfaces) > 1: + # print("need to allocate holes") + for hf in newfaces[1:]: + if pa.ContainsPoly(hf, self.polyarea.points): + # print("add", hf, "to", pa.poly) + pa.holes.append(hf) + elif pa2.ContainsPoly(hf, self.polyarea.points): + # print("add", hf, "to", pa2.poly) + pa2.holes.append(hf) + else: + print("whoops, hole in neither poly!") + self.inneroffsets = [Offset(pa, newt, self.vspeed), \ + Offset(pa2, newt, self.vspeed)] + else: + # A hole was split. New faces just replace the split one. + pa.poly = newfaces[0] + pa.holes = newfaces[0:findex] + [newface0, newface1] + \ + newfaces[findex + 1:] + else: + # A join + (_, findex, othfindex, newface0) = splitjoin + if findex == 0 or othfindex == 0: + # Outer poly was joined to one hole. + pa.poly = newface0 + pa.holes = [f for f in newfaces if f is not None] + else: + # Two holes were joined + pa.poly = newfaces[0] + pa.holes = [f for f in newfaces if f is not None] + \ + [newface0] + self.inneroffsets = [Offset(pa, newt, self.vspeed)] + if pa2: + self.inneroffsets.append(Offset(pa2, newt, self.vspeed)) + if nexttarget > TOL: + for o in self.inneroffsets: + o.Build(nexttarget) + + def FaceAtSpokeEnds(self, f, t): + """Return a new face that is at the spoke ends of face f at time t. + + Also merges any adjacent approximately equal vertices into one vertex, + so returned list may be smaller than len(f). + Also sets the destindex fields of the spokes to the vertex they + will now end at. + + Args: + f: list of Spoke - one of self.faces + t: float - time in this offset + Returns: + list of int - indices into self.polyarea.points + (which has been extended with new ones) + """ + + newface = [] + points = self.polyarea.points + for i in range(0, len(f)): + s = f[i] + vcoords = s.EndPoint(t, points, self.vspeed) + v = points.AddPoint(vcoords) + if newface: + if v == newface[-1]: + s.destindex = len(newface) - 1 + elif i == len(f) - 1 and v == newface[0]: + s.destindex = 0 + else: + newface.append(v) + s.destindex = len(newface) - 1 + else: + newface.append(v) + s.destindex = 0 + s.dest = v + return newface + + def MakeNewFaces(self, t): + """For each face in this offset, make new face extending spokes + to time t. + + Args: + t: double - time + Returns: + list of list of int - list of new faces + """ + + ans = [] + for f in self.facespokes: + newf = self.FaceAtSpokeEnds(f, t) + if len(newf) > 2: + ans.append(newf) + return ans + + def SplitJoinFaces(self, newfaces, ev): + """Use event ev to split or join faces. + + Given ev, an edge event, use the ev spoke to split the + other spoke's inner edge. + If the ev spoke's face and other's face are the same, this splits the + face into two; if the faces are different, it joins them into one. + We have just made faces at the end of the spokes. + We have to remove the edge going from the other spoke to its + next spoke, and replace it with two edges, going to and from + the event spoke's destination. + General situation: + __ s ____ + c\ b\ | /a /e + \ \|/ / + f----------------g + / d \ + o/ \h + + where sd is the event spoke and of is the "other spoke", + hg is a spoke, and cf, fg. ge, ad, and db are edges in + the new inside face. + What we are to do is to split fg into two edges, with the + joining point attached where b,s,a join. + There are a bunch of special cases: + - one of split fg edges might have zero length because end points + are already coincident or nearly coincident. + - maybe c==b or e==a + + Args: + newfaces: list of list of int - the new faces + ev: OffsetEvent - an edge event + Side Effects: + faces in newfaces that are involved in split or join are + set to None + Returns: one of: + ('split', int, list of int, list of int) - int is the index in + newfaces of the face that was split, two lists are the + split faces + ('join', int, int, list of int) - two ints are the indices in + newfaces of the faces that were joined, and the list is + the joined face + """ + + # print("SplitJoinFaces", newfaces, ev) + spoke = ev.spoke + other = ev.other + findex = spoke.face + othfindex = other.face + newface = newfaces[findex] + othface = newfaces[othfindex] + nnf = len(newface) + nonf = len(othface) + d = spoke.destindex + f = other.destindex + c = (f - 1) % nonf + g = (f + 1) % nonf + e = (f + 2) % nonf + a = (d - 1) % nnf + b = (d + 1) % nnf + # print("newface=", newface) + # if findex != othfindex: print("othface=", othface) + # print("d=", d, "f=", f, "c=", c, "g=", g, "e=", e, "a=", a, "b=", b) + newface0 = [] + newface1 = [] + # The two new faces put spoke si's dest on edge between + # pi's dest and qi (edge after pi)'s dest in original face. + # These are indices in the original face; the current dest face + # may have fewer elements because of merging successive points + if findex == othfindex: + # Case where splitting one new face into two. + # The new new faces are: + # [d, g, e, ..., a] and [d, b, ..., c, f] + # (except we actually want the vertex numbers at those positions) + newface0 = [newface[d]] + i = g + while i != d: + newface0.append(newface[i]) + i = (i + 1) % nnf + newface1 = [newface[d]] + i = b + while i != f: + newface1.append(newface[i]) + i = (i + 1) % nnf + newface1.append(newface[f]) + # print("newface0=", newface0, "newface1=", newface1) + # now the destindex values for the spokes are messed up + # but I don't think we need them again + newfaces[findex] = None + return ('split', findex, newface0, newface1) + else: + # Case where joining two faces into one. + # The new face is splicing d's face between + # f and g in other face (or the reverse of all of that). + newface0 = [othface[i] for i in range(0, f + 1)] + newface0.append(newface[d]) + i = b + while i != d: + newface0.append(newface[i]) + i = (i + 1) % nnf + newface0.append(newface[d]) + if g != 0: + newface0.extend([othface[i] for i in range(g, nonf)]) + # print("newface0=", newface0) + newfaces[findex] = None + newfaces[othfindex] = None + return ('join', findex, othfindex, newface0) + + def InnerPolyAreas(self): + """Return the interior of the offset (and contained offsets) as + PolyAreas. + + Returns: + geom.PolyAreas + """ + + ans = geom.PolyAreas() + ans.points = self.polyarea.points + _AddInnerAreas(self, ans) + return ans + + def MaxAmount(self): + """Returns the maximum offset amount possible. + Returns: + float - maximum amount + """ + + # Need to do Build on a copy of points + # so don't add points that won't be used when + # really do a Build with a smaller amount + test_points = geom.Points() + test_points.AddPoints(self.polyarea.points) + save_points = self.polyarea.points + self.polyarea.points = test_points + self.Build() + max_amount = self._MaxTime() + self.polyarea.points = save_points + return max_amount + + def _MaxTime(self): + if self.inneroffsets: + return max([o._MaxTime() for o in self.inneroffsets]) + else: + return self.timesofar + self.endtime + + +def _AddInnerAreas(off, polyareas): + """Add the innermost areas of offset off to polyareas. + + Assume that polyareas is already using the proper shared points. + + Arguments: + off: Offset + polyareas: geom.PolyAreas + Side Effects: + Any non-zero-area faces in the very inside of off are + added to polyareas. + """ + + if off.inneroffsets: + for o in off.inneroffsets: + _AddInnerAreas(o, polyareas) + else: + newpa = geom.PolyArea(polyareas.points) + for i, f in enumerate(off.facespokes): + newface = off.FaceAtSpokeEnds(f, off.endtime) + area = abs(geom.SignedArea(newface, polyareas.points)) + if area < AREATOL: + if i == 0: + break + else: + continue + if i == 0: + newpa.poly = newface + newpa.data = off.polyarea.data + else: + newpa.holes.append(newface) + if newpa.poly: + polyareas.polyareas.append(newpa) |