MODULE Octf;

(*
   Facet-Edge data structure. Created on 1998 by Luis Arturo Perez Lozada,
   (see notice of copyright at the end of this file), inspired in the im-
   plementation of Quad-Edge data structure (by J. Stolfi and R. Marcone)

   See :
   "Primitives for the Manipulation of Three-Dimensional Subdivisions" by
   D. Dobkin and J. Laszlo, Algorithmica 1989.

   Last Modification: 03-08-00 by stolfi:

   
*)


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

FROM Stdio IMPORT stdout;

TYPE 
  MarkBit = BOOLEAN;

REVEAL 
  Node = PublicNode BRANDED OBJECT 
    END;

REVEAL     
  Edge = PublicEdge BRANDED
  OBJECT
    marks: ARRAY [0..1] OF  MarkBit;
    (* Used by NumberEdges, ex: marks[SpinBit(a)] associated to Edge *)

  OVERRIDES  
    init := EdgeInit;
  END;

REVEAL 
  Face = PublicFace BRANDED
  OBJECT
    marks: ARRAY [0..1] OF  MarkBit; 
    (* Used by NumberFacets, ex: marks[DualBit(a)] associated to Face *)

  OVERRIDES
    init := FaceInit;
  END;
 
REVEAL 
  FacetEdge = PublicFacetEdge BRANDED
  OBJECT 
    marks: ARRAY [0..1] OF  MarkBit;

    (* Used by EnumVertices, EnumFacetEdges *)
    (* ex: marks[OrientationBit(a)] associarted to Vertices *)

    (* Element index [r] refers to pair a[r] = Pair{n,r}: *)
    (* and r is RBits the Rotations Bits, r = {0,1,2,3}   *)

    fenext: ARRAY [0..3] OF FacetEdge;     (* FacetEdge of Fnext(a[r]) *)
    bnext:  ARRAY [0..3] OF SRBits;        (* SRBits of Fnext(a[r]) *)
    ind : CARDINAL := 0;                   (* index of FacetEdge used by
                                              procedures NumberFacets, Number-
					      Edges and NumberEdgesForDegree *)
  OVERRIDES  
    init := OctInit;
  END;

  (* Description for computational of bits for Spin and Srot
    If "a = spin^s(srot^r(a0))", where "a0" is the reference pair
    of "a", "r" in [0..3], and "s" in [0..1], then "a.bits = (2r + s)".
    Therefore,

    | a.spin.bits = XOR(a.bits, 1)
    | If sbit is 0
    |   a.srot.bits = (a.bits + 2) MOD 8 
    | otherwise
    |   a.srot.bits = (a.bits + 6) MOD 8

 The pair "Fnext(a)" is given by "pair(a.facetedge.fenext[rbit],
 a.facetedge.bnext[rbit])", provided that the "spin" bit of "a" is 0.
 Otherwise "Fnext(a)" is computed by 
 the formula: "Fnext(a) = Spin(Clock(Fnext(Clock(Spin(a)))))".

 rbit is obtained by "rbit = Word.And( Word.RightShift( a.bits , 1 ), 3 ) *)
    
PROCEDURE OctInit(n: FacetEdge): FacetEdge =    
  BEGIN 
    n.marks[0] := FALSE;
    n.marks[1] := FALSE;
            
    FOR i := 0 TO 3 DO n.fenext[i] := n END;
    n.bnext[0] := 0;
    n.bnext[1] := 2;
    n.bnext[2] := 4;
    n.bnext[3] := 6;
    n.edge := MakeEdge();
    n.face := MakeFace();
    RETURN n
  END OctInit;

PROCEDURE EdgeInit(e: Edge): Edge =    
  BEGIN
    e.marks[0] := FALSE;
    e.marks[1] := FALSE;
  RETURN e;
  END EdgeInit;

PROCEDURE FaceInit(f: Face): Face =    
  BEGIN 
    f.marks[0] := FALSE;
    f.marks[1] := FALSE;
  RETURN f;
  END FaceInit;


PROCEDURE MakeEdge(): Edge =
  BEGIN
    WITH e = NEW(Edge).init() DO RETURN e END;
  END MakeEdge;

PROCEDURE MakeFace(): Face =
  BEGIN
    WITH f = NEW(Face).init() DO RETURN f END;
  END MakeFace;

PROCEDURE MakeFacetEdge (): Pair =
 VAR a : Pair;
 BEGIN
  WITH n = NEW(FacetEdge).init() DO
    a := Pair{facetedge := n, bits := 0};
    a.facetedge.edge.pa := a;
    a.facetedge.face.pa := a;
    RETURN a;
  END;
 END MakeFacetEdge;

(* ====== Updating relations ====== *)

PROCEDURE SetFace(a: Pair; n: Face) =
  BEGIN
    WITH c = a.facetedge.face DO
      c := n;
    END;
  END SetFace;

PROCEDURE SetEdge(a: Pair; n: Edge) =
  BEGIN
    WITH c = a.facetedge.edge DO
      c := n;
    END;
  END SetEdge;

PROCEDURE SetEdgeAll(a: Pair; n: Edge) =
  VAR t: Pair := a;
  BEGIN
    REPEAT
      SetEdge(t,n);
      t := Fnext(t);
    UNTIL (t = a);
  END SetEdgeAll;


PROCEDURE SetFaceAll(a: Pair; n: Face) =
  VAR t: Pair := a;
  BEGIN
    REPEAT
      SetFace(t,n);
      t := Enext(t);
    UNTIL (t = a);
  END SetFaceAll;

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

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

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

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

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

PROCEDURE SpinBit(a: Pair): SBit =
  BEGIN RETURN Word.And(a.bits, 1) END SpinBit;
  
PROCEDURE DualBit(a: Pair): DBit =
  BEGIN RETURN Word.And(Word.RightShift(a.bits, 1), 1) END DualBit;

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

PROCEDURE GetPairNum(a: Pair): PairNum =
  BEGIN RETURN Word.LeftShift(a.facetedge.num, 3) + a.bits END GetPairNum;

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

PROCEDURE Spin(a: Pair): Pair =
  BEGIN
    RETURN Pair{facetedge := a.facetedge, bits := Word.Xor( a.bits, 1 )};
  END Spin;
  
PROCEDURE Clock(a: Pair): Pair =
  BEGIN
    RETURN Pair{facetedge := a.facetedge, bits := Word.And( (a.bits + 4), 7 )};
  END Clock;
  
PROCEDURE Tors(a: Pair): Pair =
  BEGIN
    RETURN Pair{facetedge := a.facetedge,
               bits := Word.And(a.bits + 6 + Word.LeftShift( 
                       Word.And(a.bits, 1), 2), 7) };    
  END Tors;


PROCEDURE Fnext(a: Pair) : Pair =

  PROCEDURE Fnext_(a: Pair): Pair =
    VAR r : SRBits;
    BEGIN
      WITH s = Word.And( Word.RightShift(  a.bits, 1 ), 3 ) DO 
        r := a.facetedge.bnext[s]; 
        a.facetedge := a.facetedge.fenext[s];
        a.bits := r;  
        RETURN a;
      END;
    END Fnext_;

  BEGIN
    IF ( Word.And( a.bits, 1 ) = 0 ) THEN RETURN Fnext_(a);
    ELSE RETURN Spin(Clock(Fnext_(Clock(a))));
    END;
  END Fnext;

PROCEDURE Sdual(a: Pair): Pair =
  BEGIN RETURN Spin(Srot(a)) END Sdual;

PROCEDURE Enext(a: Pair): Pair =
  BEGIN RETURN Srot(Fnext(Tors(a)))  END Enext;

PROCEDURE Enext_1(a: Pair): Pair =
  BEGIN RETURN Clock(Srot(Fnext(Srot(a))))  END Enext_1;

PROCEDURE Fnext_1(a: Pair): Pair =
  BEGIN RETURN Clock(Fnext(Clock(a)))  END Fnext_1;

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

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

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

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

PROCEDURE SpliceEdges(a,b: Pair) =
  BEGIN
    SpliceFacets(Sdual(a), Sdual(b));
  END SpliceEdges;
 
PROCEDURE DeleteFacetEdge (a : Pair) =
  BEGIN
    SpliceFacets(a, Fnext_1(a));
    SpliceFacets(Clock(a), Fnext_1(Clock(a))); 
  END DeleteFacetEdge;

PROCEDURE Meld(a,b : Pair) =
  (* Delete the first pair facetedge argument *)
  VAR firsta,a_1 : Pair;
  BEGIN
    firsta := a;
    REPEAT
      a_1 := Fnext_1(a);
      IF Fnext(a) # b THEN
        SpliceFacets(a, Fnext_1(b));
        <* ASSERT Fnext(a) = b *>
      END; 
      DeleteFacetEdge(a);
      IF SpinBit(a_1) = SpinBit(b) THEN
        <* ASSERT Fnext(a_1) = b *>
      END;
      a := Enext(a);
      b := Enext(b);
    UNTIL a = firsta;
  END Meld;

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

PROCEDURE Onext(s: Pair): Pair =
  BEGIN
    IF DualBit(s)= 0 THEN 
      RETURN Clock(Fnext(Enext_1(s)));
    ELSE
      RETURN Enext_1(s); 
    END;  
  END Onext;

PROCEDURE Onext_1(s: Pair): Pair =
  (* Procedure equivalent to function Onext but change Fnext by Fnext_1 *) 
  BEGIN
    RETURN Clock(Fnext_1(Enext_1(s)));
  END Onext_1;

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

PROCEDURE EnumFacetEdges(a: Pair; visit: VisitProc;
                         facetedges: BOOLEAN := FALSE) =
  CONST IniStackSize = 1024;

  VAR festack := NEW(REF ARRAY OF FacetEdge, IniStackSize);
      bstack  := NEW(REF ARRAY OF SRBits, IniStackSize);
      top: CARDINAL;

  PROCEDURE DoubleStack() =
    BEGIN
      WITH 
        sz = NUMBER(festack^), szNew = 2*sz,
        festackNew = NEW(REF ARRAY OF FacetEdge, szNew),
        bstackNew = NEW(REF ARRAY OF SRBits, szNew)
      DO
        SUBARRAY(festackNew^, 0, sz) := festack^; festack := festackNew;
        SUBARRAY(bstackNew^, 0, sz) := bstack^; bstack := bstackNew;
      END
    END DoubleStack;

  PROCEDURE VisitandMark(c: Pair)=
    (* If facetedge(c) is unmarked: visit, mark, and stack it. *)
    BEGIN
      IF NOT c.facetedge.marks[DualBit(c)] THEN 
        visit(c);
        IF NOT facetedges THEN visit(Clock(c)) END;
        c.facetedge.marks[DualBit(c)] := TRUE; (* c.facetedge marked *) 
        IF top >= NUMBER(festack^) THEN DoubleStack() END;
        festack[top] := c.facetedge; bstack[top] := c.bits;
        top := top + 1;
      END;
    END VisitandMark;

  VAR seen: CARDINAL; (* # of quads whose childeren were looked at *)
  BEGIN
    <* ASSERT NOT a.facetedge.marks[DualBit(a)] *> 
    top := 0;
    seen := 0;
    TRY
      VisitandMark(a);
      WHILE seen < top DO 
        WITH b = Pair{facetedge := festack[seen], bits := bstack[seen]} DO
          VisitandMark(Fnext_1(b));
          VisitandMark(Enext_1(b));
          VisitandMark(Enext(b));
          VisitandMark(Fnext(Enext_1(b)));	  
        END;
        seen := seen + 1
      END;
    FINALLY
      (* Erase all marks *)
      WHILE top > 0 DO
        top := top - 1;
        WITH b = Pair{facetedge := festack[top], bits := bstack[top]} DO
          b.facetedge.marks[DualBit(b)] := FALSE;
        END
      END
    END
  END EnumFacetEdges;

PROCEDURE NumberEdges(READONLY a: ARRAY OF Pair): REF ARRAY OF Edge =
  (* Note: I change the SpinBit by DualBit for purpouses of gluing octahedra. *)

  CONST
    IniStackSize = 1024;
  VAR
    estack := NEW(REF ARRAY OF Edge, IniStackSize);
    stack := NEW(REF ARRAY OF Pair, IniStackSize);
    top: CARDINAL;         (* top for stack "estack" *)
    nstack: CARDINAL;      (* top for stack "nstack" *)

  PROCEDURE DoubleStack() =
    BEGIN
      WITH 
        sz = NUMBER(stack^), szNew = 2*sz,
        estackNew = NEW(REF ARRAY OF Edge, szNew),
        stackNew = NEW(REF ARRAY OF Pair, szNew)
      DO
        SUBARRAY(estackNew^, 0, sz) := estack^; estack := estackNew;
        SUBARRAY(stackNew^, 0, sz) := stack^; stack := stackNew;
      END
    END DoubleStack;

  PROCEDURE VisitAndMark(t: Pair) =
    (* If t is unmarked: visit, mark, and stack it. *)
    VAR tn : Pair;
    BEGIN
      WITH ind = t.facetedge.ind, tnum1 = t.facetedge.edge.num DO
        IF ind < nstack AND stack[ind].facetedge = t.facetedge THEN
                  (* FacetEdge is already marked, not do nothing *)
        ELSE      (* FacetEdge not marked, do *)
          (* If component edge of pair "t" not are marked, then  *)
          IF NOT t.facetedge.edge.marks[DualBit(t)] THEN
             (* mark t.facetedge.edge and alls components edges
                of pairs adjacents to "t" *)
            tn := t;
            REPEAT
              tn.facetedge.edge.marks[DualBit(tn)] := TRUE;
              tn := Fnext(tn);
            UNTIL (tn = t);
         
            IF top >= NUMBER(estack^) THEN DoubleStack() END;
            estack[top] := t.facetedge.edge;
            estack[top].pa := t; 
            tnum1 := top;
            top := top + 1;
          END;
 
         IF nstack = NUMBER(stack^) THEN DoubleStack() END;
         ind := nstack;
         stack[nstack] := t;
         INC(nstack);
        END;
      END; 
    END VisitAndMark;

  VAR seen: CARDINAL := 0; (* # of facetedges whose childeren were looked at *)
  BEGIN
    seen := 0;
    top := 0;
    nstack := 0;
    FOR i := 0 TO LAST(a) DO
      VisitAndMark (a[i])
    END;
    WHILE seen < nstack DO
      WITH s = stack[seen] DO
        VisitAndMark(Fnext(s));
        VisitAndMark(Fnext(Enext_1(s)));
        VisitAndMark(Enext(s));
        VisitAndMark(Fnext_1(s)); 
      END;
      seen := seen + 1
    END;

    (* Erase all marks *)
    WHILE nstack > 0 DO
      nstack := nstack - 1;
      WITH b = stack[nstack] DO
        VAR atn: Pair := b;
        BEGIN
          REPEAT
            atn.facetedge.edge.marks[DualBit(atn)] := FALSE;
            atn:= Fnext(atn);
          UNTIL (atn = b);
        END;
      END
    END;

    WITH
      r = NEW(REF ARRAY OF Edge, top)
    DO
      r^ := SUBARRAY(estack^, 0, top);
      RETURN r
    END;
  END NumberEdges;

PROCEDURE NumberEdgesForDegree(READONLY a: ARRAY OF Pair): REF ARRAY OF Pair =
  (* Computes the degree of a vertex, i.e. the number of edges incidents to this
     vertex.
  *)
  CONST
    IniStackSize = 20;
  VAR
    estack := NEW(REF ARRAY OF Pair, IniStackSize);
    stack  := NEW(REF ARRAY OF Pair, IniStackSize);
    top: CARDINAL;         (* top for stack "estack" *)
    nstack: CARDINAL;      (* top for stack "nstack" *)

  PROCEDURE DoubleStack() =
    BEGIN
      WITH 
        sz = NUMBER(stack^), szNew = 2*sz,
        estackNew = NEW(REF ARRAY OF Pair, szNew),
        stackNew = NEW(REF ARRAY OF Pair, szNew)
      DO
        SUBARRAY(estackNew^, 0, sz) := estack^; estack := estackNew;
        SUBARRAY(stackNew^, 0, sz) := stack^; stack := stackNew;
      END
    END DoubleStack;

  PROCEDURE VisitAndMark(t: Pair) =
    (* If t is unmarked: visit, mark, and stack it. *)
    VAR tn : Pair;
    BEGIN
      WITH tind = t.facetedge.ind, tdg = t.facetedge.edge.dg DO
        IF tind < nstack AND stack[tind].facetedge = t.facetedge THEN
                  (* FacetEdge is already marked, not do nothing *)
        ELSE      (* FacetEdge not marked, do *)
          (* If component edge of pair "t" not are marked, then  *)
          IF NOT t.facetedge.marks[DualBit(t)] THEN
             (* mark t.facetedge.edge and alls components edges
                of pairs adjacents to "t" *)
            tn := t;
            REPEAT
              tn.facetedge.marks[DualBit(tn)] := TRUE;
              tn := Fnext(tn);
            UNTIL (tn = t);
          
            IF top >= NUMBER(estack^) THEN DoubleStack() END;
            estack[top] := t;
            tdg := top;
            top := top + 1;
          END;  
        
         IF nstack = NUMBER(stack^) THEN DoubleStack() END;
         tind := nstack;
         stack[nstack] := t;
         INC(nstack);
        END;
      END; 
    END VisitAndMark;

  VAR seen: CARDINAL := 0;
  BEGIN
    seen := 0;
    top := 0;
    nstack := 0;
    
    FOR i := 0 TO LAST(a) DO
      VisitAndMark(a[i]);
    END;
    
    WHILE seen < nstack DO
      WITH s = stack[seen] DO
        IF DualBit(s) = 0 THEN
          VisitAndMark(Onext(s));
          VisitAndMark(Onext(Fnext(s)));
        ELSE
          VisitAndMark(Clock(Onext(s)));
          VisitAndMark(Clock(Onext(Fnext(s))));
        END;
      END;
      seen := seen + 1
    END;
    
    (* Erase all marks *)
    WHILE nstack > 0 DO
      nstack := nstack - 1;
      WITH b = stack[nstack] DO
        VAR atn: Pair := b;
        BEGIN
          REPEAT
            atn.facetedge.marks[DualBit(atn)] := FALSE;
            atn:= Fnext(atn);
          UNTIL (atn = b);
        END
      END
    END;
    
    WITH
      r = NEW(REF ARRAY OF Pair, top)
    DO
      r^ := SUBARRAY(estack^, 0, top);
      RETURN r
    END;
  END NumberEdgesForDegree;

PROCEDURE NumberFacets(READONLY a: ARRAY OF Pair): REF ARRAY OF Face =
  CONST IniStackSize = 1024;

  VAR fstack := NEW(REF ARRAY OF Face, IniStackSize);
      stack := NEW(REF ARRAY OF Pair, IniStackSize);
      top: CARDINAL;        (* top for stack "festack" *)
      nstack: CARDINAL;     (* top for stack "stack"   *)

  PROCEDURE DoubleStack() =
    BEGIN
      WITH 
        sz = NUMBER(stack^), szNew = 2*sz,
        fstackNew = NEW(REF ARRAY OF Face, szNew),
        stackNew = NEW(REF ARRAY OF Pair, szNew)
      DO
        SUBARRAY(fstackNew^, 0, sz) := fstack^; fstack := fstackNew;
        SUBARRAY(stackNew^, 0, sz) := stack^; stack := stackNew;
      END
    END DoubleStack;

  PROCEDURE VisitAndMark (t: Pair) =
    (* If t is unmarked: visit, mark, and stack it. *)
    VAR tn : Pair;
        
    BEGIN
      WITH ind = t.facetedge.ind, tnum1 = t.facetedge.face.num DO
        IF ind < nstack AND stack[ind].facetedge = t.facetedge THEN
                  (* FacetEdge is already marked, not do nothing *)
        ELSE      (* FacetEdge no marked, do *)
          (* if component face of pair "t" not are marked, then *)
          IF NOT t.facetedge.face.marks[DualBit(t)] THEN
             (* mark t.facetedge.face and alls components facets 
                of pairs adjacents to "t" *)
            tn := t;
            REPEAT
              tn.facetedge.face.marks[DualBit(tn)] := TRUE;
              tn := Enext_1(tn);
            UNTIL (tn = t);
         
            IF top >= NUMBER(fstack^) THEN DoubleStack() END;
            fstack[top] := t.facetedge.face;
            fstack[top].pa := t; 
            tnum1 := top; 
            top := top + 1;  
          END;
 
         IF nstack = NUMBER(stack^) THEN DoubleStack() END;
         ind := nstack;
         stack[nstack] := t;
         INC(nstack);
        END;
      END; 

    END VisitAndMark;

  VAR seen: CARDINAL := 0; (* # of facetedges whose childeren were looked at *)
  BEGIN
    seen := 0;
    top := 0;
    nstack := 0;
    FOR i := 0 TO LAST(a) DO
      VisitAndMark (a[i])
    END;
    WHILE seen < nstack DO
      WITH s = stack[seen] DO
        VisitAndMark(Fnext(s));
        VisitAndMark(Fnext(Enext_1(s)));
        VisitAndMark(Enext(s));
        VisitAndMark(Fnext_1(s));
      END;
      seen := seen + 1
    END;

    (* Erase all marks *)
    WHILE nstack > 0 DO
      nstack := nstack - 1;
      WITH b = stack[nstack] DO
        VAR atn: Pair := b;
        BEGIN
          REPEAT
            atn.facetedge.face.marks[DualBit(atn)] := FALSE;
            atn:= Enext_1(atn);
          UNTIL (atn = b);
        END;
      END
    END;
    WITH
      r = NEW(REF ARRAY OF Face, top)
    DO
      r^ := SUBARRAY(fstack^, 0, top);
      RETURN r
    END;
  END NumberFacets;

<* UNUSED *>
PROCEDURE PrintMark(a: Pair) =
  (* Print the mark of facetedge.edge *)
  <* FATAL Wr.Failure, Thread.Alerted *>
  VAR an: Pair := a;
  BEGIN
   REPEAT
    PrintPair(stdout, an);
    Wr.PutText(stdout, " **" &
         Fmt.Bool(an.facetedge.marks[SpinBit(an)]) & "\n");
    an := Fnext(an);
   UNTIL(an = a) 
   END PrintMark;

PROCEDURE NumberFacetEdges(READONLY a: ARRAY OF Pair): REF ARRAY OF Pair =
  CONST IniStackSize = 1024;
  VAR stack := NEW(REF ARRAY OF Pair, IniStackSize);
      nstack: CARDINAL := 0;

  PROCEDURE DoubleStack() =
    BEGIN
      WITH 
        sz = NUMBER(stack^), szNew = 2*sz,
        stackNew = NEW(REF ARRAY OF Pair, szNew)
      DO
        SUBARRAY(stackNew^, 0, sz) := stack^; stack := stackNew;
      END
    END DoubleStack;
      
  (* An facetedge "n" is "marked" if stack[n.num].facetedge = n.
     i.e whether stacked *)
 
  PROCEDURE VisitAndMark(t: Pair) =
    (* If t is unmarked: visit, mark, and stack it. *)
    BEGIN
      WITH num = t.facetedge.num DO
        IF num < nstack AND stack[num].facetedge = t.facetedge THEN
          (* FacetEdge is already marked, not do nothing *)
        ELSE
          IF nstack = NUMBER(stack^) THEN DoubleStack() END;
          num := nstack;
          stack[nstack] := t;
          INC(nstack);
        END
      END
  END VisitAndMark;
  

  VAR seen: CARDINAL := 0; (* # of facetedges whose childeren were looked at *)
  BEGIN
    nstack := 0;
    seen := 0;
    FOR i := 0 TO LAST(a) DO
      VisitAndMark (a[i]) 
    END; 
    WHILE seen < nstack DO
      WITH s = stack[seen] DO
        VisitAndMark(Fnext_1(s));
        VisitAndMark(Enext_1(s));
        VisitAndMark(Enext(s));
        VisitAndMark(Fnext(Enext_1(s)));
      END;
      seen := seen + 1
    END;
    WITH
      r = NEW(REF ARRAY OF Pair, nstack)
    DO
      r^ := SUBARRAY(stack^, 0, nstack);
      RETURN r
    END;
  END NumberFacetEdges;

PROCEDURE SetFnext(a,b: Pair) =
  BEGIN
    IF Fnext(a) # b THEN SpliceFacets(a, Fnext_1(b)) END
  END SetFnext;

PROCEDURE SetEnext(a, b: Pair) =
  BEGIN
    IF Enext(a) # b THEN SpliceEdges(a, Enext_1(b)) END
  END SetEnext;

PROCEDURE PrintPair(
  wr: Wr.T;
  a: Pair;
  feWidth: CARDINAL := 1;
  nl : BOOLEAN := FALSE;
  ) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    IF nl THEN  
      Wr.PutText(wr,Fmt.Pad(Fmt.Int(a.facetedge.num), feWidth) & ":" &
      Fmt.Int(SrotBits(a)) & ":" & Fmt.Int(SpinBit(a))  & "\n");
    ELSE
      Wr.PutText(wr,Fmt.Pad(Fmt.Int(a.facetedge.num), feWidth) & ":" &
      Fmt.Int(SrotBits(a)) & ":" & Fmt.Int(SpinBit(a)));
    END
  END PrintPair;

PROCEDURE ReadPair(rd: Rd.T; READONLY map: ARRAY OF FacetEdge): Pair =
  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 Pair{facetedge := map[m], bits := 2 * r + s}
  END ReadPair;

PROCEDURE PrintFacetEdge(wr: Wr.T; n: FacetEdge; feWidth: CARDINAL := 1) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    FOR i := 0 TO 3 DO 
      IF i > 0 THEN Wr.PutChar(wr, ' ') END;
      PrintPair(wr, Pair{facetedge := n.fenext[i], bits := n.bnext[i]},
                feWidth);
    END
  END PrintFacetEdge;
  
PROCEDURE ReadFacetEdge(rd: Rd.T;
                        n: FacetEdge; READONLY map: ARRAY OF FacetEdge) =
  BEGIN
    FOR i := 0 TO 3 DO 
      SetFnext(Pair{facetedge := n, bits := 2 * i}, ReadPair(rd, map))
    END
  END ReadFacetEdge;


BEGIN
END Octf.
 
(**************************************************************************)
(*                                                                        *)
(* 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.                        *)
(*                                                                        *)
(**************************************************************************)