# ##### 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 #####
#
"""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)