MODULE SPTriang;

FROM SPQuad IMPORT 
  MakeEdge, Splice, NumberArcs,
  Sym, Rot, Onext, Oprev, Dnext, Dprev, Lnext, Lprev, Rnext, Rprev,
  Odata, Ddata, Ldata, Rdata, 
  SetOdata, SetDdata, SetLdata, SetRdata, SetEdgeNum,
  NullArc, ToText;
  
IMPORT Wr, Rd, Thread, FileWr, Fmt, PSPlot, HLR3, LR4, LR3, LR3x3, LR3Extras, OSError;
IMPORT FGet, FPut, NGet, NPut, FileFmt;

CONST TriangFileVersion = "1997-12-15";

PROCEDURE NumberVertices(READONLY arc: ARRAY OF Arc): REF ARRAY OF Arc =
  VAR r: REF ARRAY OF Arc := NEW(REF ARRAY OF Arc, 1);
      nVertices: CARDINAL := 0;
      
  PROCEDURE VisitedOrg(a: Arc): BOOLEAN =
    (*
      TRUE if the Org of "a" has been visited already. *)
    BEGIN
      WITH v = Org(a) DO
        RETURN v # NIL
        AND v.num < nVertices
        AND Org(r[v.num]) = v
      END
    END VisitedOrg;
  
  PROCEDURE VisitOrg(a: Arc) =
    BEGIN
      IF nVertices = NUMBER(r^) THEN 
        WITH t = NEW(REF ARRAY OF Arc, 2*nVertices) DO
          SUBARRAY(t^, 0, nVertices) := r^;
          r := t
        END
      END;
      WITH v = Org(a) DO
        v.num := nVertices;
        VAR b: Arc := a;
        BEGIN
          REPEAT
            <* ASSERT Org(b) = v *>
            b := Onext(b)
          UNTIL b = a
        END
      END;
      r[nVertices] := a;
      nVertices := nVertices + 1
    END VisitOrg;
  
  BEGIN
    FOR i := 0 TO LAST(arc) DO
      IF NOT VisitedOrg(arc[i]) THEN VisitOrg(arc[i]) END
    END;
    WITH t = NEW(REF ARRAY OF Arc, nVertices) DO
      t^ := SUBARRAY(r^, 0, nVertices);
      RETURN t
    END
  END NumberVertices;
  
PROCEDURE CollectVertices(e: Arc; VAR out: ARRAY OF Arc) =
  BEGIN
    WITH 
      NSITES = NUMBER(out)
    DO
      FOR i := 0 TO NSITES-1 DO out[i] := NullArc END;
      
      PROCEDURE DoCollect(f: Arc) =
        BEGIN
          WITH
            no = Org(f).num
          DO
            out[no] := f;
            VAR g: Arc := f;
            BEGIN
              REPEAT
                <* ASSERT Org(g) = Org(f) *>
                IF out[Dest(g).num] = NullArc THEN DoCollect(Sym(g)) END;
                g := Onext(g)
              UNTIL g = f
            END
          END;
        END DoCollect;
      BEGIN
        DoCollect(e);
      END 
    END
  END CollectVertices;
  
PROCEDURE NumberTriangles(READONLY arc: ARRAY OF Arc): REF ARRAY OF Arc =
  VAR r: REF ARRAY OF Arc := NEW(REF ARRAY OF Arc, 2);
      nTriangles: CARDINAL := 0;
      
  PROCEDURE VisitedLeft(a: Arc): BOOLEAN =
    (*
      TRUE if the Left face of "a" has been visited already. *)
    BEGIN
      WITH t = Left(a) DO
        RETURN t # NIL
        AND t.num < nTriangles
        AND Left(r[t.num]) = t
      END
    END VisitedLeft;
  
  PROCEDURE VisitLeft(a: Arc) =
    BEGIN
      IF nTriangles = NUMBER(r^) THEN 
        WITH t = NEW(REF ARRAY OF Arc, 2*nTriangles) DO
          SUBARRAY(t^, 0, nTriangles) := r^;
          r := t
        END
      END;
      WITH t = NEW(REF Triangle) DO
        t.num := nTriangles;
        VAR b: Arc := a;
        BEGIN
          REPEAT
            SetLdata(b, t);
            b := Lnext(b)
          UNTIL b = a
        END
      END;
      r[nTriangles] := a;
      nTriangles := nTriangles + 1
    END VisitLeft;
  
  BEGIN
    FOR i := 0 TO LAST(arc) DO
      IF NOT VisitedLeft(arc[i]) THEN VisitLeft(arc[i]) END
    END;
    WITH t = NEW(REF ARRAY OF Arc, nTriangles) DO
      t^ := SUBARRAY(r^, 0, nTriangles);
      RETURN t
    END
  END NumberTriangles;
  
PROCEDURE CollectTriangles(e: Arc; VAR side: ARRAY OF Arc) =
  BEGIN
    WITH 
      NTRIANGLES = NUMBER(side) 
    DO
      FOR i := 0 TO NTRIANGLES-1 DO side[i] := NullArc END;
      
      PROCEDURE DoCollect(f: Arc) =
        BEGIN
          WITH
            ft = Left(f),
            no = ft.num
          DO
            side[no] := f;
            VAR g: Arc := f;
            BEGIN
              REPEAT
                <* ASSERT Left(g) = Left(f) *>
                WITH
                  gt = Right(g)
                DO
                  IF side[gt.num] = NullArc THEN DoCollect(Sym(g)) END
                END;
                g := Lnext(g)
              UNTIL g = f
            END
          END;
        END DoCollect;
      
      BEGIN
        DoCollect(e);
      END 
    END
  END CollectTriangles;
  
PROCEDURE ComputeTriangleMatrices(VAR tri: T) = 
  BEGIN
    tri.b2c := NEW(REF ARRAY OF LR3x3.T, NUMBER(tri.side^));
    tri.c2b := NEW(REF ARRAY OF LR3x3.T, NUMBER(tri.side^));
    WITH
      side = tri.side^,
      b2c = tri.b2c^,
      c2b = tri.c2b^
    DO
      FOR i := 0 TO LAST(side) DO 
        WITH 
          e = side[i],
          u = Dest(Onext(e)).c,
          v = Org(e).c,
          w = Dest(e).c
        DO
          b2c[i] := LR3x3.T{
             LR3.T{u[0], v[0], w[0]},
             LR3.T{u[1], v[1], w[1]},
             LR3.T{u[2], v[2], w[2]}
          };
          c2b[i] := LR3x3.Inv(b2c[i])
        END
      END
    END
  END ComputeTriangleMatrices;
  
PROCEDURE BuildTables(e: Arc): T =
  VAR tri: T;
  BEGIN
    tri.arc := NumberArcs(e);
    tri.out := NumberVertices(tri.arc^);
    tri.side := NumberTriangles(tri.arc^);
    ComputeTriangleMatrices(tri);
    RETURN tri
  END BuildTables;
  
PROCEDURE Degree(e: Arc): CARDINAL =
  VAR t: Arc := e;
      n: CARDINAL := 0;
  BEGIN
    REPEAT
      INC(n); t := Onext(t)
    UNTIL t = e;
    RETURN n
  END Degree;

PROCEDURE Org (e: Arc): REF Site =
  BEGIN 
    RETURN Odata(e)
  END Org;
  
PROCEDURE Dest (e: Arc): REF Site =
  BEGIN
    RETURN Ddata(e)
  END Dest;

PROCEDURE Left (e: Arc): REF Triangle =
  BEGIN 
    RETURN Ldata(e)
  END Left;
  
PROCEDURE Right (e: Arc): REF Triangle =
  BEGIN
    RETURN Rdata(e)
  END Right;

PROCEDURE BuildTetrahedron(READONLY site: ARRAY [0..3] OF REF Site): Arc =
  VAR e1, e2, e: Arc;
  BEGIN
    e1 := MakeEdge();
    SetOdata(e1, site[0]);  
    SetDdata(e1, site[1]);
    e2 := MakeEdge();
    SetOdata(e2, site[1]);  
    SetDdata(e2, site[2]);
    Splice (Sym(e1), e2);
    e := Connect(e2, e1);   
    IF PositiveSide(Org(e1).c, Org(e2).c, Dest(e2).c, site[3].c) THEN 
      Pyramid(e, site[3]); 
      Consist(e);
      Consist(Rnext(e));
      Consist(Rprev(e));
    ELSE
      Pyramid(Sym(e), site[3]);
      Consist(e);
      Consist(Lnext(e));
      Consist(Lprev(e));
    END;
    RETURN e
  END BuildTetrahedron;
  
PROCEDURE Consist(e : Arc) =
  VAR t := e;
  BEGIN
    <* ASSERT e # Onext(e) *>
    <* ASSERT Oprev(e) # Onext(e) *>
    <* ASSERT e # Dnext(e) *>
    <* ASSERT Dprev(e) # Dnext(e) *>
    <* ASSERT Org(e) # Dest(e) *>
    REPEAT
      WITH
        a = Org(t)^, 
        b = Dest(t)^,
        u = Dest(Onext(t))^,
        v = Dest(Oprev(t))^
      DO
        <* ASSERT PositiveSide(a.c, u.c, b.c, v.c) *>
      END;
      t := Onext(t)
    UNTIL t = e
  END Consist;

PROCEDURE Pyramid(e: Arc; x: REF Site) =
 VAR base: Arc;
     pri: REF Site;
     
  BEGIN
    base := MakeEdge();
    pri := Org(e);
    SetOdata(base, pri);
    SetDdata(base, x); 
    Splice(base, e);
    REPEAT 
      base := Connect(e, Sym(base));
      e := Oprev(base);
    UNTIL Dest(e) = pri;
  END Pyramid;

PROCEDURE Swap(e: Arc) = 
  VAR a, b : Arc;
  BEGIN
    a:= Oprev(e);
    b:= Oprev(Sym(e));
    Splice(e, a); Splice(Sym(e), b);
    Splice(e, Lnext(a)); Splice(Sym(e), Lnext(b));
    SetOdata(e, Dest(a));
    SetDdata(e, Dest(b));
  END Swap;

PROCEDURE DeleteEdge(e: Arc) = 
  BEGIN
    Splice(e, Oprev(e));
    Splice(Sym(e), Oprev(Sym(e)));
  END DeleteEdge;
  
PROCEDURE RemoveSite(b: Arc) =
  (* Dest(b) = no a ser retirado *)
  VAR pri: REF Site;
      e: Arc;
  BEGIN
    pri:= Org(b);
    REPEAT
      e := Oprev(b);
      DeleteEdge(b);
      b := Lnext(e)
    UNTIL Dest(e) = pri
  END RemoveSite;
  
PROCEDURE Retriangulate(
    e: Arc; 
    deg: INTEGER;
    VAR tr: ARRAY OF Arc;
    VAR nTr: CARDINAL;
  ) = 
  
  PROCEDURE RemoveCorner(a: Arc): Arc =
    (*
      Removes the corner "Org(a)", "Dest(a)", "Dest(Lnext(a))"
      from the face "Left(a)", if that triangle is a Delaunay
      triangle; and returns the new edge that is 
      a side of the new triangle. Otherwise returns NullArc. *)
    VAR r: Arc;
    BEGIN
      WITH 
        b = Lnext(a)
      DO
        (* Check whether "Org(a)", "Dest(a)", "Dest(b)" is a valid face; that is, *)
        (* whether all other vertices are on the correct side of plane: *) 
        r := Lnext(Lnext(b));
        <* ASSERT r # a *>
        WHILE PositiveSide(Org(r).c, Org(a).c, Dest(a).c, Dest(b).c) DO  
          r := Lnext(r);
          IF r = a THEN RETURN Connect(b, a) END
        END;
        RETURN NullArc
      END
    END RemoveCorner;
  
  BEGIN
    nTr := 0;
    WHILE deg > 3 DO
      WITH
        g = RemoveCorner(e)
      DO
        IF g # NullArc THEN
          tr[nTr] := g; INC(nTr);
          DEC(deg); e := Sym(g)
        ELSE
          e := Lnext(e)
        END
      END
    END;
    tr[nTr] := e; INC(nTr);
  END Retriangulate;   

PROCEDURE SamePoint(READONLY p, q: Coords): BOOLEAN =
  BEGIN
    RETURN p[0] = q[0] AND p[1] = q[1] AND p[2] = q[2]
  END SamePoint;
   
PROCEDURE InsertSite(READONLY s: REF Site; e: Arc): Arc =
  VAR t, res: Arc;
     pri: REF Site;
  BEGIN
    e:= Locate(s.c, e);
    IF SamePoint(s.c, Org(e).c) 
    OR SamePoint(s.c, Dest(e).c)  THEN 
      RETURN e
    END;
    WITH
      a = Org(e)^,
      b = Dest(e)^,
      c = Dest(Onext(e))^
    DO
      <* ASSERT PositiveSide(a.c, b.c, c.c, s.c) *>
    END;
    pri := Dest(e);
    Pyramid(e, s);
    LOOP 
      t := Dnext(e);
      IF PositiveSide(Org(e).c, Org(t).c, Dest(e).c, s.c) THEN 
        Swap(e);
        e := t;
      ELSIF  Org(e) = pri  THEN 
        res := Sym(Onext(e)); EXIT
      ELSE 
        e := Lprev(Onext(e))
      END
    END;
    Consist(res);
    t := res;
    REPEAT
      Consist(Sym(t));
      t := Onext(t)
    UNTIL t = res;
    RETURN res
  END InsertSite; 

PROCEDURE BuildInc(READONLY site: ARRAY OF REF Site): Arc =
  VAR e: Arc;
  BEGIN
    <* ASSERT NUMBER(site) >= 4 *>
    e := BuildTetrahedron(SUBARRAY(site, 0, 4));
    FOR i := 4 TO LAST(site) DO
      e := InsertSite(site[i], e)
    END;
    RETURN e;
  END BuildInc;

PROCEDURE Locate(READONLY p: Coords; e: Arc ): Arc =
  BEGIN
    WITH q = FourCenter(e) DO 
    LOOP 
      WITH 
        a = Org(e).c,
        b = Dest(e).c,
        c = Dest(Onext(e)).c
      DO
        IF SamePoint(p, a) THEN RETURN e
        ELSIF SamePoint(p, b) THEN RETURN e
        ELSIF SamePoint(p, c) THEN RETURN e
        ELSIF PositiveSide(q, b, a, p) THEN e := Sym(e)
        ELSIF NOT PositiveSide(q, c, a, p) THEN e := Onext(e)
        ELSIF NOT PositiveSide(q, b, c, p) THEN e := Dprev(e) 
        ELSE RETURN e 
        END
      END;
    END
   END
  END Locate;
  
PROCEDURE FourCenter(e: Arc): Coords =
  VAR c: Coords := Org(e).c;
  BEGIN
    c := LR3.Add(c, Dest(e).c);
    c := LR3.Add(c, Dest(Onext(e)).c);
    c := LR3.Add(c, Dest(Oprev(e)).c);
    RETURN LR3.Scale(0.25d0, c)
  END FourCenter;
  
PROCEDURE Connect (a, b: Arc): Arc =
  (* Connects Dest(a) to Org(b) across the face Left(a) = Left(b). *)
  VAR e: Arc := MakeEdge();
  BEGIN
    SetOdata(e, Dest(a));
    SetDdata(e, Org(b));
    Splice(e, Lnext(a));
    Splice(Sym(e), b);
    RETURN e
  END Connect;
  
<*UNUSED*>
PROCEDURE Disconnect(a: Arc) =
  BEGIN
    IF a # Onext(a) THEN Splice(a, Oprev(a)) END;
    a := Sym(a);
    IF a # Onext(a) THEN Splice(a, Oprev(a)) END;
  END Disconnect;
  
PROCEDURE PositiveSide(READONLY a, b, c, d: Coords): BOOLEAN =
  (* TRUE iff "d" is on the positive side of the plane "a b c" *)
  BEGIN
    WITH
      ba = LR3.Sub(b, a),
      ca = LR3.Sub(c, a),
      da = LR3.Sub(d, a),
      det = LR3Extras.Det(ba, ca, da)
    DO
      RETURN det > 0.0d0
    END
  END PositiveSide;

PROCEDURE PerspDraw(
    READONLY arc: ARRAY OF Arc; (* Result of NumberArcs *)
    name: TEXT; 
    obs, upp: HLR3.Point;
  ) =

  VAR f := NEW(PSPlot.PSFile).open(name & ".ps") ;
      Map := HLR3.PerspMap(
        obs := obs,
        foc := HLR3.Point{c := LR4.T{1.0d0,0.0d0,0.0d0,0.0d0}},
        upp := upp
      );
      
  TYPE ProjectedPoint = RECORD x, y: LONGREAL  END;
  
  PROCEDURE Project(READONLY s: Site): ProjectedPoint =
    BEGIN
      WITH 
        q = HLR3.Point{c := LR4.T{1.0d0, s.c[0], s.c[1], s.c[2]}},
        g = HLR3.MapPoint(q, Map)
      DO 
        RETURN ProjectedPoint{g.c[1]/g.c[0], g.c[2]/g.c[0]}
      END
    END Project;

  PROCEDURE PlotDelaunayEdge(e:Arc) =
    BEGIN 
      WITH 
        op = Project(Org(e)^), 
        dp = Project(Dest(e)^) 
      DO
        f.segment(op.x, op.y, dp.x, dp.y);
      END;
    END PlotDelaunayEdge;
    
  BEGIN
    f.beginPage();
    f.beginDrawing();
    f.setScale(PSPlot.Axis.X, -1.5D0, 1.5D0);
    f.setScale(PSPlot.Axis.Y, -1.5D0, 1.5D0);
    f.caption("Delaunay triangulation on sphere");
    f.setLineWidth(0.1);
    f.setLineColor(PSPlot.Black);
    
    FOR i := 0 TO LAST(arc) BY 2 DO
      IF arc[i] # NullArc THEN PlotDelaunayEdge(arc[i]) END
    END;
    
    f.frame();
    f.endDrawing();
    f.endPage();
    f.close();
  END PerspDraw;

PROCEDURE Check(
    READONLY arc: ARRAY OF Arc;
  ) =

  PROCEDURE CheckDelaunayEdge(e:Arc) =
    BEGIN
      WITH 
        a = Org(e)^, 
        b = Dest(e)^,
        u = Dest(Onext(e))^,
        v = Dest(Oprev(e))^
      DO
        <* ASSERT PositiveSide(a.c, u.c, b.c, v.c) *>
      END;
    END CheckDelaunayEdge;

  BEGIN
    FOR i := 0 TO LAST(arc) BY 2 DO
      IF arc[i] # NullArc THEN CheckDelaunayEdge(arc[i]) END
    END;
  END Check;

PROCEDURE OutputX3D(
    READONLY side: ARRAY OF Arc; (* Result of NumberTriangles *)
    name: TEXT;
  ) =

  <* FATAL Wr.Failure, Thread.Alerted, OSError.E *>

  VAR f: Wr.T := FileWr.Open(name & ".poly") ;

  PROCEDURE OutputDelaunayTriangle(e:Arc) =
    BEGIN
      WITH 
        a = Org(e)^, 
        b = Dest(e)^,
        c = Dest(Onext(e))^
      DO
        Wr.PutText(f, "\n");
        Wr.PutText(f, "0.5 0.5 0 # color\n");
        OutputX3DCoords(f, a.c);
        Wr.PutText(f, "\n");
        OutputX3DCoords(f, b.c);
        Wr.PutText(f, "\n");
        OutputX3DCoords(f, c.c);
        Wr.PutText(f, "\n");
      END;
    END OutputDelaunayTriangle;

  BEGIN
    FOR i := 0 TO LAST(side) DO
      IF side[i] # NullArc THEN OutputDelaunayTriangle(side[i]) END;
    END;
    Wr.Close(f);
  END OutputX3D;

PROCEDURE OutputX3DCoords(f: Wr.T; READONLY p: Coords) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    FOR i := 0 TO 2 DO
      IF i # 0 THEN Wr.PutText(f, "  ") END;
      Wr.PutText(f, Fmt.Pad(Fmt.LongReal(p[i] * 1000.0d0, Fmt.Style.Fix, 2), 8));
    END;
  END OutputX3DCoords;

PROCEDURE Print(READONLY t: T; name: TEXT) =

  <* FATAL Wr.Failure, Thread.Alerted, OSError.E *>

  VAR f: Wr.T := FileWr.Open(name & ".txt") ;
      
  PROCEDURE PrintDelaunayEdge(e:Arc) =
    BEGIN
      WITH 
        op = Org(e)^, 
        dp = Dest(e)^ 
      DO
        Wr.PutText(f, Fmt.Pad(ToText(e), 5));
        Wr.PutText(f, " ");
        Wr.PutText(f,  Fmt.Pad( ToText(Onext(e)), 5));
        Wr.PutText(f, " ");
        Wr.PutText(f,  Fmt.Pad(ToText(Oprev(e)), 5));
        Wr.PutText(f, "  ");
        Wr.PutText(f, Fmt.Pad(Fmt.Int(op.num), 4));
        Wr.PutText(f, " - ");
        Wr.PutText(f, Fmt.Pad(Fmt.Int(dp.num), 4));
        Wr.PutText(f, "   ");
        WITH tr = NARROW(Ldata(e), REF Triangle) DO 
          IF tr # NIL THEN 
            Wr.PutText(f, Fmt.Pad(Fmt.Int(tr.num), 4))
          ELSE
            Wr.PutText(f, Fmt.Pad("NIL", 4))
          END;
        END;
        Wr.PutText(f, " - ");
        WITH tr = NARROW(Rdata(e), REF Triangle) DO 
          IF tr # NIL THEN 
            Wr.PutText(f, Fmt.Pad(Fmt.Int(tr.num), 4))
          ELSE
            Wr.PutText(f, Fmt.Pad("NIL", 4))
          END;
        END;
        Wr.PutText(f, "\n");
      END;
    END PrintDelaunayEdge;

  BEGIN
    Wr.PutText(f, "# " & Fmt.Int(NUMBER(t.out^)) & " Vertices:\n");
    FOR i := 0 TO LAST(t.out^) DO 
      IF t.out[i] # NullArc THEN
        WITH v = Org(t.out[i]) DO 
          Wr.PutText(f, Fmt.Pad(Fmt.Int(v.num), 4));
          Wr.PutText(f, "  ");
          PrintCoords(f, v.c);
          Wr.PutText(f, "\n");
        END
      END;
    END;
  
    Wr.PutText(f, "# " & Fmt.Int(NUMBER(t.arc^) DIV 2) & " Edges:\n");
    FOR i := 0 TO LAST(t.arc^) BY 2 DO
      IF t.arc[i] # NullArc THEN 
        PrintDelaunayEdge(t.arc[i])
      END;
    END;
    
    (* Should print triangles, too *)
    
    Wr.Close(f);
  END Print;
  
PROCEDURE Write(wr: Wr.T; READONLY t: T) =

  <* FATAL Wr.Failure, Thread.Alerted *>
  
  PROCEDURE WriteCoords(READONLY p: Coords) =
    BEGIN
      FOR i := 0 TO 2 DO
        IF i # 0 THEN FPut.Char(wr, ' ') END;
        Wr.PutText(wr, Fmt.LongReal(p[i]));
      END;
    END WriteCoords;

  PROCEDURE WriteArc(a: Arc) =
    BEGIN
      Wr.PutText(wr, ToText(a));
    END WriteArc;
    
  PROCEDURE WriteOrg(a: Arc) =
    BEGIN
      FPut.Char(wr, 'v'); FPut.Int(wr, Org(a).num);
    END WriteOrg;
    
  PROCEDURE WriteLeft(a: Arc) =
    BEGIN
      FPut.Char(wr, 'f'); FPut.Int(wr, Left(a).num);
    END WriteLeft;
    
  PROCEDURE WriteEdgeData(e:Arc) =
    BEGIN
      WriteArc(e); FPut.Char(wr, ' ');
      FPut.Char(wr, ' ');

      WriteArc(Onext(e)); FPut.Char(wr, ' ');
      WriteArc(Onext(Sym(e))); FPut.Char(wr, ' ');
      FPut.Char(wr, ' ');

      WriteOrg(e); FPut.Char(wr, ' ');
      WriteOrg(Sym(e)); FPut.Char(wr, ' ');
      WriteLeft(e); FPut.Char(wr, ' ');
      WriteLeft(Sym(e)); 
      FPut.EOL(wr);
    END WriteEdgeData;

  BEGIN
    FileFmt.WriteHeader(wr, "triangulation", TriangFileVersion);
    WITH
      NV = NUMBER(t.out^),
      NE = NUMBER(t.arc^) DIV 2,
      NF = NUMBER(t.side^)
    DO
      NPut.Int(wr, "vertices", NV);  FPut.EOL(wr);
      NPut.Int(wr, "edges", NE);  FPut.EOL(wr);
      NPut.Int(wr, "faces", NF); FPut.EOL(wr);
      
      Wr.PutText(wr, "vertices:"); FPut.EOL(wr);
      FOR i := 0 TO NV-1 DO 
        <* ASSERT t.out[i] # NullArc *>
        WITH a = t.out[i], v = Org(a) DO 
          FPut.Int(wr, i); FPut.Char(wr, ' ');
          WriteArc(a); FPut.Char(wr, ' ');
          WriteCoords(v.c); FPut.EOL(wr);
        END;
      END;

      Wr.PutText(wr, "edges:"); FPut.EOL(wr);
      FOR i := 0 TO LAST(t.arc^) BY 2 DO
        <* ASSERT t.arc[i] # NullArc *>
        FPut.Int(wr, i DIV 2); FPut.Char(wr, ' ');
        WriteEdgeData(t.arc[i])
      END;

      Wr.PutText(wr, "faces:"); FPut.EOL(wr);
      FOR i := 0 TO LAST(t.side^) DO
        <* ASSERT t.side[i] # NullArc *>
        FPut.Int(wr, i); FPut.Char(wr, ' ');
        WriteArc(t.side[i]); FPut.EOL(wr);
      END;

    END;
    FileFmt.WriteFooter(wr, "triangulation");
    Wr.Flush(wr);
  END Write;
  
PROCEDURE Read(rd: Rd.T): T =
  VAR NV, NE, NF: CARDINAL;
  VAR t: T;
  BEGIN
    FileFmt.ReadHeader(rd, "triangulation", TriangFileVersion);
    NV := NGet.Int(rd, "vertices"); FGet.Skip(rd);
    NE := NGet.Int(rd, "edges"); FGet.Skip(rd);
    NF := NGet.Int(rd, "faces"); FGet.Skip(rd);
    t.out := NEW(REF ARRAY OF Arc, NV);
    t.arc := NEW(REF ARRAY OF Arc, 2*NE);
    t.side := NEW(REF ARRAY OF Arc, NF);

    WITH
      site = NEW(REF ARRAY OF REF Site, NV)^,
      edge = NEW(REF ARRAY OF Arc, NE)^,
      face = NEW(REF ARRAY OF REF Triangle, NF)^
    DO
      
      (* Create all vertices, edges, and faces: *)
      FOR i := 0 TO NV-1 DO 
        WITH v = site[i]  DO v := NEW(REF Site); v.num := i END
      END;
      FOR i := 0 TO NE-1 DO
        WITH a = MakeEdge() DO 
          SetEdgeNum(a, i);
          edge[i] := a;
        END
      END;
      FOR i := 0 TO NF-1 DO 
        WITH f = face[i]  DO f := NEW(REF Triangle); f.num := i END
      END;

      PROCEDURE GetArc(): Arc =
        VAR qnum, rnum: CARDINAL; a: Arc;
        BEGIN
          FGet.Skip(rd);
          qnum := FGet.Int(rd);
          FGet.Match(rd, ":");
          rnum := FGet.Int(rd);
          <* ASSERT rnum >= 0 AND rnum <= 3 *>
          a := edge[qnum];
          WHILE a.bits # rnum DO a := Rot(a) END;
          RETURN a
        END GetArc;

      PROCEDURE GetVertex(): REF Site =
        VAR vnum: CARDINAL;
        BEGIN
          FGet.Skip(rd);
          FGet.Match(rd, "v");
          vnum := FGet.Int(rd);
          RETURN site[vnum]
        END GetVertex;

      PROCEDURE GetFace(): REF Triangle =
        VAR fnum: CARDINAL;
        BEGIN
          FGet.Skip(rd);
          FGet.Match(rd, "f");
          fnum := FGet.Int(rd);
          RETURN face[fnum]
        END GetFace;
        
      PROCEDURE ReadCoords(VAR c: LR3.T) =
        BEGIN
          FOR i := 0 TO 2 DO
            FGet.Skip(rd); c[i] := FGet.LongReal(rd);
          END;
        END ReadCoords;

      BEGIN

        (* Read vertex coordinates: *)
        FGet.Match(rd, "vertices:"); FGet.EOL(rd);
        FOR i := 0 TO NV-1 DO 
          WITH j = FGet.Int(rd) DO <* ASSERT j = i *> END;
          FGet.Skip(rd);
          t.out[i] := GetArc();
          FGet.Skip(rd);
          WITH v = site[i] DO
            ReadCoords(v.c); FGet.EOL(rd);
          END;
        END;

        (* Read edges pointers: *)
        FGet.Match(rd, "edges:"); FGet.EOL(rd);
        FOR i := 0 TO NE-1 DO
          WITH j = FGet.Int(rd) DO <* ASSERT j = i *> END;
          FGet.Skip(rd);
          WITH 
            a = GetArc(),
            b = GetArc(), 
            c = GetArc(),
            o = GetVertex(),
            d = GetVertex(),
            l = GetFace(),
            r = GetFace()
          DO
            t.arc[2*i] := a; 
            t.arc[2*i+1] := Sym(a);
            Splice(a, Oprev(b));
            IF b # Sym(a) THEN Splice(Sym(a), Oprev(c)) END;
            SetOdata(a, o);
            SetDdata(a, d);
            SetLdata(a, l);
            SetRdata(a, r);
          END;
          FGet.EOL(rd);
        END;
        
        (* Read face data: *)
        FGet.Match(rd, "faces:"); FGet.EOL(rd);
        FOR i := 0 TO NF-1 DO 
          WITH j = FGet.Int(rd) DO <* ASSERT j = i *> END;
          FGet.Skip(rd);
          t.side[i] := GetArc();
          FGet.EOL(rd);
        END;
      END
    END;
    FileFmt.ReadFooter(rd, "triangulation");
    ComputeTriangleMatrices(t);
    RETURN t
  END Read;

PROCEDURE PrintCoords(
    f: Wr.T; 
    READONLY p: Coords; 
    lp: TEXT := "( ";
    rp: TEXT := " )";
    cm: TEXT := ", ";
  ) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    FOR i := 0 TO 2 DO
      IF i = 0 THEN Wr.PutText(f, lp) ELSE Wr.PutText(f, cm) END;
      Wr.PutText(f, Fmt.Pad(Fmt.LongReal(p[i], Fmt.Style.Fix, 5), 8));
    END;
    Wr.PutText(f, rp)
  END PrintCoords;

BEGIN
END SPTriang.