-- To do: debug ShowRegion - sometimes shows complement of (infinite) region instead of region

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

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

-- To do: Repaint and Reset commands

-- To do: Make of it a "Geometrical Viewer module"

-- To do: check locking scheme

{RETURN [(e.vNext).sym]};

{RETURN [e.fNext]};

{RETURN [(e.sym).fNext]};

{RETURN [e.vNext]};

{RETURN [e.dest]};

{RETURN [(e.sym).dest]};

{RETURN [e.leftVV]};

{RETURN [(e.sym).leftVV]};

BEGIN
-- Computes the Voronoi vertex that corresponds to the Delaunay face passing through p, q and r (i.e., their circumcenter). If p, q and r do not form a convex angle, they are assumed to be on the convex hull; in that case, the infinity point perpendicular to pq is returned.
IF Homo.Convex [p, q, r]
THEN {v ← Homo.CircumCenter [p, q, r]}
ELSE {dir: Homo.Vector ← Homo.Sub[q, p];
v ← IF dir.w < 0.0 THEN [dir.y, -dir.x, 0.0] ELSE [-dir.y, dir.x, 0.0]}
END;

BEGIN
-- inserts an edge e (and its symmetrical) between the Delaunay vertices p and q.
-- ep must be the edge out of p preceding immediately e c.c.w.,
-- and eq must be the edge out of q immediately preceding e.sym c.c.w.
out.PutF["Adding edge (%g, %g)", int[p.no], int[q.no]];
IF ep # NIL THEN out.PutF[" after (%g, %g) at p ", int[Org[ep].no], int[Dest[ep].no]];
IF eq # NIL THEN out.PutF[" after (%g, %g) at q", int[Org[eq].no], int[Dest[eq].no]];
out.PutF["...\n"];
e ← NEW [Edge ← [dest: q, fNext: NIL, vNext: NIL, sym: NIL, leftVV: lv]];
e.sym ← NEW [Edge ← [dest: p, fNext: NIL, vNext: NIL, sym: e, leftVV: rv]];
lv.anEdge ← e;
rv.anEdge ← e.sym;
IF ep # NIL THEN
{e.vNext ← ep.vNext; e.sym.fNext ← ep;
ep.vNext.sym.fNext ← e; ep.vNext ← e}
ELSE
{e.vNext ← e; e.sym.fNext ← e; p.anEdge ← e};
IF eq # NIL THEN
{e.sym.vNext ← eq.vNext; e.fNext ← eq;
eq.vNext.sym.fNext ← e.sym; eq.vNext ← e.sym}
ELSE
{e.sym.vNext ← e.sym; e.fNext ← e.sym; q.anEdge ← e.sym}
END;

BEGIN
-- Deletes the edge e (and its symmetrical).
p: REF DVertex = Org[e];
q: REF DVertex = Dest[e];
ep: REF Edge = VPrev[e];
eq: REF Edge = VPrev[e.sym];
out.PutF["Deleting edge (%g, %g)...\n", int[p.no], int[q.no]];
IF ep # e THEN
{ep.vNext ← e.vNext; e.vNext.sym.fNext ← ep;
IF p.anEdge = e THEN p.anEdge ← ep}
ELSE
p.anEdge ← NIL;
IF eq # e THEN
{eq.vNext ← e.sym.vNext; e.sym.vNext.sym.fNext ← eq;
IF q.anEdge = e.sym THEN q.anEdge ← eq}
ELSE
q.anEdge ← NIL
END;

BEGIN
-- Returns the face of the Voronoi containing point p
-- Temporary version - uses exhaustive enumeration
pts: LIST OF REF DVertex ← Diag.points;
pa: REF DVertex;
hp: HPoint ← [p.x, p.y, 1.0];
minD: REAL;
paD: REAL;
np ← pts.first; pts ← pts.rest;
minD ← Homo.DistSq[np.dp, hp];
WHILE pts # NIL DO
pa ← pts.first;
IF Homo.FinPt[pa.dp] AND minD > (paD ← Homo.DistSq[pa.dp, hp]) THEN
{minD ← paD; np ← pa};
pts ← pts.rest
ENDLOOP
END;

BEGIN
-- Adds a new finite vData point newP to the Delaunay diagram Diag
-- Temporary version - produces a triangulation in any case.
f0: REF DVertex ← DumbLocateInVoronoi [newP, Diag];
p: REF DVertex = NEW [DVertex ←
[dp: [newP.x, newP.y, 1.0],
no: (nDV ← nDV+1),
anEdge: NIL]];
q, f, fn: REF DVertex;
ef: REF Edge ← f0.anEdge;
ep, eq, en: REF Edge;
out.PutF ["Inserting point at [%g, %g]...\n", real[newP.x], real[newP.y]];
ShowPoint [Diag, p.dp, thin];
ShowPoint [Diag, f0.dp, bold];
IF p.dp = f0.dp THEN
{out.PutF["Duplicate point - ignored\n"];
RETURN};
-- locate some neighbor f of f0 following p c.c.w. --
WHILE Homo.Convex [p.dp, f0.dp, Dest[ef].dp] DO ef ← VNext[ef] ENDLOOP;
-- locate the new neighbors of p (in c.c.w. order) and delete diagonal edges between them --
q ← f0;
DO
out.PutF["Will connect to vertex %g.\n", int[q.no]];
-- find the first neighbor of q that immediately precedes p c.c.w. around q --
DO
ef ← VPrev[ef];
IF Homo.Convex [p.dp, q.dp, Dest[ef].dp] THEN EXIT
ENDLOOP;
f ← Dest [ef];
-- delete diagonal edges from q --
DO
-- check whether ef is gonna stay in the final Delaunay --
en ← VPrev[ef]; fn ← Dest[en];
IF NOT Homo.Convex[p.dp, q.dp, fn.dp] OR
NOT Homo.Convex[fn.dp, f.dp, p.dp] OR
Homo.ShouldConnect24[p.dp, q.dp, fn.dp, f.dp] THEN EXIT;
ShowEdge [Diag, ef, boldGray];
DeleteEdge [ef];
ef ← en; f ← fn
ENDLOOP;
-- found the next neighbor of p: it is f --
-- fix LeftVV of edge ef
ef.leftVV ← NEW [VVertex ←
[vp: ComputeVoronoiVertex [p.dp, q.dp, f.dp],
anEdge: ef
]];
ShowEdge [Diag, ef, bold];
q ← f; ef ← ef.sym; IF q = f0 THEN EXIT
ENDLOOP;
-- add p and the Delaunay edges from it to the Delaunay face just created
ep ← NIL; ef ← VPrev[ef];
DO
-- connect p and q (ef is the edge immediately preceding (q,p) c.c.w around q,
-- and ep is the one immediately preceding (p,q) c.c.w. around p)--
eq ← ConnectPoints [q, p, ef, ep, LeftVV[ef], RightVV[VNext[ef]]]; ep ← eq.sym;
ShowEdge [Diag, ep, bold];
q ← Dest [ef]; ef ← FNext [ef]; IF q = f0 THEN EXIT
ENDLOOP;
-- add p to the list of points --
Diag.points ← CONS [p, Diag.points]
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 Edge;

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 EdgeList
];
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: EdgeList, 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;

BEGIN
vData: ViewerData = NARROW [Diag.viewer.data];
delau:BOOL = vData.showDelaunay;
p: REF DVertex = Org[e];
q: REF DVertex = Dest[e];
thing: REF GraphicThing;
IF (IF delau THEN Homo.FinPt[p.dp] OR Homo.FinPt [q.dp]
ELSE Homo.FinPt[p.dp] AND Homo.FinPt [q.dp]) THEN
{thing ← NEW [GraphicThing ←
[style: style,
obj: NEW [LineSegment ←
IF delau THEN [org: Org[e].dp, dest: Dest[e].dp]
ELSE [org: LeftVV [e].vp, dest: RightVV [e].vp]]
]];
vData.extraThings ← CONS [thing, vData.extraThings];
ViewerOps.PaintViewer [Diag.viewer, client, FALSE, thing];
Rest[]}
END;

BEGIN
vData: ViewerData = NARROW [Diag.viewer.data];
delau:BOOL = vData.showDelaunay;
e0, ea: REF Edge;
prevInf: BOOL ← TRUE;
vs: Polygon ← NIL;
thing: REF GraphicThing;
IF NOT Homo.FinPt[p.dp] THEN RETURN;
thing ← NEW [GraphicThing ← [style: style, obj: NIL]];
e0 ← ea ← p.anEdge;
DO
IF Homo.FinPt [Dest[ea].dp] THEN
{IF prevInf THEN vs ← CONS [RightVV[ea].vp, vs];
vs ← CONS [LeftVV[ea].vp, vs];
prevInf ← FALSE}
ELSE
{prevInf ← TRUE};
ea ← VNext[ea]; IF ea = e0 THEN EXIT
ENDLOOP;
thing.obj ← vs;
vData.extraThings ← CONS [thing, vData.extraThings];
ViewerOps.PaintViewer [Diag.viewer, client, FALSE, thing];
Rest[]
END;

BEGIN
vData: ViewerData = NARROW [Diag.viewer.data];
thing: REF GraphicThing =
NEW [GraphicThing ← [style: style, obj: NEW [Homo.Point ← pt]]];
vData.extraThings ← CONS [thing, vData.extraThings];
ViewerOps.PaintViewer [Diag.viewer, client, FALSE, thing];
Rest[]
END;

BEGIN
InsertInVoronoi [newP, Diag]
END;

BEGIN
[] ← SortRegions [Diag, TRUE]
END;

BEGIN
-- Maps a point from user coordinates to Viewer coordinates
OPEN scales;
xx ← x*scale + orgX;
yy ← y*scale + orgY
END;

BEGIN
-- Maps a point from Viewer coordinates to user coordinates
OPEN scales;
x ← (xx - orgX)/scale;
y ← (yy - orgY)/scale
END;

BEGIN
-- Maps a point from homogeneous user coordinates to Viewer coordinates
OPEN scales;
s: REAL;
IF p.w # 0.0 THEN s ← scale/p.w
ELSE {IF p.x = 0.0 AND p.y = 0.0 THEN p.x ← p.y ← 1.0;
s ← 20.0*scrSize/MAX [ABS [p.x], ABS [p.y]]};
xx ← p.x*s + orgX;
yy ← p.y*s + orgY
END;

BEGIN OPEN Random, Real;
p.x ← Float[Next[]]/LAST[INT];
p.y ← Float[Next[]]/LAST[INT];
END;

BEGIN OPEN Graphics, GraphicsColor;

END;

BEGIN
-- parameters: [self: Viewer]
vData: ViewerData ← NARROW[self.data];
vData.live ← FALSE;
TellMother[]
END;

BEGIN
-- parameters: [self: Viewer]
out.PutF["Called InitMe...\n"]
END;

BEGIN
-- parameters: [viewer: Viewer, clientData: REF ANY, redButton: BOOL]
vData: ViewerData ← NARROW [viewer.data];
out.PutF["Called ShowDualOfMe...\n"];
vData.showDelaunay ← NOT vData.showDelaunay;
vData.action ← $Repaint;
TellMother[]
END;

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

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

BEGIN
-- parameters [self: Viewer, input: LIST OF REF ANY]
-- N.B. Called at Process.priorityForeground!
-- we have conspired to make the leading item in the list
-- an atom that tells us what to do with the rest of the list
IF input # NIL AND input.rest # NIL AND
(input.first = $LeftDown OR input.first = $CenterDown OR
input.first = $RightDown) THEN
{vData: ViewerData = NARROW[self.data];
coord: TIPTables.TIPScreenCoords ← NARROW[input.rest.first, TIPTables.TIPScreenCoords];
vData.action ← input.first;
[vData.pt.x, vData.pt.y] ← UnmapPoint [coord.mouseX, coord.mouseY, vData.scales];
TellMother []}
END;

BEGIN
-- creates the menu
-- the leftmost 4 entries are in standard order
m ← Menus.CreateMenu[];
Menus.AppendMenuEntry[menu: m, name: "Close", proc: ViewerMenus.Close];
Menus.AppendMenuEntry[menu: m, name: "Grow", proc: ViewerMenus.Grow];
Menus.AppendMenuEntry[menu: m, name: "<-->", proc: ViewerMenus.Move];
Menus.AppendMenuEntry[menu: m, name: "Quit", proc: ViewerMenus.Destroy, fork: TRUE];
Menus.AppendMenuEntry[menu: m, name: "Dual", proc: ShowDualOfMe];
Menus.AppendMenuEntry[menu: m, name: "Sort", proc: SortMe];
Menus.AppendMenuEntry[menu: m, name: "Throw", proc: ThrowMe]
END;

BEGIN
-- returns a new Delaunay diagram with only three dummy vData points
-- at the Mercedes infinities.
sqrt3: REAL = SqRt[3.0];
pi1: REF DVertex = NEW [DVertex ← [dp: [0.0, sqrt3/3.0, 0.0], no: -1, anEdge: NIL]];
pi2: REF DVertex = NEW [DVertex ← [dp: [-0.5, -sqrt3/6.0, 0.0], no: -2, anEdge: NIL]];
pi3: REF DVertex = NEW [DVertex ← [dp: [0.5, -sqrt3/6.0, 0.0], no: -3, anEdge: NIL]];
v123: REF VVertex = NEW [VVertex ← [vp: [0.0, 0.0, 1.0], anEdge: NIL]];
v321: REF VVertex = NEW [VVertex ← [vp: [0.0, 0.0, 0.0], anEdge: NIL]];
e12: REF Edge = ConnectPoints [pi1, pi2, NIL, NIL, v123, v321];
e23: REF Edge = ConnectPoints [pi2, pi3, e12.sym, NIL, v123, v321];
e31: REF Edge = ConnectPoints [pi3, pi1, e23.sym, e12, v123, v321];
nDV ← 0;
RETURN [NEW [DelaunayDiagram ←
[points: LIST [pi1, pi2, pi3],
viewer: NIL]]]
END;

BEGIN
-- creates the viewer
tipTable: TIPUser.TIPTable ← TIPUser.InstantiateNewTIPTable["COG.tip"];
viewerClass: ViewerClasses.ViewerClass ← NEW
[ViewerClasses.ViewerClassRec ←
[paint: PaintMe,   -- called whenever the Viewer should repaint
  notify: TipMe,   -- TIP input events
  modify: NIL,    -- InputFocus changes reported through here
  destroy: DestroyMe,   -- called before Viewer structures freed on destroy op
  copy: NIL,   -- copy data to new Viewer
  set: NIL,   -- set the viewer contents
  get: NIL,   -- get the viewer contents
  init: InitMe,   -- called on creation or reset to init data
  save: NIL,   -- requests client to write contents to disk
  scroll: NIL,   -- document scrolling
  icon: document,  -- picture to display when small
  tipTable: tipTable,  -- could be moved into Viewer instance if needed
  cursor: crossHairsCircle -- standard cursor when mouse is in viewer
  ]];
menu: Menus.Menu ← MakeMenu[];
Diag: REF DelaunayDiagram;
vData: ViewerData;

out.PutF["Registering Viewer Class...\n"];
ViewerOps.RegisterViewerClass[$VoroPlotter, viewerClass];
out.PutF["Initializing the Voronoi...\n"];
Diag ← MakeDelaunay[];
vData ← NEW[ViewerDataRec ← [Diag: Diag]];
out.PutF["Creating Viewer...\n"];
viewer ← Diag.viewer ← ViewerOps.CreateViewer
[flavor: $VoroPlotter,
info:
[menu: menu,
name: "Voronoi Diagram",
column: left,
iconic: FALSE,
data: vData],
paint: TRUE]
END;

BEGIN
ENABLE ABORTED => GO TO bye;
vData: ViewerData ← NARROW [viewer.data];
WHILE TRUE DO
WaitStatusChanged;
IF NOT vData.live THEN GO TO bye;
vData.extraThings ← NIL;
ViewerOps.PaintViewer [viewer, client, FALSE, NIL];
IF vData.action = $LeftDown THEN
{LockAndInsertInVoronoi[vData.pt, vData.Diag]}
ELSE IF vData.action = $RightDown THEN
{out.PutF["Sorry, can't delete yet.\n"]}
ELSE IF vData.action = $Repaint THEN
{}
ELSE IF vData.action = $Sort THEN
{LockAndSortRegions [vData.Diag]}
ELSE IF vData.action = $Throw THEN
{n: NAT;
out.PutF ["How many points? "];
n ← in.GetInt[]; out.PutF["\n"];
THROUGH [1..n] DO
LockAndInsertInVoronoi [ThrowAPoint [], vData.Diag]
ENDLOOP;
};
vData.action ← NIL
ENDLOOP;

EXITS bye =>
{Close[in]; Close[out]};
END;