MODULE Map3D;
(*
   Facet-Edge data structure. Created on 1998 by Luis Arturo Perez Lozada, with
   modifications by J. Stolfi. (See notice of copyright at the end of this file.)
   Details based on the implementation of the Quad-Edge data structure 
   by J. Stolfi and R. Marcone Rosi.
*)


IMPORT Lex, Word, Wr, Rd, Fmt, FloatMode, Thread, Stdio;

FROM Stdio IMPORT stdout;

TYPE 
  MarkBit = BITS 8 FOR BOOLEAN;

REVEAL 
  FaceEdge = PublicFaceEdge BRANDED
  OBJECT 
    (* Component elements: *)
    edg: ARRAY VBit OF EFNode;      (* Component edge and face *)   

    (* Incident elements: *)
    rorg: ARRAY VEBits OF VPNode;   (* Vertex/cell data records for each handle. *)
    borg: ARRAY VEBits OF VBit;     (* Vertex/cell dualization bit. *)

    (* Topological links between fece-edges. *)
    rnext: ARRAY VEBits OF FaceEdge;  (* FaceEdge of nextF(a[r]) *)
    bnext: ARRAY VEBits OF VFEBits;   (* VFEBits of nextF(a[r]) *)
    
    emark: ARRAY VBit OF MarkBit;   (* Marks associated with edge/face component. *)
    vmark: ARRAY FEBits OF MarkBit;  (* Marks associated with vertices and cells. *)

 *)

    (* Edge and face records. *)
    (* Elements are indexed with primal/dual bit. *)

    indXX : CARDINAL := 0;             (* index of FaceEdge used by NumberFacets, etc. *)
  END;
  (*
    (* Element index "[2*d+e]" of fields "fe.fnext" and "fe.bnext" refers to 
      handle "h.Dual^d.RevE^e" where "h" is the reference handle for "fe".
      Handle{n,r}, *)
    (* for each rotations bitfield r = {0,1,2,3}.  Ditto for bnext and org. *)
    

  (*
    Meaning of bits for RevE and Srot:
    If "a = spin^s(srot^r(a0))", where "a0" is the reference handle
    of the fe underlying "a", with "r" in [0..3], and "s" in [0..1],
    then "a.bits = (2r + s)".   Therefore,

    | spin(a).bits = XOR(a.bits, 1);

    | If sbit=0 then srot(a).bits = (a.bits + 2) MOD 8 
    | otherwise      srot(a).bits = (a.bits + 6) MOD 8;

 The handle "nextF(a)" is given by "Handle{a.fe.fenext[rbits(a)],
 a.fe.bnext[rbits(a)]}", where rbits(a) = a.bits DIV 2,
 provided that the "spin" bit of "a" is 0.
 Otherwise "nextF(a)" is computed by 
 the formula: "nextF(a) = RevE(Flip(nextF(Flip(RevE(a)))))".

*)

(* ====== Facet-Edge rotation and reversal functions ====== *)

PROCEDURE Srot(a: Handle): Handle =
  BEGIN
    RETURN Handle{fe := a.fe, 
               bits := Word.And(a.bits + 2 + Word.LeftShift( 
               Word.And(a.bits, 1), 2), 7) };
  END Srot;

PROCEDURE RevE(a: Handle): Handle =
  BEGIN
    RETURN Handle{fe := a.fe, bits := Word.Xor(a.bits, 1)};
  END RevE;
  
PROCEDURE Flip(a: Handle): Handle =
  BEGIN
    RETURN Handle{fe := a.fe, bits := Word.And((a.bits + 4), 7)};
  END Flip;
  
PROCEDURE Tors(a: Handle): Handle =
  BEGIN
    RETURN Handle{fe := a.fe,
               bits := Word.And(a.bits + 6 + Word.LeftShift( 
                       Word.And(a.bits, 1), 2), 7) };    
  END Tors;

(* ====== Facet-Edge dual function ====== *)

PROCEDURE Sdual(a: Handle): Handle =
  BEGIN RETURN RevE(Srot(a)) END Sdual;

(* ====== Facet-Edge traversal functions ====== *)

PROCEDURE nextF(a: Handle) : Handle =

  PROCEDURE nextFBase(a: Handle): Handle =
    (* nextF for a handle with zero spin bit. *)
    BEGIN
      WITH rbits = Word.And(Word.RightShift(a.bits, 1), 3) DO 
        RETURN Handle{fe := a.fe.fenext[rbits], bits := a.fe.bnext[rbits]}
      END;
    END nextFBase;

  BEGIN
    IF Word.And(a.bits, 1) = 0 THEN RETURN nextFBase(a);
    ELSE RETURN RevE(Flip(nextFBase(Flip(a))));
    END;
  END nextF;

PROCEDURE nextE(a: Handle): Handle =
  BEGIN RETURN Srot(nextF(Tors(a))) END nextE;

PROCEDURE prevE(a: Handle): Handle =
  BEGIN RETURN Flip(Srot(nextF(Srot(a)))) END prevE;

PROCEDURE prevF(a: Handle): Handle =
  BEGIN RETURN Flip(nextF(Flip(a))) END prevF;

(* ====== QuadEdge emulation ========= *)

PROCEDURE PnegOnext(s: Handle): Handle =
  BEGIN RETURN Flip(prevF(prevE(s))) END PnegOnext;

PROCEDURE PnegOnext_1(s: Handle): Handle =
  BEGIN RETURN Flip(prevE(prevF(s))) END PnegOnext_1;

PROCEDURE PnegLnext(s: Handle): Handle =
  BEGIN RETURN nextE(s) END PnegLnext;

PROCEDURE PnegLnext_1(s: Handle): Handle =
  BEGIN RETURN prevE(s) END PnegLnext_1;

PROCEDURE PnegSym(s: Handle): Handle =
  BEGIN RETURN Flip(prevF(s)) END PnegSym;

(* ====== Splicing ====== *)

PROCEDURE SpliceFacets(a,b: Handle) =
  BEGIN
    <* ASSERT b # RevE(nextF(a)) *>
    IF a # b THEN
      WITH 
        ta = nextF(a),  tb = nextF(b),
        c =  Flip(ta),  d = Flip(tb), 
        tc = nextF(c),  td = nextF(d)
      DO
        IF Word.And(a.bits, 1) = 0 THEN 
          WITH
            rba = SrotBits(a) 
          DO
            a.fe.fenext[rba] := tb.fe;
            a.fe.bnext[rba] := tb.bits;
          END
        ELSE
          WITH
            ra = Flip(RevE(a)), fd = RevE(d), rbra = SrotBits(ra)
          DO
            ra.fe.fenext[rbra] := fd.fe;
            ra.fe.bnext[rbra] := fd.bits;
          END
        END;

        IF Word.And(b.bits, 1) = 0 THEN 
          WITH
            rbb = SrotBits(b)
          DO
            b.fe.fenext[rbb] := ta.fe;
            b.fe.bnext[rbb] := ta.bits;          
          END
        ELSE
          WITH
            rb = Flip(RevE(b)), fc = RevE(c), rbrb = SrotBits(rb)
          DO
            rb.fe.fenext[rbrb] := fc.fe;
            rb.fe.bnext[rbrb] := fc.bits;
          END
        END;

        IF Word.And(c.bits, 1) = 0 THEN 
          WITH
            rbc = SrotBits(c) 
          DO
            c.fe.fenext[rbc] := td.fe;
            c.fe.bnext[rbc] := td.bits;
          END
        ELSE
          WITH
            rc = Flip(RevE(c)), fb = RevE(b), rbrc = SrotBits(rc)
          DO
            rc.fe.fenext[rbrc] := fb.fe;
            rc.fe.bnext[rbrc] := fb.bits;
          END          
        END;

        IF Word.And(d.bits, 1) = 0 THEN
        WITH
          rbd = SrotBits(d)
        DO
          d.fe.fenext[rbd] := tc.fe;
          d.fe.bnext[rbd] := tc.bits;          
        END
        ELSE
          WITH
            rd = Flip(RevE(d)), fa = RevE(a), rbrd = SrotBits(rd) 
          DO
            rd.fe.fenext[rbrd] := fa.fe;
            rd.fe.bnext[rbrd] := fa.bits;
          END
        END;
      END;
    END;  
  END SpliceFacets;

PROCEDURE SpliceEdges(a,b: Handle) =
  BEGIN
    SpliceFacets(Sdual(a), Sdual(b));
  END SpliceEdges;
 
(* ====== Handy splicicing tools ====== *)

PROCEDURE Meld(a,b : Handle) =
  (* Delete the first handle fe argument *)
  VAR firsta, firstb, a1: Handle;
  BEGIN
    firsta := a; firstb := b;
    REPEAT
      <* ASSERT (a = firsta) = (b = firstb) *>
      a1 := prevF(a);
      IF nextF(a) # b THEN
        SpliceFacets(a, prevF(b));
        <* ASSERT nextF(a) = b *>
      END; 

      (* Delete the Handle "a" from its nextF ring: *)
      SpliceFacets(a, prevF(a));
      SpliceFacets(Flip(a), prevF(Flip(a))); 

      IF a1.fe # a.fe THEN <* ASSERT nextF(a1) = b *> END;
      a := nextE(a);
      b := nextE(b);
    UNTIL a = firsta;
    <* ASSERT b = firstb *>
  END Meld;

PROCEDURE SetnextF(a,b: Handle) =
  BEGIN
    IF nextF(a) # b THEN SpliceFacets(a, prevF(b)) END
  END SetnextF;

PROCEDURE SetnextE(a, b: Handle) =
  BEGIN
    IF nextE(a) # b THEN SpliceEdges(a, prevE(b)) END
  END SetnextE;

(* ====== Computing Bits ====== *)

PROCEDURE OrientationBit(a: Handle): EBit =
  BEGIN RETURN Word.RightShift(a.bits, 2) END OrientationBit;

PROCEDURE RevEBit(a: Handle): EBit =
  BEGIN RETURN Word.And(a.bits, 1) END RevEBit;
  
PROCEDURE DualBit(a: Handle): VBit =
  BEGIN RETURN Word.And(Word.RightShift(a.bits, 1), 1) END DualBit;

PROCEDURE SrotBits(a: Handle): FEBits =
  BEGIN RETURN Word.And( Word.RightShift(a.bits, 1), 3) END SrotBits;

(* ====== Counting ====== *)

PROCEDURE DegreeFaceRing(a: Handle): CARDINAL =
  VAR
    n: CARDINAL := 0;
    s: Handle := a;
  BEGIN
    REPEAT
      INC(n);
      s := nextF(s);
    UNTIL s = a;
    RETURN n
  END DegreeFaceRing;

PROCEDURE DegreeEdgeRing(a: Handle): CARDINAL =
  VAR
    n: CARDINAL := 0;
    s: Handle := a;
  BEGIN
    REPEAT
      INC(n);
      s := nextE(s);
    UNTIL s = a;
    RETURN n
  END DegreeEdgeRing;

(* ====== Creation ====== *)

PROCEDURE MakeRawFaceEdge (): Handle =
  BEGIN
    WITH n = NEW(FaceEdge) DO
      n.emark[0] := FALSE; n.emark[1] := FALSE; 
      FOR i := 0 TO 3 DO
    	n.vmark[i] := FALSE; n.fenext[i] := n; n.bnext[i] := 2*i
      END;
      RETURN Handle{fe := n, bits := 0};
    END;
  END MakeRawFaceEdge;

(* ======= Element records ======= *)

REVEAL
  ElemNode = PublicElemNode BRANDED OBJECT
    END;

(* ======= Edge/face data records ======= *)

REVEAL
  EFNode = PublicEFNode BRANDED OBJECT
      marksXXX: ARRAY [0..1] OF  MarkBit;  (* Used by NumberFacets, ex: emark[DualBit(a)] associated to Face *)
    END;
  
PROCEDURE Edg(a: Handle): EFNode =
  BEGIN
    WITH dbit = DualBit(a) DO RETURN a.fe.edg[dbit] END
  END Edg;

PROCEDURE SetEdgSingle(a: Handle; n: EFNode) =
  BEGIN
    WITH dbit = DualBit(a) DO a.fe.edg[dbit] := n END
  END SetEdgSingle;

PROCEDURE SetEdgAll(a: Handle; n: EFNode) =
  VAR t: Handle := a;
  BEGIN
    REPEAT
      SetEdgSingle(t,n);
      t := nextF(t);
    UNTIL (t = a);
  END SetEdgAll;

PROCEDURE Fce(a: Handle): EFNode =
  BEGIN
    WITH dbit = DualBit(a) DO RETURN a.fe.edg[1-dbit] END
  END Fce;

PROCEDURE SetFceSingle(a: Handle; n: EFNode) =
  BEGIN
    WITH dbit = DualBit(a) DO a.fe.edg[1-dbit] := n END
  END SetFceSingle;

PROCEDURE SetFceAll(a: Handle; n: EFNode) =
  VAR t: Handle := a;
  BEGIN
    REPEAT
      SetFceSingle(t,n);
      t := nextE(t);
    UNTIL (t = a);
  END SetFceAll;

(* ======= Vertex/cell data records ======= *)

REVEAL
  VPNode = PublicVPNode BRANDED OBJECT 
    END;

PROCEDURE Org(a: Handle): VPNode =
  BEGIN
    WITH srbits = SrotBits(a) DO RETURN a.fe.org[srbits] END
  END Org;

PROCEDURE SetOrgSingle(a: Handle; n: VPNode) =
  BEGIN
    WITH srbits = SrotBits(a) DO a.fe.org[srbits] := n END
  END SetOrgSingle;

PROCEDURE SetOrgAll(a: Handle; n: VPNode) =
  VAR c : Handle := a;
  BEGIN
    WITH nei = OutgoingEdges(a) DO
      FOR i := 0 TO LAST(nei^) DO
        WITH b = nei[i] DO
	  c := b;  
	  REPEAT
            SetOrgSingle(c,n);
	    c := nextF(c)    
          UNTIL (c = b);
	END
      END
    END
  END SetOrgAll;

PROCEDURE Dst(a: Handle): VPNode =
  BEGIN RETURN Org(Flip(a)) END Dst;

PROCEDURE SetDstSingle(a: Handle; n: VPNode) =
  BEGIN SetOrgSingle(Flip(a), n) END SetDstSingle;

PROCEDURE SetDstAll(a: Handle; n: VPNode) =
  BEGIN SetOrgAll(Flip(a), n) END SetDstAll;

PROCEDURE Ppos(a: Handle): VPNode =
  BEGIN RETURN Org(Tors(a)) END Ppos;

PROCEDURE SetPposSingle(a: Handle; n: VPNode) =
  BEGIN SetOrgSingle(Tors(a), n) END SetPposSingle;

PROCEDURE SetPposAll(a: Handle; n: VPNode) =
  BEGIN SetOrgAll(Tors(a), n) END SetPposAll;

PROCEDURE Pneg(a: Handle): VPNode =
  BEGIN RETURN Org(Srot(a)) END Pneg;

PROCEDURE SetPnegSingle(a: Handle; n: VPNode) =
  BEGIN SetOrgSingle(Srot(a), n) END SetPnegSingle;

PROCEDURE SetPnegAll(a : Handle; n: VPNode) =
  BEGIN SetOrgAll(Srot(a), n) END SetPnegAll;

(* ================ TRAVESAL ===================== *)

CONST
   InitHandleStackSize = 1024;
   InitEdgeStackSize = 1024;
   InitVertexStackSize = 1024;
   InitOutStackSize = 20;

PROCEDURE EnsureStackSpace(VAR stack: REF ARRAY OF Handle; VAR top: CARDINAL) =
  (* If the stack is full, allocates a new one, copying the old contents to it. *)
  BEGIN
    IF top >= NUMBER(stack^) THEN
      WITH 
    	sz = NUMBER(stack^), szNew = 2*sz,
    	stackNew = NEW(REF ARRAY OF Handle, szNew)
      DO
    	SUBARRAY(stackNew^, 0, sz) := stack^; stack := stackNew;
      END
    END
  END EnsureStackSpace;

PROCEDURE EnumHandles(READONLY a: ARRAY OF Handle; visit: VisitProc; fes: BOOLEAN := FALSE) =
  VAR
    stack := NEW(REF ARRAY OF Handle, InitHandleStackSize);
    np: CARDINAL;

  (* The fe record of handle "t" is marked by setting "t.fe.emark[DualBit(t)]". *)
  (* Note that only the primal or dual half of the record is marked, but each mark *)
  (* applies to both t and Flip(t). *)
  PROCEDURE VisitAndMark(t: Handle)=
    (* If fe(t) is unmarked: visit, mark, and stack it. *)
    BEGIN
      IF NOT t.fe.emark[DualBit(t)] THEN 
        visit(t);
        IF NOT fes THEN visit(Flip(t)) END;
        t.fe.emark[DualBit(t)] := TRUE; 
        EnsureStackSpace(stack, np);
        stack[np] := t; INC(np);
      END;
    END VisitAndMark;

  VAR seen: CARDINAL; (* # of quads whose childeren were looked at *)
  BEGIN
    IF NUMBER(a) > 0 THEN <* ASSERT NOT a[0].fe.emark[DualBit(a[0])] *> END;
    np := 0;
    seen := 0;
    TRY
      FOR i := 0 TO LAST(a) DO VisitAndMark(a[i]) END; 
      WHILE seen < np DO 
        WITH b = stack[seen] DO
          VisitAndMark(nextE(b));
          VisitAndMark(nextF(b));
        END;
        seen := seen + 1
      END;
    FINALLY
      (* Erase all marks *)
      FOR i := 0 TO np-1 DO
        WITH b = stack[i] DO
          b.fe.emark[DualBit(b)] := FALSE;
        END
      END
    END
  END EnumHandles;

(* ====== Numbering ====== *)

PROCEDURE NumberFaceEdges(READONLY a: ARRAY OF Handle): REF ARRAY OF Handle =

  VAR stack := NEW(REF ARRAY OF Handle, InitHandleStackSize);
      nfe: CARDINAL := 0;

  (* The fe record of handle "t" is marked by setting "t.fe.num" to the index "k" *)
  (* such that stack[k] = t. *)
  (* Note that the mark applies to all eight handles of that facet-edge. *)
  PROCEDURE VisitAndMark(t: Handle) =
    (* If t is unmarked: visit, mark, and stack it. *)
    BEGIN
      WITH num = t.fe.num DO
        IF NOT (num < nfe AND stack[num].fe = t.fe) THEN
          EnsureStackSpace(stack, nfe);
          num := nfe;
          stack[nfe] := t; INC(nfe);
        END
      END
  END VisitAndMark;

  VAR seen: CARDINAL := 0; (* # of fes whose childeren were looked at *)
  BEGIN
    nfe := 0;
    seen := 0;
    FOR i := 0 TO LAST(a) DO VisitAndMark(a[i]) END; 
    WHILE seen < nfe DO
      WITH s = stack[seen] DO
        VisitAndMark(nextF(s));
        VisitAndMark(nextE(s));
      END;
      seen := seen + 1
    END;
    WITH r = NEW(REF ARRAY OF Handle, nfe) DO
      r^ := SUBARRAY(stack^, 0, nfe);
      RETURN r
    END;
  END NumberFaceEdges;

PROCEDURE NumberEdgesOrFacets(READONLY a: ARRAY OF Handle; facets: BOOLEAN := FALSE): REF ARRAY OF Handle =
  VAR
    stack := NEW(REF ARRAY OF Handle, InitHandleStackSize);
    ne: CARDINAL; 

  (* The edge "e" is marked by setting t.fe.mark[DualBit(t)] for every handle t *)
  (* with that same edge.  *)
    
  PROCEDURE VisitAndMark(t: Handle) =
    (* If t is unmarked: visit, mark, and stack it. *)
    VAR tn: Handle;
    BEGIN
      IF NOT t.fe.emark[DualBit(t)] THEN
        tn := t;
        REPEAT
          tn.fe.emark[DualBit(tn)] := TRUE;
          WITH e = Edg(tn) DO IF e # NIL THEN e.num := ne END END;
          tn := nextF(tn);
        UNTIL (tn = t);
        EnsureStackSpace(stack, ne);
        stack[ne] := t; INC(ne);
      END; 
    END VisitAndMark;

  VAR seen: CARDINAL := 0; sn: Handle;
  BEGIN
    seen := 0;
    ne := 0;
    FOR i := 0 TO LAST(a) DO
      IF facets THEN VisitAndMark(Sdual(a[i])) ELSE VisitAndMark(a[i]) END
    END;
    WHILE seen < ne DO
      WITH s = stack[seen] DO
        sn := s;
        REPEAT
          VisitAndMark(nextE(sn));
          sn := nextF(sn);
        UNTIL (sn = s);
      END;
      seen := seen + 1
    END;

    (* Erase all marks *)
    FOR i := 0 TO ne-1 DO
      WITH b = stack[i] DO
        VAR bn: Handle := b;
        BEGIN
          REPEAT
            bn.fe.emark[DualBit(bn)] := FALSE;
            bn:= nextF(bn);
          UNTIL (bn = b);
        END;
      END
    END;

    WITH r = NEW(REF ARRAY OF Handle, ne) DO
      r^ := SUBARRAY(stack^, 0, ne);
      IF facets THEN FOR i := 0 TO ne-1 DO r[i] := Sdual(r[i]) END END;
      RETURN r
    END;
  END NumberEdgesOrFacets;

PROCEDURE NumberFaceEdges(READONLY a: ARRAY OF Handle): REF ARRAY OF Handle;
  (*
    Assigns distinct serial numbers to all fes reachable from "a" by
    nextE/nextF chains. Stores those numbers in the "num" field
    of each face-edge record. Returns a vector with one reachable handle 
    from each face-edge (not necessarily the reference one).  *)

PROCEDURE NumberEdges(READONLY a: ARRAY OF Handle): REF ARRAY OF Handle;
  (*
    Enumerates all edges that are reachable from "a" by
    nextE/nextF chains. Assigns distinct serial numbers to those edges and stores
    them in the ".num" field of the corresponding EFNode records, if they exist. 
    In any case, returns a representative handle "b" for each edge "e" such that 
    Edge(b) = e (resp. for each facet f such that Face(b) = f). *)

PROCEDURE NumberVerticesOrCells(READONLY a: ARRAY OF Handle; cells: BOOLEAN := FALSE): REF ARRAY OF Handle; 
  (*
    Enumerates all vertices (resp. cells, if "cells = TRUE") that are
    reachable from "a" by nextE/nextF chains. Assigns distinct serial numbers to
    those vertices (resp. cells) and stores them in the ".num" field
    of the corresponding VPNode records, if they exist.  In any case, 
    returns a representative handle "b" for of each vertex v such that Org(b) = v
    (resp. for each cell p such that NegC(b) = p). *)


(* ====== Low-level splicicing tools ====== *)

(* 
  The perators SpliceFaces" and "SpliceEdges" modify the "NextF" and
  "NextE" pointers, respectively.
  
  These procedures only guarantee preservation of the locally testable
  axioms, excluding `Bipolar'. The user is responsible for restoring
  the validity of the latter. *)

PROCEDURE SpliceFaces(a, b: Handle);
  (* Merges or splits the face-rings of handles "a" and "b", so that
    "NextF(a) = b1" and "NextF(b)=a1", where "a1" and "b1" are the 
    original sucessors of "a" and "b", resepctively. *)

PROCEDURE SpliceEdges(a, b: Handle);
  (* Merges or splits the edge-rings of handles, so that
    "NextF(a) = b1" and "NextF(b)=a1", where "a1" and "b1" are the 
    original sucessors of "a" and "b", resepctively. *)

PROCEDURE SetNextF(a,b: Handle);
  (*
    If "NextF(a) # b", performs "SpliceFacets(a, PrevF(b))". After this call
    "NextF(a)" will be equal to "b". Valid whenever "SpliceFacets(a,b)" is valid. *)

PROCEDURE SetnextE(a,b: Handle);
  (*
    If "nextE(a) # b", performs "SpliceEdges(a, prevE(b))".
    After this call, "nextE(a)" will be equal to "b".
    Valid whenever "SpliceEdges(a,b)" is valid. *)


PROCEDURE FaceSides(a: Handle): REF ARRAY OF Handle =
  VAR an : Handle := a; n: CARDINAL := 0;
  BEGIN
    an := a; REPEAT an := nextE(an); INC(n) UNTIL an = a;
    WITH rs = NEW(REF ARRAY OF Handle, n), s = rs^ DO
      FOR i := 0 TO n-1 DO
        s[i] := an; an := nextE(an);
      END;
      RETURN rs
    END
  END FaceSides;

PROCEDURE EdgeWings(a: Handle): REF ARRAY OF Handle =
  BEGIN
    WITH rs = FaceSides(Sdual(a)), s = rs^ DO
      FOR i := 0 TO LAST(s) DO s[i] := Sdual(s[i]) END;
      RETURN rs
    END
  END EdgeWings;

PROCEDURE OutgoingEdges(a: Handle): REF ARRAY OF Handle =

  VAR
    stack := NEW(REF ARRAY OF Handle, InitOutStackSize);
    nev: CARDINAL; (* Number of edges out of vertex. *)

  BEGIN
    DoListOutgoingEdges(a, stack, nev, -1);
    (* Erase all marks *)
    FOR i := 0 TO nev-1 DO
      WITH b = stack[i] DO
        b.fe.vmark[SrotBits(b)] := FALSE;
      END
    END;
    
    WITH r = NEW(REF ARRAY OF Handle, nev) DO
      r^ := SUBARRAY(stack^, 0, nev);
      RETURN r
    END;
  END OutgoingEdges;

PROCEDURE BoundingFaces(a: Handle): REF ARRAY OF Handle =
  BEGIN
    WITH r = OutgoingEdges(Sdual(a)) DO
      FOR i := 0 TO LAST(r^) DO r[i] := Sdual(r[i]) END;
      RETURN r
    END
  END BoundingFaces;

PROCEDURE DoListOutgoingEdges(a: Handle; VAR stack: REF ARRAY OF Handle; VAR top: CARDINAL; vnum: INTEGER) =
  (*
     Enumerates all handles with same origin as "a", sets the respective marks
     (indexed by SrotBits), and pushes one respresentative handle for each edge
     out of that vertex on the given "stack". A handle t is marked by setting
     t.fe.vmark[SrotBits(t)].  If vnum >= 0, sets the "num" field of all vertex records
     seen to "vnum".
  *)

  PROCEDURE VisitAndMark(t: Handle) =
    (* If t is unmarked: visit, mark, and stack it. *)
    VAR tn: Handle;
    BEGIN
      IF NOT t.fe.vmark[SrotBits(t)] THEN
        tn := t;
        REPEAT
          IF vnum >= 0 THEN
            WITH v = Org(tn) DO IF v # NIL THEN v.num := vnum END END
          END;
          tn.fe.vmark[SrotBits(tn)] := TRUE;
          tn := nextF(tn)
        UNTIL tn = t;
        EnsureStackSpace(stack, top);
        stack[top] := t; INC(top);
      END
    END VisitAndMark;

  VAR seen: CARDINAL; sn: Handle;
  BEGIN
    seen := top;
    VisitAndMark(a);
    WHILE seen < top DO
      WITH s = stack[seen] DO
        sn := s;
        REPEAT
          VisitAndMark(Flip(prevE(sn)));
          sn := nextF(sn)
        UNTIL sn = s;
      END;
      seen := seen + 1
    END;
  END DoListOutgoingEdges;

PROCEDURE NumberVerticesOrCells(READONLY a: ARRAY OF Handle; cells: BOOLEAN := FALSE): REF ARRAY OF Handle =
  VAR
    estack := NEW(REF ARRAY OF Handle, InitEdgeStackSize);   (* One representative handle out of each directed edge. *)
    ne: CARDINAL := 0;

    vstack := NEW(REF ARRAY OF Handle, InitVertexStackSize); (* One representative handle out of each vertex. *)
    nv: CARDINAL := 0;

  (* A vertex v is marked by setting t.fe.vmark[SrotBits(t)] for all handles t out of v. *)
  PROCEDURE VisitAndMark(t: Handle) =
    (* If t is unmarked: visit, mark, and stack it. *)
    BEGIN
      IF NOT t.fe.vmark[SrotBits(t)] THEN
        DoListOutgoingEdges(t, estack, ne, nv);
        EnsureStackSpace(vstack, nv);
        vstack[nv] := t; INC(nv);
      END
    END VisitAndMark;

  VAR seen: CARDINAL := 0;
  BEGIN
    seen := 0;
    ne := 0;
    nv := 0;
    
    FOR i := 0 TO LAST(a) DO
      IF cells THEN VisitAndMark(Sdual(a[i])) ELSE VisitAndMark(a[i]) END
    END;
    
    WHILE seen < ne DO
      WITH s = estack[seen] DO
        VisitAndMark(Flip(s))
      END;
      seen := seen + 1
    END;
    
    (* Erase all marks *)
    FOR i := 0 TO ne-1 DO
      WITH b = estack[i] DO
        VAR bn: Handle := b;
        BEGIN
          REPEAT
            bn.fe.vmark[SrotBits(bn)] := FALSE;
            bn:= nextF(bn);
          UNTIL (bn = b);
        END
      END
    END;
    
    WITH r = NEW(REF ARRAY OF Handle, nv) DO
      r^ := SUBARRAY(vstack^, 0, nv);
      IF cells THEN FOR i := 0 TO nv-1 DO r[i] := Sdual(r[i]) END END;
      RETURN r
    END;
  END NumberVerticesOrCells;

(* ====== Handle numbering ====== *)

PROCEDURE GetHandleNum(a: Handle): HandleNum =
  BEGIN RETURN Word.LeftShift(a.fe.num, 3) + a.bits END GetHandleNum;

(* ================= Traversal's functions of QuadEdge ========= *)

<* UNUSED *>
PROCEDURE PrintVMark(a: Handle) =
  (* Print the mark of fe.edge *)
  <* FATAL Wr.Failure, Thread.Alerted *>
  VAR an: Handle := a;
  BEGIN
   REPEAT
    Wr.PutText(stdout, FmtHandle(an, 1));
    Wr.PutText(stdout, " **" & Fmt.Bool(an.fe.vmark[SrotBits(an)]) & "\n");
    an := nextF(an);
   UNTIL(an = a) 
   END PrintVMark;

PROCEDURE FmtHandle(a: Handle;  feWidth: CARDINAL := 1): TEXT =
  BEGIN
    RETURN
      Fmt.Pad(Fmt.Int(a.fe.num), feWidth) & ":" &
      Fmt.Int(SrotBits(a)) & ":" & Fmt.Int(RevEBit(a));
  END FmtHandle;

PROCEDURE ReadHandle(rd: Rd.T; READONLY map: ARRAY OF FaceEdge): Handle =
  VAR m, r, s: CARDINAL;
  <* FATAL Rd.Failure, Rd.EndOfFile, Thread.Alerted, FloatMode.Trap,
     Lex.Error *>
  BEGIN
    Lex.Skip(rd);
    m := Lex.Int(rd);
    <* ASSERT Rd.GetChar(rd) = ':' *>
    r := Lex.Int(rd);
    <* ASSERT Rd.GetChar(rd) = ':' *>
    s := Lex.Int(rd);
    <* ASSERT m < NUMBER(map) *>
    <* ASSERT r < 4 AND s < 2 *>
    RETURN Handle{fe := map[m], bits := 2 * r + s}
  END ReadHandle;

PROCEDURE PrintFaceEdge(wr: Wr.T; n: FaceEdge; feWidth: CARDINAL := 1) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    FOR i := 0 TO 3 DO 
      IF i > 0 THEN Wr.PutChar(wr, ' ') END;
      Wr.PutText(wr, FmtHandle(Handle{fe := n.fenext[i], bits := n.bnext[i]}, feWidth));
    END
  END PrintFaceEdge;
  
PROCEDURE ReadFaceEdge(rd: Rd.T;  n: FaceEdge; READONLY map: ARRAY OF FaceEdge) =
  BEGIN
    FOR i := 0 TO 3 DO 
      SetnextF(Handle{fe := n, bits := 2 * i}, ReadHandle(rd, map))
    END
  END ReadFaceEdge;


BEGIN
END Map3D.
 
(**************************************************************************)
(*                                                                        *)
(* Copyright (C) 1998 Universidade Estadual de Campinas (UNICAMP)         *)
(*                                                                        *)
(* Authors:                                                               *)
(*   L. P. Lozada  & J. Stolfi  - UNICAMP <lozada@dcc.unicamp.br>         *)
(*                                                                        *)
(* This file can be freely used, distributed, and modified, provided that *)
(* this copyright and authorship notice is included in every copy or      *)
(* derived version.                                                       *)
(*                                                                        *)
(* DISCLAIMER: This software is offered ``as is'', without any guarantee  *)
(*   as to fitness for any particular purpose.  Neither the copyright     *)
(*   holder nor the authors or their employers can be held responsible    *)
(*   for any damages that may result from its use.                        *)
(*                                                                        *)
(* Last edited on 2001-05-06 23:59:36 by stolfi                           *)
(**************************************************************************)