-- COGVorRecImpl.mesa: Voronoi diagram by recursive split and merge
-- last modified by Stolfi, October 18, 1982 3:54 pm
-- To do: rethink sequence of painting/erasing/cleaning/showing etc..
-- To run: run COGAll; run COGVoronoiImpl; run COGVorTestImpl; run COGVorRecImpl
DIRECTORY
Rope USING [ROPE],
IO USING [PutF, int, GetChar],
List USING [CompareProc, UniqueSort, Length, Comparison],
Real USING [CompareREAL],
GraphicsColor USING [green, yellow, red, black],
COGDebug USING [in, out],
COGCart USING [Point, Vector],
COGHomo USING
[Point, CCW, InCircle, LineThrough, SameSide, PtPtDir, CircumCenter],
COGVoronoi USING
[Delaunay, DVertex, VVertex, VVertexCoords, Dest, Org, VertexRec,
VertexNo, DelaunayRec, MakeDelaunayObject, CheckDelaunay, CheckFaceConvexity],
COGDiagram USING
[DEdge, Sym, ONext, OPrev, LNext, RPrev, FixLeft, FixRight, AddVertex,
DeleteEdge, MakeSphere, MakeBridge, PutDiagram, PutORing, PutE,
EdgeSide, CloseFace, SetLeft, MergeOrSplitVertices],
COGDrawing USING
[Drawing, Color, Object, Size, PainterProc,
MouseProc, MenuProc, AddMouseAction, AddMenuAction, Make, Add,
DrawDot, DoPaint, RemoveAll, PutProp, GetProp, Remove],
COGVorTest USING [theDrawing, CleanTheDrawing, ShowEdge];
COGVorRecImpl: CEDAR PROGRAM
IMPORTS
IO, List, Real, COGHomo, COGDrawing, COGDebug, COGDiagram,
COGVoronoi, COGVorTest =
BEGIN
OPEN
Rope, GraphicsColor, IO,
Homo: COGHomo, Cart: COGCart, Draw: COGDrawing, COGVoronoi, Bug: COGDebug,
Diag: COGDiagram, COGVorTest;
THE POINT BAG - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-- The points are collected in a "point bag" (a list of points) until the user clicks the "Build" menu entry. The point bag is also an object of the main Drawing, and is displayed as a set of unnumbered dots.
PointBag: TYPE = LIST OF REF ANY; -- actually REF Cart.Point
PointBagPainter: Draw.PainterProc =
BEGIN
-- parameters: [dr: Drawing, context, sf: Cart.ScaleFactors, parms: ObjectParms]
bag: PointBag ← NARROW [parms.data];
p: REF Cart.Point;
WHILE bag # NIL DO
p ← NARROW [bag.first];
Draw.DrawDot [context, sf, [p.x, p.y, 1.0], parms.size, parms.color];
bag ← bag.rest
ENDLOOP
END;
AddPoint: Draw.MouseProc =
BEGIN
-- Called with write rights to the drawing
-- Parameters: [dr: Draw.Drawing, eventData: REF ANY, x, y: REAL, event: Draw.MouseEvent]
-- The eventData is NIL. GetProp[dr, $BagObj] is a PointBag object where the point is to be inserted
bagObj: Draw.Object = NARROW[Draw.GetProp[dr, $BagObj]];
bag: PointBag = NARROW[bagObj.parms.data];
p: REF Cart.Point = NEW[Cart.Point ← [x, y]];
bagObj.parms.data ← CONS [p, bag];
Draw.DoPaint[dr, bagObj]
END;
-- CONSTRUCTION OF THE DELAUNAY DIAGRAM - - - - - - - - - - - - - - - - - - - - - - - - -
BuildTheDelaunay: Draw.MenuProc =
BEGIN
-- Called with write access to the drawing
-- Parameters: [dr: Drawing, menuData: REF ANY, button: MouseButton]
-- The menuData is NIL. The current PointBag object is GetProp[dr, $BagObj].
bagObj: Draw.Object = NARROW [Draw.GetProp[dr, $BagObj]];
bag: PointBag ← NARROW [bagObj.parms.data];
nPoints: NAT; -- number of data points (excluding repetitions)
dg: Diag.DEdge;
dlObj: Draw.Object;
LeftToRight: List.CompareProc =
{p: REF Cart.Point = NARROW [ref1];
q: REF Cart.Point = NARROW [ref2];
xc: List.Comparison = Real.CompareREAL [p.x, q.x];
RETURN [IF xc = equal THEN Real.CompareREAL [p.y, q.y] ELSE xc]};
-- sort data points from left to right, removing duplicated points
bag ← List.UniqueSort [bag, LeftToRight];
nPoints ← List.Length[bag];
IF nPoints = 0 THEN RETURN;
-- erase previous diagram and repaint data points
Draw.RemoveAll[dr];
Draw.DoPaint[dr, bagObj];
-- put back the (emptied) point bag object, for future runs
bagObj.parms.data ← NIL;
Draw.Add [dr, bagObj, 3];
-- build a new Delaunay diagram by divide and conquer, and add it to the Drawing
[dg, , , dlObj] ← Voronoi[bag, 1, nPoints];
CheckDelaunay[dg];
Draw.PutProp [dr, $DelObj, dlObj]
END;
-- VORONOI CONSTRUCTION BY DIVIDE-AND-CONQUER - - - - - - - - - - - - - - - - - - - - -
ConvexQ:
PROC [p, q, r, s: Homo.Point]
RETURNS [b:
BOOLEAN] =
BEGIN
-- returns TRUE if p, q, r, s form a convex quadrilateral in that order (does NOT check for CCW)
OPEN Homo;
RETURN [ NOT (Homo.SameSide[p, r, LineThrough[q, s]] OR
Homo.SameSide[q, s, LineThrough[p, r]])
]
END;
ValidL:
PROC [lc, bl: Diag.DEdge]
RETURNS [
BOOLEAN] =
BEGIN
RETURN
[Homo.CCW [Dest[lc].pt, Dest[bl].pt, Org[bl].pt]
]
END;
ValidR:
PROC [rc, br: Diag.DEdge]
RETURNS [
BOOLEAN] =
BEGIN
RETURN
[Homo.CCW [Org[br].pt, Dest[br].pt, Dest[rc].pt]
]
END;
MakeNextVertex:
PROC [bag: PointBag, no: VertexNo]
RETURNS [p: DVertex, restBag: PointBag] =
-- Makes a new DVertex out of the first point in the bag, giving it the specified serial number.
BEGIN
cp: REF Cart.Point ← NARROW [bag.first];
p ← NEW [VertexRec ← [pt: [cp.x, cp.y, 1], no: no]];
restBag ← bag.rest
END;
Voronoi:
PROC [bag: PointBag, i, j:
NAT]
RETURNS [ld, rd: Diag.DEdge, restBag: PointBag, dlObj: Draw.Object] =
-- Computes the Voronoi diagram for the first n = j-i+1 points of the given bag. These points should be sorted by increasing x-coordinate, with ties broken by increasing y-coordinate. The returned restBag will contain the unused part of the bag. The Delaunay vertices will be numbered starting from i.
-- Returns also two edges ld and rd on the convex hull of the resulting Delaunay. The origin of ld is the leftmost data point, and it goes counterclockwise around the hull. The origin of rd is the rightmost point, and it goes clockwise around the hull. They will be dummy iff n=1.
BEGIN
OPEN Diag;
n: NAT = j - i + 1; -- number of vertices to consider
dlau: Delaunay;
-- Simple cases - direct construction
IF n <= 3
THEN
BEGIN
-- Simple cases - direct construction
[ld, rd, restBag] ← SmallVoronoi [bag, i, j]
END
ELSE
BEGIN
k: NAT = (i + j)/2; -- last vertex of first half
lcand, rcand, basel, baser, ldi, rdi, fl, fr, nl, nr, firstCrossEdge, t: DEdge;
dir: Cart.Vector; -- (temp) direction along some edge
vc: Homo.Point; -- (temp) coordinates of some Voronoi vertex
vv: VVertex; -- (temp) some Voronoi Vertex
lObj, rObj: Draw.Object;
-- General case - construction by divide and conquer
-- recursively compute the Voronoi of the left and right halves
[ld, ldi, restBag, lObj] ← Voronoi[bag, i, k];
[rdi, rd, restBag, rObj] ← Voronoi[restBag, k+1, j];
CleanTheDrawing[repaint: TRUE]; -- erase all debugging edge objects
-- compute the lower common tangent
DO
WHILE Homo.CCW [Org[ldi].pt, Dest[ldi].pt, Org[rdi].pt]
DO
ldi ← LNext[ldi];
ENDLOOP;
IF Homo.CCW [Dest[rdi].pt, Org[rdi].pt, Org[ldi].pt]
THEN
rdi ← RPrev[rdi]
ELSE EXIT -- now the edge from left to right is the "lower" common tangent
ENDLOOP;
-- add the left-right edge, with Voronoi vertex at infinity
dir ← Homo.PtPtDir [Org[ldi].pt, Org[rdi].pt];
vv ← NEW [VertexRec ← [pt: [dir.y,-dir.x, 0.0], no: 0]];
firstCrossEdge ← MakeBridge[Org[rdi], Org[ldi], vv];
Bug.out.PutF["\n>> firstCrossEdge = "]; PutE[firstCrossEdge];
Bug.out.PutF["from %g to %g", int[Org[firstCrossEdge].no], int[Dest[firstCrossEdge].no]];
-- PASS 1
-- The loop below removes all edges that have to be deleted from the left and right Delaunay, and builds a separate diagram with all the edges that have to be added between them.
basel ← firstCrossEdge; baser ← Sym[basel];
-- initialize lcand and rcand in preparation for the merge
lcand ← ONext[ldi]; rcand ← OPrev[rdi];
DO -- this is the merge pass
Bug.out.PutF["\n >> lcand is "]; PutE[lcand];
Bug.out.PutF["Dest = %g", int[Dest[lcand].no]];
-- find place to connect to on the left, using InCircle test
WHILE ValidL[lcand, basel] AND (ONext[lcand] # lcand)
AND ConvexQ[Org[lcand].pt, Org[rcand].pt, Dest[lcand].pt, Dest[ONext[lcand]].pt]
AND Homo.InCircle
[Dest[lcand].pt, Dest[ONext[lcand]].pt, Org[lcand].pt, Org[rcand].pt] DO
t ← ONext[lcand];
[] ← DeleteEdge[lcand]; ShowEdge[lcand, yellow, 2]; -- delete skipped edge
lcand ← t; -- update lcand
Bug.out.PutF[" >> now lc is "]; PutE[lcand];
Bug.out.PutF["Dest = %g", int[Dest[lcand].no]]
ENDLOOP;
Bug.out.PutF["\n >> rcand is "]; PutE[lcand];
Bug.out.PutF["Dest = %g", int[Dest[rcand].no]];
-- find place to connect to on the right
WHILE ValidR[rcand, baser] AND (rcand # OPrev[rcand])
AND ConvexQ[Org[lcand].pt, Org[rcand].pt, Dest[OPrev[rcand]].pt, Dest[rcand].pt]
AND Homo.InCircle
[Org[rcand].pt, Dest[OPrev[rcand]].pt, Dest[rcand].pt, Org[lcand].pt] DO
t ← OPrev[rcand];
[] ← DeleteEdge[rcand]; ShowEdge[rcand, yellow, 2]; -- delete skipped edge
rcand ← t; -- update rcand
Bug.out.PutF[" >> now rc is "]; PutE[lcand];
Bug.out.PutF["Dest = %g", int[Dest[rcand].no]]
ENDLOOP;
-- must choose the best candidate, if there is a choice
-- If last edge was the upper outer tangent, done!
IF (NOT ValidL[lcand, basel]) AND (NOT ValidR[rcand, baser]) THEN EXIT;
IF (NOT ValidL[lcand, basel])
OR ((ValidR[rcand, baser])
AND Homo.InCircle
[Dest[lcand].pt, Org[lcand].pt, Org[rcand].pt, Dest[rcand].pt]) THEN
BEGIN
-- connect to the right
-- compute Voronoi vertex
vc ← Homo.CircumCenter [Org[basel].pt, Dest[rcand].pt, Dest[basel].pt];
vv ← NEW [VertexRec ← [pt: vc, no: 0]];
FixRight[basel, vv]; FixLeft[rcand, vv];
ShowEdge[basel, black, 0]; ShowEdge[rcand, red, 2];
-- add new edge to connecting diagram and advance rcand
baser ← AddVertex[Dest[rcand], basel, right]; basel ← Sym[baser];
Bug.out.PutF["\n >> connected to right by edge "]; PutE[baser];
Bug.out.PutF["from %g to %g", int[Org[baser].no], int[Dest[baser].no]];
rcand ← LNext[rcand];
END
ELSE
BEGIN
-- connect to the left
-- compute Voronoi vertex
vc ← Homo.CircumCenter [Org[baser].pt, Dest[baser].pt, Dest[lcand].pt];
vv ← NEW [VertexRec ← [pt: vc, no: 0]];
FixLeft[baser, vv]; FixRight[lcand, vv];
ShowEdge[baser, black, 0]; ShowEdge[lcand, red, 2];
-- add new edge to connecting diagram and advance rcand
basel ← AddVertex[Dest[lcand], baser, left]; baser ← Sym[basel];
Bug.out.PutF["\n >> connected to left by edge "]; PutE[basel];
Bug.out.PutF["from %g to %g", int[Org[basel].no], int[Dest[basel].no]];
lcand ← RPrev[lcand];
END
ENDLOOP;
-- fix upper (infinite) Voronoi vertex of last edge
dir ← Homo.PtPtDir [Org[basel].pt, Org[baser].pt];
vv ← NEW [VertexRec ← [pt: [dir.y, -dir.x, 0.0], no: 0]];
FixRight [basel, vv];
ShowEdge[basel, black, 0];
PutDiagram[basel]; -- DEBUG --
-- PASS 2
-- This pass glues together the connecting edges and left and right Delaunays. The sequence of vertices to be connected is now implicit in the connecting diagram. The Voronoi vertices are already OK.
basel ← firstCrossEdge; baser ← Sym[basel];
lcand ← ONext[ldi]; rcand ← OPrev[rdi];
fl ← basel; fr ← baser;
nl ← OPrev[basel]; nr ← ONext[baser];
DO
IF nr = fr THEN
BEGIN -- stitch left vertex
t ← RPrev[lcand];
MergeOrSplitVertices[baser, OPrev[lcand]];
Bug.out.PutF["\nVorRec: stitched at left - vertex %g = ", int[Org[lcand].no]];
PutORing[lcand];
lcand ← t;
IF nl = fl THEN {MergeOrSplitVertices[fl, rcand]; EXIT};
basel ← nl; baser ← Sym[basel];
fr ← baser;
END
ELSE IF nl = fl THEN
BEGIN -- stitch right vertex
t ← LNext[rcand];
MergeOrSplitVertices[fl, rcand];
Bug.out.PutF["\nVorRec: stitched at right - vertex %g = ", int[Org[rcand].no]];
PutORing[rcand];
rcand ← t;
baser ← nr; basel ← Sym[baser];
fl ← basel;
END;
nl ← OPrev[basel]; nr ← ONext[baser];
ENDLOOP;
PutDiagram[firstCrossEdge];
IF Org[ld] = Dest[firstCrossEdge] THEN ld ← Sym[firstCrossEdge];
IF Org[rd] = Org[firstCrossEdge] THEN rd ← firstCrossEdge;
CheckFaceConvexity[ld];
Draw.Remove [lObj];
Draw.Remove [rObj];
END;
Build new Delaunay diagram and add it to theDrawing
dlau ← NEW [DelaunayRec ←
[dg: ld,
nDV: n,
nVV: 0, -- provisory
minDV: 0, limDV: n]];
dlObj ← MakeDelaunayObject [dlau];
Draw.Add [theDrawing, dlObj, 3];
-- end of Voronoi
END;
SmallVoronoi:
PROC [bag: PointBag, i, j:
NAT]
RETURNS [ld, rd: Diag.DEdge, restBag: PointBag] =
-- Computes the Voronoi diagram for the first n = j - i + 1 (1, 2 or 3) points in the given bag. These points should be sorted by increasing x with ties broken by increasing y. The returned restBag will contain the unused part of bag. The Delaunay vertices will be numbered sequentially starting from i.
-- Returns also two edges ld and rd on the convex hull of the resulting Delaunay. The origin of ld is the leftmost data point, and it goes counterclockwise around the hull. The origin of rd is the rightmost point, and it goes clockwise around the hull. They will be dummy iff n=1.
BEGIN
OPEN Diag;
n: NAT = j - i + 1; -- number of points to consider
IF n = 1 THEN-- one point
BEGIN
pi: DVertex;
[pi, restBag] ← MakeNextVertex[bag, i];
ld ← rd ← MakeSphere[pi, NIL];
RETURN
END;
IF n = 2 THEN -- two points
BEGIN
pi, pj: DVertex;
dij: Cart.Vector;
vl, vr: VVertex;
[pi, restBag] ← MakeNextVertex[bag, i];
[pj, restBag] ← MakeNextVertex[restBag, j];
-- build one-edge delaunay
ld ← MakeBridge[pi, pj, NIL];
rd ← Sym[ld];
-- compute voronoi vertices at infinity
dij ← Homo.PtPtDir [pi.pt, pj.pt];
vl ← NEW [VertexRec ← [pt: [-dij.y, dij.x, 0.0], no: 0]];
vr ← NEW [VertexRec ← [pt: [dij.y, -dij.x, 0.0], no: 0]];
FixLeft[ld, vl]; FixRight[ld, vr];
RETURN
END;
IF n = 3 THEN -- three points
BEGIN
e1, e2, e3, e: DEdge;
pi, pk, pj: DVertex;
dir: Cart.Vector;
vc: Homo.Point;
vv: VVertex;
[pi, restBag] ← MakeNextVertex[bag, i];
[pk, restBag] ← MakeNextVertex[restBag, i+1];
[pj, restBag] ← MakeNextVertex[restBag, i+2];
e1 ← MakeBridge[pi, pk, NIL];
e2 ← AddVertex[pj, e1, left]; -- or right, indifferent
IF Homo.CCW[pi.pt, pj.pt, pk.pt] THEN
{e3 ← CloseFace[e2, e1, NIL, right];
ld ← Sym[e3]; rd ← e3}
ELSE
{ld ← e1; rd ← Sym[e2];
IF Homo.CCW[pi.pt, pk.pt, pj.pt] THEN
{e3 ← CloseFace[e2, e1, NIL, left]}};
-- compute Voronoi vertex at center of triangle (if any)
IF ONext[ld] # ld THEN
{-- a proper triangle - compute center vertex
vc ← VVertexCoords[Dest[ONext[ld]].pt, Org[ld].pt, Dest[ld].pt];
SetLeft[ld, NEW [VertexRec ← [pt: vc, no: 0]]]};
-- compute Voronoi vertices at infinity
e ← ld;
DO -- this loop should work even if the triangle is degenerate
dir ← Homo.PtPtDir[Org[e].pt, Dest[e].pt];
vv ← NEW [VertexRec ← [pt: [dir.y, -dir.x, 0.0], no: 0]];
FixRight[e, vv];
e ← RPrev[e];
IF e = ld THEN EXIT
ENDLOOP;
RETURN
END;
END;
theBagObj: Draw.Object;
Bug.out.PutF["\nVorRec: Hello! Say something: "]; [] ← Bug.in.GetChar[];
theBagObj ← Draw.Make [PointBagPainter, [data: NIL, color: green, size: 2]];
Draw.PutProp[theDrawing, $BagObj, theBagObj];
Draw.AddMouseAction [theDrawing, [button: red], AddPoint, NIL, write];
Draw.AddMenuAction [theDrawing, "Build", BuildTheDelaunay, NIL, write];
Bug.out.PutF["\nVorRec: Enjoy yourself. Bye!"]
END.