# SPDX-License-Identifier: GPL-2.0-or-later """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) if newfaces: 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): r"""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, True) 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)