-- COGVoronoiImpl.mesa: Voronoi/Delaunay Diagrams
-- last modified by Stolfi - October 16, 1982 7:30 pm
-- To do: make CheckDelaunay and CheckConvexFaces more demanding
-- To do: fix RegionPainterin voronoi style - take care of infinite vertices!
-- To run: run COGAll; run COGVoronoiImpl
DIRECTORY
Rope USING [ROPE],
COGDebug USING [Error, out],
IO USING [PutF, PutFR, int],
Real USING [SqRt],
Graphics USING [Context, SetColor, DrawRope, DrawArea, DrawStroke, Path, NewPath, black],
GraphicsColor USING [red],
COGCart USING [ScaleFactors, Vector, Length, Mul],
COGHomo USING
[Point, CircumCenter, Vector, Add, Sub, FinPt, InfPt, Way, CCW, VecDir,
InCircle],
COGDrawing USING
[PainterProc, Color, Size, DrawSegment, GetProp, ObjectParms, MoveTo, LineTo,
SetCP, DrawDot],
COGDiagram USING
[DEdge, MakeSphere, Sym, Dual, Org, Dest, Left, Right, VisitProc, No, EdgeRec, ONext,
PutE, PutLRing, Traverse, MakeBridge, AddVertex, CloseFace, LNext, LPrev, RNext],
COGVoronoi;
COGVoronoiImpl: CEDAR PROGRAM
IMPORTS
IO, COGDebug, COGHomo, Real, Graphics, COGDrawing, COGDiagram, COGCart
EXPORTS
COGVoronoi =
BEGIN
OPEN
Real, Rope, IO, Homo: COGHomo, Cart: COGCart, Draw: COGDrawing, Diag: COGDiagram,
Bug: COGDebug, Gr: Graphics, COGVoronoi, GrC: GraphicsColor;
-- PROCEDURES FOR CREATING AND MODIFYING DELAUNAY DIAGRAMS - - - - - - - - - - -
Org: PUBLIC PROC [e: Diag.DEdge] RETURNS [d: DVertex] =
-- Origin DVertex of e.
{RETURN [NARROW[Diag.Org[e]]]};
Dest: PUBLIC PROC [e: Diag.DEdge] RETURNS [d: DVertex] =
-- Destination DVertex of e.
{RETURN [NARROW[Diag.Dest[e]]]};
Left: PUBLIC PROC [e: Diag.DEdge] RETURNS [v: VVertex] =
-- Voronoi vertex at left of e.
{RETURN [NARROW[Diag.Left[e]]]};
Right: PUBLIC PROC [e: Diag.DEdge] RETURNS [v: VVertex] =
-- Voronoi vertex at right of e.
{RETURN [NARROW[Diag.Right[e]]]};
TrivialDelaunay: PUBLIC PROC [pt: Homo.Point] RETURNS [dlau: Delaunay] =
BEGIN
d: DVertex = NEW [VertexRec ← [pt: pt, no: 0]];
v: VVertex = NEW [VertexRec ← [pt: [0, 0, 0], no: 0]];
Bug.out.PutF["\nVoronoi.TrivialDelaunay"];
RETURN [NEW [DelaunayRec ←
[dg: Diag.MakeSphere [v: d, f: v],
nDV: 1, nVV: 1,
limDV: 1, limVV: 1]
]]
END;
InfiniteTriangle: PUBLIC PROC RETURNS [dlau: Delaunay] = TRUSTED
BEGIN
sqrt3: REAL = SqRt[3.0];
pi0: DVertex = NEW [VertexRec ← [pt: [0.0, sqrt3/3.0, 0.0], no: 0]];
pi1: DVertex = NEW [VertexRec ← [pt: [-0.5, -sqrt3/6.0, 0.0], no: 1]];
pi2: DVertex = NEW [VertexRec ← [pt: [0.5, -sqrt3/6.0, 0.0], no: 2]];
v0: VVertex = NEW [VertexRec ← [pt: [0, 0, 0], no: 0]];
v1: VVertex = NEW [VertexRec ← [pt: [0, 0, 1], no: 1]];
e01: Diag.DEdge = Diag.MakeBridge [pi0, pi1, v0];
e12: Diag.DEdge = Diag.AddVertex [pi2, e01, left];
e20: Diag.DEdge = Diag.CloseFace [e12, e01, v1, left];
Bug.out.PutF["\nVoronoi.InfiniteTriangle: edges are "]; Diag.PutLRing[e01];
RETURN [NEW [DelaunayRec ←
[dg: e12,
nDV: 3, nVV: 2,
limDV: 3, limVV: 2]
]]
END;
VVertexCoords: PUBLIC PROC [p, q, r: Homo.Point] RETURNS [vPt: Homo.Point] =
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.CCW [p, q, r]
THEN {vPt ← Homo.CircumCenter [p, q, r]}
ELSE {dir: Homo.Vector ← Homo.Sub[q, p];
vPt ← IF dir.w < 0.0 THEN [dir.y, -dir.x, 0.0] ELSE [-dir.y, dir.x, 0.0]}
END;
-- PROCEDURES FOR DRAWING DELAUNAY DIAGRAMS - - - - - - - - - - - - - - - - - - - - -
DelaunayPainter: PUBLIC Draw.PainterProc =
BEGIN
dlau: Delaunay = NARROW [parms.data];
voronoi: BOOL ← Draw.GetProp[dr, $Voronoi] # NIL;
enum: BOOL ← Draw.GetProp[dr, $ENum] # NIL;
vnum: BOOL ← Draw.GetProp[dr, $VNum] # NIL;
prm: Draw.ObjectParms = parms;

DEdgeVisit: Diag.VisitProc =
{-- Bug.out.PutF["\nVoronoi.Painter: Visited edge "]; Diag.PutE[e]; -- -- DEBUG --
DrawEdge [context, sf, e, voronoi, enum, prm.color, prm.size]};

DVertexVisit: Diag.VisitProc =
{d: DVertex = NARROW [Diag.Org[e]];
DrawVertex [context, sf, d, vnum, prm.color, prm.size]};

[] ← Diag.Traverse [dlau.dg, DVertexVisit, NIL, vertices];
[] ← Diag.Traverse [dlau.dg, DEdgeVisit, NIL, oneWay]
END;
VertexPainter: PUBLIC Draw.PainterProc =
BEGIN
v: REF VertexRec = NARROW [parms.data];
vnum: BOOL ← Draw.GetProp[dr, $VNum] # NIL;
DrawVertex [context, sf, v, vnum, parms.color, parms.size]
END;
RegionPainter: PUBLIC Draw.PainterProc = TRUSTED
BEGIN
erec: REF Diag.EdgeRec = NARROW [parms.data];
eix: NAT [0..4) = parms.style;
e: Diag.DEdge = [erec, eix];
voronoi: BOOL ← Draw.GetProp [dr, $Voronoi] # NIL;
vnum: BOOL ← Draw.GetProp [dr, $VNum] # NIL;
path: Graphics.Path ← Graphics.NewPath[10];
IF voronoi THEN
{et: Diag.DEdge ← e;
fst, break: BOOLTRUE;
IF NOT Homo.FinPt[Org[e].pt] THEN RETURN;
DO
IF Homo.FinPt[Dest[et].pt] THEN
{IF NOT (Homo.FinPt[Left[et].pt] OR Homo.FinPt[Right[et].pt]) THEN
{-- two-point diagram - should paint correct half-plane
-- ??? --}
ELSE
{IF break THEN
{IF fst THEN Draw.MoveTo [path, sf, Right[et].pt]
ELSE Draw.LineTo [path, sf, Right[et].pt]};
Draw.LineTo [path, sf, Left[et].pt];
fst ← break ← FALSE}
}
ELSE
{break ← TRUE};
et ← Diag.ONext [et];
IF et = e THEN EXIT
ENDLOOP;
Graphics.SetColor [context, parms.color];
Graphics.DrawArea [context, path];
Graphics.DrawStroke [context, path]}
ELSE
{Draw.DrawDot [context, sf, Org[e].pt, parms.size, parms.color]};
IF vnum THEN
{Draw.SetCP [context, sf, Org[e].pt];
Gr.SetColor[context, GrC.red];
Gr.DrawRope[context, PutFR["%g", int[Org[e].no]]]}
END;
EdgePainter: PUBLIC Draw.PainterProc =
BEGIN
erec: REF Diag.EdgeRec = NARROW [parms.data];
eix: NAT [0..4) = parms.style;
voronoi: BOOL ← Draw.GetProp [dr, $Voronoi] # NIL;
enum: BOOL ← Draw.GetProp [dr, $ENum] # NIL;
DrawEdge [context, sf, [erec, eix], voronoi, enum, parms.color, parms.size]
END;
DrawVertex: PROC
[context: Gr.Context, sf: Cart.ScaleFactors, v: REF VertexRec,
drawNo: BOOLTRUE,
color: Draw.Color ← Gr.black, side: Draw.Size ← 0] = TRUSTED
BEGIN
Draw.DrawDot [context, sf, v.pt, side, color];
IF drawNo THEN
{Draw.SetCP[context, sf, v.pt];
Gr.SetColor[context, color];
Gr.DrawRope[context, PutFR["%g", int[v.no]]]}
END;
FiniteMidpoint: PROC [p, q: Homo.Point, eps: REAL] RETURNS [Homo.Point] =
-- If both p and q are finite, returns the midpoint of p and q. If one of them is infinite, returns a point eps away from the other, in the proper direction.
BEGIN
IF Homo.FinPt[p] AND Homo.FinPt[q] THEN
{RETURN [Homo.Way[0.5, p, q]]}
ELSE
{dir: Cart.Vector ← Homo.VecDir[Homo.Sub[q, p]];
dir ← Cart.Mul[eps/Cart.Length[dir], dir];
RETURN [IF Homo.FinPt[p] THEN
Homo.Add[p, [dir.x, dir.y, 1.0]]
ELSE
Homo.Sub [q, [dir.x, dir.y, 1.0]]]
}
END;
DrawEdge: PROC
[context: Gr.Context, sf: Cart.ScaleFactors, e: Diag.DEdge,
voronoi: BOOLFALSE, drawNo: BOOLTRUE,
color: Draw.Color ← Gr.black, width: Draw.Size ← 0] = TRUSTED
-- Draws e as a Delaunay or Voronoi edge (depending on voronoi) in the indicated context and with the given scale factors. If drawNo is TRUE, also draws the edge number by its midpoint.
BEGIN
org: DVertex = NARROW [Diag.Org[e]];
dest: DVertex = NARROW [Diag.Dest[e]];
mid: Homo.Point;
IF NOT voronoi AND (Homo.FinPt[org.pt] OR Homo.FinPt[dest.pt]) THEN
{mid ← FiniteMidpoint[org.pt, dest.pt, 0.418]; -- 0.418 here is just a guess
Draw.DrawSegment [context, sf, org.pt, dest.pt, width, color]}
ELSE IF voronoi AND (Homo.FinPt[org.pt] AND Homo.FinPt[dest.pt]) THEN
{lv: VVertex = NARROW [Diag.Left[e]];
rv: VVertex = NARROW [Diag.Right[e]];
IF Homo.InfPt[lv.pt] AND Homo.InfPt[rv.pt] THEN
{-- the diagram has only two (finite) data points
mid ← Homo.Way[0.5, org.pt, dest.pt];
Draw.DrawSegment [context, sf, lv.pt, mid, width, color];
Draw.DrawSegment [context, sf, mid, rv.pt, width, color]}
ELSE
{mid ← FiniteMidpoint[lv.pt, rv.pt, 0.418]; -- 0.418 here is just a guess
Draw.DrawSegment [context, sf, lv.pt, rv.pt, width, color]}
}
ELSE
{RETURN};
IF drawNo THEN
{Draw.SetCP[context, sf, mid];
Gr.SetColor[context, color];
Gr.DrawRope[context, PutFR["%g", int[Diag.No[e]]]]}
END;
RenumberVertices: PUBLIC PROC [dlau: Delaunay] =
BEGIN
doingDV: BOOLTRUE;
SetDVNo: Diag.VisitProc =
{dv: DVertex = NARROW [Org[e]];
dv.no ← dlau.limDV; dlau.limDV ← dlau.limDV + 1};
SetVVNo: Diag.VisitProc =
{dv: DVertex = NARROW [Org[e]];
dv.no ← dlau.limDV; dlau.limDV ← dlau.limDV + 1};
dlau.limDV ← dlau.minDV;
[] ← Diag.Traverse [dlau.dg, SetDVNo, NIL, vertices];
IF dlau.nDV # 0 THEN
{IF dlau.limDV - dlau.minDV # dlau.nDV THEN Bug.Error["lost DVertices"]}
ELSE
{dlau.nDV ← dlau.limDV - dlau.minDV};
dlau.limVV ← dlau.minVV; doingDV ← FALSE;
[] ← Diag.Traverse [Diag.Dual[dlau.dg], SetDVNo, NIL, vertices];
IF dlau.nVV # 0 THEN
{IF dlau.limVV - dlau.minVV # dlau.nVV THEN Bug.Error["lost VVertices"]}
ELSE
{dlau.nVV ← dlau.limVV - dlau.minVV}
END;
CheckDelaunay: PUBLIC PROC [dg: Diag.DEdge] =
BEGIN
DoCheckEdge: Diag.VisitProc =
{IF NOT (Homo.FinPt[Org[e].pt] AND Homo.FinPt[Dest[e].pt]) THEN
RETURN
ELSE
{enl: Diag.DEdge = Diag.LNext[e];
enr: Diag.DEdge = Diag.RNext[e];
conL: BOOL = Homo.CCW[Org[e].pt, Dest[e].pt, Dest[enl].pt];
conR: BOOL = Homo.CCW[Dest[e].pt, Org[e].pt, Org[enr].pt];
conU: BOOL = Homo.CCW[Org[enr].pt, Dest[e].pt, Dest[enl].pt];
conD: BOOL = Homo.CCW[Dest[enl].pt, Org[e].pt, Org[enr].pt];
IF conL AND conR
AND NOT Homo.InCircle [Dest[enl].pt, Org[e].pt, Org[enr].pt, Dest[e].pt] THEN
{Bug.out.PutF ["\nCheckDelaunay: edge "]; Diag.PutE[e];
Bug.out.PutF [" shouldn't be there"]};
IF conU AND conD
AND Homo.InCircle [Dest[e].pt, Dest[enl].pt, Org[e].pt, Org[enr].pt] THEN
{Bug.out.PutF ["\nCheckDelaunay: the partner of edge "]; Diag.PutE[e];
Bug.out.PutF [" should be there"]}
}
};
[] ← Diag.Traverse [dg: dg, Visit: DoCheckEdge, choice: all]
END;
CheckFaceConvexity: PUBLIC PROC [dg: Diag.DEdge, convex: BOOLTRUE] =
BEGIN
DoCheckFace: Diag.VisitProc =
{et: Diag.DEdge ← Diag.Dual[e]; -- Org[e] = Left[et] is the Delaunay face being visited
ea: Diag.DEdge ← Diag.LPrev[et];
DO
IF convex # Homo.CCW [Org[ea].pt, Org[et].pt, Dest[et].pt] THEN
{Bug.out.PutF ["\nCheckFaceConvexity: wrong face "]; Diag.PutLRing[et];
RETURN};
ea ← et; et ← Diag.LNext[et];
IF et = Diag.Dual[e] THEN EXIT
ENDLOOP
};
[] ← Diag.Traverse [dg: Diag.Dual[dg], Visit: DoCheckFace, choice: vertices]
END;
END...