[Cyan]<CedarChest6.1>Combiner>ConvexEuclideanGraphsImpl.mesa!1

ConvexEuclideanGraphsImpl.mesa

Last Edited by: Arnon, June 20, 1985 4:37:28 pm PDT

DIRECTORY

Rope,

Convert,

IO,

Imager,

ImagerPath,

RealFns,

List,

RatNums,

EuclideanGraphs,

ConvexEuclideanGraphs;

ConvexEuclideanGraphsImpl: CEDAR PROGRAM

IMPORTS IO, Rope, Convert, ImagerPath, RealFns, List, RatNums, EuclideanGraphs

EXPORTS ConvexEuclideanGraphs

= BEGIN OPEN RN: RatNums, EG: EuclideanGraphs, ConvexEuclideanGraphs;

************** Hull Computation **************

ConvexHull: PUBLIC PROC [L: PointList] RETURNS [Hull: PointList] ~ {

Convex Hull of up to 30 points.

Graham's algorithm is used (Info. Proc. Letters, I, 1972, 132-133). Modification: there is no need to do his Step 4, i.e. to check for points with equal polar coordinates. If there are exactly two such points, they will be removed in the course of step 5. Not clear what may happen if there are three or more such points.

xtemp, ytemp: RN.RatNum;

n: CARDINAL ← 1;

I: CARDINAL;

numLoops: CARDINAL;

nRat: RN.RatNum;

interiorP: Point;

PointArray: TYPE = ARRAY [1..30] OF Point;

pointArray: PointArray;

tempPoint: Point;

flip: BOOL;

in, out: IO.STREAM; -- viewer stream for debugging verbosity

Compactify: PROC [j: CARDINAL] ~ {

Compactify pointArray to remove the jth element

n ← n-1;

FOR jC: CARDINAL IN [j..n] DO

pointArray[jC] ← pointArray[jC+1];

ENDLOOP;

};

MyMOD: PROC [K, N: CARDINAL] RETURNS [CARDINAL] ~ {

IF K MOD N = 0 THEN RETURN[N] ELSE RETURN[K MOD N];

};

[in, out] ← ViewerIO.CreateViewerStreams[IO.PutFR["ConvexHull Log"]];

Load into array

n ← 0;

WHILE L # NIL DO

n ← n + 1;

pointArray[n] ← [ L.first.x, L.first.y ];

out.PutF["\n pointArray[%g] = ( %g , %g )", IO.card[n], IO.rope[RN.RatNumToRatRope[pointArray[n].x, FALSE, FALSE]], IO.rope[RN.RatNumToRatRope[pointArray[n].y, FALSE, FALSE]] ];

L ← L.rest;

ENDLOOP;

Step 1 - Find interior point P - search for noncollinear set of three points, compute com of that triangle

I ← 1;

WHILE I <= n-2 DO

IF EG.PointLineTest[pointArray[I], pointArray[I+1], pointArray[I+2]]#ONLINE THEN EXIT;

I ← I+1;

ENDLOOP;

IF EG.PointLineTest[pointArray[I], pointArray[I+1], pointArray[I+2]]=ONLINE THEN ERROR; -- noncollinear set of three points not found

xtemp ← pointArray[I].x;

ytemp ← pointArray[I].y;

xtemp ← RN.RatNumAdd[xtemp, pointArray[I+1].x];

ytemp ← RN.RatNumAdd[ytemp, pointArray[I+1].y] ;

xtemp ← RN.RatNumAdd[xtemp, pointArray[I+2].x];

ytemp ← RN.RatNumAdd[ytemp, pointArray[I+2].y];

nRat ← RN.RatNumFromSmallCards[1, 3, 1 ];

interiorP ← [x: RN.RatNumDivide[ xtemp , nRat ], y: RN.RatNumDivide[ ytemp , nRat ] ];

out.PutF["\n I, interiorP = %g , ( %g , %g )", IO.card[I], IO.rope[RN.RatNumToRatRope[interiorP.x, FALSE, FALSE]], IO.rope[RN.RatNumToRatRope[interiorP.y, FALSE, FALSE]] ];

Step 2 - Translate points by coordinates of interiorP

FOR I IN [1..n] DO

pointArray[I] ← [x: RN.RatNumSubtract[ pointArray[I].x, interiorP.x], y: RN.RatNumSubtract[ pointArray[I].y , interiorP.y] ];

ENDLOOP;

Step 3 - Bubble sort translated points according to polar angles of vectors from interiorP

- what happens to points equal to interiorP, and multiple occurrences of points?

FOR I DECREASING IN [1..n-1] DO

FOR J:CARDINAL IN [1..I] DO

SELECT RN.RatNumCompare[pointArray[J].y, RN.MakeRatNumZero[]] FROM

ratGreater => SELECT RN.RatNumCompare[pointArray[J+1].y, RN.MakeRatNumZero[]] FROM

ratGreater =>

flip ← RealFns.ArcTan[ RN.RatNumToREAL[pointArray[J].y], RN.RatNumToREAL[pointArray[J].x] ]

RealFns.ArcTan[ RN.RatNumToREAL[pointArray[J+1].y], RN.RatNumToREAL[pointArray[J+1].x] ];

ratLess =>

flip ← FALSE;

ratEqual => SELECT RN.RatNumPositive[pointArray[J+1].x ] FROM

TRUE => flip ← TRUE; -- J+1 on pos x -axis, hence less than J

FALSE => flip ← FALSE; -- J+1 on neg x -axis, hence greater

ENDCASE => ERROR;

ratLess => SELECT RN.RatNumCompare[pointArray[J+1].y, RN.MakeRatNumZero[]] FROM

ratGreater =>

flip ← TRUE;

ratLess =>

flip ← RealFns.ArcTan[ RN.RatNumToREAL[pointArray[J].y], RN.RatNumToREAL[pointArray[J].x] ]

RealFns.ArcTan[ RN.RatNumToREAL[pointArray[J+1].y], RN.RatNumToREAL[pointArray[J+1].x] ];

ratEqual => flip ← TRUE; -- whether on positive or neg x-axis, J+1 less than J

ENDCASE => ERROR;

ratEqual => SELECT RN.RatNumCompare[pointArray[J].x, RN.MakeRatNumZero[] ] FROM

ratGreater => flip ← FALSE; -- J on pos x -axis, hence less than J+1

ratLess => -- J on neg x-axis

SELECT RN.RatNumCompare[pointArray[J+1].y, RN.MakeRatNumZero[]] FROM

ratGreater =>

flip ← TRUE;

ratLess =>

flip ← FALSE;

ratEqual => SELECT RN.RatNumPositive[pointArray[J+1].x ] FROM

TRUE => flip ← TRUE; -- J+1 on pos x -axis, hence less than J

FALSE => flip ← FALSE; -- J+1 on neg x -axis, hence equal

ENDCASE => ERROR;

ratEqual => ERROR; -- J = interiorP; not yet handled

ENDCASE => ERROR;

out.PutF["\n I, J, flip = %g , %g, %g", IO.card[I],IO.card[J], IO.bool[flip] ];

IF flip THEN {

tempPoint ← pointArray[J];

pointArray[J] ← pointArray[J+1];

pointArray[J+1] ← tempPoint;

};

ENDLOOP;

out.PutF["\n New pointArray" ];

FOR I IN [1..n] DO

out.PutF["\n pointArray[%g] = ( %g , %g ), arctan = %g", IO.card[I], IO.rope[RN.RatNumToRatRope[pointArray[I].x, FALSE, FALSE]], IO.rope[RN.RatNumToRatRope[pointArray[I].y, FALSE, FALSE]],
IO.real[RealFns.ArcTan[ RN.RatNumToREAL[pointArray[I].y], RN.RatNumToREAL[pointArray[I].x] ] ] ];

ENDLOOP;

Step 5 - Walk around applying "left of line" test

I ← 1;

numLoops ← 2 * n; -- doesn't matter if we continue walking around longer than necessary (n may decrease in the course of this loop)

FOR J:CARDINAL IN [1.. numLoops] DO

WHILE EG.PointEqual[pointArray[MyMOD[I,n]], pointArray[MyMOD[I+1,n] ]] DO

Compactify[MyMOD[I+1,n]]; -- check for and remove duplicate points

ENDLOOP;

WHILE EG.PointLineTest[pointArray[MyMOD[I,n]], pointArray[MyMOD[I+1,n]], pointArray[MyMOD[I+2,n]] ] = RIGHTOFLINE DO

Compactify[MyMOD[I+1,n]];

I ← MyMOD[I-1,n];

ENDLOOP;

I ← MyMOD[I+1,n];

ENDLOOP;

Return untranslated points of hull.

Hull ← NIL;

FOR I DECREASING IN [1..n] DO

Hull ← CONS [ [x: RN.RatNumAdd[ pointArray[I].x, interiorP.x], y: RN.RatNumAdd[ pointArray[I].y, interiorP.y] ], Hull];

ENDLOOP;

};

************** Graph Structure **************

ThreeOrMoreVertices: PUBLIC PROC [w: EG.EuclideanGraph] RETURNS [BOOL] = {

hasThreeVertices: BOOL;

hasThreeVertices ← w # NIL; -- test for 0 vertices

IF hasThreeVertices THEN hasThreeVertices ← w.adjacentVertices # NIL; -- test for 1 vertex only

IF hasThreeVertices THEN hasThreeVertices ← w.adjacentVertices.first.vertex # w.adjacentVertices.rest.first.vertex; -- test for 2 vertices only

RETURN[ hasThreeVertices ];

};

AddEdge: PUBLIC PROC [root, adjacentVertex: Vertex, exterior: BOOL ← FALSE, IPEdge: BOOL ← FALSE, clientData: REF ← NIL, visitedValue: BOOL ← FALSE] ~ {

cEGEdgeData: CEGEdgeData ← NEW[CEGEdgeDataRec ← [exteriorEdge: exterior, intersectionPolygonEdge: IPEdge, clientData: clientData]];

EG.AddEdge[root, adjacentVertex, cEGEdgeData, visitedValue];

};

InsertVertexInEdges: PUBLIC PROC [v, w: Vertex, p: Point, vertexData: REF ← NIL] RETURNS [newVertex: Vertex] ~ {

vTow, wTov: Adjacency;

vTowData, wTovData: CEGEdgeData;

z: Vertex ← EG.CreateVertex[p, vertexData, v.visited];

[vTow, wTov] ← EG.FindAdjacency[v, w];

vTowData ← NARROW[vTow.data];

wTovData ← NARROW[wTov.data];

AddEdge[z, w, vTowData.exteriorEdge, vTowData.intersectionPolygonEdge, vTowData.clientData, vTow.visited];

EG.DeleteEdge[v, w];

AddEdge[v, z, vTowData.exteriorEdge, vTowData.intersectionPolygonEdge, vTowData.clientData, vTow.visited];

AddEdge[z, v, wTovData.exteriorEdge, wTovData.intersectionPolygonEdge, wTovData.clientData, wTov.visited];

EG.DeleteEdge[w, v];

AddEdge[w, z, wTovData.exteriorEdge, wTovData.intersectionPolygonEdge, wTovData.clientData, wTov.visited];

RETURN[z];

};

JoinVertices: PUBLIC PROC [start, end: Vertex, startToEndData, endToStartData: REF ← NIL, startToEndIPEdge, endToStartIPEdge: BOOL ← FALSE, startToEndExterior, endToStartExterior: BOOL ← FALSE, startToEndVisited, endToStartVisited: BOOL ← FALSE] ~ {

startToEndCEGData: CEGEdgeData ← NEW[CEGEdgeDataRec ← [exteriorEdge: startToEndExterior, intersectionPolygonEdge: startToEndIPEdge, clientData: startToEndData]];

endToStartCEGData: CEGEdgeData ← NEW[CEGEdgeDataRec ← [exteriorEdge: endToStartExterior, intersectionPolygonEdge: endToStartIPEdge, clientData: endToStartData]];

EG.AddEdge[start, end, startToEndCEGData, startToEndVisited];

EG.AddEdge[end, start, endToStartCEGData, endToStartVisited];

};

************** Graph Traversal **************

StepVertex: PUBLIC PROC[ v1, w1: Vertex, outline: BOOL] RETURNS [ v2, w2: Vertex] = {

RETURN [w1, NextVertex[ v1, w1, outline] ]

};

NextVertex: PUBLIC PROC[ v, w: Vertex, outline: BOOL] RETURNS [ z: Vertex] = {

IF outline THEN z ← NextOutlineVertex[v,w] ELSE z ← NextPolygonVertex[v,w];

RETURN [z]

};

PreviousVertex: PUBLIC PROC[ v, w: Vertex, outline: BOOL] RETURNS [ z: Vertex] = {

IF outline THEN z ← PreviousOutlineVertex[v,w] ELSE z ← PreviousPolygonVertex[v,w];

RETURN [z]

};

NextOutlineVertex: PUBLIC PROC[ v, w: Vertex] RETURNS [ z: Vertex] = {

A: AdjacencyList ← w.adjacentVertices;

z ← A.first.vertex; A ← A.rest;

WHILE A.first.vertex # v DO

z ← A.first.vertex; A ← A.rest;

ENDLOOP;

RETURN [z]

};

NextPolygonVertex: PUBLIC PROC[ v, w: Vertex] RETURNS [ Vertex] = {

A: AdjacencyList ← w.adjacentVertices;

WHILE A.first.vertex # v DO

A ← A.rest

ENDLOOP;

RETURN [A.rest.first.vertex]

};

SpecialPreviousOutlineVertex: PUBLIC PROC [w: Vertex] RETURNS [v:Vertex] ~ {

Assuming that w is a vertex on the outline of some Figure containing at least two vertices, find the vertex v which precedes w in the counterclockwise ordering of vertices (we think of this counterclockwise ordering as being defined with respect to the INTERIOR of the figure on whose outline w lies).

ahl: AdjacencyList ← w.adjacentVertices;

WHILE NOT GetEdgeExterior[ahl.first] DO ahl ← ahl.rest ENDLOOP;

v ← ahl.first.vertex;

};

PreviousOutlineVertex: PUBLIC PROC [v, w: Vertex] RETURNS [z: Vertex] ~ {

Assuming that [v,w] is an edge on the outline of some Figure containing at least two vertices, find the vertex z which precedes w in the counterclockwise ordering of vertices (we think of this counterclockwise ordering as being defined with respect to the INTERIOR of the figure on whose outline w lies).

A: AdjacencyList ← v.adjacentVertices;

WHILE A.first.vertex # w DO

A ← A.rest

ENDLOOP;

RETURN [A.rest.first.vertex]

};

PreviousPolygonVertex: PUBLIC PROC [v, w: Vertex] RETURNS [z:Vertex] ~ {

[v, w] is an edge on a polygon containing at least two vertices. Find the vertex z which precedes v in the counterclockwise ordering of vertices (we think of this counterclockwise ordering as being defined with respect to the INTERIOR of the polygon on which [v, w] lies).

A: AdjacencyList ← v.adjacentVertices;

z ← A.first.vertex; A ← A.rest;

WHILE A.first.vertex # w DO

z ← A.first.vertex; A ← A.rest;

ENDLOOP;

RETURN [z]

};

************** Edge Fields **************

GetEdgeClientData: PUBLIC PROC [startToEnd: Adjacency] RETURNS [data: REF ANY] = {

cEGEdgeData: CEGEdgeData ← NARROW[startToEnd.data];

RETURN[cEGEdgeData.clientData];

};

SetEdgeClientData: PUBLIC PROC [startToEnd: Adjacency, data: REF ANY] = {

cEGEdgeData: CEGEdgeData ← NARROW[startToEnd.data];

cEGEdgeData.clientData ← data;

};

SetPolygonClientData: PUBLIC PROC [start, end: Vertex, data: REF, setVisited: BOOL ← FALSE] = {

startToEnd, endToStart: Adjacency;

firstStart: Vertex ← start;

[ startToEnd, endToStart ] ← EG.FindAdjacency[ start, end ];

SetEdgeClientData[startToEnd, data];

IF setVisited THEN startToEnd.visited ← TRUE;

WHILE end # firstStart DO

[start, end] ← StepVertex[start, end, FALSE];

[ startToEnd, endToStart ] ← EG.FindAdjacency[ start, end ];

SetEdgeClientData[startToEnd, data];

IF setVisited THEN startToEnd.visited ← TRUE;

ENDLOOP;

};

GetEdgeExterior: PUBLIC PROC [startToEnd: Adjacency] RETURNS [value: BOOL] ~ {

cEGEdgeData: CEGEdgeData ← NARROW[startToEnd.data];

RETURN[cEGEdgeData.exteriorEdge];

};

SetEdgeExterior: PUBLIC PROC [startToEnd: Adjacency, value: BOOL] ~ {

cEGEdgeData: CEGEdgeData ← NARROW[startToEnd.data];

cEGEdgeData.exteriorEdge ← value;

};

GetEdgeIPEdge: PUBLIC PROC [startToEnd: Adjacency] RETURNS [value: BOOL] ~ {

cEGEdgeData: CEGEdgeData ← NARROW[startToEnd.data];

RETURN[cEGEdgeData.intersectionPolygonEdge];

};

SetEdgeIPEdge: PUBLIC PROC [startToEnd: Adjacency, value: BOOL] ~ {

cEGEdgeData: CEGEdgeData ← NARROW[startToEnd.data];

cEGEdgeData.intersectionPolygonEdge ← value;

};

SetPolygonIPEdges: PUBLIC PROC [start, end: Vertex, value: BOOL] ~ {

startToEnd, endToStart: Adjacency;

firstStart: Vertex ← start;

[ startToEnd, endToStart ] ← EG.FindAdjacency[ start, end ];

SetEdgeIPEdge[startToEnd, value];

WHILE end # firstStart DO

[start, end] ← StepVertex[start,end, FALSE]; -- polygon step

[ startToEnd, endToStart ] ← EG.FindAdjacency[ start, end ];

SetEdgeIPEdge[startToEnd, value];

ENDLOOP;

};

ClearGraphIPEdges: PUBLIC PROC[v, w: Vertex] = {

startToEnd, endToStart: Adjacency;

firstIntersection: Vertex ← v;

outline: BOOL ← FALSE;

Clear visited and intersectionPolygonEdge fields of edges of this polygon

[ startToEnd, endToStart ] ← EG.FindAdjacency[ v, w ];

startToEnd.visited ← FALSE;

SetEdgeIPEdge[startToEnd, FALSE];

WHILE w # firstIntersection DO

[v, w] ← StepVertex[v,w, outline];

[ startToEnd, endToStart ] ← EG.FindAdjacency[ v, w ];

startToEnd.visited ← FALSE;

SetEdgeIPEdge[startToEnd, FALSE];

ENDLOOP;

Go around edges again; check each to see if an unprocessed polygon lies across it, if so, recursive call on the flipped edge.

[v, w] ← StepVertex[v,w, outline];

[ startToEnd, endToStart ] ← EG.FindAdjacency[ v, w ];

IF endToStart.visited THEN ClearGraphIPEdges[w, v];

WHILE w # firstIntersection DO

[v, w] ← StepVertex[v,w, outline];

[ startToEnd, endToStart ] ← EG.FindAdjacency[ v, w ];

IF endToStart.visited THEN ClearGraphIPEdges[w, v];

ENDLOOP;

};

IPVertexToEdge: PUBLIC PROC [intersectionPolyVertex: Vertex] RETURNS [intersectionPolyNext: Vertex] ~ {

adjacencies: AdjacencyList ← intersectionPolyVertex.adjacentVertices;

WHILE NOT GetEdgeIPEdge[adjacencies.first] DO

adjacencies ← adjacencies.rest;

ENDLOOP;

RETURN[adjacencies.first.vertex];

};

IPVertex: PUBLIC PROC [v: Vertex] RETURNS [BOOL] ~ {

adjacencies: AdjacencyList ← v.adjacentVertices;

adjHead: AdjacencyList ← v.adjacentVertices;

IF adjacencies = NIL THEN RETURN[FALSE];

IF GetEdgeIPEdge[adjacencies.first] THEN RETURN[TRUE];

adjacencies ← adjacencies.rest;

WHILE adjacencies # adjHead DO

IF GetEdgeIPEdge[adjacencies.first] THEN RETURN[TRUE];

adjacencies← adjacencies.rest;

ENDLOOP;

RETURN[FALSE];

};

ClearEdgeVisitedFields: PUBLIC PROC [v, w: Vertex] = {

[v,w] is a counterclockwise oriented edge of a polygon contained in the outline of the current figure. Clear the visited fields of all edges of all polygons interior to a Vertex, assumed by this particular ClearEdgeVisitedFields to be an Vertex. Uses Depth First Search.

vTow, wTov: Adjacency;

firstv: Vertex ← v;

outline: BOOL ← FALSE;

Clear visited fields of edges of this polygon

[ vTow, wTov ] ← EG.FindAdjacency[ v, w ];

vTow.visited ← FALSE;

WHILE w # firstv DO

[v, w] ← StepVertex[v,w, outline];

[ vTow, wTov ] ← EG.FindAdjacency[ v, w ];

vTow.visited ← FALSE;

ENDLOOP;

Go around edges again; check (opposite orientation of) each to see if it's an interior edge (i.e. not on the outline), and if the polygon to which it belongs is uncleared, if so, recursive call on that polygon.

[v, w] ← StepVertex[v, w, outline];

[ vTow, wTov ] ← EG.FindAdjacency[ v, w ];

IF wTov.visited THEN IF NOT GetEdgeExterior[wTov] THEN

ClearEdgeVisitedFields[w, v ]

ELSE

wTov.visited ← FALSE;

WHILE w # firstv DO

[v, w] ← StepVertex[v,w, outline];

[ vTow, wTov ] ← EG.FindAdjacency[ v, w ];

IF wTov.visited THEN IF NOT GetEdgeExterior[wTov] THEN

ClearEdgeVisitedFields[w, v]

ELSE

wTov.visited ← FALSE;

ENDLOOP;

};

ClearGraphEdgeFields: PUBLIC PROC [v: Vertex] ~ {

Clears IPEdge and visited fields of all edges in the Graph pointed to by v.

Subr: PROC [v: Vertex, visitedValue: BOOL] ~ {

vAdjList: AdjacencyList ← v.adjacentVertices;

firstVertex, nextVertex: Vertex;

vToNext, nextToV: Adjacency;

Go around adjacencies clearing fields.

firstVertex ← nextVertex ← vAdjList.first.vertex;

[vToNext, nextToV] ← EG.FindAdjacency[v, nextVertex] ;

vToNext.visited ← FALSE;

SetEdgeIPEdge[vToNext, FALSE];

vAdjList ← vAdjList.rest;

nextVertex ← vAdjList.first.vertex;

WHILE nextVertex#firstVertex DO

[vToNext, nextToV] ← EG.FindAdjacency[v, nextVertex] ;

vToNext.visited ← FALSE;

SetEdgeIPEdge[vToNext, FALSE];

vAdjList ← vAdjList.rest;

nextVertex ← vAdjList.first.vertex;

ENDLOOP;

v.visited ← visitedValue; -- now we can say that we've been here.

Go around adjacencies again, making recursive calls for unvisited vertices.

IF (nextVertex.visited # visitedValue) THEN Subr[nextVertex, visitedValue];

vAdjList ← vAdjList.rest;

nextVertex ← vAdjList.first.vertex;

WHILE nextVertex#firstVertex DO

IF (nextVertex.visited # visitedValue) THEN Subr[nextVertex, visitedValue];

vAdjList ← vAdjList.rest;

nextVertex ← vAdjList.first.vertex;

ENDLOOP;

};

IF v = NIL OR v.adjacentVertices = NIL THEN RETURN ELSE Subr[v, NOT v.visited];

};

SetOutwardOutlineEdgeFields: PUBLIC PROC [start, end: Vertex, exteriorEdgeValue, intersectionPolygonEdgeValue: BOOL, clientDataValue: REF ← NIL] ~ {

cEGEdgeData: CEGEdgeData;

startToEnd, endToStart: Adjacency;

firstStart: Vertex ← start;

[ startToEnd, endToStart ] ← EG.FindAdjacency[ start, end ];

cEGEdgeData ← NARROW[endToStart.data];

cEGEdgeData.exteriorEdge ← exteriorEdgeValue;

cEGEdgeData.intersectionPolygonEdge ← intersectionPolygonEdgeValue;

cEGEdgeData.clientData ← clientDataValue;

WHILE end # firstStart DO

[start, end] ← StepVertex[start,end, TRUE]; -- outline step

[ startToEnd, endToStart ] ← EG.FindAdjacency[ start, end ];

cEGEdgeData ← NARROW[endToStart.data];

cEGEdgeData.exteriorEdge ← exteriorEdgeValue;

cEGEdgeData.intersectionPolygonEdge ← intersectionPolygonEdgeValue;

cEGEdgeData.clientData ← clientDataValue;

ENDLOOP;

};

SetPolygonFields: PUBLIC PROC [start, end: Vertex, exteriorEdgeValue, intersectionPolygonEdgeValue: BOOL ← FALSE, clientDataValue: REF ← NIL] ~ {

cEGEdgeData: CEGEdgeData;

startToEnd, endToStart: Adjacency;

firstStart: Vertex ← start;

[ startToEnd, endToStart ] ← EG.FindAdjacency[start, end];

cEGEdgeData ← NARROW[startToEnd.data];

cEGEdgeData.exteriorEdge ← exteriorEdgeValue;

cEGEdgeData.intersectionPolygonEdge ← intersectionPolygonEdgeValue;

cEGEdgeData.clientData ← clientDataValue;

WHILE end # firstStart DO

[start, end] ← StepVertex[start, end, FALSE]; -- polygon step

[ startToEnd, endToStart ] ← EG.FindAdjacency[ start, end ];

cEGEdgeData ← NARROW[startToEnd.data];

cEGEdgeData.exteriorEdge ← exteriorEdgeValue;

cEGEdgeData.intersectionPolygonEdge ← intersectionPolygonEdgeValue;

cEGEdgeData.clientData ← clientDataValue;

ENDLOOP;

};

************** Convex (internal) polygon operations **************

CenterOfMass: PUBLIC PROC[ v, w : Vertex, outline: BOOL ] RETURNS [Point] = {

[v,w] is a counterclockwise oriented edge on the polygon

vfirst: Vertex ← v;

xt: RN.RatNum ← v.coordinates.x;

yt: RN.RatNum ← v.coordinates.y;

n: INTEGER ← 1;

nRat: RN.RatNum;

WHILE w # vfirst DO

xt ← RN.RatNumAdd[ xt , w.coordinates.x ];

yt ← RN.RatNumAdd[ yt , w.coordinates.y ];

n ← n + 1;

[v,w] ← StepVertex[ v, w, outline ];

ENDLOOP;

nRat ← RN.RatNumFromSmallCards[1, n, 1 ];

RETURN[ [x: RN.RatNumDivide[ xt , nRat ], y: RN.RatNumDivide[ yt , nRat ] ] ];

};

CountPolygonVertices: PUBLIC PROC [start, end: Vertex, outline: BOOL ← FALSE] RETURNS [numVertices: CARDINAL] ~ {

firstStart: Vertex ← start;

numVertices ← 1;

WHILE end # firstStart DO

numVertices ← numVertices + 1;

[start, end] ← StepVertex[start, end, outline];

ENDLOOP;

};

AddVertexToPolygon: PUBLIC PROC [entryPrior, entryNext: Vertex, coordinates: Point, clientData: REF ← NIL] RETURNS [exitPrior, exitNext: Vertex] = {

We assume that all vertices, and all edges, in the previously existing polygon have the same visited values, which are propagated to vertices and edges in the new structure.

newVertex: Vertex ← EG.CreateVertex[coordinates, NIL, FALSE];

numVertices: CARDINAL;

priorToNext, nextToPrior: Adjacency;

edgeVisitedValue: BOOL;

IF entryNext = NIL THEN RETURN[NIL, newVertex]; -- no vertices previously

newVertex.visited ← entryNext.visited;

IF entryPrior = NIL THEN { -- one vertex previously

AddEdge[newVertex, entryNext, FALSE, FALSE, clientData, FALSE];

AddEdge[entryNext, newVertex, TRUE, FALSE, NIL, FALSE];

RETURN[newVertex, entryNext];

};

[priorToNext, nextToPrior] ← EG.FindAdjacency[entryPrior, entryNext];

edgeVisitedValue ← priorToNext.visited;

numVertices ← CountPolygonVertices[entryPrior, entryNext]; -- do before SeparateVertices

EG.SeparateVertices[entryPrior, entryNext];

AddEdge[entryPrior, newVertex, FALSE, FALSE, clientData, edgeVisitedValue];

AddEdge[newVertex, entryPrior, TRUE, FALSE, NIL, edgeVisitedValue];

AddEdge[newVertex, entryNext, FALSE, FALSE, clientData, edgeVisitedValue];

AddEdge[entryNext, newVertex, TRUE, FALSE, NIL, edgeVisitedValue];

IF numVertices = 2 THEN {

AddEdge[entryNext, entryPrior, FALSE, FALSE, clientData, edgeVisitedValue];

AddEdge[entryPrior, entryNext, TRUE, FALSE, NIL, edgeVisitedValue];

};

RETURN[newVertex, entryNext];

};

PointSectorTest: PUBLIC PROC [rightvertex, leftvertex, centerOfMass, test: Point] RETURNS[PointSectorRelation ← EQUALCOM] = {

Determine if the test point lies in the sector of a polygon determined by leftvertex, rightvertex, centerOfMass A sector is defined as the (interior of the) triangle having as vertices the center of mass of p, the ith vertex of p (rightvertex as you face out from the center of mass) and the i+1st vertex of p (leftvertex).

pointEdgeRelation: EG.PointDirectedEdgeRelation ← EG.PointDirectedEdgeTest[ centerOfMass, rightvertex, test];

SELECT pointEdgeRelation FROM

BEFORESTART => RETURN[ OUTSIDESUBTEND ];

EQUALSTART => RETURN[ EQUALCOM ];

BETWEENSTARTEND => RETURN[ RIGHTSPOKEINTERIOR ];

EQUALEND => RETURN[ EQUALRIGHTVERTEX ];

AFTEREND => RETURN[ RIGHTRAYINTERIOR ];

RIGHTOFEDGE => RETURN[ OUTSIDESUBTEND ];

LEFTOFEDGE => {

pointEdgeRelation ← EG.PointDirectedEdgeTest[ centerOfMass, leftvertex, test];

SELECT pointEdgeRelation FROM

BEFORESTART => RETURN[ OUTSIDESUBTEND ];

BETWEENSTARTEND => RETURN[ LEFTSPOKEINTERIOR ];

EQUALEND => RETURN[ EQUALLEFTVERTEX ];

AFTEREND => RETURN[ LEFTRAYINTERIOR ];

LEFTOFEDGE => RETURN[ OUTSIDESUBTEND ];

RIGHTOFEDGE => {

pointLineRelation: EG.PointLineRelation ← EG.PointLineTest[ rightvertex, leftvertex, test];

SELECT pointLineRelation FROM

ONLINE => RETURN[POLYGONEDGEINTERIOR];

LEFTOFLINE => RETURN[SECTORINTERIOR];

RIGHTOFLINE => RETURN[SUBTENDINTERIOR];

ENDCASE;

};

ENDCASE => ERROR; -- EQUALSTART shouldn't occur

};

ENDCASE;

};

FindSubtendForPoint: PUBLIC PROC[v, w: Vertex, centerOfMass, test: Point, outline: BOOL ] RETURNS [status: PointSectorRelation, t, u: Vertex] = {

[v,w] is a counterclockwise oriented edge on the polygon. [t, u] is a counterclockwise oriented edge on the polygon such that test lies in the closure of the subtend defined by [t, u] and the polygon's center of mass. If test = centerOfMass, then status = EQUALCOM and [t, u] ← [v, w].

pointSectorRelation: PointSectorRelation ← PointSectorTest[ v.coordinates, w.coordinates, centerOfMass, test];

WHILE pointSectorRelation = OUTSIDESUBTEND DO

[v, w] ← StepVertex [ v, w, outline];

pointSectorRelation ← PointSectorTest[ v.coordinates, w.coordinates, centerOfMass, test];

ENDLOOP;

RETURN[ pointSectorRelation, v, w ];

};

VertexExternalToPolygon: PUBLIC PROC [firstCOM: Point, firstStart, firstEnd, secondStart, secondEnd: Vertex] RETURNS [found: BOOL, newFirstStart, newFirstEnd, newSecondStart, newSecondEnd: Vertex ← NIL] ~ {

pointSectorStatus: PointSectorRelation;

initSecondVertex: Vertex ← secondStart;

done: BOOL ← FALSE;

[pointSectorStatus, newFirstStart, newFirstEnd] ← FindSubtendForPoint[ firstStart, firstEnd, firstCOM, secondStart.coordinates, FALSE ]; -- [v5,w5] defines a particular sector of this polygon

WHILE NOT done DO

IF pointSectorStatus = RIGHTRAYINTERIOR OR pointSectorStatus = LEFTRAYINTERIOR OR pointSectorStatus = SUBTENDINTERIOR THEN

RETURN[TRUE, newFirstStart, newFirstEnd, secondStart, secondEnd];

[secondStart, secondEnd] ← StepVertex[secondStart, secondEnd, FALSE];

IF secondStart = initSecondVertex THEN

done ← TRUE

ELSE

[pointSectorStatus, newFirstStart, newFirstEnd] ← FindSubtendForPoint[newFirstStart, newFirstEnd, firstCOM, secondStart.coordinates, FALSE ];

ENDLOOP;

RETURN[FALSE, newFirstStart, newFirstEnd, newSecondStart, newSecondEnd];

};

ConvexPolygon: PUBLIC PROC[ v, w : Vertex, outline: BOOL ] RETURNS [BOOL] = {

[v,w] is a counterclockwise oriented edge on the polygon, must have at least 2 (3?) vertices.

vfirst: Vertex ← v;

z: Vertex ← NextVertex[v, w, outline];

IF z = v OR z = w THEN RETURN[TRUE];

SELECT EG.PointLineTest[v.coordinates, w.coordinates, z.coordinates] FROM -- check first angle

RIGHTOFLINE => RETURN[FALSE];

LEFTOFLINE => NULL;

ONLINE => NULL; -- 1/7/86 allow a 180deg angle

ENDCASE;

v ← w; [w, z] ← StepVertex[ w, z, outline ]; -- step past

WHILE v # vfirst DO

SELECT EG.PointLineTest[v.coordinates, w.coordinates, z.coordinates] FROM -- check first angle

RIGHTOFLINE => RETURN[FALSE];

LEFTOFLINE => NULL;

ONLINE => NULL; -- 1/7/86 allow a 180deg angle

ENDCASE;

v ← w; [w, z] ← StepVertex[ w, z, outline ];

ENDLOOP;

RETURN[TRUE];

};

ConvexPolygonCompare: PUBLIC PROC[ innerCOM: Point, innerStart, innerEnd: Vertex, innerOutline: BOOL, outerStart, outerEnd: Vertex, checkExternal: BOOL ← TRUE] RETURNS [status: ConvexPolygonCompareStatus, innerPrevious, innerCommon, innerNext, outerPrevious, outerCommon, outerNext: Vertex ← NIL] = {

external: BOOL ← FALSE;

innerStartFirst: Vertex ← innerStart;

outerStartFirst: Vertex ← outerStart;

pointSectorRelation: PointSectorRelation;

pointLineRelation: EG.PointLineRelation;

edgeCompareStatus: EG.BiasedEdgeCompareStatus;

intersection: Point;

polygonStep: BOOL ← FALSE;

done: BOOL ← FALSE;

CheckEdgeContact: PROC ~ INLINE {

Determine which side of [innerStart, innerEnd] outerStart is on.

pointLineRelation ← EG.PointLineTest[innerStart.coordinates, innerEnd.coordinates, outerStart.coordinates];

SELECT pointLineRelation FROM

ONLINE => edgeCompareStatus ← Disjoint; -- we assume outerStart is disjoint from the edge.

LEFTOFLINE => edgeCompareStatus ← Disjoint; --

General claim: assuming that any edge of outer polygon that contacts an edge of inner polygon starts at an outerStart that is outside the inner polygon, then that outerStart must lie to the right of the inner edge that is contacted. Proof: outerStart is external to inner polygon, hence to the right of line of support through any inner edge contacted by any outer edge starting at outerStart.

RIGHTOFLINE => {

edgeCompareStatus ← EG.BiasedEdgeCompare[outerStart, outerEnd, innerStart, innerEnd];

SELECT edgeCompareStatus FROM

OuterEndEqInnerStart => {

outerCommon ← outerEnd;

outerEnd ← NextVertex[outerStart, outerCommon, polygonStep];

innerCommon ← innerStart;

innerStart ← PreviousVertex[innerCommon, innerEnd, innerOutline];

};

OuterEndInInner => {

outerCommon ← outerEnd;

innerCommon ← InsertVertexInEdges[innerStart, innerEnd, outerEnd.coordinates];

outerEnd ← NextVertex[outerStart, outerCommon, polygonStep]; -- this statement needs to come after previous one

};

OuterEndEqInnerEnd => {

outerCommon ← outerEnd;

outerEnd ← NextVertex[outerStart, outerCommon, polygonStep];

innerCommon ← innerEnd;

innerEnd ← NextVertex[innerStart, innerCommon, innerOutline];

};

InnerStartInOuter => {

innerCommon ← innerStart;

outerCommon ← InsertVertexInEdges[outerStart, outerEnd, innerStart.coordinates];

innerStart ← PreviousVertex[innerCommon, innerEnd, innerOutline]; -- this statement needs to come after previous one

};

ProperIntersection => {

intersection ← EG.IntersectLines[innerStart.coordinates, innerEnd.coordinates, outerStart.coordinates, outerEnd.coordinates];

innerCommon ← InsertVertexInEdges[innerStart, innerEnd, intersection];

outerCommon ← InsertVertexInEdges[outerStart, outerEnd, intersection];

};

InnerEndInOuter => {

innerCommon ← innerEnd;

outerCommon ← InsertVertexInEdges[outerStart, outerEnd, innerEnd.coordinates];

innerEnd ← NextVertex[innerStart, innerCommon, innerOutline]; -- this statement needs to come after previous one

};

Disjoint => NULL; -- this one falls through

ENDCASE;

};

ENDCASE;

}; -- CheckEdgeContact

pointSectorRelation ← PointSectorTest[ innerStart.coordinates, innerEnd.coordinates, innerCOM, outerEnd.coordinates ];

WHILE NOT done DO

SELECT pointSectorRelation FROM

RIGHTRAYINTERIOR, LEFTRAYINTERIOR, SUBTENDINTERIOR => NULL;

ENDCASE => { -- outerEnd either inside (closed) inner polygon, or R or L of this sector

CheckEdgeContact[];

SELECT edgeCompareStatus FROM

Disjoint => NULL;

ENDCASE => RETURN[Intersection, innerStart, innerCommon, innerEnd, outerStart, outerCommon, outerEnd];

WHILE pointSectorRelation = OUTSIDESUBTEND DO

[innerStart, innerEnd] ← StepVertex[innerStart, innerEnd, innerOutline]; -- step inner polygon to catch up with outer polygon edge

CheckEdgeContact[];

SELECT edgeCompareStatus FROM

Disjoint => NULL;

ENDCASE => RETURN[Intersection, innerStart, innerCommon, innerEnd, outerStart, outerCommon, outerEnd];

pointSectorRelation ← PointSectorTest[ innerStart.coordinates, innerEnd.coordinates, innerCOM, outerEnd.coordinates ];

ENDLOOP;

};

If we arrive here, then outerEnd is in RIGHTRAYINTERIOR, LEFTRAYINTERIOR, SUBTENDINTERIOR of current subtend of inner polygon, and there is no contact of [outerStart, outerEnd] with inner polygon.

IF checkExternal THEN IF EG.PointLineTest[ outerStart.coordinates, outerEnd.coordinates, innerStartFirst.coordinates] = RIGHTOFLINE THEN external ← TRUE; -- if inner polygon is external to outer polygon, then any vertex of the inner polygon (we arbitrarily picked innerStartFirst) must lie 'external' to at least one edge of the outer polygon. This test will eventually uncover such an outer edge if one exists.

[outerStart, outerEnd] ← StepVertex[outerStart, outerEnd, polygonStep];

IF outerStart = outerStartFirst THEN

done ← TRUE

ELSE

pointSectorRelation ← PointSectorTest[ innerStart.coordinates, innerEnd.coordinates, innerCOM, outerEnd.coordinates ];

ENDLOOP;

IF external THEN

RETURN[External, innerStart, innerCommon, innerEnd, outerStart, outerCommon, outerEnd]

ELSE

RETURN [Encloses, innerStart, innerCommon, innerEnd, outerStart, outerCommon, outerEnd];

};

************** Graph Searching **************

FindInternalPolygonForPoint: PUBLIC PROC [v, w: Vertex, test: Point] RETURNS [pointSectorStatus: PointSectorRelation, y, z: Vertex] ~ {

polygonCOM: Point;

v5: Vertex ← v;

w5: Vertex ← w;

temp: Vertex;

polygonCOM ← CenterOfMass[v5, w5, FALSE]; -- COM of current internal polygon.

[pointSectorStatus, v5, w5] ← FindSubtendForPoint[ v5, w5, polygonCOM, test, FALSE];

WHILE pointSectorStatus = RIGHTRAYINTERIOR OR pointSectorStatus = LEFTRAYINTERIOR OR pointSectorStatus = SUBTENDINTERIOR DO -- flip flop polygons until you find the one that contains test

temp ← v5; v5 ← w5; w5 ← temp; -- flip v5 and w5 to move to next polygon

polygonCOM ← CenterOfMass[ v5, w5, FALSE];

[pointSectorStatus, v5, w5] ← FindSubtendForPoint[ v5, w5, polygonCOM, test, FALSE ];

ENDLOOP;

RETURN[ pointSectorStatus, v5, w5];

};

SearchInternalPolygonsForPoint: PUBLIC PROC [v, w: Vertex, test: Point] RETURNS [found: BOOL, y, z: Vertex] ~ {

polygonCOM: Point;

pointSectorStatus: PointSectorRelation;

v5: Vertex ← v;

w5: Vertex ← w;

temp: Vertex;

v5Tow5, w5Tov5: Adjacency;

[v5Tow5, w5Tov5] ← EG.FindAdjacency[v5, w5];

IF GetEdgeExterior[v5Tow5] THEN ERROR;

polygonCOM ← CenterOfMass[ v5, w5, FALSE]; -- COM of current internal polygon.

[pointSectorStatus, v5, w5] ← FindSubtendForPoint[ v5, w5, polygonCOM, test, FALSE ]; -- [v5,w5] defines a particular sector of this polygon

WHILE pointSectorStatus = RIGHTRAYINTERIOR OR pointSectorStatus = LEFTRAYINTERIOR OR pointSectorStatus = SUBTENDINTERIOR DO -- flip flop polygons until you find the one that contains test

temp ← v5; v5 ← w5; w5 ← temp; -- flip v5 and w5 to move to next polygon

[ v5Tow5, w5Tov5] ← EG.FindAdjacency[ v5, w5];

IF GetEdgeExterior[v5Tow5] THEN RETURN[FALSE, v, w];

polygonCOM ← CenterOfMass[ v5, w5, FALSE];

[pointSectorStatus, v5, w5] ← FindSubtendForPoint[ v5, w5, polygonCOM, test, FALSE ];

ENDLOOP;

RETURN[TRUE, v5, w5];

};

************** Graph Outline **************

ConvexPolygonOnOutline: PUBLIC PROC[prew2, w2, prew1, w1: Vertex] RETURNS [BOOL] = {

We assume that the vertices are given to us in counterclockwise order (with respect to figure interior) along the outline of the figure.

x, y, z: Vertex;

SELECT EG.PointLineTest[w2.coordinates, prew2.coordinates, w1.coordinates] FROM

RIGHTOFLINE => RETURN[FALSE];

LEFTOFLINE => NULL;

ONLINE => NULL; -- 180deg angle ok

ENDCASE;

SELECT EG.PointLineTest[prew2.coordinates, w1.coordinates, prew1.coordinates] FROM

RIGHTOFLINE => RETURN[FALSE];

LEFTOFLINE => NULL;

ONLINE => NULL;

ENDCASE;

x ← w1; y ← prew1; z ← PreviousOutlineVertex[y, x];

WHILE x # w2 DO

SELECT EG.PointLineTest[x.coordinates, y.coordinates, z.coordinates] FROM

RIGHTOFLINE => RETURN[FALSE];

LEFTOFLINE => NULL;

ONLINE => NULL;

ENDCASE;

x ← y; y ← z; z ← PreviousOutlineVertex[y, x];

ENDLOOP;

RETURN[TRUE];

};

GrowToConvexOutline: PUBLIC PROC[w2, v, w1: Vertex, verbosity: Verbosity ← NIL] RETURNS [v3, w3: Vertex] = {

w2, v, and w1 are assumed to be in counterclockwise order (with respect to the interior of the outline) on an outline. [v3,w3] is a counterclockwise-oriented (with respect to the interior of the figure) edge on the new outline.

It is assumed that [w2, v, w1] is a nontrivial triangle, i.e. that the outline is not in the neighborhood of [w2, v, w1] have a convex outline, in particular that w1 is to the right of the line from w2 through v.

The triangle [w2, v, w1] is the initial convex polygon. WHILE the figure does not (in the neighborhood of w2, v, w1) have a convex outline, the algorithm enlarges this polygon, breaking it off (i.e. inserting an edge) and starting a new polygon as necessary.

ok: BOOL ← FALSE;

prew1, forw1, prew2, forw2: Vertex; -- in the counterclockwise order on the outline (with respect to interior of figure), we have (prew2, w2, forw2, prew1, w1, forw1)

pointLineRelation: EG.PointLineRelation;

prew2 ← PreviousOutlineVertex[w2, v];

forw2 ← v;

prew1 ← v;

forw1 ← NextOutlineVertex[v, w1];

Invariant assumed to hold at the beginning of this loop: w1 and w2 have the property that an edge joining them will "close off", i.e. form, a convex polygon "on the side of [w1,w2] towards the interior of the figure".

WHILE NOT ok DO

ok ← TRUE;

Move w2 backwards (clockwise) along the outline until edges [w1, w2] added to the outline will give a "local convex extension" of the outline, i.e. the outline will satisfy convexity check from before w2 up to and including w1. For each successive w2 value, we check whether adding edges [w1, prew2] will leave a convex polygon; if not, we put in edges [w1, w2] (by loop invariant, this will yield a convex polygon).

pointLineRelation ← EG.PointLineTest[w1.coordinates, w2.coordinates, prew2.coordinates];

WHILE pointLineRelation = LEFTOFLINE DO

ok ← FALSE;

IF NOT ConvexPolygonOnOutline[prew2, w2, prew1, w1] THEN {

JoinVertices[w1, w2, NIL, NIL, FALSE, FALSE, TRUE, FALSE]; -- [w1, w2] is an exterior edge

SetPolygonFields[w2, w1, FALSE, FALSE, NIL]; -- go around newly formed internal polygon setting leftRegionExterior = FALSE and leftClientData= NIL.

prew1 ← w2; -- take account of presence of new edge joining w1 and w2

IF verbosity#NIL THEN verbosity.out.PutF["\nGrowToConvexOutline has joined vertices #%g and #%g\n", IO.int[w1.index], IO.int[w2.index]];

};

forw2 ← w2; w2 ← prew2; prew2 ← PreviousOutlineVertex[w2, forw2];

pointLineRelation ← EG.PointLineTest[w1.coordinates, w2.coordinates, prew2.coordinates];

ENDLOOP;

Move w1 forwards (counterclockwise) along the outline until edges [w2, w1] added to the outline will give a "local convex extension" of the outline, i.e. the outline will satisfy convexity check from before w2 up through and beyond w1. For each successive w1 value, we check whether adding edges [w2, forw1] will leave a convex polygon; if not, we put in edges [w2, w1] (by loop invariant, this will yield a convex polygon).

pointLineRelation ← EG.PointLineTest[w2.coordinates, w1.coordinates, forw1.coordinates];

WHILE pointLineRelation = RIGHTOFLINE DO

ok ← FALSE;

IF NOT ConvexPolygonOnOutline[ w2, forw2, w1, forw1] THEN {

JoinVertices[w2, w1, NIL, NIL, FALSE, FALSE, FALSE, TRUE]; -- [w1, w2] is an exterior edge

SetPolygonFields[w2, w1, FALSE, FALSE, NIL]; -- go around newly formed internal polygon setting leftRegionExterior = FALSE and leftClientData= NIL.

forw2 ← w1; -- take account of presence of new edge joining w1 and w2

IF verbosity#NIL THEN verbosity.out.PutF["\nGrowToConvexOutline has joined vertices #%g and #%g\n", IO.int[w2.index], IO.int[w1.index]];

};

prew1 ← w1; w1 ← forw1; forw1 ← NextOutlineVertex[prew1, w1];

pointLineRelation ← EG.PointLineTest[w2.coordinates, w1.coordinates, forw1.coordinates];

ENDLOOP;

JoinVertices[w2, w1, NIL, NIL, FALSE, FALSE, FALSE, TRUE]; -- [w1, w2] is an exterior edge

SetPolygonFields[w2, w1, FALSE, FALSE, NIL]; -- go around newly formed internal polygon setting leftRegionExterior = FALSE and leftClientData= NIL.

IF verbosity#NIL THEN verbosity.out.PutF["\nGrowToConvexOutline has joined vertices #%g and #%g\n", IO.int[w2.index], IO.int[w1.index]];

RETURN [w2, w1]

};

ConvexifyOutline: PUBLIC PROC [outlineStart, outlineEnd: Vertex, verbosity: Verbosity ← NIL] RETURNS [newOutlineStart, newOutlineEnd: Vertex ← NIL] ~ {

firstOutlineStart: Vertex ← outlineStart;

outlineNext: Vertex ← NextVertex[outlineStart, outlineEnd, TRUE];

SELECT EG.PointLineTest[outlineStart.coordinates, outlineEnd.coordinates, outlineNext.coordinates] FROM

LEFTOFLINE, ONLINE => {

[outlineStart, outlineEnd] ← StepVertex[outlineStart, outlineEnd, TRUE];

IF outlineStart = firstOutlineStart THEN RETURN[outlineStart, outlineEnd];

};

RIGHTOFLINE => {

[outlineStart, outlineEnd] ← GrowToConvexOutline[outlineStart, outlineEnd, outlineNext, verbosity];

firstOutlineStart ← outlineStart;

};

ENDCASE;

ENDLOOP;

};

************** Region Extraction **************

InternalPolygons: PUBLIC PROC [v, w: Vertex] RETURNS [RegionList] = {

vTow, wTov: Adjacency;

firstv: Vertex;

outline: BOOL ← FALSE;

polygonOutline: Imager.Trajectory;

clientData: REF;

vec: Imager.VEC;

region: Region;

outList: RegionList;

firstv ← v;

vec ← EG.ImagerVecFromPoint[v.coordinates];

polygonOutline ← ImagerPath.MoveTo[vec];

Extract clientData

[ vTow, wTov ] ← EG.FindAdjacency[ v, w ];

clientData ← GetEdgeClientData[vTow]; -- only need to set state once

Set visited fields of edges of this polygon, and build polygonOutline

vec ← EG.ImagerVecFromPoint[w.coordinates];

polygonOutline ← ImagerPath.LineTo[polygonOutline, vec ];

vTow.visited ← TRUE;

WHILE w # firstv DO

[v, w] ← StepVertex[v,w, outline];

vec ← EG.ImagerVecFromPoint[w.coordinates];

polygonOutline ← ImagerPath.LineTo[polygonOutline, vec ];

[ vTow, wTov ] ← EG.FindAdjacency[ v, w ];

vTow.visited ← TRUE;

ENDLOOP;

region ← NEW[RegionRec ← [outline: polygonOutline, clientData: clientData] ];

outList ← LIST[region];

Go around edges again; check (opposite orientation of) each to see if it's an interior edge (i.e. not on the outline), and if the polygon to which it belongs is unprocessed, if so, recursive call on that polygon.

[v, w] ← StepVertex[v, w, outline];

[ vTow, wTov ] ← EG.FindAdjacency[ v, w ];

IF NOT GetEdgeExterior[wTov] THEN

IF NOT wTov.visited THEN TRUSTED {

outList ← LOOPHOLE[List.Append[LOOPHOLE[outList], LOOPHOLE[InternalPolygons[w, v]]]];

};

WHILE w # firstv DO

[v, w] ← StepVertex[v,w, outline];

[ vTow, wTov ] ← EG.FindAdjacency[ v, w ];

IF NOT GetEdgeExterior[wTov] THEN

IF NOT wTov.visited THEN TRUSTED {

outList ← LOOPHOLE[List.Append[LOOPHOLE[outList], LOOPHOLE[InternalPolygons[w, v]]]];

};

ENDLOOP;

RETURN[outList];

};

RegionBoundaryEdge: PROC [start, end: Vertex, isA: IsAProc] RETURNS [BOOL] ~ {

XOR: PROC[a, b: BOOLEAN] RETURNS[BOOLEAN] = {

IF (a AND b) THEN RETURN[FALSE];

IF (a OR b) THEN RETURN[TRUE];

RETURN[FALSE];

};

startToEnd, endToStart: Adjacency;

[startToEnd, endToStart] ← EG.FindAdjacency[start, end];

RETURN[XOR[isA[GetEdgeClientData[startToEnd]], isA[GetEdgeClientData[endToStart]]]];

};

NextRegionBoundaryEdge: PROC[start, end: Vertex, isA: IsAProc] RETURNS [next: Vertex] = {

A: AdjacencyList ← end.adjacentVertices;

WHILE A.first.vertex # start DO

A ← A.rest

ENDLOOP;

A ← A.rest; next ← A.first.vertex;

WHILE NOT RegionBoundaryEdge[end, next, isA] DO

A ← A.rest; next ← A.first.vertex;

ENDLOOP;

RETURN[next];

};

MaximalRegions: PUBLIC PROC [outlineVertex: Vertex, isA: IsAProc, verbosity: Verbosity ← NIL] RETURNS [outRegionList: RegionList] ~ {

MaximalRegionsSubr: PROC [v: Vertex, vertexVisitedValue: BOOLEAN, isA: IsAProc, inRegionList: RegionList, verbosity: Verbosity ← NIL] RETURNS [outRegionList: RegionList] ~ {

Basically a Depth First Search, pausing to build an entire region each time a region boundary edge is discovered.

vAdjList: AdjacencyList ← v.adjacentVertices;

firstVertex, nextVertex: Vertex;

vToNext, nextToV: Adjacency;

outRegionList ← inRegionList;

IF verbosity#NIL THEN verbosity.out.PutF["\nMaximalRegionsSubr entered with argument #%g\n", IO.int[v.index] ];

Go around adjacencies, looking for region boundary edges.

firstVertex ← nextVertex ← vAdjList.first.vertex;

[vToNext, nextToV] ← EG.FindAdjacency[v, nextVertex] ;

IF RegionBoundaryEdge[v, nextVertex, isA] AND isA[GetEdgeClientData[vToNext]] AND NOT vToNext.visited THEN {

IF verbosity#NIL THEN verbosity.out.PutF["\nMaximalRegionsSubr finds region boundary edge [%g, %g]\n", IO.int[v.index], IO.int[nextVertex.index]];

outRegionList ← CONS[RegionBuilder[v, nextVertex, FALSE, vertexVisitedValue, isA, verbosity], outRegionList];

};

vAdjList ← vAdjList.rest;

nextVertex ← vAdjList.first.vertex;

WHILE nextVertex#firstVertex DO

[vToNext, nextToV] ← EG.FindAdjacency[v, nextVertex] ;

IF RegionBoundaryEdge[v, nextVertex, isA] AND isA[GetEdgeClientData[vToNext]] AND NOT vToNext.visited THEN {

IF verbosity#NIL THEN verbosity.out.PutF["\nMaximalRegionsSubr finds region boundary edge [%g, %g]\n", IO.int[v.index], IO.int[nextVertex.index]];

outRegionList ← CONS[RegionBuilder[v, nextVertex, FALSE, vertexVisitedValue, isA, verbosity], outRegionList];

};

vAdjList ← vAdjList.rest;

nextVertex ← vAdjList.first.vertex;

ENDLOOP;

v.visited ← vertexVisitedValue; -- now we can say that we've been here.

Go around adjacencies again, making recursive calls for unvisited vertices.

IF (nextVertex.visited # vertexVisitedValue) THEN outRegionList ← MaximalRegionsSubr[nextVertex, vertexVisitedValue, isA, outRegionList, verbosity];

vAdjList ← vAdjList.rest;

nextVertex ← vAdjList.first.vertex;

WHILE nextVertex#firstVertex DO

IF (nextVertex.visited # vertexVisitedValue) THEN outRegionList ← MaximalRegionsSubr[nextVertex, vertexVisitedValue, isA, outRegionList, verbosity];

vAdjList ← vAdjList.rest;

nextVertex ← vAdjList.first.vertex;

ENDLOOP;

RETURN[outRegionList];

};

IF NOT ThreeOrMoreVertices[outlineVertex] THEN RETURN[NIL];

outRegionList ← MaximalRegionsSubr[outlineVertex, NOT outlineVertex.visited, isA, NIL, verbosity];

ClearGraphEdgeFields[outlineVertex];

RETURN[outRegionList];

};

RegionBuilder: PUBLIC PROC [outlineStart, outlineEnd: Vertex, setOutlineVerticesVisited: BOOL, vertexVisitedValue: BOOLEAN, isA: IsAProc, verbosity: Verbosity ← NIL] RETURNS [region: Region] ~ {

SearchInternalPolygons: PROC [start, end: Vertex, vertexVisitedValue: BOOLEAN, isA: IsAProc, inHoles: RegionList] RETURNS [outHoles: RegionList] ~ {

firstStart: Vertex ← start;

startToEnd, endToStart: Adjacency;

nextStartToNextEnd, nextEndTonextStart: Adjacency;

polygonStep: BOOL ← FALSE;

outHoles ← inHoles;

IF verbosity#NIL THEN verbosity.out.PutF["\nSearchInternalPolygons entered with edge [%g, %g]\n", IO.int[start.index], IO.int[end.index]];

Set edge visited fields, and vertex visited fields of non-outline vertices.

[ startToEnd, endToStart ] ← EG.FindAdjacency[start, end];

WHILE end # firstStart DO

[start, end] ← StepVertex[start, end, polygonStep];

[nextStartToNextEnd, nextEndTonextStart] ← EG.FindAdjacency[start, end];

nextStartToNextEnd.visited ← TRUE;

IF NOT IPVertex[start] THEN start.visited ← vertexVisitedValue;

startToEnd ← nextStartToNextEnd;

endToStart ← nextEndTonextStart;

ENDLOOP;

[start, end] ← StepVertex[start, end, polygonStep];

[nextStartToNextEnd, nextEndTonextStart] ← EG.FindAdjacency[start, end];

nextStartToNextEnd.visited ← TRUE;

IF NOT IPVertex[start] THEN start.visited ← vertexVisitedValue;

Go around edges again, checking opposite facing edges of non-outline edges. If such an edge is not yet visited, then if a region boundary edge, call RegionBuilder recursively to build a new hole, else call SearchInternalPolygons recursively.

[ startToEnd, endToStart ] ← EG.FindAdjacency[start, end];

IF NOT GetEdgeIPEdge[startToEnd] THEN

IF NOT endToStart.visited THEN

IF RegionBoundaryEdge[end, start, isA] THEN {

IF verbosity#NIL THEN verbosity.out.PutF["\nSearchInternalPolygons finds boundary edge [%g, %g]\n", IO.int[end.index], IO.int[start.index]];

outHoles ← CONS[ RegionBuilder[end, start, TRUE, vertexVisitedValue, isA, verbosity], outHoles]

}

ELSE

outHoles ← SearchInternalPolygons[end, start, vertexVisitedValue, isA, outHoles];

WHILE end # firstStart DO

[start, end] ← StepVertex[start, end, polygonStep];

[startToEnd, endToStart] ← EG.FindAdjacency[start, end];

IF NOT GetEdgeIPEdge[startToEnd] THEN

IF NOT endToStart.visited THEN

IF RegionBoundaryEdge[end, start, isA] THEN {

IF verbosity#NIL THEN verbosity.out.PutF["\nSearchInternalPolygons finds boundary edge [%g, %g]\n", IO.int[end.index], IO.int[start.index]];

outHoles ← CONS[ RegionBuilder[end, start, TRUE, vertexVisitedValue, isA, verbosity], outHoles];

}

ELSE

outHoles ← SearchInternalPolygons[end, start, vertexVisitedValue, isA, outHoles];

ENDLOOP;

RETURN[outHoles];

};

polygonOutline: Imager.Trajectory ← NIL;

vec: Imager.VEC;

clientData: REF;

startToEnd, endToStart: Adjacency;

firstStart: Vertex ← outlineStart;

tempVertex: Vertex;

holes: RegionList;

firstTrajPoint: EG.Point;

IF verbosity#NIL THEN verbosity.out.PutF["\nRegionBuilder entered with region boundary edge [%g, %g]\n", IO.int[outlineStart.index], IO.int[outlineEnd.index]];

Trace outline, setting edge IPEdge and visited fields, setting outline vertex visited fields if called for by setOutlineVerticesVisited, and building outline trajectory. Redundant vertices are omitted from outline trajectory.

[ startToEnd, endToStart ] ← EG.FindAdjacency[outlineStart, outlineEnd];

SetEdgeIPEdge[startToEnd, TRUE];

startToEnd.visited ← TRUE;

IF setOutlineVerticesVisited THEN outlineEnd.visited ← vertexVisitedValue;

tempVertex ← NextRegionBoundaryEdge[outlineStart, outlineEnd, isA];

IF EG.PointLineTest[outlineStart.coordinates, tempVertex.coordinates, outlineEnd.coordinates] # ONLINE THEN { -- don't output this if collinear with prev and next

firstTrajPoint ← outlineEnd.coordinates;

vec ← EG.ImagerVecFromPoint[firstTrajPoint];

polygonOutline ← ImagerPath.MoveTo[vec];

};

WHILE outlineEnd # firstStart DO

outlineStart ← outlineEnd; outlineEnd ← tempVertex;

[ startToEnd, endToStart ] ← EG.FindAdjacency[outlineStart, outlineEnd];

SetEdgeIPEdge[startToEnd, TRUE];

startToEnd.visited ← TRUE;

IF setOutlineVerticesVisited THEN outlineEnd.visited ← vertexVisitedValue;

tempVertex ← NextRegionBoundaryEdge[outlineStart, outlineEnd, isA];

IF EG.PointLineTest[outlineStart.coordinates, tempVertex.coordinates, outlineEnd.coordinates] # ONLINE THEN {

IF polygonOutline # NIL THEN { -- don't output this if collinear with prev and next

vec ← EG.ImagerVecFromPoint[outlineEnd.coordinates];

polygonOutline ← ImagerPath.LineTo[polygonOutline, vec];

}

ELSE {

firstTrajPoint ← outlineEnd.coordinates;

vec ← EG.ImagerVecFromPoint[firstTrajPoint];

polygonOutline ← ImagerPath.MoveTo[vec];

};

ENDLOOP;

Get back to initial outlineStart, outlineEnd

outlineStart ← outlineEnd; outlineEnd ← tempVertex;

vec ← EG.ImagerVecFromPoint[firstTrajPoint];

polygonOutline ← ImagerPath.LineTo[polygonOutline, vec]; -- close trajectory

Client data of outline - holes.

clientData ← GetEdgeClientData[startToEnd];

Search for holes

holes ← SearchInternalPolygons[outlineStart, outlineEnd, vertexVisitedValue, isA, NIL];

region ← NEW[RegionRec ← [outline: polygonOutline, clientData: clientData, holes: holes] ];

};

**************** Input and Output ****************

ropeFromEdgeData: EG.RopeFromEdgeData ~ {

cEGEdgeData: CEGEdgeData ← NARROW[edgeData];

IF cEGEdgeData = NIL THEN RETURN["NIL"];

out ← Rope.Cat[Convert.RopeFromBool[cEGEdgeData.exteriorEdge], ",", Convert.RopeFromBool[cEGEdgeData.intersectionPolygonEdge], "," ];

out ← Rope.Concat[out, clientProc1[cEGEdgeData.clientData] ];

};

DumpGraph: PUBLIC PROC [v: ConvexEuclideanGraph, ropeFromClientData: RopeFromClientData, filename: Rope.ROPE] ~ {

EG.DumpGraph[v, ropeFromEdgeData, filename, ropeFromClientData];

};

edgeDataFromRope: EG.EdgeDataFromRope ~ {

dataStream: IO.STREAM ← IO.RIS[in];

char: CHAR;

exteriorEdge, intersectionPolygonEdge: BOOLEAN;

clientData: REF;

DataFromRopeFail: PUBLIC ERROR [subclass: ATOM ← $Unspecified] = CODE;

[]← dataStream.SkipWhitespace[]; char ← dataStream.PeekChar[];

IF char = 'N THEN {

[] ← dataStream.GetChar[]; -- toss N

IF dataStream.GetChar[] # 'I THEN DataFromRopeFail[$IExpected];

IF dataStream.GetChar[] # 'L THEN DataFromRopeFail[$LExpected];

RETURN[NIL]

};

IF char # 'T AND char # 'F THEN DataFromRopeFail[$TorFExpected];

exteriorEdge ← IO.GetBool[dataStream];

IF dataStream.GetChar[] # ', THEN DataFromRopeFail[$commaExpected];

char ← dataStream.PeekChar[];

IF char # 'T AND char # 'F THEN DataFromRopeFail[$TorFExpected];

intersectionPolygonEdge ← IO.GetBool[dataStream];

IF dataStream.GetChar[] # ', THEN DataFromRopeFail[$commaExpected];

char ← dataStream.PeekChar[];

IF char = 'N THEN {

[] ← dataStream.GetChar[]; -- toss N

IF dataStream.GetChar[] # 'I THEN DataFromRopeFail[$IExpected];

IF dataStream.GetChar[] # 'L THEN DataFromRopeFail[$LExpected];

clientData ← NIL;

}

ELSE {

index: INT ← IO.GetIndex[dataStream];

rope: Rope.ROPE ← Rope.Substr[in, index];

clientData ← clientProc1[rope];

};

RETURN[NEW[CEGEdgeDataRec ← [exteriorEdge: exteriorEdge, intersectionPolygonEdge: intersectionPolygonEdge, clientData: clientData] ] ];

};

ReadIOGraphFromFile: PUBLIC PROC [filename: Rope.ROPE, clientDataFromRope: ClientDataFromRope, MessageStream: IO.STREAM] RETURNS [iOGraph: EG.IOGraph] ~ {

RETURN[EG.ReadIOGraphFromFile[filename, edgeDataFromRope, MessageStream, clientDataFromRope] ];

};

GraphFromIOGraph: PUBLIC PROC [iOGraph: EG.IOGraph] RETURNS [v: ConvexEuclideanGraph, numberVertices: CARDINAL] ~ {

Convert an IOGraph to a Vertex. The particular Vertex v returned will be chosen so as to be on the outline of the figure. numberVertices is the number of vertices encountered in the converted IOGraph.

vertexList, lastCell: VertexList;

scratchVertexList: VertexList;

scratchIOGraph: EG.IOGraph;

vertex: Vertex;

adjacentVertex: Vertex;

aList, aListLast: AdjacencyList;

adjacency : Adjacency;

ioAList: EG.IOAdjacencyList;

iOAdjacency: EG.IOAdjacency;

Create initial list of vertex indices

vertex ← NEW[EG.VertexRec ← [coordinates: iOGraph.first.coordinates,

index: iOGraph.first.index, adjacentVertices: NIL, data: NIL, visited: FALSE] ];

lastCell ← vertexList ← LIST[vertex];

numberVertices ← 1;

scratchIOGraph ← iOGraph.rest;

WHILE scratchIOGraph # NIL DO

vertex ← NEW[EG.VertexRec ← [coordinates: scratchIOGraph.first.coordinates,

index: scratchIOGraph.first.index, adjacentVertices: NIL, data: NIL, visited: FALSE] ];

lastCell.rest ← LIST[vertex]; lastCell ← lastCell.rest;

numberVertices ← numberVertices + 1;

scratchIOGraph ← scratchIOGraph.rest;

ENDLOOP;

scratchVertexList ← vertexList; -- Need to keep vertexList intact for searching

Fill in links to adjacent vertices

WHILE scratchVertexList # NIL DO

vertex ← scratchVertexList.first;

ioAList ← iOGraph.first.adjacentVertices;

aList ← NIL; -- initialize

WHILE ioAList # NIL DO -- assumes adjacent vertices occur in reverse clockwise order; this loop will restore correct order.

iOAdjacency ← ioAList.first;

adjacentVertex ← EG.SearchVertexList[vertexList, iOAdjacency.vertex ]; -- these calls assume call by value, i.e. that vertexList always points to the head of the list

adjacency ← NEW[EG.AdjacencyRec ← [vertex: adjacentVertex, data: iOAdjacency.data, visited: FALSE] ];

IF GetEdgeExterior[adjacency] THEN v ← vertex; -- identify vertex on outline of figure

aList ← CONS[adjacency, aList]; -- adjacencies in ioAList are in backwards (counterclockwise) order; using CONS here puts them back in clockwise order.

IF aList.rest = NIL THEN aListLast ← aList; -- save pointer for setting circularity

ioAList ← ioAList.rest;

ENDLOOP;

aListLast.rest ← aList; -- make list circular

vertex.adjacentVertices ← aList; -- attach to current vertex

scratchVertexList ← scratchVertexList.rest;

iOGraph ← iOGraph.rest; -- advance scratchVertexList and iOGraph in lock step

ENDLOOP;

RETURN[ v, numberVertices ];

};

GraphFromFile: PUBLIC PROC [filename: Rope.ROPE, clientDataFromRope: ClientDataFromRope, MessageStream: IO.STREAM] RETURNS [v: ConvexEuclideanGraph, numberVertices: CARDINAL] ~ {

[v, numberVertices] ← GraphFromIOGraph[ReadIOGraphFromFile[filename, clientDataFromRope, MessageStream] ];

};

END.