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