DIRECTORY
  Rope USING [ROPE],
  IO USING [STREAM, CreateTTYStreams, PutF, int, real, Close, GetInt],
  Random USING [Init, Next],
  List USING [Comparison, Sort, Length, equal, less, greater, IsAListOfRefAny],
  Real USING [SqRt, LargestNumber, Float],
  COGCart,
  COGHomo,
  COGVoronoi,
  COGViewers USING [];
COGVoronoiImpl: CEDAR PROGRAM
  IMPORTS IO, Random, COGHomo, Real, List, COGViewers =
BEGIN
  OPEN Real, Rope, IO, Homo: COGHomo, Cart: COGCart, Vw: COGViewers, Vor: COGVoronoi;
SortMe: Menus.MenuProc = TRUSTED
BEGIN
?? rewrite!
  -- parameters: [viewer: Viewer, clientData: REF ANY, redButton: BOOL]
  vData: ViewerData _ NARROW [viewer.data];
  out.PutF["Called SortMe...\n"];
  vData.action _ $Sort;
  TellMother[]
END;

SortRegions: INTERNAL PROC [Diag: REF DelaunayDiagram, trace: BOOL]
  RETURNS [sortedPts: LIST OF REF DVertex] = TRUSTED
BEGIN
  -- Sorts the finite vertices of a Delaunay diagram in such a way that, for any two
  --   adjacent points, the leftmost one occurs first. When this rule (and its transitive
  --   closure) are insufficient to decide the ordering of two points, the one with lowest
  --   y occurs first.
  -- The final ordering is given by the 'no' field of each vertex and by the list sortedPts.
  --   also, p.anEdge will point to the (finite) neighbor of p that follows p and has maximum y.
  -- Infinite data points will get p.no = -1, -2, -3, ... and will be omitted from the sorted list.
  -- The number of finite points is returned in n. maximum y.
  pts, finPts, stack, free, t: LIST OF REF DVertex;
  pBq, pBprev: BOOL;
  p, q: REF DVertex;
  order, infOrder: NAT;
  e0, ea, ePrev: REF DEdge;

  pts _ Diag.points;
  out.PutF ["SortRegions - first pass...\n"];
  WHILE pts # NIL DO
    p _ pts.first; pts _ pts.rest;
    p.nBlockers _ 0;
    -- Counts blocking neighbors of p
    -- Also puts in eTop and in p.anEdge the edge blocked by p with Dest of maximum y
    -- Also finds first vertex in the ordering and initializes the Qyu to it.
    e0 _ ea _ p.anEdge;
    ePrev _ VPrev [e0]; pBprev _ Homo.Precedes [p.dp, Dest[ePrev].dp];
    out.PutF ["Looking at point %g\n", int[p.no]];
    DO
       q _ Dest[ea];
       pBq _ Homo.Precedes [p.dp, q.dp];
       IF NOT pBq THEN
         {IF pBprev THEN p.anEdge _ ePrev;
          p.nBlockers _ p.nBlockers + 1};
       ePrev _ ea; pBprev _ pBq; ea _ VNext[ea];
       IF ea = e0 THEN EXIT
      ENDLOOP;
     -- delete p from points list and put in free list (or prime the stack if p is first in ordering)
     IF p.nBlockers = 0 THEN
       {stack _ LIST [p];
        out.PutF["Point %g is first.\n", int[p.no]]}
    ENDLOOP;
   
   -- topological sorting: removes a vertex from the stack,
   --  decrements ref.counts of blocked neighbors,
   --  pushes freed neighbors into stack (in proper y order)
   out.PutF ["SortRegions - second pass...\n"];
   order _ 0; sortedPts _ finPts _ NIL;
   infOrder _ -1;
   WHILE stack # NIL DO
     p _ stack.first;
     out.PutF["Looking at point %g\n", int[p.no]];
     IF trace THEN ShowRegion [Diag, p, thin];
     t _ stack; stack _ stack.rest;
     IF Homo.FinPt[p.dp] THEN
       {p.no _ order; order _ order + 1;
         IF sortedPts = NIL THEN sortedPts _ t ELSE finPts.rest _ t;
        finPts _ t; t.rest _ NIL}
     ELSE
       {p.no _ infOrder; infOrder _ infOrder - 1};
     e0 _ ea _ p.anEdge;
     WHILE Homo.Precedes [p.dp, (q _ Dest[ea]).dp] DO
       IF trace THEN ShowEdge [Diag, ea, bold];
       q.nBlockers _ q.nBlockers - 1;
       IF q.nBlockers = 0 THEN
         {stack _ CONS [q, stack]};
       ea _ VPrev[ea]; -- note: cw order here
       IF ea = e0 THEN EXIT
      ENDLOOP;
    ENDLOOP;
   Diag.points _ pts;
   out.PutF ["Exit SortRegions.\n"]
END;
YTest: TYPE = RECORD
     [testY: REAL,
      up, dn: REF -- to YTest or EdgeTest or DVertex
      ];
EdgeTest: TYPE = RECORD
     [edge: REF DEdge,
      left, right: REF -- to YTest or EdgeTest or DVertex
      ];
TreeRoot: TYPE = REF ANY; -- root of a Lee-Preparata tree: EdgeTest, YTest, or DVertex
BuildLeePreparataStructure: INTERNAL PROC[Diag: REF DelaunayDiagram, trace: BOOL]
  RETURNS [leep: REF] = TRUSTED
BEGIN
  -- Builds the Lee-Preparata tree structure for point location
  
  -- Each chain is at first collected as list of edges; these are sorted by ascending y coordinate, and a chain tree is built from them. Since every edge is represented in only one chain (of highest possible level), there may be one or more gaps in the ranges of y-values covered by a given chain. A DummyEdge is temporarily inserted in each of those gaps (but not in the final structure).  
  
  DummyEdge: TYPE = RECORD
     [topY: REAL,
      onLeft: BOOL -- TRUE if gap corresponds to edge(s) of a left ancestor of the chain
      ];

  -- Once the Voronoi regions have been sorted from left to right, each edge can be assigned to a specific chain based on the ranks of the two regions on each side of it. If the regions are numbered in the range [0..n) (n>1), the root chain separates regions [0..k) and [k..n), where k = n DIV 2 is the chain index. Each group of regions is split in half by a child of chain k, and so forth; in generla, if at some step the set of regions is [fst..lim) (with lim > fst+1), the root chain will separate regions [fst..r) from [r..lim), where r = (fst+lim) DIV 2. The chain indices will lie therefore in the range (0..n).
    -- Given two adjacent regions i and j, the edge between them will belong to a unique chain whose index, LCA[i, j], is computed by the following procedure.
  ChainIndex: PROC [i, j: NAT] RETURNS [ix: NAT] = TRUSTED
BEGIN
  -- Returns the index of the monotone chain that separates region i from region j (i#j)
  --   and has highest possible level in the Lee-Preparata tree.  
  fst: NAT _ 0;
  lim: NAT _ n;
  DO
    -- assert: i and j IN [fst..lim)
    ix _ (fst + lim)/2;
    IF i < ix THEN
      IF j < ix THEN
        lim _ ix
      ELSE
        RETURN
    ELSE
      IF j < ix THEN
        RETURN
      ELSE
        fst _ ix
   ENDLOOP
END;

  Descendants: PROC [ix: NAT] RETURNS [fst, lim: NAT] = TRUSTED
BEGIN
  -- Returns fst and lim such that LCA[i, j] = ix iff i IN [fst..ix) and j IN [ix..lim).
  av: NAT;
  fst _ 0; lim _ n;
  DO
    -- assert: ix IN [fst..lim)
    av _ (fst + lim)/2;
    IF ix < av THEN lim _ av
    ELSE IF ix > av THEN fst _ av
    ELSE RETURN
   ENDLOOP
END;
  -- The folowing records are temporary data structures.
   
  ChainTable: TYPE = RECORD
     [elist: SEQUENCE nLists: NAT OF DEdgeList
      ];
  TipsTable: TYPE = RECORD
     [tip: SEQUENCE nRegs: NAT OF REF VVertex
      ];
  RegionTable: TYPE = RECORD
     [dvertex: SEQUENCE nRegs: NAT OF REF DVertex
      ];

  n: NAT; -- number of finite data points in the diagram
  edges: REF ChainTable; -- edges in each Lee-Preparata chain
  tips:  REF TipsTable;  -- upper Voronoi vertices of each region
  regions: REF RegionTable; -- upper Voronoi vertices of each region
  
  SetUpTables: INTERNAL PROC = TRUSTED
BEGIN
  -- Sorts the regions from left to right and sets the chain table in such a way that
  --   an edge between region i and region j is in edges[ChainIndex[i, j]].
  -- Returns also in regions.dvertex[k] a pointer to region k, and in
  --   regions.tip[k] the its point with maximum y.
  pts: DVertexList _ SortRegions [Diag, FALSE];
  p, q: REF DVertex;
  e, e0: REF DEdge;
  highest: BOOL;
  ix: NAT;
  n _ List.Length[List.IsAListOfRefAny[pts].list];
  edges _ NEW [ChainTable[n]];
  tips _ NEW [TipsTable[n]];
  regions _ NEW [RegionTable[n]];
  -- enumerate the edges of the delaunay
  WHILE pts # NIL DO
    p _ pts.first; pts _ pts.rest;
    e _ e0 _ p.anEdge;
    regions[p.no] _ p;
    WHILE e # e0 DO
      q _ Dest[e];
      IF p.no < q.no THEN
        {ix _ ChainIndex [p.no, q.no];
         IF Homo.Below[tips[p.no].vp, LeftVV[e].vp] THEN
           {tips[p.no] _ LeftVV[e]};
         edges[ix] _ CONS [e, edges[ix]]};
      e _ VNext[e]
     ENDLOOP
   ENDLOOP
END;
 
  CookYCoord: PROC [p: Homo.Point] RETURNS [y: REAL] = TRUSTED
BEGIN
  -- Returns the cartesian y-coordinate of p, or plus/minus LasgestNumber (with correct sign) if p is infinite.

  RETURN
     [IF Homo.InfPt[p] THEN
        IF p.y > 0.0 THEN LargestNumber
        ELSE IF p.y < 0.0 THEN - LargestNumber
        ELSE 0.0
      ELSE p.y/p.w]
END;

  SortChainEdges: INTERNAL PROC [le: DEdgeList, fst, lim: NAT, leftTip: REF VVertex]
    RETURNS [sle: LIST OF REF ANY] = TRUSTED
BEGIN OPEN List;
  -- Sorts the edges of a chain by descending Y, filling the gaps with dummy edges. The chain is assumed to bisect the regions [fst..lim); these parameters are used to figure out whether gaps lie on the left or right boundary of the current region set.
  -- A Voronoi vertex lies on the left boundary iff it is adjacent to some region with sequence number lim or more.
  -- The leftTip parameter gives the Voronoi vertex with maximum y on the regions at the left of the chain. It is used as a separator when the chain is empty.
    
  DecreasingY: INTERNAL PROC [x: REF ANY, y: REF ANY] RETURNS [Comparison] = TRUSTED
     {xe: REF DEdge = NARROW [x];
      ye: REF DEdge = NARROW [y];
      xlv: Homo.Point = LeftVV[xe].vp;
      ylv: Homo.Point = LeftVV[ye].vp;
      RETURN
        [IF Homo.Below[ylv, xlv] THEN less ELSE
         IF Homo.Below[xlv, ylv] THEN greater ELSE equal]
       };
         
  OnLeftBoundary: INTERNAL PROC [v: REF VVertex] RETURNS [is: BOOL] = TRUSTED
     {e0: REF DEdge = v.anEdge;
      e: REF DEdge _ e0;
      DO
        IF Dest[e].no > lim THEN RETURN [TRUE];
        e _ FNext[e];
        IF e = e0 THEN RETURN [FALSE]
       ENDLOOP
      };
  
  tle: LIST OF REF ANY _ List.Sort [IsAListOfRefAny[le].list, DecreasingY];
  t: LIST OF REF ANY;
  e: REF DEdge;
  de: REF DummyEdge;
  lastV: Homo.Point _ [0.0, 1.0, 0.0]; -- y = plus infinity
  eTop: Homo.Point;
  
  sle _ NIL;
  IF tle = NIL THEN
    {-- no proper edges in chain - put a dummy edge from top to y separator
      de _ NEW [DummyEdge _ 
              [topY: CookYCoord[lastV],
               onLeft: TRUE]];
           sle _ CONS[de, sle];
           lastV _ leftTip.vp}
  ELSE
    {WHILE tle # NIL DO
        t _ tle; e _ NARROW [t.first]; tle _ t.rest;
        eTop _ LeftVV[e].vp;
        IF Homo.Below [eTop, lastV] THEN
          {-- Found a gap
           de _ NEW [DummyEdge _ 
              [topY: CookYCoord[lastV],
               onLeft: OnLeftBoundary[LeftVV[e]]]];
           sle _ CONS[de, sle];
           lastV _ eTop};
        t.rest _ sle; sle _ t;
        lastV _ RightVV[e].vp;
       ENDLOOP};
  -- Put last dummy edge if needed
  IF Homo.Below [[0.0, -1.0, 0.0], lastV] THEN
     {de _ NEW [DummyEdge _ 
              [topY: CookYCoord[lastV],
               onLeft: FALSE]];
           sle _ CONS[de, sle]}
END;

  Error: ERROR [what: ROPE] = CODE;

  BuildPartialTreeFromChain: INTERNAL PROC
        [sle: LIST OF REF ANY, ne: NAT, left, right: TreeRoot]
    RETURNS
        [root: TreeRoot, chTop: REAL, rest: LIST OF REF ANY] = TRUSTED
BEGIN
  -- Builds a search tree for the first ne edges (including dummy ones) in the sorted
  --   list sle. The parameters left and right give the roots of the childern chains on
  --   either side. Returns also the y coordinate chTop of the upper end of the partial chain,
  --   and the list of unused edges.
  IF ne = 1 THEN
    {WITH sle.first SELECT FROM
        de: REF DummyEdge => 
          {RETURN
              [root: IF de.onLeft THEN right ELSE left,
               chTop: de.topY,
               rest: sle.rest]};
        e: REF DEdge =>
           {RETURN
              [root: NEW [EdgeTest _ [edge: e, left: left, right: right]],
               chTop: CookYCoord[LeftVV[e].vp],
               rest: sle.rest]};
       ENDCASE => ERROR Error["UFO in chain."]}
   ELSE
     {up, dn: TreeRoot;
      dnTop: REAL;
      [dn, dnTop, sle] _ BuildPartialTreeFromChain [sle, ne/2, left, right];
      [up, chTop, rest] _ BuildPartialTreeFromChain [sle, ne - ne/2, left, right];
      root _ NEW [YTest _ [testY: dnTop, up: up, dn: dn]];
      RETURN}
END;

  TreeData: TYPE = RECORD [root: TreeRoot, tip: REF VVertex];
 
  BuildTreeOfChains: INTERNAL PROC [fst, lim: NAT]
    RETURNS [rslt: TreeData] = TRUSTED
BEGIN
  -- Builds a Lee-Preparata tree for regions [fst..lim).

  IF fst+1 = lim THEN
    {RETURN [[regions[fst], tips[fst]]]}
  ELSE
    {rtix: NAT = (fst + lim)/2;
     left: TreeData = BuildTreeOfChains [fst, rtix];
     right: TreeData = BuildTreeOfChains [rtix, lim];
     sle: LIST OF REF ANY = SortChainEdges [edges[rtix], fst, lim, left.tip];
     ne: NAT = List.Length[sle];
     RETURN [[root: BuildPartialTreeFromChain[sle, ne, left.root, right.root].root,
                 tip: IF Homo.Below[right.tip.vp, left.tip.vp] THEN left.tip ELSE right.tip]]}
END;

  SetUpTables[];
  RETURN [BuildTreeOfChains [0, n].root]
END;
 
  [in,out] _ CreateTTYStreams["Voronoi.TTY"];
  [] _ Random.Init[,-1];

  out.PutF["Hello!\n"];
  Process.Detach[FORK Mother[MakeViewer[]]]
END.



Ί-- COGVorLeeP.mesa: Lee/Preparata Point location algorithm
-- last modified by Stolfi August 23, 1982 5:37 pm
-- To do: Update for COGVoronoi, COGDraw
-- To do: Throw: repaint diagram after each point
-- To do: Throw: use a better range (now is astronomical)
-- To do: check locking scheme
-- compile COGVoronoiImpl

-- to run:
--   run RandomImpl
--   run COGHomoImpl
--   run COGCartImpl
--   run COGVoronoiImpl

-- LEFT-TO-RIGHT SORT OF VORONOI REGIONS - - - - - - - - - - - - - - - - - - - - - - - - - 
-- THE LEE-PREPARATA ALGORITHM - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-- The Lee-Preparata structure is a binary search tree, whose leaves are the Voronoi regions and whose internal nodes are monotone chains of Voronoi edges. Each chain is by itself a binary search tree: the nodes of the latter are of two kinds, called YTests and EdgeTests.
-- YTest nodes branch depending on whether the query point q is above or below the testY field. When an EdgeTest is encountered, q.y is known to be between those of the edge endpoints; that being the case, q.x can be tested against the edge.
-- A node of either kind has two pointers, which can refer to Voronoi regions (ultimate leaves), to its children in the same chain C, or to the roots of the chains that are children of C in the global tree. The latter is always the case for EdgeTests, and sometimes also for the YTests. The latter case occurs when the edge of C that covers q.y belongs also to some ancestor C' of C. At this point, we know on which side of C' the query point q lies, and therefore we do not have to test it again against that edge; we can proceed immediately to the root of the appropriate descendant chain C'' of C.
Κ(–
"Mesa" style˜IprocšΟcΠbc(™:š
2™2KšœΟbœ ™(KšœŸœ)™1KšœŸœ1™9KšœŸœ™—Kš
™KšžX™aKšΟk	œ œ œ œ
 œ œ> œ œK œP œ˜ΈKš
Ÿœ œ
 œ œ
 œ+˜UKš œ œ
 œ@˜[šŸœ ˜ Kš
 œ
ŸœFœ œY œ
˜Σ
—K˜Kš
[™[šΟnœ œ
 œ œ œ œ
 œ
 œ
 œ ˜xKšˆ
 œSœXœYœœ[œ_œbœ<œ œ
 œ
 œ œ
 œ œ œN œ œ
 œ="œRœJœ–
 œF œ
 œ œ œ œu œ	 œ
 œ œaœ œ œ œ? œ	9œ0œ:œS œ œ	 œ
 œO œ œF œ œ3 œ
 œ
 œ œ/ œ œR œ) œ œ œI œ œ œ%œ œ	 œ
 œ œ œ< œ
˜·—š
\™\Kš
™Kš;ž
EžJž ™ρ
Kš
Ψ™Ψ—KšŸœ œ œ œ œ
#œ˜eKšŸœ œ œ
 œ œ
#œ˜qKš	Ÿœ œ œ
 œ<˜Vš‘œ œ
 œ œ œ œ œ ˜qKš œ>œ†œ œ œ
 œ œ
Fœιœš
˜ν
šœΠbn
œ œ œ œ œ ˜:Kš( œWœAœ œ
 œ œ!œ œ œ œ œ œ	 œ œ œ œ	 œ œ œ
 œ
˜―—šœ‘œ œ œ œ œ ˜@Kš œWœ œ œœ œ	 œ œ
 œ	 œ œ
 œ œ
 œ
˜™—Kš5œ7œ œ œ œ	 œ
 œ! œ œ œ œ
 œ
 œ! œ œ œ œ
 œ
 œ œ/œ	 œ
%œ	 œ
)œ œ)œ˜Όšœ‘œ œ
 œ ˜&Kš6 œTœJœDœ2œ( œ œ œ œ œ? œ œ œ'œ œ œ
 œV œ œ œ
 œ1 œ) œ; œ* œ œ
 œ
˜Χ—š
œ‘
œ œ œ œ ˜@Kš œnœ œ œ œ	 œ œ œ
 œ œ œ œ
 œ
˜œ—šœ‘œ œ
 œ œ œ œ œ
 œ
 œ
 œ ˜‚
Kš’
 œ
 œ	ϋ
œrœ
œ‘œ œ
 œ œ
 œ œ
 œ œ œ œ	 œ œ	 œZ œ
 œ œ œ
 œ œ	 œ‘œ œ
 œ œ
 œ œ œ œ œ œ	 œ œ
 œ œ! œ œ
 œ œ	 œ œ
 œ
 œ
 œ; œ
 œ
 œ
 œ œ œ3œ œ œ œ
 œGœ œN œ œ, œ œ œ
 œ œ? œ œœ œ| œj œ!œ œ& œ œN œ œ œ
˜―—Kš	œŸœ œ œ œ
˜$šœ‘œ œ
 œ œ
 œ
 œ
 œ œ œ! œ œ
 œ
 œ
 œ ˜½
Kš2 œTœVœ]œ#œ œ œ œ œ
 œ
 œ œ œ œ œR œ œ œ‹
 œ œ œ& œ―
 œ1 œ œ
˜¨—Kš	œŸœ œ œ œ
˜>šœ‘œ œ
 œ œ œ ˜[Kš" œ7œ œ
 œ œ  œ œ‡
 œ
 œ
 œ
 œ> œ œ_ œ' œ
 œ œ
˜·—Kšœ œ! œ
˜?—Kšœs œ œ˜”
J˜—…—2κAΜ