-- 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).

-- 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 = (0+n) DIV 2 is the chain index. Each group of regions is split in half by a child of the chain k, and so forth; in general, a group of regions [fst..lim) obtained in this process (lim>fst+1) is split into two subgroups [fst..k) and [k..lim) by the chain whose index is k = (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 is ChainIndex[i, j]

-- Returns the index of the monotone chain that separates region i from region j (i#j) and has highest possible level in the Lee-Preparata tree. Assumes the regions have been renumbered from 0 up.

BEGIN
fst: NAT ← 0;
lim: NAT ← n;
DO
-- assert: i and j IN [fst..lim)
ix ← (fst + lim)/2;
IF i < ix THEN
IF j < ix THEN
lim ← ix
ELSE
RETURN
ELSE
IF j < ix THEN
RETURN
ELSE
fst ← ix
ENDLOOP
END;

-- Returns fst and lim such that LCA[i, j] = ix iff i IN [fst..ix) and j IN [ix..lim).

BEGIN
-- Returns fst and lim such that LCA[i, j] = ix iff i IN [fst..ix) and j IN [ix..lim).
av: NAT;
fst ← 0; lim ← n;
DO
-- Assert: ix IN [fst..lim)
av ← (fst + lim)/2;
IF ix < av THEN lim ← av
ELSE IF ix > av THEN fst ← av
ELSE RETURN
ENDLOOP
END;

-- Sorts the regions from left to right and sets the chain table in such a way that an edge between region i and region j is in elists[ChainIndex[i, j]].

-- Returns also in regions[k] n edge out of region k, and in tips[k] the its point with maximum y.

-- Finally, for every voronoi vertex v saves in maxchix[v.no] the maximum chain index of all edges for which it is the upper endpoint.

BEGIN

ThrowEdge: Diag.VisitProc = TRUSTED

n ← SortRegions [dlau, TRUE];
elists ← NEW [ChainTable[n]];
tips ← NEW [TipsTable[n]];
regions ← NEW [RegionTable[n]];
maxix ← NEW [MaxIxTable [dlau.limVV]]; -- better they be numbered from 0 to nVV
[] ← Diag.Traverse [dlau.dg, ThrowEdge, NIL, all];
END;

BEGIN
-- Returns the cartesian y-coordinate of p, or plus/minus LasgestNumber (with correct sign) if p is infinite.

RETURN
[IF Homo.InfPt[p] THEN
IF p.y > 0.0 THEN LargestNumber
ELSE IF p.y < 0.0 THEN - LargestNumber
ELSE 0.0
ELSE p.y/p.w]
END;

-- Sorts the edges of a chain by increasing Y, filling the gaps with dummy edges. The chain is assumed to bisect the regions [fst..lim); these parameters are used to figure out whether gaps lie on the left or right boundary of the current region set.

-- A Voronoi vertex lies on the left boundary iff it is adjacent to some region with sequence number lim or more.

-- The leftTip parameter gives the Voronoi vertex with maximum y on the regions at the left of the chain. It is used as a separator when the chain is empty.

BEGIN OPEN List;

DecreasingY: CompareProc = TRUSTED

OnLeftBoundary: PROC [v: VVertex] RETURNS [is: BOOL] = TRUSTED
{RETURN[maxix[v.no] < lim]};
tle: LIST OF REF ANY ← List.Sort [IsAListOfRefAny[le].list, DecreasingY];
t: LIST OF REF ANY;
e: Diag.DEdge;
de: REF DummyEdge;
lastV: VVertex ← NIL; -- Bottom vertex of chain processed so far (NIL at beginning)
eTop: VVertex;

sle ← NIL;
IF tle = NIL THEN
{-- no proper edges in chain - put a dummy edge from top to y separator
de ← NEW [DummyEdge ←
[topV: lastV,
onLeft: TRUE]];
sle ← CONS[de, sle];
lastV ← leftTip}
ELSE
{WHILE tle # NIL DO
t ← tle; e ← NARROW [t.first, REF Diag.DEdge]^; tle ← t.rest;
eTop ← Left[e];
IF Homo.Below [eTop.pt, IF lastV = NIL THEN [0.0, 1.0, 0.0] ELSE lastV.pt] THEN
{-- Found a gap
de ← NEW [DummyEdge ←
[topV: lastV,
onLeft: OnLeftBoundary[Left[e]]]];
sle ← CONS[de, sle];
lastV ← eTop};
t.rest ← sle; sle ← t;
lastV ← Right[e];
ENDLOOP};
-- Put last dummy edge if needed
IF Homo.Below [[0.0, -1.0, 0.0], lastV.pt] THEN
{de ← NEW [DummyEdge ←
[topV: lastV,
onLeft: OnLeftBoundary[lastV]]];
sle ← CONS[de, sle]}
END;

-- Builds a search tree for the first ne edges (including dummy ones) in the sorted list sle. The parameters left and right give the roots of the childern chains on either side.

-- Returns also the vertex chTop at the upper end of the partial chain, and the list rest of unused edges.

BEGIN
IF ne = 1 THEN
{WITH sle.first SELECT FROM
de: REF DummyEdge =>
{RETURN
[root: IF de.onLeft THEN right ELSE left,
chTop: de.topV,
rest: sle.rest]};
re: REF Diag.DEdge =>
{RETURN
[root: NEW [EdgeTest ← [edge: re^, left: left, right: right]],
chTop: Left[re^],
rest: sle.rest]};
ENDCASE => ERROR Bug.Error["UFO in chain."]}
ELSE
{up, dn: Tlee;
dnTop, upTop: VVertex;
vt: REF VTest;
[dn, dnTop, sle] ← BuildPartialTreeFromChain [sle, ne/2, left, right];
[up, upTop, sle] ← BuildPartialTreeFromChain [sle, ne - ne/2, left, right];
vt ← NEW [VTest ← [testV: dnTop, up: up, dn: dn]];
RETURN[vt, upTop, sle]}
END;

BEGIN
-- Builds a Lee-Preparata tree for regions [fst..lim).

IF fst+1 = lim THEN
{RETURN [[Org[regions[fst]], tips[fst]]]}
ELSE
{rtix: NAT = (fst + lim)/2;
left: TreeData = BuildTreeOfChains [fst, rtix];
right: TreeData = BuildTreeOfChains [rtix, lim];
sle: LIST OF REF ANY = SortChainEdges [elists[rtix], fst, lim, left.tip];
ne: NAT = List.Length[sle];
rootChain: Tlee = BuildPartialTreeFromChain[sle, ne, left.root, right.root].root;
IF trace THEN
{CleanTheDrawing[repaint: TRUE];
FOR no: VertexNo IN [fst..lim) DO ShowRegion [regions[no], yellow, 2] ENDLOOP;
ShowChain [rootChain, left.root, right.root, green, 2]};
RETURN [[root: rootChain,
tip: IF Homo.Below[right.tip.pt, left.tip.pt] THEN left.tip ELSE right.tip]]}
END;