What's wrong with this pseudocode for Forster-Overfelt's version of the Greiner-Horman polygon clipping algorithm?

Question

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?

Answer 1

In this particular case it appears that the entry/exit flags will be set correctly by your pseudo-code. It is only the walking of the tree at the end which is incorrect. Since this tree walking is the same as the original Greiner-Horman it is not really explained well in the paper. It is the last paragraph of section 2 of the paper. If you have read the paper then you know the notation used to represent the various flags that are set by your pseudo-code, as in Figure 2 of the paper. There are four intersections which I label as

I1=(5,6), I2=(4,6), I3=(3,6) and I4=(2,6).  

This is from traversing the subject polygon, but you could also walk the clipping polygon. Working this out with pencil and paper instead of with your pseudo-code I get the following flags. For notation the Ss are the corners as listed in the subjpoly array and the Cs are the corners of the clippoly array and the I’s are the inserted intersections.

The corners and intersections of the subject polygon:

S1   S2   S3  I1  I2   I3  I4  S4 
              en  ex   en  ex

And clipping:

C1   I3   C2  I2   C3   I1   C4   C5   I4
     ex       en        ex             en

The only tricky ones are the (on,on) intersections of the I2 and I3 intersections of the subject polygon (seeing that the I’s to either side are ‘on’). These are labeled by examining the neighbor, or clipping, polygon and setting to the opposite. Since the clipping polygon has only simple (out,in) and (in,out) for (prev,next) since these intersections are surround by simple corners, there is no doubt what the subject polygon labels are for these intersections. I see you did this correctly in your pseudo-code. Checking we can see that all the intersections are consistent between subject and clipping, i.e. always opposite for this simple case. I think your pseudo code would correctly produce this but you should check.

Next is the traversal. In your case you started with the I1 intersection on the subject polygon. Since it is an “enter” intersection, we stay on the subject. The traversal tree looks like:

I1 (subject) -> I2 -> I2 (clip) -> C3 -> I1 -> I1 (subject)

Since I2 on the subject polygon an “exit” we switch from subject to clipping polygon and follow it back to I1 which completes the loop. I see that is just what you have.

Now to continue it is not obvious how to get to the next intersection but you seen to have done the correct thing by traversing the clipping polygon through C4 and C5 to find I4. Now I4 is an “entry” for the clipping polygon so the traversal is

 I4 (clip) -> C1 -> I3 -> I3 (subject) ->  I4 -> I4 (clip)

And that uses up all of the corners and intersections of the clipping polygon and it appears to have worked. I’m not sure why your second polygon is off. It appears that you have I4=(1,6) labeled as an “exit” for the clipping polygon and jumped to the neighbor, or subject, polygon.