-- 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: BOOL ← TRUE;
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:
BOOL ←
TRUE,
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:
BOOL ←
FALSE, drawNo:
BOOL ←
TRUE,
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: BOOL ← TRUE;
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:
BOOL ←
TRUE] =
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...