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;

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

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] =
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] =
BEGIN OPEN Diag;
 n: NAT = j - i + 1;  -- number of vertices to consider
dlau: Delaunay;
 IF n <= 3 THEN
  BEGIN
    [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;
-- 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]];


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

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



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

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.
-- CONSTRUCTION OF THE DELAUNAY DIAGRAM - - - - - - - - - - - - - - - - - - - - - - - - - 
-- VORONOI CONSTRUCTION BY DIVIDE-AND-CONQUER - - - - - - - - - - - - - - - - - - - - -
-- Makes a new DVertex out of the first point in the bag, giving it the specified serial number.
-- 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.
-- Simple cases - direct construction
-- Simple cases - direct construction
-- General case - construction by divide and conquer
-- 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.

-- 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.
Build new Delaunay diagram and add it to theDrawing
-- 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.
�Êä–
"Mesa" style˜IprocšÏcÐbc�.™Cš
�4™4KšœÏbœ=™E—Kš�ž�J™SKš Ïk	œ œ œ œ
 œ œ7 œ  œ) œ œ œ
 œH œ¨
 œØ
 œ¿
 œ)˜ëKš
Ÿ
œ œ
 œ œ œU˜‚
Kš œ œ œr˜š
K˜š
�_™_Kš
�Ð
™Ð
—Kšœ
 œ œ
 œ
 œ
 œ�˜<˜#Kš œ�Qœ œ œ œ œ
 œ	 œn œ
 œ
˜²—šŸœ˜Kš œ�ž�œ�Zœ�ž	�	ž�ž�*œ œ/ œ œ œ- œ& œ
˜Ò—Kš
�Z™ZšŸœ˜!Kš> œ�ž�œ�Eœ�ž�ž�ž�ž�œ œ0 œ! œ�1œ-Ÿœ œ œ œ œF œ œ œ œ
�CœL œ
 œ
 œ�2œ7�<œ œ#�Qœg œ
˜ã	—K˜Kš
�W™WšÏnœ œ œ œ˜=Kš œ�aœ œ	 œ œ) œN œ
˜�—š¡œ œ œ œ˜5Kš œ œ œ/ œ
˜O—š¡œ œ œ œ˜5Kš œ œ œ. œ
˜N—š¡œ œ œ"˜\Kš
�`™`Kš
 œ œ œ œC œ
˜�
—š¡œ œ œ œ>˜oKš�.ž
�ž
�
ž
�ž�vž�%ž�7ž
�
™­Kš�ž�ž�=ž�Yž�Wž
�™˜š œ
 œ˜Kšœ œ�"œ˜GKš
�%™%šœ
 œ œ ˜Kš
�%™%—Kš
œ0˜0šœ œ
 œ ˜Kšœ œ�œb�$œ�-œ
�œ˜¸Kš
�4™4Kš�@œ~ œ�$œ
�$ œ œ œ- œ œ œ œ, œ œ
 œ
�Aœ
 œ
˜œKš�<œ4 œÿ
˜òK˜Kš�ž™	Kš
�²
™²
Kšfœ2�;œ) œ
�œf�=œ œ œ œS œ] œP�œ�œg œm�)œ œ œ œR œ] œO�œ�œg œ�8œ�3œ œ œ œ œ œ
 œ œ œ œ! œ[ œ œ�œ�œY œ�
�8œü
 œ œ œ�œ�œY œ�
�8œú
 œ œ�4œ8 œv�˜üK™Kš
ž	™	Kš
�Ç
™Ç
Kš)œ”
 œ œ	 œ œ
�œÌ
 œ	 œ# œ? œ œ
 œ	 œ œ
�œù
 œ+ œ œ  œ œ œU˜Î—Kšœ œ
˜Kš
�3™3Kšœ œ7�
œ_�˜¿
—Kš œ
˜—š¡œ œ œ œ*˜_Kš�.ž
�-ž�_ž�!ž�Dž
�
™¯Kš
�ž�ž�=ž�!œ�6ž�Wž
�™˜š œ
 œ˜Kšœ œ�˜3Kš
 œ �
œ œ[ œ œ œ
˜“
Kš œ œ
�œ œ¢
�œ œ�(œ0 œ; œ] œ œ
˜›Kš< œ œ
�œ œ› œ%�œ œ œ œ œ, œ& œ œ œ# œ�9œ œ œ	�-œ] œ&�(œ œ
�<œ< œd œ œ
 œ œ œ œ
˜Ù	—Kš œ
˜—Kš	œ—
 œ‡
 œG œ> œ˜²—…—2°CÈ