(* 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;
    
(* === INIT METHODS === *)

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 NodeInit(e: Node; num: CARDINAL; exists: BOOLEAN; color: Color.T) =
  BEGIN
    e.num := num;
    e.exists := exists;
    e.color := color
  END NodeInit;

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

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

(* === ELEMENT CREATION === *)
  
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 MakeVertex(): Vertex =
  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
      VertexInit(e, 0, TRUE, FALSE, color, coords, radius);
      RETURN e
    END;
  END MakeVertex;

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
      FaceInit(e, 0, TRUE, color, transp, patch);
      RETURN e
    END;  
  END MakeFace;
 
(* === ARC PROPERTIES === *)

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 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;
  
(* === EDGE PROPERTIES: === *)

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

PROCEDURE SetEdgeExists(e: Edge; exists: BOOLEAN) =
  BEGIN
    e.exists := exists
  END SetEdgeExists;

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

PROCEDURE SetEdgeSpring(e: Edge; spring: BOOLEAN) =
  BEGIN
    e.spring := spring
  END SetEdgeSpring;
  
PROCEDURE EdgeColor(e: Edge): Color.T =
  BEGIN
    RETURN e.color
  END EdgeColor; 

PROCEDURE SetEdgeColor(e: Edge; color: Color.T) =
  (* Set edge color *)    
  BEGIN
    e.color := color
  END SetEdgeColor;
  
PROCEDURE EdgeRadius(e: Edge): REAL =
  BEGIN
    RETURN e.radius
  END EdgeRadius;
  
PROCEDURE SetEdgeRadius(e: Edge; radius: REAL) =
  BEGIN
    e.radius := radius
  END SetEdgeRadius;
 
(* === NODE PROPERTIES === *)

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

PROCEDURE SetNodeExists(n: Node; b: BOOLEAN) =
  BEGIN
    n.exists := b
  END SetNodeExists;

PROCEDURE NodeColor(n: Node): Color.T =
  BEGIN
    RETURN n.color
  END NodeColor; 

PROCEDURE SetNodeColor(n: Node; color: Color.T) =
  BEGIN
    n.color := color
  END SetNodeColor;
  
PROCEDURE NNode(n: Node): CARDINAL =
  BEGIN
    RETURN n.num
  END NNode;

PROCEDURE SetNNode(n: Node; num: CARDINAL) =
  BEGIN
    n.num := num
  END SetNNode;

(* === VERTEX PROPERTIES === *)

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

PROCEDURE SetVertexFixed(v: Vertex; fixed: BOOLEAN) =
  BEGIN
    v.fixed := fixed
  END SetVertexFixed;

PROCEDURE VertexRadius(v: Vertex): REAL =
  BEGIN
    RETURN v.radius
  END VertexRadius; 

PROCEDURE SetVertexRadius(v: Vertex; radius: REAL) =
  BEGIN
    v.radius := radius
  END SetVertexRadius;
  
PROCEDURE VertexCoords(v: Vertex): R3.T =
  BEGIN
    RETURN v.coords
  END VertexCoords;

PROCEDURE SetVertexCoords(v: Vertex; coords: R3.T) =
  BEGIN
    v.coords := coords
  END SetVertexCoords; 

(* === FACE PROPERTIES === *)

PROCEDURE NPatch(f: Face): CARDINAL =
  BEGIN
    RETURN f.patch
  END NPatch;
  
PROCEDURE SetNPatch(f: Face; patch: CARDINAL) =
  BEGIN
    f.patch := patch
  END SetNPatch;

PROCEDURE FaceTransp(f: Face): Color.T =
  BEGIN
    RETURN f.transp
  END FaceTransp;

PROCEDURE SetFaceTransp(f: Face; color: Color.T) =
  BEGIN
    f.transp := color
  END SetFaceTransp;

(* === CONSTRUCTION TOOLS === *)

PROCEDURE Connect(a, b: Arc): Arc =
  VAR e: Arc;
  BEGIN
    e := Make();
    Splice(e, a);
    Splice(Sym(e), b);
    RETURN e;
  END Connect;
  
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.edge, edgeColor);
        SetEdgeRadius(a.edge, 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;
 
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;

(* === GLOBAL PROCEDURES === *)

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

(* === GEOMETRIC TOOLS === *)

PROCEDURE InitCoords(
    coins: Random.T; 
    VAR coords: Coords; 
    radius: REAL := 1.0
  ) =

  PROCEDURE RR(): REAL =
    BEGIN
      RETURN radius * (2.0*Random.Real(coins) - 1.0)
    END RR;
    
  BEGIN
    FOR i := 0 TO LAST(coords) DO
      coords[i] := R3.T{RR(), RR(), RR()}
    END;
  END InitCoords;

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(1.0d0/FLOAT(N, LONGREAL), B)
    END Barycenter;

  PROCEDURE MeanDist(): LONGREAL =
    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.Dot(vc, vc); INC(N) END;
        END;
      END;
      RETURN 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: LONGREAL) =
    BEGIN
      FOR i := 0 TO LAST(coords) DO 
        WITH v = coords[i] DO
          v := R3.Scale(1.0d0/d, v);
        END;
      END;
    END Reduce;

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

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

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.Dir(normal)
      END
    END
  END FaceNormal;

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.Dir(normal)
    END
  END VertexNormal;
  
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; 
  
(* === INPUT/OUTPUT === *)

CONST
  Boole = ARRAY BOOLEAN OF CHAR {'F', 'T'};
  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 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 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 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;

BEGIN
END Triang.