(* Triang data structure. *)
(* Last edited by Rober Marcone Rosi on Mon Oct 18 1993 *)

MODULE Triang;

IMPORT 
  Oct, R3, Color, FileStream, Random, Math, Thread, 
  Convert, Rd, Wr, Fmt, Lex;
FROM Oct IMPORT 
  RBits, RotBits, Sym, Onext, Lnext, Rprev, Rot, Tor, Flip, Oprev, 
  Enum, DumpDoc, LoadDoc, SetEdgeNo, GetEdgeNo, EnumVertices, Splice;
  
REVEAL
  Edge = PublicEdge BRANDED OBJECT
    org: ARRAY RBits OF Node;
    exists: BOOLEAN;
    spring: BOOLEAN;
    color: Color.T;
    radius: REAL;
  OVERRIDES
    init := EdgeInit;
  END;

REVEAL 
  Node = PublicNode BRANDED OBJECT 
      num: CARDINAL;
      exists: BOOLEAN;
      color: Color.T;
    END;

REVEAL     
  Vertex = PublicVertex BRANDED OBJECT 
      fixed: BOOLEAN;
      coords: R3.T;
      radius: REAL;
    END;

REVEAL 
  Face = PublicFace BRANDED OBJECT
      transp: Color.T;
      patch: CARDINAL;
    END;
    
PROCEDURE EdgeInit(e: Edge; no: Oct.EdgeNo := 0): Edge=
  BEGIN
    EVAL NARROW(e,Oct.Edge).init(no); 
    e.org[0] := NIL; (* For now *)
    e.org[1] := NIL; (* For now *)
    e.org[2] := NIL; (* For now *)
    e.org[3] := NIL; (* For now *)
    e.exists := TRUE; (* For now *)
    e.spring := TRUE; (* For now *)  
    e.color :=  Color.T{1.0, 1.0, 1.0}; (* For now *)
    e.radius := 0.0; (* For now *)
    RETURN  e;
  END EdgeInit;

PROCEDURE Make(no: Oct.EdgeNo := 0): Arc =
  BEGIN
    WITH e = NEW(Edge).init(no) DO
      RETURN Arc{edge := e, bits := 0};
    END;
  END Make;

PROCEDURE Org(e: Arc): Node =
  BEGIN
    WITH c = NARROW(e.edge, Edge).org[RotBits(e)] DO
      RETURN c;
    END;
  END Org;

PROCEDURE NOrg(e: Arc): CARDINAL =
  BEGIN
    WITH c = Org(e).num DO
      RETURN c;
    END;
  END NOrg;

PROCEDURE Degree(a: Arc): CARDINAL =
  VAR an: Arc := a;
      deg: CARDINAL := 0;
  BEGIN
    REPEAT
      an := Onext(an);
      INC(deg);
    UNTIL (an = a);
    RETURN deg
  END Degree;

PROCEDURE SetOrg(e: Arc; n: Node) =
  BEGIN
    WITH c = NARROW(e.edge, Edge).org[RotBits(e)] DO
      c := n
    END;
  END SetOrg;

PROCEDURE SetAllOrgs(a: Arc; n: Node) =
  VAR t: Arc := a;
  BEGIN
    REPEAT
      SetOrg(t, n);
      t := Onext(t)
    UNTIL t = a
  END SetAllOrgs;

PROCEDURE Left(e: Arc): Node =
  BEGIN
    RETURN Org(Tor(e))
  END Left;
   
PROCEDURE NLeft(e: Arc): CARDINAL =
  BEGIN
    WITH c = Left(e).num DO
      RETURN c;
    END;
  END NLeft;

PROCEDURE SetLeft(e: Arc; n: Node) =
  BEGIN
    SetOrg(Tor(e), n)
  END SetLeft;

PROCEDURE SetAllLefts(a: Arc; n: Node) =
  VAR t: Arc := a;
  BEGIN
    REPEAT
      SetLeft(t, n);
      t := Lnext(t)
    UNTIL t = a
  END SetAllLefts;
  
PROCEDURE NPatch(f: Face): CARDINAL =
  BEGIN
    RETURN f.patch
  END NPatch;
  
PROCEDURE SetNPatch(f: Face; patch: CARDINAL) =
  BEGIN
    f.patch := patch
  END SetNPatch;

PROCEDURE SetNodeColor(f: Node; color: Color.T) =
  BEGIN
    f.color := color
  END SetNodeColor;
 
PROCEDURE SetNodeExists(n: Node; b: BOOLEAN) =
  BEGIN
    WITH c = n.exists DO
      c := b
    END;
  END SetNodeExists;

PROCEDURE NodeExists(n: Node): BOOLEAN =
  BEGIN
    WITH c = n.exists DO
      RETURN c;
    END;
  END NodeExists; 

PROCEDURE SetVertexRadius(v: Vertex; radius: REAL) =
  BEGIN
    v.radius := radius
  END SetVertexRadius;
  
PROCEDURE SetFaceTransp(f: Face; color: Color.T) =
  BEGIN
    f.transp := color
  END SetFaceTransp;

PROCEDURE SetEdgeExists(e: Arc; b: BOOLEAN) =
  BEGIN
    WITH c = NARROW(e.edge, Edge).exists DO
      c := b
    END;
  END SetEdgeExists;

PROCEDURE SetEdgeSpring(e: Arc; b: BOOLEAN) =
  BEGIN
    WITH c = NARROW(e.edge, Edge).spring DO
      c := b
    END;
  END SetEdgeSpring;

PROCEDURE SetEdgeColor(e: Arc; color: Color.T) =
  (* Set edge color *)    
  BEGIN
    WITH c = NARROW(e.edge, Edge).color DO
      c := color
    END;
  END SetEdgeColor;
  
PROCEDURE SetEdgeRadius(e: Arc; radius: REAL) =
  BEGIN
    WITH ee = NARROW(e.edge, Edge) DO
      ee.radius := radius
    END;
  END SetEdgeRadius;
 
PROCEDURE SetVertexFixed(v: Vertex; b: BOOLEAN) =
  BEGIN
    v.fixed := b
  END SetVertexFixed;

PROCEDURE SetVertexCoords(v: Vertex; coords: R3.T) =
  BEGIN
    v.coords := coords
  END SetVertexCoords; 
  
PROCEDURE SetCoords(READONLY t: Topology; coords: Coords) =
  BEGIN
    FOR i := 0 TO LAST(coords) DO
      WITH
        v = t.vertex[i]
      DO
        SetVertexCoords(v, coords[i])
      END
    END;
  END SetCoords;

PROCEDURE GetVertexCoords(v: Vertex): R3.T =
  BEGIN
    WITH c = v.coords DO
      RETURN c
    END;
  END GetVertexCoords;
  
PROCEDURE GetCoords(READONLY t: Topology; VAR coords: Coords) =
  BEGIN
    FOR i := 0 TO LAST(coords) DO
      WITH
        v = t.vertex[i]
      DO
        coords[i] := GetVertexCoords(v)
      END;
    END;
  END GetCoords;

PROCEDURE EdgeExists(e: Edge): BOOLEAN =
  BEGIN
    WITH c = e.exists DO
      RETURN c;
    END;
  END EdgeExists; 

PROCEDURE EdgeSpring(e: Edge): BOOLEAN =
  BEGIN
    WITH c = e.spring DO
      RETURN c;
    END;
  END EdgeSpring; 

PROCEDURE VertexFixed(v: Vertex): BOOLEAN =
  BEGIN
    WITH c = v.fixed DO
      RETURN c;
    END;
  END VertexFixed; 

PROCEDURE InitVariable(READONLY t: Topology; VAR variable: ARRAY OF BOOLEAN) =
  BEGIN
    FOR i := 0 TO t.NV-1 DO
      WITH v = t.vertex[i] DO 
        (* ASSERT i = NARROW(v, Vertex).num *)
        variable[i] := (NOT VertexFixed(v)) AND NodeExists(v)
      END;
    END;
  END InitVariable;

PROCEDURE InitNode(e: Node; num: CARDINAL; exists: BOOLEAN; color: Color.T) =
  BEGIN
    e.num := num;
    e.exists := exists;
    e.color := color
  END InitNode;

PROCEDURE InitVertex(
    e: Vertex; 
    num: CARDINAL; 
    exists, fixed: BOOLEAN; 
    color: Color.T; 
    coords: R3.T; 
    radius: REAL
  ) =
  BEGIN
    InitNode(e, num, exists, color);
    e.fixed := fixed;
    e.coords := coords;
    e.radius := radius 
  END InitVertex;

PROCEDURE MakeVertex(): Node =
  BEGIN
    WITH 
      e = NEW(Vertex),
      color = Color.T{1.0, 1.0, 1.0},
      coords = R3.T{Random.Real(), Random.Real(), Random.Real()},
      radius = 0.0
    DO
      InitVertex(e, 0, TRUE, FALSE, color, coords, radius);
      RETURN e
    END;
  END MakeVertex;

PROCEDURE InitFace(
    e: Face; 
    num: CARDINAL; 
    exists: BOOLEAN; 
    color: Color.T; 
    transp: Color.T;
    patch: CARDINAL
  ) =
  BEGIN
    InitNode(e, num, exists, color);
    e.transp := transp;
    e.patch := patch
  END InitFace;

PROCEDURE MakeFace(): Face =
  BEGIN
    WITH
      e = NEW(Face),
      color = Color.T{1.0, 1.0, 1.0},
      transp = Color.T{0.8, 0.8, 0.8},
      patch = 0
    DO
      InitFace(e, 0, TRUE, color, transp, patch);
      RETURN e
    END;  
  END MakeFace;
 
PROCEDURE Connect(a, b: Arc): Arc =
  VAR e: Arc;
  BEGIN
    e := Make();
    Splice(e, a);
    Splice(Sym(e), b);
    RETURN e;
  END Connect;
  
PROCEDURE Glue(a, b: Arc; n: CARDINAL): Arc =
  VAR ta, tb, pa: Arc;
  BEGIN
    (* Sanity check: *)
    <* ASSERT n >= 1 *>
    ta := a; tb := b;
    FOR i := 2 TO n DO 
      ta := Lnext(ta); tb := Rprev(tb);
      <* ASSERT ta # a *>
      <* ASSERT tb # b *>
    END;
    
    Splice(a, Oprev(b));
    SetAllOrgs(a, Org(a));
    
    FOR i := 1 TO n DO
      ta := Lnext(a); tb := Rprev(b);
      Splice(ta, Oprev(tb));
      SetAllOrgs(ta, Org(ta));
      SetLeft(a, Left(b));
      (* Disconnect b: *)
      Splice(b, Oprev(b));
      Splice(Sym(b), Oprev(Sym(b)));
      pa := a; 
      a := ta; b := tb
    END;
    
    RETURN pa
  END Glue;

<*UNUSED*>
PROCEDURE Cos(x: REAL): REAL =
  BEGIN
    RETURN FLOAT(Math.cos(FLOAT(x, LONGREAL)))
  END Cos;

<*UNUSED*>
PROCEDURE Sin(x: REAL): REAL =
  BEGIN
    RETURN FLOAT(Math.sin(FLOAT(x, LONGREAL)))
  END Sin;

PROCEDURE R3Print(wr: Wr.T; w: R3.T) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    Wr.PutText(wr,
      "( " & Fmt.Real(w[0]) & " "
           & Fmt.Real(w[1]) & " "
           & Fmt.Real(w[2]) & " )"
    );
  END R3Print;

PROCEDURE FaceNormal(a: Arc; READONLY coords: Coords): R3.T =
  BEGIN
    WITH
      b = Onext(a),
      vo = coords[NOrg(a)],
      va = coords[NOrg(Sym(a))],
      vb = coords[NOrg(Sym(b))],
      normal =  R3.cross(R3.sub(va,vo), R3.sub(vb,vo))
    DO
      IF normal = R3.Zero THEN 
        RETURN R3.Zero
      ELSE
        RETURN R3.reduce(normal)
      END
    END
  END FaceNormal;


PROCEDURE NormalizeCoords(READONLY t: Topology; VAR coords: Coords) =

  PROCEDURE Barycenter(): R3.T =
    VAR B: R3.T := R3.T{0.0, ..};
        N: CARDINAL := 0;
    BEGIN
      FOR i := 0 TO LAST(coords) DO 
        WITH v = t.vertex[i], vc = coords[i] DO
          IF NodeExists(v) THEN B := R3.add(B, vc); INC(N) END;
        END;
      END;
      RETURN R3.scale( B, 1.0 /FLOAT(N) )
    END Barycenter;

  PROCEDURE MeanDist(): REAL =
    VAR S: LONGREAL := 0.0d0;
        N: CARDINAL := 0;
    BEGIN
      FOR i := 0 TO LAST(coords) DO 
        WITH v = t.vertex[i], vc = coords[i] DO
          IF NodeExists(v) THEN S := S + R3.ddot(vc, vc); INC(N) END;
        END;
      END;
      RETURN FLOAT(Math.sqrt(S/FLOAT(N,LONGREAL) ) )
    END MeanDist;

  PROCEDURE Displace(d: R3.T) =
    BEGIN
      FOR i := 0 TO LAST(coords) DO 
        WITH v = coords[i] DO
          v := R3.add(v, d);
        END;
      END;
    END Displace;

  PROCEDURE Reduce(d: REAL) =
    BEGIN
      FOR i := 0 TO LAST(coords) DO 
        WITH v = coords[i] DO
          v := R3.scale(v, 1.0 / d);
        END;
      END;
    END Reduce;

  BEGIN
    WITH b = Barycenter() DO
      Displace(R3.neg(b))
    END;

    WITH s = MeanDist() DO
      Reduce(s)
    END;
  END NormalizeCoords;

PROCEDURE CountEdges(a: Arc): CARDINAL =
  VAR ne: CARDINAL := 0;
  PROCEDURE Visit(<*UNUSED*> c: Arc) =
    BEGIN
      INC(ne, 1); 
    END Visit;
    
  BEGIN
    Enum(a, Visit, edges := TRUE);
    RETURN ne; 
  END CountEdges;

PROCEDURE NumberVertices(a: Arc): CARDINAL = 
  VAR nv: CARDINAL := 0;
  PROCEDURE Visit(c: Arc) =
    VAR cn: Arc;
    BEGIN
      cn := c;
      REPEAT
        NARROW(cn.edge, Edge).org[RotBits(cn)].num := nv;
        cn := Onext(cn);
      UNTIL (cn = c);
      INC(nv,1);
    END Visit;
    
  BEGIN
    EnumVertices(a, Visit);
    RETURN nv
  END NumberVertices;

PROCEDURE InitCoords(VAR coords: Coords) =
  BEGIN
    FOR i := 0 TO LAST(coords) DO
      WITH
        r = Random.Real(), s = Random.Real(), t =  Random.Real()
      DO
        coords[i] := R3.T{r, s, t}
      END;
    END;
  END InitCoords;

CONST boole = ARRAY BOOLEAN OF CHAR {'F', 'T'};

PROCEDURE MakeTopology(a: Arc): Topology =
  VAR NV: CARDINAL := NumberVertices(a); 
      NF: CARDINAL := NumberVertices(Rot(a));
      NE: CARDINAL := CountEdges(a);    
      adj := NEW(REF ARRAY OF ARRAY OF BOOLEAN, NV, NV); 
      vertex := NEW(REF ARRAY OF Vertex, NV);
      face := NEW(REF ARRAY OF Face, NF);
      edge := NEW(REF ARRAY OF Arc, NE);
      out := NEW(REF ARRAY OF Arc, NV); 
      side := NEW(REF ARRAY OF Arc, NF); 

      t := Topology{ NV , NF , NE, adj, vertex, face, edge, out, side};
      ne: CARDINAL := 0;
 
  PROCEDURE Visit(c: Arc) =
    VAR cn: Arc;
    BEGIN
      WITH
        cv = Org(c),  
        i = cv.num
      DO
        t.vertex[i] := cv; 
        t.out[i] := c;
        cn := c;
        REPEAT
          WITH
            j = NARROW(cn.edge, Edge).org[RotBits(Sym(cn))].num
          DO
            IF t.adj[i,j] = FALSE THEN 
              t.adj[i,j] := TRUE;
              t.adj[j,i] := TRUE;
              t.edge[ne] := cn;
              SetEdgeNo(t.edge[ne].edge, ne);
              INC(ne, 1);
            END;
            cn := Onext(cn);
          END; 
        UNTIL (cn = c);
      END;  
    END Visit;
  
  PROCEDURE VisitF(c: Arc) = 
    BEGIN
      WITH
        cf = Org(c),  
        i = cf.num 
      DO
        t.side[i] := c;
        t.face[i] := cf 
      END;
    END VisitF;
  
  BEGIN
    EnumVertices(a, Visit);
    t.NE := ne;  
    EnumVertices(Rot(a), VisitF);
    (* 
      Este c'odigo supoe que nao existem lacos nem
      arestas paralelas *)
    RETURN t
  END MakeTopology;
  
PROCEDURE AddSpokes(a: Arc): Arc =
  VAR e, b, c, d: Arc;
  BEGIN
    WITH v = MakeVertex() DO 
      e := Make();
      Splice(e, a); 
      SetOrg(Sym(e), v);
      SetOrg(e, Org(a));

      b := a;
      d := e;
      c := Lnext(a);
      WHILE c # e DO
        WITH
          r = Connect(c, Sym(d)),
          f = MakeFace()
        DO
          SetOrg(r, Org(c));
          SetOrg(Sym(r), v);
          SetLeft(r, f);
          SetLeft(Sym(d), f);
          SetLeft(b, f);
          b := c;
          d := r;
          c := Lnext(c)
        END
      END;

      WITH
        f = MakeFace()
      DO
        SetLeft(e, f);
        SetLeft(Sym(d), f);
        SetLeft(b, f);
      END;

      RETURN e
    END
  END AddSpokes;

PROCEDURE MakeGrid(order: CARDINAL): ARRAY [0..7] OF Arc =
      
  VAR EdgeCount: CARDINAL := 0;
  
  PROCEDURE MakeCell(s: Arc): Arc =
    (* 
      Adds four triangles nested in the corner between "s" 
      (on the left side, vertical, upwards) and "Oprev(s)" 
      (on the bottom, horizontal, to the right).
      Returns the right edge of the new cell, directed upwards. *)
    VAR b, c, t: Arc;
        v: Vertex;
    BEGIN
      v := MakeVertex();
      t := Sym(Oprev(s));

      b := Make(EdgeCount); INC(EdgeCount);
      Splice(Oprev(t), b);
      SetOrg(b, Org(t));
      SetOrg(Sym(b), v);
      
      c := Make(EdgeCount); INC(EdgeCount);
      Splice(c, Sym(b));
      Splice(Sym(s), Sym(c));
      SetOrg(c, v);
      SetOrg(Sym(c), Org(Sym(s)));

      EVAL AddSpokes(b);

      RETURN b
    END MakeCell;
  
  PROCEDURE MakeCellRow(a: Arc): Arc =
  (* 
    Adds a new row of cells, bounded on the left by arc "a"
    (vertical, upwards).  Returns the right edge of the 
    rightmost cell (vertical, directed upwards). *)
    BEGIN
      FOR col := 0 TO order-1 DO
        a := MakeCell(a);
      END;
      RETURN a
    END MakeCellRow;

  VAR b, t: Arc;
      ca: ARRAY [0..7] OF Arc;
  BEGIN
    (* Create bottom row of edges: *)  
    FOR col := 0 TO order-1 DO
      WITH c = Make(EdgeCount) DO
        INC(EdgeCount);
        IF col = 0 THEN
          WITH v = MakeVertex() DO
            SetOrg(c, v);
          END;
          t := c;
        ELSE
          Splice(Sym(b), c);
          SetOrg(c, Org(Sym(b)));
        END;
        b := c;
        WITH w = MakeVertex() DO
          SetOrg(Sym(b), w)
        END;
      END
    END;

    ca[1] := Flip(t);
    ca[2] := Sym(b);
    FOR row := 0 TO order-1 DO
      (* "t" is the first edge on the top of the last row of cells *)
      (* "b" is the last edge on the top of the last row of cells *)
      (* Both are pointing from left to right *)
      (* Add another row of cells, beginning just above "t": *)
      WITH a = Make(EdgeCount) DO
        INC(EdgeCount);
        IF row = 0 THEN ca[0] := a END;
        IF row = order-1 THEN ca[7] := Flip(Sym(a)) END;

        Splice(a, t);
        SetOrg(a, Org(t));
        WITH u = MakeVertex() DO
          SetOrg(Sym(a), u)
        END;
        WITH f = MakeCellRow(a) DO
          t := Oprev(Sym(a));
          b := Sym(Onext(Sym(f)));

          IF row = 0 THEN ca[3] := Flip(f) END;
          IF row = order-1 THEN ca[4] := Sym(f) END;
        END
      END
    END;
    ca[6] := t;
    ca[5] := Flip(Sym(b));
    RETURN ca
  END MakeGrid;
  
PROCEDURE SetQuarterGridProperties(
    c: Arc; 
    order: CARDINAL; 
    vertexColor: Color.T;
    vertexRadius: REAL;
    edgeColor: Color.T;
    edgeRadius: REAL;
    facePatch: CARDINAL;
    faceColor: Color.T;
    faceTransp: Color.T;
  ) =
   
  PROCEDURE SetCornerVertex(v: Vertex) =
    BEGIN
      SetNodeColor(v, vertexColor);
      SetVertexRadius(v, vertexRadius);
    END SetCornerVertex;
  
  PROCEDURE SetQuarterEdge(a: Arc) =
    BEGIN
      FOR i := 0 TO order DO 
        IF i # 0 THEN 
          WITH v = Org(a) DO
            SetNodeColor(v, edgeColor);
            SetVertexRadius(v, edgeRadius);
          END;
        END;
        SetEdgeColor(a, edgeColor);
        SetEdgeRadius(a, edgeRadius);
        a := Onext(Onext(Sym(a)));
        IF i MOD 2 # 0 THEN a := Onext(Onext(a)) END;
      END;
    END SetQuarterEdge;
    
  PROCEDURE SetLeftTriangle(t: Face) =
    BEGIN
      SetNPatch(t, facePatch);
      SetNodeColor(t, faceColor);
      SetFaceTransp(t, faceTransp);
    END SetLeftTriangle;
  
  PROCEDURE SetRowPatches(d: Arc; row: CARDINAL) =
    BEGIN
      d := Sym(d);
      SetLeftTriangle(Left(d));
      d := Onext(d);
      FOR col := 1 TO order-row-1 DO 
        d := Onext(d);
        SetLeftTriangle(Left(d));
        d := Sym(d);
        SetLeftTriangle(Left(d));
        d := Onext(Flip(d));
      END
    END SetRowPatches;
    
  BEGIN
    SetCornerVertex(Org(c));
    SetQuarterEdge(Oprev(c));
    FOR row := 0 TO order-1 DO 
      SetRowPatches(c, row);
      c := Oprev(Sym(c))
    END
  END SetQuarterGridProperties;
 
CONST
  Digits = SET OF CHAR{'0'..'9'};
  AlphaChars = SET OF CHAR{'\t', '\n', '\r', '\f', ' ', 'A'..'Z', 'a'..'z'};

PROCEDURE ReadTopTri(shape: TEXT): Topology =
  <* FATAL Rd.Failure, Rd.EndOfFile, Thread.Alerted, Convert.Failed, Lex.Error *>
  VAR NV: INTEGER := 0;
      NF: CARDINAL := 0;
      NE: CARDINAL := 0;

  BEGIN
    WITH 
      name = shape & ".top",
      rd = FileStream.OpenRead(name)
    DO 
      EVAL LoadDoc(rd);
      Lex.Skip(rd, cs := AlphaChars);
      NV := Lex.Int(rd, cs := Digits);
      Lex.Skip(rd, cs := AlphaChars);
      NF := Lex.Int(rd, cs := Digits);
      Lex.Skip(rd, cs := AlphaChars);
      NE := Lex.Int(rd, cs := Digits);
      Lex.Skip(rd);
      WITH
        map = NEW(REF ARRAY OF Oct.Edge, NE),
        adj = NEW(REF ARRAY OF ARRAY OF BOOLEAN, NV, NV),
        vertex = NEW(REF ARRAY OF Vertex, NV),
        face = NEW(REF ARRAY OF Face, NF),
        edge = NEW(REF ARRAY OF Arc, NE),
        out = NEW(REF ARRAY OF Arc, NV),
        side = NEW(REF ARRAY OF Arc, NF),
        t = Topology{ NV , NF , NE, adj, vertex, face, edge, out, side}
      DO

        FOR i := 0 TO t.NE-1 DO 
          t.edge[i] := Make(i);
          map[i] := t.edge[i].edge
        END;
        FOR i := 0 TO NV-1 DO t.vertex[i] := MakeVertex(); END;
        FOR i := 0 TO NF-1 DO t.face[i] := MakeFace() END;

        EVAL LoadDoc(rd);
        FOR j := 0 TO t.NV-1 DO 
          Lex.Skip(rd);
          WITH 
            nv = Lex.Int(rd, cs := Digits),
            v = t.vertex[nv]
          DO 
            v.num := nv;
            Lex.Skip(rd); 
            v.exists := Lex.Bool(rd, cs := Lex.BoolChars);
            Lex.Skip(rd);
            v.fixed := Lex.Bool(rd);
            Lex.Skip(rd);
            v.coords[0] := Lex.Real(rd);
            Lex.Skip(rd);
            v.coords[1] := Lex.Real(rd);
            Lex.Skip(rd);
            v.coords[2] := Lex.Real(rd);
            Lex.Skip(rd);
            v.color.r := Lex.Real(rd);
            Lex.Skip(rd);
            v.color.g := Lex.Real(rd);
            Lex.Skip(rd);
            v.color.b := Lex.Real(rd); 
            Lex.Skip(rd);
            v.radius := Lex.Real(rd);
          END;      
        END;
        Lex.Skip(rd);
        EVAL LoadDoc(rd);
        FOR j := 0 TO t.NF-1 DO 
          Lex.Skip(rd);
          WITH 
            nf = Lex.Int(rd, cs := Digits),
            f = t.face[nf]
          DO
            f.num := nf;
            Lex.Skip(rd);
            f.exists := Lex.Bool(rd);
            Lex.Skip(rd);
            f.color.r := Lex.Real(rd);
            Lex.Skip(rd);
            f.color.g := Lex.Real(rd);
            Lex.Skip(rd);
            f.color.b := Lex.Real(rd);
            Lex.Skip(rd);
            f.transp.r := Lex.Real(rd);
            Lex.Skip(rd);
            f.transp.g := Lex.Real(rd);
            Lex.Skip(rd);
            f.transp.b := Lex.Real(rd);
            Lex.Skip(rd);
            f.patch := Lex.Int(rd);
          END;
        END;
        Lex.Skip(rd);
        EVAL LoadDoc(rd);
        FOR k := 0 TO t.NE-1 DO 
          Lex.Skip(rd);
          WITH 
            ne = Lex.Int(rd, cs := Digits),
            e = NARROW(t.edge[ne].edge, Edge)
          DO
            <* ASSERT t.edge[ne].bits = 0 *>
            Oct.ReadEdge(rd, e, map^);
            FOR l := 0 TO 3 DO 
              Lex.Skip(rd);
              WITH 
                n = Lex.Int(rd, cs := Digits),
                vf = e.org[l]
              DO 
                IF Rd.GetChar(rd) = 'v' THEN 
                  vf := t.vertex[n];
                  t.out[n] := Arc{edge := e, bits := 2*l};
                ELSE
                  vf := t.face[n];
                  t.side[n] := Arc{edge := e, bits := 2*l};                  
                END;
              END;
            END;
            Lex.Skip(rd);
            e.exists := Lex.Bool(rd);
            Lex.Skip(rd);
            e.spring := Lex.Bool(rd);
            Lex.Skip(rd);
            e.color.r := Lex.Real(rd);
            Lex.Skip(rd);
            e.color.g := Lex.Real(rd);
            Lex.Skip(rd);
            e.color.b := Lex.Real(rd);
            Lex.Skip(rd);
            e.radius := Lex.Real(rd); 
          END;
        END;
        CheckOutAndSide(t);
        ComputeAdjacencyMatrix(t);
        Rd.Close(rd); 
        RETURN t;    
      END;
    END;  
  END ReadTopTri;
  
PROCEDURE CheckOutAndSide(READONLY t: Topology) =
  BEGIN
    FOR i := 0 TO t.NV-1 DO 
      WITH
        v = t.vertex[i],
        e = t.out[i]
      DO
        <* ASSERT v.num = i *>
        <* ASSERT Org(e) = v *>
      END
    END;
    FOR i := 0 TO t.NF-1 DO
      WITH
        f = t.face[i],
        e = t.side[i]
      DO
        <* ASSERT f.num = i *>
        <* ASSERT Org(e) = f *>
      END
    END
  END CheckOutAndSide;
  
PROCEDURE ComputeAdjacencyMatrix(READONLY t: Topology) =
  VAR p: Arc;
  BEGIN
    FOR i := 0 TO t.NV-1 DO 
      WITH ta = t.adj[i] DO 
        FOR j := 0 TO t.NV-1 DO 
          ta[j] := FALSE
        END
      END
    END;
    FOR i := 0 TO t.NV-1 DO 
      WITH
        v = t.vertex[i],
        e = t.out[i]
      DO
        <* ASSERT Org(e) = v *>
        p := e;
        REPEAT
          WITH
            j = Org(Sym(p)).num
          DO
            <* ASSERT Org(Sym(p)) = t.vertex[j] *>
            t.adj[i,j] := TRUE;
            t.adj[j,i] := TRUE;
            p := Onext(p);
          END; 
        UNTIL p = e;
      END
    END
  END ComputeAdjacencyMatrix;

PROCEDURE PrintTopTri(t: Topology; shape: TEXT) =
  VAR comment: TEXT;
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    WITH 
      name = shape & ".top",
      wr = FileStream.OpenWrite(name)
    DO 
      comment := "Description file of the topology of " & shape;
      DumpDoc(wr, comment);
      Wr.PutText(wr, "|\n");
      Wr.PutText(wr, "vertices "); Wr.PutText(wr, Fmt.Int(t.NV) & "\n");  
      Wr.PutText(wr, "faces "); Wr.PutText(wr, Fmt.Int(t.NF) & "\n"); 
      Wr.PutText(wr, "edges "); Wr.PutText(wr, Fmt.Int(t.NE) & "\n");
      Wr.PutText(wr, "|\n");
      comment := "vertices:";
      DumpDoc(wr, comment);
      Wr.PutText(wr, "|\n");
      FOR i := 0 TO t.NV-1 DO
        WITH v = t.vertex[i] DO 
          Wr.PutText(wr, Fmt.Int(v.num) & " "); 
          Wr.PutText(wr, Fmt.Char(boole[v.exists]) & " ");
          Wr.PutText(wr, Fmt.Char(boole[v.fixed]) & " ");
          Wr.PutText(wr, Fmt.Real(v.coords[0]) & " ");
          Wr.PutText(wr, Fmt.Real(v.coords[1]) & " ");
          Wr.PutText(wr, Fmt.Real(v.coords[2]) & " ");
          Wr.PutText(wr, Fmt.Real(v.color.r, precision := 4) & " ");  
          Wr.PutText(wr, Fmt.Real(v.color.g, precision := 4) & " ");  
          Wr.PutText(wr, Fmt.Real(v.color.b, precision := 4) & " "); 
          Wr.PutText(wr, Fmt.Real(v.radius, precision := 4) & "\n")
        END       
      END;
      Wr.PutText(wr, "|\n");
      comment := "Faces:";
      DumpDoc(wr, comment);
      Wr.PutText(wr, "|\n");
      FOR i := 0 TO t.NF-1 DO
        WITH f = t.face[i] DO 
          Wr.PutText(wr, Fmt.Int(f.num) & " "); 
          Wr.PutText(wr, Fmt.Char(boole[f.exists]) & " ");
          Wr.PutText(wr, Fmt.Real(f.color.r, precision := 4) & " ");  
          Wr.PutText(wr, Fmt.Real(f.color.g, precision := 4) & " ");  
          Wr.PutText(wr, Fmt.Real(f.color.b, precision := 4) & " "); 
          Wr.PutText(wr, Fmt.Real(f.transp.r, precision := 4) & " ");  
          Wr.PutText(wr, Fmt.Real(f.transp.g, precision := 4) & " ");  
          Wr.PutText(wr, Fmt.Real(f.transp.b, precision := 4) & " ");  
          Wr.PutText(wr, Fmt.Int(f.patch) & "\n")
        END
      END;
      Wr.PutText(wr, "|\n");
      comment := "Edges:";
      DumpDoc(wr, comment);
      Wr.PutText(wr, "|\n");
      FOR i := 0 TO t.NE-1 DO
        WITH e = t.edge[i] DO 
          Wr.PutText(wr, Fmt.Int(GetEdgeNo(e.edge)) & " "); 
          Oct.PrintEdge(wr, e.edge);
          Wr.PutText(wr, " ");
          FOR j := 0 TO 3 DO
            WITH n = NARROW(e.edge, Edge).org[j] DO 
              TYPECASE n OF 
              | NULL => <* ASSERT FALSE *>
              | Vertex(v) => Wr.PutText(wr, Fmt.Int(v.num) & "v ");
              | Face(f) => Wr.PutText(wr, Fmt.Int(f.num) & "f ");
              ELSE <* ASSERT FALSE *>
              END;
            END;
          END;
          Wr.PutText(wr, Fmt.Char(boole[NARROW(e.edge, Edge).exists]) & " ");
          Wr.PutText(wr, Fmt.Char(boole[NARROW(e.edge, Edge).spring]) & " ");
          Wr.PutText(wr, Fmt.Real(NARROW(e.edge, Edge).color.r, precision := 4) & " ");  
          Wr.PutText(wr, Fmt.Real(NARROW(e.edge, Edge).color.g, precision := 4) & " ");  
          Wr.PutText(wr, Fmt.Real(NARROW(e.edge, Edge).color.b, precision := 4) & " "); 
          Wr.PutText(wr, Fmt.Real(NARROW(e.edge, Edge).radius, precision := 4));
          Wr.PutText(wr, "\n");
        END;
      END;       
      Wr.Close(wr); 
    END;     
  END PrintTopTri;
  
PROCEDURE PutRay(
    wr: Wr.T;
    READONLY t: Topology; 
    READONLY coords: Coords;
  ) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  
  PROCEDURE PutRayCoord(c: REAL) =
    BEGIN
      Wr.PutText(wr, Fmt.Real(c, 4, Fmt.Style.Flo))
    END PutRayCoord;

  PROCEDURE PutRayPoint(p: R3.T) =
    BEGIN
      Wr.PutText(wr, "<");
      PutRayCoord(p[1]); 
      Wr.PutText(wr, ",");
      PutRayCoord(p[2]); 
      Wr.PutText(wr, ",");
      PutRayCoord(-p[0]);
      Wr.PutText(wr, ">");
    END PutRayPoint;
  
  PROCEDURE PutRayColor(cr: Color.T) =
    BEGIN
      Wr.PutText(wr, "rgb <");
      Wr.PutText(wr, Fmt.Real(cr.r) & "," & Fmt.Real(cr.g) & "," & Fmt.Real(cr.b));
      Wr.PutText(wr, ">");
    END PutRayColor;
  
  PROCEDURE PutRayColorTransp(cr: Color.T; tr: REAL) =
    BEGIN
      Wr.PutText(wr, "rgbf <");
      Wr.PutText(wr, Fmt.Real(cr.r) & "," & Fmt.Real(cr.g) & "," & Fmt.Real(cr.b));
      Wr.PutText(wr, "," & Fmt.Real(tr));
      Wr.PutText(wr, ">");
    END PutRayColorTransp;
  
  PROCEDURE PutRayCylinder(READONLY o, d: R3.T; radius: REAL; READONLY cr: Color.T) =
    BEGIN
      Wr.PutText(wr, "  cylinder {\n");
      Wr.PutText(wr, "    "); PutRayPoint(o); Wr.PutText(wr, ",\n");
      Wr.PutText(wr, "    "); PutRayPoint(d); Wr.PutText(wr, ",\n");
      Wr.PutText(wr, "    " & Fmt.Real(radius) & "\n");
      Wr.PutText(wr, "    open\n");
      Wr.PutText(wr, "    pigment { color "); PutRayColor(cr); Wr.PutText(wr, " }\n");
      Wr.PutText(wr, "  }\n");
      Wr.PutText(wr, "\n");
      Wr.Flush(wr);
    END PutRayCylinder;

  PROCEDURE PutRaySphere(READONLY p: R3.T; radius: REAL; READONLY cr: Color.T) =
    BEGIN
      Wr.PutText(wr, "  sphere {\n");
      Wr.PutText(wr, "    "); PutRayPoint(p);
        Wr.PutText(wr, ", " & Fmt.Real(radius) & "\n");
      Wr.PutText(wr, "    pigment { color "); PutRayColor(cr); Wr.PutText(wr, " }\n");
      Wr.PutText(wr, "  } \n");
      Wr.PutText(wr, "\n");
      Wr.Flush(wr);
    END PutRaySphere;

  <*UNUSED*>
  PROCEDURE PutRayTriangle(READONLY a, b, c: R3.T; READONLY cr: Color.T; tr: REAL) =
    BEGIN
      Wr.PutText(wr, "  triangle {");
      Wr.PutText(wr, "    "); PutRayPoint(a); Wr.PutText(wr, ",\n");
      Wr.PutText(wr, "    "); PutRayPoint(b); Wr.PutText(wr, ",\n");
      Wr.PutText(wr, "    "); PutRayPoint(c); Wr.PutText(wr, "\n");
      Wr.PutText(wr, "    pigment { color "); 
        PutRayColorTransp(cr, tr); 
        Wr.PutText(wr, " }\n");
      Wr.PutText(wr, "  } \n");
      Wr.PutText(wr, "\n");
      Wr.Flush(wr);
    END PutRayTriangle;
    
  PROCEDURE PutRaySmoothTriangle(
      READONLY a, an: R3.T;
      b, bn: R3.T;
      c, cn: R3.T; 
      READONLY cr: Color.T; 
      tr: REAL
    ) =
    BEGIN
      Wr.PutText(wr, "  smooth_triangle {\n");
      Wr.PutText(wr, "    "); PutRayPoint(a); Wr.PutText(wr, ", ");
        PutRayPoint(an); Wr.PutText(wr, ",\n");
      Wr.PutText(wr, "    "); PutRayPoint(b); Wr.PutText(wr, ", ");
        PutRayPoint(bn); Wr.PutText(wr, ",\n");
      Wr.PutText(wr, "    "); PutRayPoint(c); Wr.PutText(wr, ", ");
        PutRayPoint(cn); Wr.PutText(wr, "\n");
      Wr.PutText(wr, "    pigment { color "); 
        PutRayColorTransp(cr, tr); 
        Wr.PutText(wr, " }\n");
      Wr.PutText(wr, "  } \n");
      Wr.PutText(wr, "\n");
      Wr.Flush(wr);
    END PutRaySmoothTriangle;
    
  BEGIN
    FOR i := 0 TO t.NE-1 DO 
      WITH 
        e = t.edge[i],
        er = NARROW(e.edge, Edge)
      DO 
        IF EdgeExists(er) AND er.radius > 0.0 THEN 
          WITH o = NOrg(e), d = NOrg(Sym(e)) DO    
            PutRayCylinder(coords[o], coords[d], er.radius, er.color)
          END
        END;
      END;
    END;
    
    FOR i := 0 TO t.NV-1 DO 
      WITH v = t.vertex[i] DO 
        IF NodeExists(v) AND v.radius > 0.0 THEN 
          PutRaySphere(coords[i], v.radius, v.color)
        END
      END
    END;

    WITH
      vn = ComputeAllVertexNormals(t, coords)^ 
    DO
      FOR i := 0 TO t.NF-1 DO 
        WITH f = t.face[i] DO
          IF NodeExists(f) THEN 
            WITH 
              ea = t.side[i],  a = NOrg(Tor(ea)),
              eb = Onext(ea),  b = NOrg(Tor(eb)),
              ec = Onext(eb),  c = NOrg(Tor(ec)),
              t3 = f.transp,
              transp = (t3.r + t3.g + t3.b) / 3.0
            DO    
              PutRaySmoothTriangle(
                coords[a], vn[a],
                coords[b], vn[b],
                coords[c], vn[c],
                f.color, transp
              )
            END
          END
        END
      END;
    END;

    Wr.Flush(wr);

  END PutRay;
  
PROCEDURE ComputeAllVertexNormals(
    READONLY t: Topology; 
    READONLY coords: Coords;
  ): REF ARRAY OF R3.T =
  BEGIN
    WITH
      rvn = NEW(REF ARRAY OF R3.T, t.NV),
      vn = rvn^
    DO
      FOR i := 0 TO t.NV-1 DO 
        vn[i] := VertexNormal(t.out[i], coords)
      END;
      RETURN rvn
    END;
  END ComputeAllVertexNormals; 
  
PROCEDURE VertexNormal(a: Arc; READONLY coords: Coords): R3.T =
  VAR an, ao: Arc;
      normal:= R3.T{0.0, ..};
  BEGIN
    ao := Oprev(a);     
    an := a;
    WITH
      uc = coords[NOrg(a)]
    DO
      REPEAT
        WITH
          aov = Org(Sym(ao)), aon = aov.num, (* to get data value of Dest(ao) *)  
          anv = Org(Sym(an)), ann = anv.num  (* to get data value of Dest(an) *)
        DO
          IF NodeExists(aov) AND NodeExists(anv) THEN
            normal := R3.add(normal, 
              R3.cross(
                R3.sub(coords[aon], uc), 
                R3.sub(coords[ann], uc)
              )
            )
          END
        END;
        ao := an;
        an := Onext(an);
      UNTIL (an = a);
    END;
    IF normal = R3.Zero THEN 
      RETURN R3.Zero
    ELSE
      RETURN R3.reduce(normal)
    END
  END VertexNormal;

PROCEDURE PutWire(
    wr: Wr.T; 
    READONLY t: Topology; 
    READONLY coords: Coords;
    all: BOOLEAN;
  ) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  
  PROCEDURE PutWireColor(cr: Color.T) =
    BEGIN
      Wr.PutText(wr, Fmt.Int(ROUND(255.0*cr.r)));
      Wr.PutText(wr, " ");
      Wr.PutText(wr, Fmt.Int(ROUND(255.0*cr.g)));
      Wr.PutText(wr, " ");
      Wr.PutText(wr, Fmt.Int(ROUND(255.0*cr.b)));
    END PutWireColor;
    
  PROCEDURE PutWirePoint(p: R3.T) =
    BEGIN
      Wr.PutText(wr, Fmt.Real(p[0]*1000.0));
      Wr.PutText(wr, " ");
      Wr.PutText(wr, Fmt.Real(p[1]*1000.0));
      Wr.PutText(wr, " ");
      Wr.PutText(wr, Fmt.Real(p[2]*1000.0));
    END PutWirePoint;
    
  BEGIN
    FOR i := 0 TO t.NF-1 DO 
      IF all OR NodeExists(t.face[i]) THEN 
        WITH
          c = t.side[i],
          k = NOrg(Tor(c)),
          d = Onext(c), l = NOrg(Tor(d)),
          e = Onext(d), m = NOrg(Tor(e)),
          color = Org(c).color 
        DO
          Wr.PutText(wr, "# face " & Fmt.Int(i) & "\n");
          Wr.PutText(wr, "  "); PutWireColor(color); Wr.PutText(wr, "# color \n\n");
          Wr.PutText(wr, "  "); PutWirePoint(coords[k]); Wr.PutText(wr, "\n");
          Wr.PutText(wr, "  "); PutWirePoint(coords[l]); Wr.PutText(wr, "\n");
          Wr.PutText(wr, "  "); PutWirePoint(coords[m]); Wr.PutText(wr, "\n");
          Wr.PutText(wr, "\n");
          Wr.Flush(wr);
        END;
      END;
    END;
  END PutWire;

  
PROCEDURE WriteAnima(
    wr1, wr2: Wr.T; 
    READONLY top: Topology; 
    READONLY c: Coords;
    all: BOOLEAN;
  ) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  
  PROCEDURE WriteAnimaCoord(x: REAL) =
    BEGIN
      Wr.PutText(wr2, Fmt.Real(x*1000.0))
    END WriteAnimaCoord;
    
  PROCEDURE WriteAnimaPoint(READONLY p: R3.T) =
    BEGIN
      WriteAnimaCoord(p[0]);
      Wr.PutText(wr2, " ");
      WriteAnimaCoord(p[1]);
      Wr.PutText(wr2, " ");
      WriteAnimaCoord(p[2]);
    END WriteAnimaPoint;
  
  PROCEDURE WriteAnimaSu() =  
    BEGIN
      Wr.PutText(wr1, "begin surface");
      Wr.PutText(wr1, "\n");
      Wr.PutText(wr1, "triangles = " & Fmt.Int(top.NF));
      Wr.PutText(wr1, "\n");

      FOR i := 0 TO top.NF-1 DO 
        IF all OR top.face[i].exists THEN 
          WITH
            ec = Rot(top.side[i]),
            ed = Lnext(ec), 
            ee = Lnext(ed), 
            k = Org(ec).num,
            l = Org(ed).num,
            m = Org(ee).num
          DO
            IF ed # ec AND ee # ec AND Lnext(ee) = ec THEN
              (* Face is a triangle, we can show it: *)
              Wr.PutText(wr1, Fmt.Int(k) & " ");
              Wr.PutText(wr1, Fmt.Int(l) & " ");
              Wr.PutText(wr1, Fmt.Int(m) & " ");
              Wr.PutText(wr1, "   0"); (* One type material *)
              Wr.PutText(wr1, "   0 0 0");           
              Wr.PutText(wr1, "\n");
              Wr.Flush(wr1);
            END
          END;
        END;
      END;
      Wr.PutText(wr1, "end surface");
      Wr.PutText(wr1, "\n");
    END WriteAnimaSu;
  
  PROCEDURE WriteAnimaSt() =  
    BEGIN
      Wr.PutText(wr2, "begin state");
      Wr.PutText(wr2, "\n");
      Wr.PutText(wr2, "vertices = " & Fmt.Int(top.NV));
      Wr.PutText(wr2, "\n");
      Wr.PutText(wr2, "t = 0");
      Wr.PutText(wr2, "\n");

      FOR i := 0 TO top.NV-1 DO 
        IF all OR (top.vertex[i].exists AND top.vertex[i].radius > 0.0) THEN 
          WITH
            v = top.vertex[i],
            ec = top.out[i],
            ed = Lnext(ec), 
            ee = Lnext(ed)
          DO
            IF ed # ec AND ee # ec AND Lnext(ee) = ec THEN
              (* Vertex belong to a triangle, we can show it: *)
              Wr.PutText(wr2, Fmt.Int(v.num)); 
              Wr.PutText(wr2, "  "); WriteAnimaPoint(c[v.num]);
              Wr.PutText(wr2, "   0 0 0");           
              Wr.PutText(wr2, "\n");
              Wr.Flush(wr2);
            END
          END;
        END;
      END;
      Wr.PutText(wr2, "end state");
      Wr.PutText(wr2, "\n");
    END WriteAnimaSt;
  BEGIN
    WriteAnimaSu();
    WriteAnimaSt()
  END WriteAnima;

BEGIN

END Triang.