What's wrong with this pseudocode for Forster-Overfelt's version of the Greiner-Horman polygon clipping algorithm?
The Problem
I'm trying to understand and implement the Forster-Overfelt version of the Greiner-Horman polygon clipping algorithm. I've read the other Stackoverflow post about clarifying this algorithm, but I still can't seem to get it to work.
I know there's something wrong bc it produces the wrong intersection of two polygons even for a simple example that doesnt have degenerates:
subjpoly = [(0,0),(6,0),(6,6),(0,6),(0,0)]
clippoly = [(1,4),(3,8),(5,4),(5,10),(1,10),(1,4)]
which produces an intersection of:
[ [(5.0, 6.0), (4.0, 6.0), (5, 4), (5.0, 6.0)],
[(1.0, 6.0), (2.0, 6.0), (4.0, 6.0)] ]
Visualized it looks like this:
So this question is not about a particular piece of code or language syntax, but about understanding the algorithm and putting it into pseudocode. Thanks in advance!
The Algorithm
The algorithm presented here is based directly off and should mimic the one described in section 4.2 in the Forster-Oberfelt paper, but obviously there is something I'm missing that's giving me wrong results.
Part 1:
Start by looping both subj and clip and marking each vertex location as "in", "out", or "on" the other polygon.
for s in subj.iter():
s.loc = testLocation(s, clip)
for c in clip.iter():
c.loc = testLocation(c, subj)
Part 2:
Proceed to loop the intersection points of the subj polygon
for s in subj.iter():
if s.intersect:
Subpart 2.1:
Process each subj intersection by either unmarking them as an intersection or marking them as entry or exit, and do the same for the neighbour intersection point. NOTE: the algorithm explained in the article only explains how to mark the main subject polygon, but never says how to mark the neighbour, so here I'm just assuming both are marked using the same set of rules.
mark(s)
mark(s.neighbour)
Where the mark() processing rules is defined as:
def mark(s):
curlocs = (s.prev.loc,s.next.loc)
neighlocs = (s.neighbour.prev.loc,s.neighbour.next.loc)
# in in
if curlocs == ("in","in"):
if neighlocs == ("in","in")\
or neighlocs == ("out","out")\
or neighlocs == ("on","on"):
s.intersect = False
else:
s.entry = True
# out out
elif curlocs == ("out","out"):
if neighlocs == ("in","in")\
or neighlocs == ("out","out")\
or neighlocs == ("on","on"):
s.intersect = False
else:
s.entry = False
# on on
elif curlocs == ("on","on"):
if neighlocs == ("in","in")\
or neighlocs == ("out","out")\
or neighlocs == ("on","on"):
s.intersect = False
else:
# label opposite of neighbour
# NOTE: this is not specified in the article,
# but one cannot take the opposite of the neighbour's entry flag
# if the neighbour hasn't been marked yet,
# thus the decision to mark the neighbour first
mark(s.neighbour)
s.entry = not s.neighbour
# partial exit
elif curlocs == ("in","on")\
or curlocs == ("on","out"):
s.entry = False
# partial entry
elif curlocs == ("on","in")\
or curlocs == ("out","on"):
s.entry = True
# normal exit
elif curlocs == ("in","out"):
s.entry = False
# normal entry
elif curlocs == ("out","in"):
s.entry = True
Subpart 2.2:
Finally make sure curr and neighbour dont have same entry or exit flags; if they do disable their intersection flag and change the location flags.
if s.entry and s.neighbour.entry:
s.intersect = False
s.neighbour.intersect = False
s.loc = "in"
elif not s.entry and not s.neighbour.entry:
s.intersect = False
s.neighbour.intersect = False
s.loc = "out"
Bonus Question
A bonus-question is how to make this algorithm support both union and intersection operations since the original Greiner algorithm's support for union was by simply inverting the initial entry/exit flag, but this Forster algorithm doesn't use such a flag?
One more comment about unions instead of intersections. The main idea is that a union operation will proceed in the opposite direction as compared to the intersection operation. That is, if one is suppose to move backwards along a polygon for the intersection operation then for the union operation one would move forward and vice versa.
Now onto the algorithm: first, let's start with the outline of the algorithm. The algorithm I have here will only create one polygon per intersection operation, so you will have to adapt it to create more than one.
'''
The following is an adaptation of the above Greiner-Hormann* algorithm to deal
with degenerate cases. The adaptation was briefly described by Liu et al.**
*Greiner, G. and Hormann K., Efficient Clipping of Arbitrary Polygons, ACM
Trans. on Graphics, 17(2), 1998, pp.71-83
**Liu, Y. K., Wang X. Q., Bao S. Z., Gombosi M., and Zalik B, An Algorithm for
Polygon Clipping and for Determining Polygon Intersections and Unions, Comp. &
Geo, 33, pp. 589-598, 2007
'''
def clip(subject, constraint):
subject, constraint = inside_outside(subject, constraint) #label vertices as inside or outside
subject, constraint = poly_inters(subject, constraint) #find intersections
subject, constraint = label(subject, constraint) #label intersections and entry or exit and possibly remove
flag = True #loop flag
#set our current location to the first point in subject
current = subject.first
#loop through our polygon until we have found the first intersection
while flag:
current = current.next
#Either an intersection has been found or no intersections found
if current.intersect or current.pt == subject.first.pt:
flag = False
#reset our flag for the new loop
flag = True
#If a point lies outside of c and there was an intersection clip s
if current.intersect:
append(clipped, current.pt)
While flag:
#Entry
if current.en:
clipped, current = forward(clipped, current)
#Exit
else:
clipped, current = backward(clipped, current)
#Check to see if we have completed a polygon
if current.pt == clipped.first.pt:
#check if the polygon intersect at a point
if clipped.num_vertices is not 1:
#remove the last vertex because it is also the first
remove(clipped, clipped.last)
#we have created our polygon so we can exit
flag = .FALSE.
#change to the neighbor polygon since we are at a new intersection
current = current.neighbor
#Check if one polygon is contained wholly within the other
elif contained(subject, constraint):
clipped = subject
elif contained(subject, constraint):
clipped = constraint
return clipped
Now we can discuss the labelling. The following code is the loop for labelling intersections as inside or outside. It does not include the logic for determining inside/outside and is only the order of operations.
#label intersections as entry or exit
def label(poly1, poly2):
#cycle through first polygon and label intersections as en or ex
current = poly1.first
for i in range(0,poly1.num_vertices):
if current.intersect:
current = intersect_cases(current)
#Make sure current is still an intersection
if current.isect:
current.neighbor = intersect_cases(current.neighbor)
#if the intersection is en/en or ex/ex
if current.en == current.neighbor.en:
current = remove_inter(current)
current = current.next #move to the next point
return poly1, poly2
And finally dealing with the various cases for labelling.
#deal with the cases
#on/on, on/out, on/in, out/on, out/out, out/in, in/on, in/out, in/in
def intersect_cases(current):
neighbor = current.neighbor
#on/on
if current.prev.intersect and current.next.intersect:
#Determine what to do based on the neighbor
#en tag is the opposite of the neighbor's en tag
if neighbor.prev.intersect and neighbor.next.intersect:
current = remove_inter(current)
current.en = True
neighbor.en = True
elif neighbor.prev.intersect and not neighbor.next.en:
current.en = False
elif neighbor.prev.intersect and neighbor.next.en:
current.en = True
elif not neighbor.prev.en and neighbor.next.intersect:
current.en = False
elif not (neighbor.prev.en or neighbor.next.en):
current = remove_inter(current)
current.en = True
neighbor.en = False
elif not neighbor.prev.en and neighbor.next.en:
current.en = False
elif neighbor.prev.en and neighbor.next.isect:
current.en = True
elif neighbor.prev.en and not neighbor.next.en:
current.en = True
elif neighbor.prev.en and neighbor.next.en:
current = remove_inter(current)
current.en = False
neighbor.en = True
#on/out
elif current.prev.intersect and not current.next.en:
current.en = False
#on/in
elif current.prev.intersect and current.next.en:
current.en = True
#out/on
elif not current.prev.en and current.next.intersect:
current.en = True
#out/out
elif not (current.prev%en or current.next.en):
if neighbor.prev%intersect and neighbor.next.intersect:
current = remove_inter(current)
neighbor.en = True
elif neighbor.prev.en == neighbor.next.en:
current = remove_inter(current)
else:
if neighbor.prev.en and not neighbor.next.en:
current.en = True
else:
current.en = False
#out/in
elif not current.prev.en and current.next.en:
current.en = True
#in/on
elif current.prev.en and current.next.intersect:
current.en = False
#in/out
elif current.prev.en and not current.next.en:
current.en = False
#in/in
elif current.prev.en and current.next.en:
if neighbor.prev.intersect and neighbor.next.intersect:
current = remove_inter(current)
neighbor.en = False
elif neighbor.prev.en == neighbor.next.en:
current = remove_inter(current)
else:
if neighbor.prev.en and not neighbor.next.en:
current.en = True
else:
current.en = False
return current
The above code has not been tested and it is not written for efficiency, rather it is written for readability and understanding.