-- To do: Update for COGVoronoi, COGDraw

-- To do: Throw: repaint diagram after each point

-- To do: Throw: use a better range (now is astronomical)

-- To do: check locking scheme

BEGIN
?? rewrite!
-- parameters: [viewer: Viewer, clientData: REF ANY, redButton: BOOL]
vData: ViewerData ← NARROW [viewer.data];
out.PutF["Called SortMe...\n"];
vData.action ← $Sort;
TellMother[]
END;

BEGIN
-- Sorts the finite vertices of a Delaunay diagram in such a way that, for any two
-- adjacent points, the leftmost one occurs first. When this rule (and its transitive
-- closure) are insufficient to decide the ordering of two points, the one with lowest
-- y occurs first.
-- The final ordering is given by the 'no' field of each vertex and by the list sortedPts.
-- also, p.anEdge will point to the (finite) neighbor of p that follows p and has maximum y.
-- Infinite data points will get p.no = -1, -2, -3, ... and will be omitted from the sorted list.
-- The number of finite points is returned in n. maximum y.
pts, finPts, stack, free, t: LIST OF REF DVertex;
pBq, pBprev: BOOL;
p, q: REF DVertex;
order, infOrder: NAT;
e0, ea, ePrev: REF DEdge;

pts ← Diag.points;
out.PutF ["SortRegions - first pass...\n"];
WHILE pts # NIL DO
p ← pts.first; pts ← pts.rest;
p.nBlockers ← 0;
-- Counts blocking neighbors of p
-- Also puts in eTop and in p.anEdge the edge blocked by p with Dest of maximum y
-- Also finds first vertex in the ordering and initializes the Qyu to it.
e0 ← ea ← p.anEdge;
ePrev ← VPrev [e0]; pBprev ← Homo.Precedes [p.dp, Dest[ePrev].dp];
out.PutF ["Looking at point %g\n", int[p.no]];
DO
q ← Dest[ea];
pBq ← Homo.Precedes [p.dp, q.dp];
IF NOT pBq THEN
{IF pBprev THEN p.anEdge ← ePrev;
p.nBlockers ← p.nBlockers + 1};
ePrev ← ea; pBprev ← pBq; ea ← VNext[ea];
IF ea = e0 THEN EXIT
ENDLOOP;
-- delete p from points list and put in free list (or prime the stack if p is first in ordering)
IF p.nBlockers = 0 THEN
{stack ← LIST [p];
out.PutF["Point %g is first.\n", int[p.no]]}
ENDLOOP;

-- topological sorting: removes a vertex from the stack,
-- decrements ref.counts of blocked neighbors,
-- pushes freed neighbors into stack (in proper y order)
out.PutF ["SortRegions - second pass...\n"];
order ← 0; sortedPts ← finPts ← NIL;
infOrder ← -1;
WHILE stack # NIL DO
p ← stack.first;
out.PutF["Looking at point %g\n", int[p.no]];
IF trace THEN ShowRegion [Diag, p, thin];
t ← stack; stack ← stack.rest;
IF Homo.FinPt[p.dp] THEN
{p.no ← order; order ← order + 1;
IF sortedPts = NIL THEN sortedPts ← t ELSE finPts.rest ← t;
finPts ← t; t.rest ← NIL}
ELSE
{p.no ← infOrder; infOrder ← infOrder - 1};
e0 ← ea ← p.anEdge;
WHILE Homo.Precedes [p.dp, (q ← Dest[ea]).dp] DO
IF trace THEN ShowEdge [Diag, ea, bold];
q.nBlockers ← q.nBlockers - 1;
IF q.nBlockers = 0 THEN
{stack ← CONS [q, stack]};
ea ← VPrev[ea]; -- note: cw order here
IF ea = e0 THEN EXIT
ENDLOOP;
ENDLOOP;
Diag.points ← pts;
out.PutF ["Exit SortRegions.\n"]
END;

-- The Lee-Preparata structure is a binary search tree, whose leaves are the Voronoi regions and whose internal nodes are monotone chains of Voronoi edges. Each chain is by itself a binary search tree: the nodes of the latter are of two kinds, called YTests and EdgeTests.

-- YTest nodes branch depending on whether the query point q is above or below the testY field. When an EdgeTest is encountered, q.y is known to be between those of the edge endpoints; that being the case, q.x can be tested against the edge.

-- A node of either kind has two pointers, which can refer to Voronoi regions (ultimate leaves), to its children in the same chain C, or to the roots of the chains that are children of C in the global tree. The latter is always the case for EdgeTests, and sometimes also for the YTests. The latter case occurs when the edge of C that covers q.y belongs also to some ancestor C' of C. At this point, we know on which side of C' the query point q lies, and therefore we do not have to test it again against that edge; we can proceed immediately to the root of the appropriate descendant chain C'' of C.

BEGIN
-- Builds the Lee-Preparata tree structure for point location

-- Each chain is at first collected as list of edges; these are sorted by ascending y coordinate, and a chain tree is built from them. Since every edge is represented in only one chain (of highest possible level), there may be one or more gaps in the ranges of y-values covered by a given chain. A DummyEdge is temporarily inserted in each of those gaps (but not in the final structure).

DummyEdge: TYPE = RECORD
[topY: REAL,
onLeft: BOOL -- TRUE if gap corresponds to edge(s) of a left ancestor of the chain
];

-- Once the Voronoi regions have been sorted from left to right, each edge can be assigned to a specific chain based on the ranks of the two regions on each side of it. If the regions are numbered in the range [0..n) (n>1), the root chain separates regions [0..k) and [k..n), where k = n DIV 2 is the chain index. Each group of regions is split in half by a child of chain k, and so forth; in generla, if at some step the set of regions is [fst..lim) (with lim > fst+1), the root chain will separate regions [fst..r) from [r..lim), where r = (fst+lim) DIV 2. The chain indices will lie therefore in the range (0..n).
-- Given two adjacent regions i and j, the edge between them will belong to a unique chain whose index, LCA[i, j], is computed by the following procedure.

ChainIndex: PROC [i, j: NAT] RETURNS [ix: NAT] = TRUSTED

Descendants: PROC [ix: NAT] RETURNS [fst, lim: NAT] = TRUSTED

-- The folowing records are temporary data structures.

ChainTable: TYPE = RECORD
[elist: SEQUENCE nLists: NAT OF DEdgeList
];
TipsTable: TYPE = RECORD
[tip: SEQUENCE nRegs: NAT OF REF VVertex
];
RegionTable: TYPE = RECORD
[dvertex: SEQUENCE nRegs: NAT OF REF DVertex
];

n: NAT; -- number of finite data points in the diagram
edges: REF ChainTable; -- edges in each Lee-Preparata chain
tips: REF TipsTable; -- upper Voronoi vertices of each region
regions: REF RegionTable; -- upper Voronoi vertices of each region

SetUpTables: INTERNAL PROC = TRUSTED

CookYCoord: PROC [p: Homo.Point] RETURNS [y: REAL] = TRUSTED

SortChainEdges: INTERNAL PROC [le: DEdgeList, fst, lim: NAT, leftTip: REF VVertex]
RETURNS [sle: LIST OF REF ANY] = TRUSTED

Error: ERROR [what: ROPE] = CODE;

BuildPartialTreeFromChain: INTERNAL PROC
[sle: LIST OF REF ANY, ne: NAT, left, right: TreeRoot]
RETURNS
[root: TreeRoot, chTop: REAL, rest: LIST OF REF ANY] = TRUSTED

TreeData: TYPE = RECORD [root: TreeRoot, tip: REF VVertex];

BuildTreeOfChains: INTERNAL PROC [fst, lim: NAT]
RETURNS [rslt: TreeData] = TRUSTED

SetUpTables[];
RETURN [BuildTreeOfChains [0, n].root]
END;