1 # ##### BEGIN GPL LICENSE BLOCK #####
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.
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.
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.
17 # ##### END GPL LICENSE BLOCK #####
21 """Creating offset polygons inside faces."""
23 __author__
= "howard.trickey@gmail.com"
28 from .triquad
import Sub2
, Add2
, Angle
, Ccw
, Normalized2
, Perp2
, Length2
, \
30 from .geom
import Points
36 """A Spoke is a line growing from an outer vertex to an inner one.
38 A Spoke is contained in an Offset (see below).
41 origin: int - index of origin point in a Points
42 dest: int - index of dest point
43 is_reflex: bool - True if spoke grows from a reflex angle
44 dir: (float, float, float) - direction vector (normalized)
45 speed: float - at time t, other end of spoke is
46 origin + t*dir. Speed is such that the wavefront
47 from the face edges moves at speed 1.
48 face: int - index of face containing this Spoke, in Offset
49 index: int - index of this Spoke in its face
50 destindex: int - index of Spoke dest in its face
53 def __init__(self
, v
, prev
, next
, face
, index
, points
):
54 """Set attribute of spoke from points making up initial angle.
56 The spoke grows from an angle inside a face along the bisector
57 of that angle. Its speed is 1/sin(.5a), where a is the angle
58 formed by (prev, v, next). That speed means that the perpendicular
59 from the end of the spoke to either of the prev->v or v->prev
60 edges will grow at speed 1.
63 v: int - index of point spoke grows from
64 prev: int - index of point before v on boundary (in CCW order)
65 next: int - index of point after v on boundary (in CCW order)
66 face: int - index of face containing this spoke, in containing offset
67 index: int - index of this spoke in its face
68 points: geom.Points - maps vertex indices to 3d coords
80 uin
= Normalized2(Sub2(vp
, prevp
))
81 uout
= Normalized2(Sub2(nextp
, vp
))
82 uavg
= Normalized2((0.5 * (uin
[0] + uout
[0]), \
83 0.5 * (uin
[1] + uout
[1])))
84 if abs(Length2(uavg
)) < TOL
:
85 # in and out vectors are reverse of each other
86 self
.dir = (uout
[0], uout
[1], 0.0)
87 self
.is_reflex
= False
90 # bisector direction is 90 degree CCW rotation of
91 # average incoming/outgoing
92 self
.dir = (-uavg
[1], uavg
[0], 0.0)
93 self
.is_reflex
= Ccw(next
, v
, prev
, points
)
94 ang
= Angle(prev
, v
, next
, points
) # in range [0, 180)
95 sin_half_ang
= math
.sin(math
.pi
* ang
/ 360.0)
96 if abs(sin_half_ang
) < TOL
:
99 self
.speed
= 1.0 / sin_half_ang
102 """Printing representation of a Spoke."""
104 return "@%d+%gt%s <%d,%d>" % (self
.origin
, \
105 self
.speed
, str(self
.dir), \
106 self
.face
, self
.index
)
108 def EndPoint(self
, t
, points
, vspeed
):
109 """Return the coordinates of the non-origin point at time t.
112 t: float - time to end of spoke
113 points: geom.Points - coordinate map
114 vspeed: float - speed in z direction
116 (float, float, float) - coords of spoke's endpoint at time t
119 p
= points
.pos
[self
.origin
]
122 return (p
[0] + v
* t
* d
[0], p
[1] + v
* t
* d
[1], p
[2] + vspeed
* t
)
124 def VertexEvent(self
, other
, points
):
125 """Intersect self with other spoke, and return the OffsetEvent, if any.
127 A vertex event is with one advancing spoke intersects an adjacent
128 adavancing spoke, forming a new vertex.
131 other: Spoke - other spoke to intersect with
134 None or OffsetEvent - if there's an intersection in the growing
135 directions of the spokes, will return the OffsetEvent for
137 if lines are collinear or parallel, return None
141 a
= vmap
[self
.origin
]
142 b
= Add2(a
, self
.dir)
143 c
= vmap
[other
.origin
]
144 d
= Add2(c
, other
.dir)
145 # find intersection of line ab with line cd
151 # lines or neither parallel nor collinear
152 si
= Perp2(v
, w
) / pp
153 ti
= Perp2(u
, w
) / pp
154 if si
>= 0 and ti
>= 0:
155 p
= LinInterp2(a
, b
, si
)
156 dist_ab
= si
* Length2(u
)
157 dist_cd
= ti
* Length2(v
)
158 time_ab
= dist_ab
/ self
.speed
159 time_cd
= dist_cd
/ other
.speed
160 time
= max(time_ab
, time_cd
)
161 return OffsetEvent(True, time
, p
, self
, other
)
164 def EdgeEvent(self
, other
, offset
):
165 """Intersect self with advancing edge and return OffsetEvent, if any.
167 An edge event is when one advancing spoke intersects an advancing
168 edge. Advancing edges start out as face edges and move perpendicular
169 to them, at a rate of 1. The endpoints of the edge are the advancing
170 spokes on either end of the edge (so the edge shrinks or grows as
171 it advances). At some time, the edge may shrink to nothing and there
172 will be no EdgeEvent after that time.
174 We represent an advancing edge by the first spoke (in CCW order
175 of face) of the pair of defining spokes.
177 At time t, end of this spoke is at
179 where o=self.origin, d=self.dir, s= self.speed.
180 The advancing edge line has this equation:
182 where oo, od, os are o, d, s for other spoke, and p is direction
183 vector parallel to advancing edge, and a is a real parameter.
184 Equating x and y of intersection point:
186 o.x + d.x*s*t = oo.x + od.x*os*t + p.x*w
187 o.y + d.y*s*t = oo.y + od.y*os*t + p.y*w
189 which can be rearranged into the form
197 other: Spoke - the edge out of this spoke's origin is the advancing
198 edge to be checked for intersection
199 offset: Offset - the containing Offset
201 None or OffsetEvent - with data about the intersection, if any
204 vmap
= offset
.polyarea
.points
.pos
205 o
= vmap
[self
.origin
]
206 oo
= vmap
[other
.origin
]
207 otherface
= offset
.facespokes
[other
.face
]
208 othernext
= otherface
[(other
.index
+ 1) % len(otherface
)]
209 oonext
= vmap
[othernext
.origin
]
210 p
= Normalized2(Sub2(oonext
, oo
))
213 b
= other
.dir[0] * other
.speed
- self
.dir[0] * self
.speed
214 e
= other
.dir[1] * other
.speed
- self
.dir[1] * self
.speed
220 t
= (d
- f
* a
/ c
) / dem
227 t
= (a
- c
* d
/ f
) / dem
234 # intersection is in backward direction along self spoke
237 # intersection on wrong side of first end of advancing line segment
239 # calculate the equivalent of w for the other end
240 aa
= o
[0] - oonext
[0]
241 dd
= o
[1] - oonext
[1]
242 bb
= othernext
.dir[0] * othernext
.speed
- self
.dir[0] * self
.speed
243 ee
= othernext
.dir[1] * othernext
.speed
- self
.dir[1] * self
.speed
247 ww
= (aa
- bb
* t
) / cc
249 ww
= (dd
- ee
* t
) / ff
254 evertex
= (o
[0] + self
.dir[0] * self
.speed
* t
, \
255 o
[1] + self
.dir[1] * self
.speed
* t
)
256 return OffsetEvent(False, t
, evertex
, self
, other
)
259 class OffsetEvent(object):
260 """An event involving a spoke during offset computation.
262 The events kinds are:
263 vertex event: the spoke intersects an adjacent spoke and makes a new
265 edge event: the spoke hits an advancing edge and splits it
268 is_vertex_event: True if this is a vertex event (else it is edge event)
269 time: float - time at which it happens (edges advance at speed 1)
270 event_vertex: (float, float) - intersection point of event
271 spoke: Spoke - the spoke that this event is for
272 other: Spoke - other spoke involved in event; if vertex event, this will
273 be an adjacent spoke that intersects; if an edge event, this is the
274 spoke whose origin's outgoing edge grows to hit this event's spoke
277 def __init__(self
, isv
, time
, evertex
, spoke
, other
):
278 """Creates and initializes attributes of an OffsetEvent."""
280 self
.is_vertex_event
= isv
282 self
.event_vertex
= evertex
287 """Printing representation of an event."""
289 if self
.is_vertex_event
:
293 return "%s t=%5f %s %s %s" % (c
, self
.time
, str(self
.event_vertex
), \
294 repr(self
.spoke
), repr(self
.other
))
297 class Offset(object):
298 """Represents an offset polygonal area, and used to construct one.
300 Currently, the polygonal area must lie approximately in the XY plane.
301 As well as growing inwards in that plane, the advancing lines also
302 move in the Z direction at the rate of vspeed.
305 polyarea: geom.PolyArea - the area we are offsetting from.
306 We share the polyarea.points, and add to it as points in
307 the offset polygonal area are computed.
308 facespokes: list of list of Spoke - each sublist is a closed face
309 (oriented CCW); the faces may mutually interfere.
310 These lists are spokes for polyarea.poly + polyarea.holes.
311 endtime: float - time when this offset hits its first
312 event (relative to beginning of this offset), or the amount
313 that takes this offset to the end of the total Build time
314 timesofar: float - sum of times taken by all containing Offsets
315 vspeed: float - speed that edges move perpendicular to offset plane
316 inneroffsets: list of Offset - the offsets that take over after this
320 def __init__(self
, polyarea
, time
, vspeed
):
321 """Set up initial state of Offset from a polyarea.
324 polyarea: geom.PolyArea
325 time: float - time so far
328 self
.polyarea
= polyarea
331 self
.timesofar
= time
333 self
.inneroffsets
= []
334 self
.InitFaceSpokes(polyarea
.poly
)
335 for f
in polyarea
.holes
:
336 self
.InitFaceSpokes(f
)
339 ans
= ["Offset: endtime=%g" % self
.endtime
]
340 for i
, face
in enumerate(self
.facespokes
):
341 ans
.append(("<%d>" % i
) + str([str(spoke
) for spoke
in face
]))
342 return '\n'.join(ans
)
344 def PrintNest(self
, indent_level
=0):
345 indent
= " " * indent_level
* 4
346 print(indent
+ "Offset timesofar=", self
.timesofar
, "endtime=",
348 print(indent
+ " polyarea=", self
.polyarea
.poly
, self
.polyarea
.holes
)
349 for o
in self
.inneroffsets
:
350 o
.PrintNest(indent_level
+ 1)
352 def InitFaceSpokes(self
, face_vertices
):
353 """Initialize the offset representation of a face from vertex list.
355 If the face has no area or too small an area, don't bother making it.
358 face_vertices: list of int - point indices for boundary of face
360 A new face (list of spokes) may be added to self.facespokes
363 n
= len(face_vertices
)
366 points
= self
.polyarea
.points
367 area
= abs(geom
.SignedArea(face_vertices
, points
))
370 findex
= len(self
.facespokes
)
371 fspokes
= [Spoke(v
, face_vertices
[(i
- 1) % n
], \
372 face_vertices
[(i
+ 1) % n
], findex
, i
, points
) \
373 for i
, v
in enumerate(face_vertices
)]
374 self
.facespokes
.append(fspokes
)
376 def NextSpokeEvents(self
, spoke
):
377 """Return the OffsetEvents that will next happen for a given spoke.
379 It might happen that some events happen essentially simultaneously,
380 and also it is convenient to separate Edge and Vertex events, so
382 But, for vertex events, only look at the event with the next Spoke,
383 as the event with the previous spoke will be accounted for when we
384 consider that previous spoke.
387 spoke: Spoke - a spoke in one of the faces of this object
389 (float, list of OffsetEvent, list of OffsetEvent) -
391 next Vertex event list and next Edge event list
394 facespokes
= self
.facespokes
[spoke
.face
]
399 # First find vertex event (only the one with next spoke)
400 next_spoke
= facespokes
[(spoke
.index
+ 1) % n
]
401 ev
= spoke
.VertexEvent(next_spoke
, self
.polyarea
.points
)
405 # Now find edge events, if this is a reflex vertex
407 prev_spoke
= facespokes
[(spoke
.index
- 1) % n
]
408 for f
in self
.facespokes
:
410 if other
== spoke
or other
== prev_spoke
:
412 ev
= spoke
.EdgeEvent(other
, self
)
414 if ev
.time
< bestt
- TOL
:
418 if abs(ev
.time
- bestt
) < TOL
:
420 return (bestt
, bestv
, beste
)
422 def Build(self
, target
=2e100
):
423 """Build the complete Offset structure or up until target time.
425 Find the next event(s), makes the appropriate inner Offsets
426 that are inside this one, and calls Build on those Offsets to continue
427 the process until only a single point is left or time reaches target.
432 for f
in self
.facespokes
:
434 (t
, ve
, ee
) = self
.NextSpokeEvents(s
)
438 if abs(t
- bestt
) < TOL
:
439 bestevs
[0].extend(ve
)
440 bestevs
[1].extend(ee
)
442 # could happen if polygon is oriented wrong
443 # or in other special cases
446 # seems to be in a loop, so quit
452 if target
< self
.endtime
:
453 self
.endtime
= target
454 newfaces
= self
.MakeNewFaces(self
.endtime
)
456 # Only vertex events.
457 # Merging of successive vertices in inset face will
458 # take care of the vertex events
459 newfaces
= self
.MakeNewFaces(self
.endtime
)
462 # First make the new faces (handles all vertex events)
463 newfaces
= self
.MakeNewFaces(self
.endtime
)
464 # Only do one edge event (handle other simultaneous edge
465 # events in subsequent recursive Build calls)
467 splitjoin
= self
.SplitJoinFaces(newfaces
, ee
[0])
468 nexttarget
= target
- self
.endtime
469 if len(newfaces
) > 0:
470 pa
= geom
.PolyArea(points
=self
.polyarea
.points
)
471 pa
.data
= self
.polyarea
.data
472 newt
= self
.timesofar
+ self
.endtime
473 pa2
= None # may make another
475 pa
.poly
= newfaces
[0]
476 pa
.holes
= newfaces
[1:]
477 elif splitjoin
[0] == 'split':
478 (_
, findex
, newface0
, newface1
) = splitjoin
480 # Outer poly of polyarea was split.
481 # Now there will be two polyareas.
482 # If there were holes, need to allocate according to
483 # which one contains the holes.
485 pa2
= geom
.PolyArea(points
=self
.polyarea
.points
)
486 pa2
.data
= self
.polyarea
.data
488 if len(newfaces
) > 1:
489 # print("need to allocate holes")
490 for hf
in newfaces
[1:]:
491 if pa
.ContainsPoly(hf
, self
.polyarea
.points
):
492 # print("add", hf, "to", pa.poly)
494 elif pa2
.ContainsPoly(hf
, self
.polyarea
.points
):
495 # print("add", hf, "to", pa2.poly)
498 print("whoops, hole in neither poly!")
499 self
.inneroffsets
= [Offset(pa
, newt
, self
.vspeed
), \
500 Offset(pa2
, newt
, self
.vspeed
)]
502 # A hole was split. New faces just replace the split one.
503 pa
.poly
= newfaces
[0]
504 pa
.holes
= newfaces
[0:findex
] + [newface0
, newface1
] + \
505 newfaces
[findex
+ 1:]
508 (_
, findex
, othfindex
, newface0
) = splitjoin
509 if findex
== 0 or othfindex
== 0:
510 # Outer poly was joined to one hole.
512 pa
.holes
= [f
for f
in newfaces
if f
is not None]
514 # Two holes were joined
515 pa
.poly
= newfaces
[0]
516 pa
.holes
= [f
for f
in newfaces
if f
is not None] + \
518 self
.inneroffsets
= [Offset(pa
, newt
, self
.vspeed
)]
520 self
.inneroffsets
.append(Offset(pa2
, newt
, self
.vspeed
))
522 for o
in self
.inneroffsets
:
525 def FaceAtSpokeEnds(self
, f
, t
):
526 """Return a new face that is at the spoke ends of face f at time t.
528 Also merges any adjacent approximately equal vertices into one vertex,
529 so returned list may be smaller than len(f).
530 Also sets the destindex fields of the spokes to the vertex they
534 f: list of Spoke - one of self.faces
535 t: float - time in this offset
537 list of int - indices into self.polyarea.points
538 (which has been extended with new ones)
542 points
= self
.polyarea
.points
543 for i
in range(0, len(f
)):
545 vcoords
= s
.EndPoint(t
, points
, self
.vspeed
)
546 v
= points
.AddPoint(vcoords
)
549 s
.destindex
= len(newface
) - 1
550 elif i
== len(f
) - 1 and v
== newface
[0]:
554 s
.destindex
= len(newface
) - 1
561 def MakeNewFaces(self
, t
):
562 """For each face in this offset, make new face extending spokes
568 list of list of int - list of new faces
572 for f
in self
.facespokes
:
573 newf
= self
.FaceAtSpokeEnds(f
, t
)
578 def SplitJoinFaces(self
, newfaces
, ev
):
579 """Use event ev to split or join faces.
581 Given ev, an edge event, use the ev spoke to split the
582 other spoke's inner edge.
583 If the ev spoke's face and other's face are the same, this splits the
584 face into two; if the faces are different, it joins them into one.
585 We have just made faces at the end of the spokes.
586 We have to remove the edge going from the other spoke to its
587 next spoke, and replace it with two edges, going to and from
588 the event spoke's destination.
597 where sd is the event spoke and of is the "other spoke",
598 hg is a spoke, and cf, fg. ge, ad, and db are edges in
600 What we are to do is to split fg into two edges, with the
601 joining point attached where b,s,a join.
602 There are a bunch of special cases:
603 - one of split fg edges might have zero length because end points
604 are already coincident or nearly coincident.
608 newfaces: list of list of int - the new faces
609 ev: OffsetEvent - an edge event
611 faces in newfaces that are involved in split or join are
614 ('split', int, list of int, list of int) - int is the index in
615 newfaces of the face that was split, two lists are the
617 ('join', int, int, list of int) - two ints are the indices in
618 newfaces of the faces that were joined, and the list is
622 # print("SplitJoinFaces", newfaces, ev)
626 othfindex
= other
.face
627 newface
= newfaces
[findex
]
628 othface
= newfaces
[othfindex
]
638 # print("newface=", newface)
639 # if findex != othfindex: print("othface=", othface)
640 # print("d=", d, "f=", f, "c=", c, "g=", g, "e=", e, "a=", a, "b=", b)
643 # The two new faces put spoke si's dest on edge between
644 # pi's dest and qi (edge after pi)'s dest in original face.
645 # These are indices in the original face; the current dest face
646 # may have fewer elements because of merging successive points
647 if findex
== othfindex
:
648 # Case where splitting one new face into two.
649 # The new new faces are:
650 # [d, g, e, ..., a] and [d, b, ..., c, f]
651 # (except we actually want the vertex numbers at those positions)
652 newface0
= [newface
[d
]]
655 newface0
.append(newface
[i
])
657 newface1
= [newface
[d
]]
660 newface1
.append(newface
[i
])
662 newface1
.append(newface
[f
])
663 # print("newface0=", newface0, "newface1=", newface1)
664 # now the destindex values for the spokes are messed up
665 # but I don't think we need them again
666 newfaces
[findex
] = None
667 return ('split', findex
, newface0
, newface1
)
669 # Case where joining two faces into one.
670 # The new face is splicing d's face between
671 # f and g in other face (or the reverse of all of that).
672 newface0
= [othface
[i
] for i
in range(0, f
+ 1)]
673 newface0
.append(newface
[d
])
676 newface0
.append(newface
[i
])
678 newface0
.append(newface
[d
])
680 newface0
.extend([othface
[i
] for i
in range(g
, nonf
)])
681 # print("newface0=", newface0)
682 newfaces
[findex
] = None
683 newfaces
[othfindex
] = None
684 return ('join', findex
, othfindex
, newface0
)
686 def InnerPolyAreas(self
):
687 """Return the interior of the offset (and contained offsets) as
694 ans
= geom
.PolyAreas()
695 ans
.points
= self
.polyarea
.points
696 _AddInnerAreas(self
, ans
)
700 """Returns the maximum offset amount possible.
702 float - maximum amount
705 # Need to do Build on a copy of points
706 # so don't add points that won't be used when
707 # really do a Build with a smaller amount
708 test_points
= geom
.Points()
709 test_points
.AddPoints(self
.polyarea
.points
, True)
710 save_points
= self
.polyarea
.points
711 self
.polyarea
.points
= test_points
713 max_amount
= self
._MaxTime
()
714 self
.polyarea
.points
= save_points
718 if self
.inneroffsets
:
719 return max([o
._MaxTime
() for o
in self
.inneroffsets
])
721 return self
.timesofar
+ self
.endtime
724 def _AddInnerAreas(off
, polyareas
):
725 """Add the innermost areas of offset off to polyareas.
727 Assume that polyareas is already using the proper shared points.
731 polyareas: geom.PolyAreas
733 Any non-zero-area faces in the very inside of off are
738 for o
in off
.inneroffsets
:
739 _AddInnerAreas(o
, polyareas
)
741 newpa
= geom
.PolyArea(polyareas
.points
)
742 for i
, f
in enumerate(off
.facespokes
):
743 newface
= off
.FaceAtSpokeEnds(f
, off
.endtime
)
744 area
= abs(geom
.SignedArea(newface
, polyareas
.points
))
752 newpa
.data
= off
.polyarea
.data
754 newpa
.holes
.append(newface
)
756 polyareas
.polyareas
.append(newpa
)