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