
<<FastIntersectImpl.mesa>>
    <<Copyright © 1985 by Xerox Corporation.  All rights reserved.>>
    <<Greene, October 3, 1986 5:15:16 pm PDT>>

DIRECTORY
Sweep;

FastIntersectImpl: CEDAR PROGRAM IMPORTS Sweep EXPORTS Sweep =

BEGIN OPEN Sweep;

LineCarrier: TYPE = REF LineCarrierRec;
LineCarrierRec: TYPE = RECORD[
carried: Line,
lineLeft, lineRight: LineCarrier _ NIL
];


    NilCopy: PUBLIC PROC [stateIn: REF ANY] RETURNS [stateOut: REF ANY] ~ {
        RETURN[stateIn]
        };

NilCombine: PUBLIC PROC [state1, state2: REF ANY] RETURNS [state: REF ANY] ~ {
    RETURN[NIL]
    };

    NilFlip: PUBLIC PROC [stateIn: REF ANY] RETURNS [stateOut: REF ANY] ~ {
        RETURN[NIL]
            };


    PowerOf2: PROC [in: CARDINAL] RETURNS [out: CARDINAL _ 32] ~ {
        UNTIL out >= in DO
            out _ out * 2;
            ENDLOOP;
        };

    ExpandIfNecessary: PROC [in: REF PointSeq] RETURNS [out: REF PointSeq] ~ INLINE {
        OPEN in;
        IF size = max THEN {
            out _ NEW[PointSeq[max*2]];
            FOR i: CARDINAL IN [0..size) DO out.data[i] _ data[i]  ENDLOOP;
            out.size _ size;
            }
        ELSE out _ in;
        };


Intersect: PUBLIC PROC [in: Graph, Copy: CopyLineProc _ NilCopy, Combine: CombineLineProc _ NilCombine, Flip: FlipLineProc _ 
NilFlip] RETURNS [out: Graph] ~ {
pointCount: CARDINAL _ 0; i: CARDINAL _ 0;
changed: BOOLEAN;

    FlipLine: PROC [line: Line] ~ {
    OPEN line;
        temp: Point _ above;  above _ below;  below _ temp;
        state _ Flip[state];
        };

    freeLC: LineCarrier _ NIL;

    NewLC: PROC [carried: Line, lineLeft, lineRight: LineCarrier _ NIL] RETURNS [alive: LineCarrier] ~ {
        IF freeLC = NIL THEN RETURN[NEW[LineCarrierRec_ [carried: carried, lineLeft: lineLeft, lineRight: lineRight]]] ELSE {
            alive _ freeLC;  freeLC _ alive.lineRight;
            alive.carried _ carried; alive.lineLeft _ lineLeft; alive.lineRight _ lineRight;
            }
        };

    TrashFreeLC: PROC [] ~ {
        UNTIL freeLC = NIL DO
            freeLC.lineLeft _ NIL; freeLC _ freeLC.lineRight;
            ENDLOOP;
        };

sortedIn: REF PointSeq _ NEW[PointSeq[in.points.max]];
indexIn: CARDINAL _ 0;
sortedOut: REF PointSeq _ NEW[PointSeq[in.points.max]];
switch: REF PointSeq;
heap: REF PointSeq _ in.points;

    Heapify: PROC [] ~ {
        FOR i: CARDINAL IN [1..heap.size) DO
            BubbleUp[i];
            ENDLOOP;
        };

    QueueInsert: PROC [p: Point] ~ {
        heap _ ExpandIfNecessary[heap];
        heap.data[heap.size] _ p; heap.size _ heap.size + 1;
        BubbleUp[heap.size - 1];
        };

    BubbleUp: PROC [pos: CARDINAL] ~ {
        pL: Point; r: Relation;
        pos _ pos + 1;  -- shift for convience
        pL _ heap.data[pos - 1];
        UNTIL pos = 1 DO
            pH: Point _ heap.data[pos/2 - 1];
            r _ ComparePoints[pL, pH];
            IF r # greater THEN {
                heap.data[pos - 1] _ pL;
                RETURN;
                }
            ELSE {
                heap.data[pos - 1] _ pH;
                pos _ pos / 2;
                };
            ENDLOOP;
        heap.data[0] _ pL;
        };

    HoleDown: PROC [pos: CARDINAL] ~ {
    OPEN heap;
    lower: CARDINAL; r: Relation;
        pos _ pos + 1;
        lower _ pos * 2;
        UNTIL lower >= size DO
            r _ ComparePoints[data[lower-1], data[lower]];
            IF r = less THEN {
                data[pos - 1] _ data[lower];
                pos _ lower + 1;
                }
            ELSE {
                data[pos - 1] _ data[lower - 1];
                pos _ lower;
                };
            lower _ pos*2;
            ENDLOOP;
        IF lower = size THEN {data[pos - 1] _ data[lower - 1]; size _ size - 1}
        ELSE {
            size _ size - 1;
            IF pos - 1 # size THEN {
                data[pos - 1] _ data[size];
                BubbleUp[pos - 1];
                };
            };
        };

    EmptyQueue: PROC [] RETURNS [BOOLEAN] ~ {
        RETURN [(indexIn = sortedIn.size) AND (heap.size = 0)]
        };

    NextEvent: PROC [] RETURNS [p: Point] ~ {
    OPEN heap;
    IF size = 0 THEN { -- heap is empty
        IF indexIn = sortedIn.size THEN ERROR ELSE {
            p _ sortedIn.data[indexIn];
            indexIn _ indexIn + 1;
            }
        }
    ELSE { -- heap is full
        WHILE size >= 3 AND ComparePoints[data[0], data[2]] = equal DO  --remove duplicates right tree
            MergePoints[data[0], data[2]];
            HoleDown[2];
            ENDLOOP;
        WHILE size >= 2 AND ComparePoints[data[0], data[1]] = equal DO  --remove duplicates left tree
            MergePoints[data[0], data[1]];
            HoleDown[1];
            ENDLOOP;
        IF indexIn = sortedIn.size THEN {
            p _ data[0];
            HoleDown[0];
            }
        ELSE {
            r: Relation _ ComparePoints[ data[0], sortedIn.data[indexIn] ];
            IF r = less THEN {
                p _ sortedIn.data[indexIn];
                indexIn _ indexIn + 1;
                }
            ELSE {
                IF r = greater THEN {
                    p _ data[0];
                    HoleDown[0]
                    }
                ELSE {
                    p _ data[0];
                    MergePoints[p, sortedIn.data[indexIn]];
                    indexIn _ indexIn + 1;
                    HoleDown[0];
                    };
                };
            };
            };
        sortedOut _ ExpandIfNecessary[sortedOut];
        sortedOut.data[sortedOut.size] _ p;
        sortedOut.size _ sortedOut.size + 1;
        };


    xOrder, neighborRight, leftMostInserted: LineCarrier _ NIL;

    OrderRemoveLines: PROC [p: Point] RETURNS [rightNeighbor: LineCarrier] ~ {
        orderScan: LineCarrier _ xOrder;
        leftNeighbor: LineCarrier;
        lineScan: Line _ p.incoming;
        oldLineScan: Line; count: INT;

        AdvanceLineCarrier: PROC [] ~ {
            rightNeighbor _ orderScan.lineRight;
            orderScan.lineRight _ freeLC;  freeLC _ orderScan;
            orderScan _ rightNeighbor;
            };
        IF p.incoming = NIL THEN RETURN[NIL];
        UNTIL lineScan = orderScan.carried DO    --bounds fault indicates algorithmic failure
            orderScan _ orderScan.lineLeft;
            ENDLOOP;
        leftNeighbor _ orderScan.lineLeft;
        WHILE leftNeighbor # NIL DO
            IF leftNeighbor.carried.below # p THEN EXIT;
            orderScan _ orderScan.lineLeft;
            leftNeighbor _ orderScan.lineLeft;
            ENDLOOP;
        UNTIL lineScan = NIL DO
            IF lineScan # orderScan.carried THEN
                {--Out of order acceptable up to slope inequality
                count _ 0;
                WHILE orderScan # NIL DO
                    IF (CompareSlopes[lineScan, orderScan.carried] # equal) OR (orderScan.carried.below # p) THEN EXIT;
                    AdvanceLineCarrier[];
                    count _ count + 1;
                    ENDLOOP;
                oldLineScan _ lineScan; lineScan _ lineScan.clockwiseAroundBelow;
                count _ count - 1;
                WHILE lineScan # NIL DO
                    IF CompareSlopes[lineScan, oldLineScan] # equal THEN EXIT;
                    lineScan _ lineScan.clockwiseAroundBelow;
                    count _ count - 1;
                    ENDLOOP;
                IF count # 0 THEN ERROR;
                }
            ELSE {
                lineScan _ lineScan.clockwiseAroundBelow;
                AdvanceLineCarrier[];
                };
            ENDLOOP;
        IF leftNeighbor # NIL THEN leftNeighbor.lineRight _ rightNeighbor;
        IF rightNeighbor # NIL THEN rightNeighbor.lineLeft _ leftNeighbor ELSE xOrder _ leftNeighbor;
        };

    OrderInsertLines: PROC [p: Point, rightNeighborIn: LineCarrier] RETURNS [leftMostInserted, rightNeighbor: LineCarrier] ~ {
        r: Relation;
        lineScan, garbage, topHalf, bottomHalf: Line;
        lc, broken, leftNeighbor: LineCarrier;
        firstTime: BOOLEAN;

        rightNeighbor _ rightNeighborIn;
        IF p.outgoing = NIL THEN RETURN[NIL, rightNeighbor];
        leftNeighbor _ IF rightNeighbor = NIL THEN xOrder ELSE rightNeighbor.lineLeft;
        UNTIL leftNeighbor = NIL DO
            r _ ComparePointLine[p, leftNeighbor.carried];
            IF r = greater THEN EXIT;
            IF r = equal THEN {
                broken _ leftNeighbor;
                leftNeighbor _ leftNeighbor.lineLeft;
                topHalf _ broken.carried;
                broken.carried _ NIL; broken.lineLeft _ broken.lineRight _ NIL;
                bottomHalf _ NEW[LineRec _ [above: p, below: topHalf.below, state: Copy[topHalf.state]]];
                RemoveLineFromPointBelow[topHalf];
                InsertLineInPointBelow[bottomHalf];
                InsertLineInPointAbove[bottomHalf];
                IF p.incoming # NIL THEN ERROR;
                topHalf.below _ p;
                InsertLineInPointBelow[topHalf];
                EXIT;
                };
            rightNeighbor _ leftNeighbor;
            leftNeighbor _ leftNeighbor.lineLeft;
            ENDLOOP;
        leftMostInserted _ rightNeighbor;
        lineScan _ p.outgoing;
        firstTime _ TRUE;
        UNTIL lineScan = NIL DO
            IF (NOT firstTime) AND (CompareSlopes[leftMostInserted.carried, lineScan] = equal) THEN
                {r _ ComparePoints[leftMostInserted.carried.below, lineScan.below];
                garbage _ lineScan;
                lineScan _ lineScan.clockwiseAroundAbove;
                IF r = equal THEN {
                    leftMostInserted.carried.state _ Combine[leftMostInserted.carried.state, garbage.state];
                    RemoveLineFromPointAbove[garbage]; RemoveLineFromPointBelow[garbage];
                    }
                ELSE {
                    IF r = greater THEN {topHalf _ leftMostInserted.carried; bottomHalf _ garbage}
                    ELSE {topHalf _ garbage; bottomHalf _ leftMostInserted.carried; leftMostInserted.carried _ topHalf};
                    topHalf.state _ Combine[topHalf.state, bottomHalf.state];
                    RemoveLineFromPointAbove[bottomHalf];
                    bottomHalf.above _ topHalf.below;
                    InsertLineInPointAbove[bottomHalf]
                    };
                }
            ELSE
                { lc _ NewLC[carried: lineScan, lineRight: leftMostInserted];
                IF leftMostInserted # NIL THEN leftMostInserted.lineLeft _ lc ELSE xOrder _ lc;
                leftMostInserted _ lc;
                firstTime _ FALSE;
                lineScan _ lineScan.clockwiseAroundAbove};
            ENDLOOP;
        leftMostInserted.lineLeft _ leftNeighbor;
        IF leftNeighbor # NIL THEN  leftNeighbor.lineRight _ leftMostInserted;
        };

    CompareLineLeft: PROC [lc: LineCarrier]  ~ {
        workLoad: LIST OF LineCarrier _ LIST[lc];
        newLc: LineCarrier;
        p: Point; s: LineRelation; s1, s2: EndPointRelation;
        ll, lr, topHalf, bottomHalf: Line;
        propagate, insertInQueue: BOOLEAN;
        r: Relation;

        BreakALine: PROC [l: Line] ~ {
            propagate _ TRUE;
            changed _ TRUE;
            topHalf _ l;
            bottomHalf  _ NEW[LineRec _ [above: p, below: l.below, state: Copy[topHalf.state]]];
            RemoveLineFromPointAbove[topHalf];
            RemoveLineFromPointBelow[topHalf];
            topHalf.below _ p;
            InsertLineInPointAbove[topHalf];
            InsertLineInPointBelow[topHalf];

            r _ ComparePoints[bottomHalf.above, bottomHalf.below];
            IF r = equal THEN {
                IF bottomHalf.above = bottomHalf.below THEN insertInQueue _ FALSE
                }
            ELSE {
                IF r = less THEN FlipLine[bottomHalf];
                InsertLineInPointAbove[bottomHalf];
                InsertLineInPointBelow[bottomHalf]
                };
            };

        UNTIL workLoad = NIL DO
            lc _ workLoad.first; workLoad _ workLoad.rest;
            propagate _ FALSE; insertInQueue _ TRUE;
            IF lc = NIL THEN LOOP;
            lr _ lc.carried;
            IF lc.lineLeft = NIL THEN LOOP;
            ll _ lc.lineLeft.carried;

            [p, s, s1, s2] _ PairIntersection[ll, lr];
            IF s = colinear THEN {
                s1 _ s2 _ interior;
                IF ComparePoints[ll.below, lr.below] = less THEN {
                    p _ lr.below;
                    s2 _ below;
                    }
                ELSE {
                    p _ ll.below;
                    s1 _ below;
                    };
                };

            IF (s = intersect) OR (s = colinear) THEN {
                IF ComparePoints[p, currentPoint] = less THEN {
                    IF s2 = interior THEN  BreakALine[lr]
                    ELSE {
                        insertInQueue _ FALSE;
                        IF s2 # below THEN ERROR;  --Remove after debug
                        };
                    IF s1 = interior THEN BreakALine[ll]
                    ELSE {
                        insertInQueue _ FALSE;
                        IF s1 # below THEN ERROR;  --Remove after debug
                        };
                    IF insertInQueue THEN QueueInsert[p];
                    };
                };

            IF s = colinear THEN {
                r _ ComparePoints[ll.above, lr.above];
                IF r = greater THEN {
                    <<Break LL>>
                    RemoveLineFromPointBelow[ll];
                    lr.state _ Combine[lr.state, ll.state];
                    ll.below _ lr.above;
                    InsertLineInPointBelow[ll];
                    lc.lineLeft.carried _ lr 
                    }
                ELSE IF r = less THEN {
                    <<Break LR>>
                    RemoveLineFromPointBelow[lr];
                    ll.state _ Combine[ll.state, lr.state];
                    lr.below _ ll.above;
                    InsertLineInPointBelow[lr];
                    }
                ELSE 
                    {ll.state _ Combine[ll.state, lr.state];
                    RemoveLineFromPointBelow[lr];
                    RemoveLineFromPointAbove[lr];
                    };
                newLc _ lc.lineLeft;
                newLc.lineRight _ lc.lineRight;
                IF lc.lineRight # NIL THEN lc.lineRight.lineLeft _ newLc ELSE xOrder _ newLc;
                lc.carried _ NIL; lc.lineLeft _ lc.lineRight _ NIL;
                workLoad _ CONS[newLc.lineRight, workLoad];
                workLoad _ CONS[newLc, workLoad];
                }
            ELSE {
                IF CompareLinesAtY[ll, lr, currentPoint] = greater THEN
                    {lc.carried _ ll;
                    lc.lineLeft.carried _ lr;
                    propagate _ TRUE};
                IF propagate THEN
                    {workLoad _ CONS[lc.lineRight, workLoad];
                    workLoad _ CONS[lc.lineLeft, workLoad];
                    };
                };
            ENDLOOP;
        };

    currentPoint: Point;

    IF in.status = invalid THEN ERROR;
    IF in.status = planar THEN RETURN[in];

<<Load the points into a priority queue>>

    Heapify[];

<<Process the points in the priority queue>>

changed _ TRUE;
    WHILE changed DO
        changed _ FALSE;
        xOrder _ NIL;
        UNTIL EmptyQueue[] DO
            currentPoint _ NextEvent[];
            neighborRight _ OrderRemoveLines[currentPoint];
            [leftMostInserted, neighborRight] _ OrderInsertLines[currentPoint, neighborRight];
            CompareLineLeft[neighborRight];
            CompareLineLeft[leftMostInserted];
            ENDLOOP;
        switch _ sortedIn; sortedIn _ sortedOut; sortedOut _ switch;
        indexIn _ 0; sortedOut.size _ 0;
        ENDLOOP;

    out _ NEW[GraphRec _ [status: planar, points: sortedIn]];
    in.status _ invalid;

    TrashFreeLC[];
};

END.