MODULE TrArScene;

IMPORT H3, R4, Wr, Fmt, Thread, PSPlot;

(* DEBUG: *)
FROM Stdio IMPORT stderr;
CONST Debug = FALSE;

VAR 
  mu: MUTEX := NEW(MUTEX);  (* Protects the following global counts: *)
  vertexCount:CARDINAL := 0;
  edgeCount:CARDINAL := 0;
  triangleCount:CARDINAL := 0;

PROCEDURE MakeVertex(
    org: REFANY; 
    READONLY wcs: H3.Point;
  ): REF Vertex =
  BEGIN
    LOCK mu DO
      WITH
        v = NEW(REF Vertex,
          num := vertexCount,
          org := org,
          mark := FALSE,
          wcs := wcs,
          val := 0.0d0
        )
      DO
        INC(vertexCount);
        RETURN v
      END
    END
  END MakeVertex;
  
PROCEDURE MakeEdge(
    org: REFANY; 
    a, b: REF Vertex;
    line: REF Line := NIL;
  ): REF Edge =
  BEGIN
    IF line = NIL THEN
      line := NEW(REF Line, ln := H3.LineFromTwoPoints(a.wcs, b.wcs));
    END;
    LOCK mu DO
      WITH
        e = NEW(REF Edge,
          num := edgeCount,
          org := org,
          mark := FALSE,
          vert := VertexPair{a, b},
          line := line
        )
      DO
        e.part := EdgePair{e, e};
        INC(edgeCount);
        RETURN e
      END
    END
  END MakeEdge;
  
PROCEDURE MakeTriangle(
    org: REFANY; 
    a: REF Edge; ad: EDir;
    b: REF Edge; bd: EDir;
    c: REF Edge; cd: EDir;
    plane: REF Plane := NIL;
    lights: LightSet := LightSet{TRUE, ..};
    rootv: [0..2] := 0;
  ): REF Triangle =
  BEGIN
    <* ASSERT a.vert[ad] = c.vert[1-cd] *>
    <* ASSERT b.vert[bd] = a.vert[1-ad] *>
    <* ASSERT c.vert[cd] = b.vert[1-bd] *>
    IF plane = NIL THEN
      WITH 
        av = a.vert[ad],
        bv = b.vert[bd],
        cv = c.vert[cd]
      DO
        plane := NEW(REF Plane, pl := H3.PlaneFromThreePoints(av.wcs, bv.wcs, cv.wcs))
      END;
    END;
    LOCK mu DO
      WITH
        ia = rootv,
        ib = (ia + 1) MOD 3,
        ic = (ib + 1) MOD 3,
        t = NEW(REF Triangle,
          num := triangleCount,
          org := org,
          mark := FALSE,
          edge := EdgeTriplet{a, b, c},
          edir := EDirTriplet{ad, bd, cd},
          plane := plane,
          lights := lights
        )
      DO
        INC(triangleCount);
        RETURN t
      END
    END
  END MakeTriangle;

PROCEDURE GatherEdges(tls: TriangleList): EdgeList =

  VAR els: EdgeList := NIL;
  
  PROCEDURE GatherEdge(e: REF Edge) =
    BEGIN
      IF NOT e.mark THEN
        els := NEW(EdgeList, e := e, rest := els);
        e.mark := TRUE;
      END
    END GatherEdge;
  
  BEGIN
    (* Clear all edge marks *)
    VAR tp: TriangleList := tls;
    BEGIN
      WHILE tp # NIL DO 
        WITH t = tp.t^ DO
          t.edge[0].mark := FALSE;
          t.edge[1].mark := FALSE;
          t.edge[2].mark := FALSE;
        END;
        tp := tp.rest
      END
    END;
    
    (* Collect edges: *)
    VAR tp: TriangleList := tls;
    BEGIN
      WHILE tp # NIL DO
        WITH t = tp.t^ DO
          GatherEdge(t.edge[0]);
          GatherEdge(t.edge[1]);
          GatherEdge(t.edge[2]);
        END;
        tp := tp.rest
      END
    END;
    RETURN els
  END GatherEdges;
    
PROCEDURE GatherVertices(els: EdgeList): VertexList =

  VAR vls: VertexList := NIL;
  
  PROCEDURE GatherVertex(v: REF Vertex) =
    BEGIN
      IF NOT v.mark THEN
        vls := NEW(VertexList, v := v, rest := vls);
        v.mark := TRUE;
      END
    END GatherVertex;
  
  BEGIN
    (* Clear all vertex marks *)
    VAR ep: EdgeList := els;
    BEGIN
      WHILE ep # NIL DO 
        WITH e = ep.e^ DO
          e.vert[0].mark := FALSE;
          e.vert[1].mark := FALSE;
        END;
        ep := ep.rest
      END
    END;
    
    (* Collect vertices: *)
    VAR ep: EdgeList := els;
    BEGIN
      WHILE ep # NIL DO
        WITH e = ep.e^ DO
          GatherVertex(e.vert[0]);
          GatherVertex(e.vert[1]);
        END;
        ep := ep.rest
      END
    END;
    RETURN vls
  END GatherVertices;

PROCEDURE CloseScene(VAR sc: T) =
  
  PROCEDURE UnmarkEdge(VAR e: Edge) = 
    BEGIN
      e.mark := FALSE; 
      e.vert[0].mark := FALSE;
      e.vert[1].mark := FALSE;
    END UnmarkEdge;
  
  PROCEDURE UnmarkTriangle(VAR t: Triangle) = 
    BEGIN
      t.mark := FALSE; 
      UnmarkEdge(t.edge[0]^);
      UnmarkEdge(t.edge[1]^);
      UnmarkEdge(t.edge[2]^);
    END UnmarkTriangle;
  
  VAR ns: T;
  
  PROCEDURE CloseVertex(v: REF Vertex) =
    BEGIN
      IF NOT v.mark THEN
        WITH l = NEW(VertexList) DO
          l.v := v; l.rest := ns.v; ns.v := l
        END;
        v.mark := TRUE
      END
    END CloseVertex;
    
  PROCEDURE CloseEdge(e: REF Edge) =
    BEGIN
      IF NOT e.mark THEN
        WITH l = NEW(EdgeList) DO
          l.e := e; l.rest := ns.e; ns.e := l
        END;
        CloseVertex(e.vert[0]);
        CloseVertex(e.vert[1]);
        e.mark := TRUE;
      END
    END CloseEdge;
  
  BEGIN
    (* Clear all marks *)
    VAR t: VertexList := sc.v;
    BEGIN
      WHILE t # NIL DO t.v.mark := FALSE; t := t.rest END
    END;
    VAR t: EdgeList := sc.e;
    BEGIN
      WHILE t # NIL DO UnmarkEdge(t.e^); t := t.rest END
    END;
    VAR t: TriangleList := sc.t;
    BEGIN
      WHILE t # NIL DO UnmarkTriangle(t.t^); t := t.rest
      END
    END;
    
    (* Collect elements (once each): *)
    ns := T{NIL, NIL, NIL};
    VAR t: VertexList := sc.v; r: VertexList;
    BEGIN
      WHILE t # NIL DO 
        r := t.rest;
        IF NOT t.v.mark THEN
          t.rest := ns.v; ns.v := t; 
          t.v.mark := TRUE;
        END;
        t := r
      END
    END;
    VAR t: EdgeList := sc.e; r: EdgeList;
    BEGIN
      WHILE t # NIL DO 
        r := t.rest;
        IF NOT t.e.mark THEN
          t.rest := ns.e; ns.e := t; 
          t.e.mark := TRUE;
          CloseVertex(t.e.vert[0]);
          CloseVertex(t.e.vert[1]);
        END;
        t := r
      END
    END;
    VAR t: TriangleList := sc.t; r: TriangleList;
    BEGIN
      WHILE t # NIL DO 
        r := t.rest;
        IF NOT t.t.mark THEN
          t.rest := ns.t; ns.t := t; 
          t.t.mark := TRUE;
          CloseEdge(t.t.edge[0]);
          CloseEdge(t.t.edge[1]);
          CloseEdge(t.t.edge[2]);
        END;
        t := r
      END
    END;
    sc := ns
  END CloseScene;

PROCEDURE ClipSceneToRegion(
    VAR in, bd, ot: T;
    READONLY plane: Planes;
    READONLY side: Sides;  
    inF, bdF, otF: BOOLEAN;
    debug: BOOLEAN := FALSE
  ) =
  VAR inP, bdP, otP: BOOLEAN;
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    (* Cut boundary and interior succesively by each plane: *)
    FOR i := 0 TO LAST(plane) DO
      IF debug THEN
        Wr.PutText(stderr, "\n");
        Wr.PutText(stderr, "ClipScene: clipping against plane " & Fmt.Int(side[i]*i) & "\n");
        Wr.PutText(stderr, "  in = "); PrintCounts(stderr, in); Wr.PutText(stderr, "\n");
        Wr.PutText(stderr, "  bd = "); PrintCounts(stderr, bd); Wr.PutText(stderr, "\n");
        Wr.PutText(stderr, "  ot = "); PrintCounts(stderr, ot); Wr.PutText(stderr, "\n");
      END;
      
      (* Let "R[i]" be the polyhedron defined by "R0 \inter reg[0..i]" *)
      (* Decide which parts of "R[i]" we want: *)
      inP := inF OR ((otF OR bdF) AND i < LAST(plane));
      bdP := bdF OR (otF AND i < LAST(plane));
      otP := otF;
      
      (* Clear "val" and "cut" fields: *)
      IF bdP OR otP THEN
        ClearVertexVals(bd.v);
        ClearLineCuts(bd.e);
        ClearPlaneCuts(bd.t)
      END;
      IF inP OR bdP OR otP THEN 
        ClearVertexVals(in.v);
        ClearLineCuts(in.e);
        ClearPlaneCuts(in.t)
      END;
      
      IF bdP OR otP THEN 
        IF debug THEN Wr.PutText(stderr, "  cutting \"bd\"...\n") END;
        EvalVertices(bd.v, plane[i], side[i]);
        AuxCutScene(bd, otP, bdP, bdP, ot, bd, bd);
        IF debug THEN
          Wr.PutText(stderr, "  result of cutting \"bd\":\n");
          Wr.PutText(stderr, "  in = "); PrintCounts(stderr, in); Wr.PutText(stderr, "\n");
          Wr.PutText(stderr, "  bd = "); PrintCounts(stderr, bd); Wr.PutText(stderr, "\n");
          Wr.PutText(stderr, "  ot = "); PrintCounts(stderr, ot); Wr.PutText(stderr, "\n");
        END;
      END;
      IF inP OR bdP OR otP THEN 
        IF debug THEN Wr.PutText(stderr, "  cutting \"in\"...\n") END;
        EvalVertices(in.v, plane[i], side[i]); 
        AuxCutScene(in, otP, bdP, inP, ot, bd, in);
        IF debug THEN
          Wr.PutText(stderr, "  result of cutting \"in\":\n");
          Wr.PutText(stderr, "  in = "); PrintCounts(stderr, in); Wr.PutText(stderr, "\n");
          Wr.PutText(stderr, "  bd = "); PrintCounts(stderr, bd); Wr.PutText(stderr, "\n");
          Wr.PutText(stderr, "  ot = "); PrintCounts(stderr, ot); Wr.PutText(stderr, "\n");
        END;
      END;
    END;
    IF debug THEN
      Wr.PutText(stderr, "\n");
      Wr.PutText(stderr, "Exit ClipScene:\n");
      Wr.PutText(stderr, "  in = "); PrintCounts(stderr, in); Wr.PutText(stderr, "\n");
      Wr.PutText(stderr, "  bd = "); PrintCounts(stderr, bd); Wr.PutText(stderr, "\n");
      Wr.PutText(stderr, "  ot = "); PrintCounts(stderr, ot); Wr.PutText(stderr, "\n");
      Wr.PutText(stderr, "\n");
    END;
  END ClipSceneToRegion;
  
PROCEDURE CutSceneByPlane(
    VAR sc: T;                 (* Scene to cut *)
    READONLY plane: H3.Plane;  (* Cutting plane *)
    side: H3.Sign := +1;       (* Multiplies "plane" *)
    READONLY null: Vertices;   (* Vertices known to be on "plane" *)
    negF, zerF, posF: BOOLEAN; (* Which parts to return *)
    VAR neg, zer, pos: T;      (* (I/O): Negative, zero, and positive parts. *)
  ) =
  BEGIN
    CloseScene(sc);
    ClearVertexVals(sc.v);
    ClearLineCuts(sc.e);
    ClearPlaneCuts(sc.t);
    EvalVertices(sc.v, plane, side, null);
    AuxCutScene(sc, negF, zerF, posF, neg, zer, pos);
  END CutSceneByPlane;

PROCEDURE ClearVertexVals(vp: VertexList) =    
  (* Clear out the marks in all supporting lines. *)
  BEGIN
    WHILE vp # NIL DO
      WITH v = vp.v^ DO v.mark := FALSE; v.val := 0.0d0 END;
      vp := vp.rest
    END
  END ClearVertexVals;

PROCEDURE ClearLineCuts(ep: EdgeList) =    
  (* Clear out the marks in all supporting lines. *)
  BEGIN
    WHILE ep # NIL DO
      WITH ln = ep.e.line^ DO ln.mark := FALSE; ln.cut := NIL END;
      ep := ep.rest
    END
  END ClearLineCuts;

PROCEDURE ClearPlaneCuts(tp: TriangleList) =
  (* Clear out the marks and cut lines in all supporting planes: *)
  BEGIN
    WHILE tp # NIL DO
      WITH pl = tp.t.plane DO pl.mark := FALSE; pl.cut := NIL END;
      tp := tp.rest
    END
  END ClearPlaneCuts;

PROCEDURE EvalVertices(
    vls: VertexList;
    READONLY pl: H3.Plane;     (* Cutting plane *)
    side: H3.Sign;             (* Multiplies "pl" *)
    READONLY null: Vertices := Vertices{}; 
  ) =
  (*
    Sets the "val" field of all vertices to their value relative to "pl",
    times the specified sign.  
  
    The "null" vertices are assumed to lie on the given plane,
    and their "val" fields are explicitly set to zero..
    The "val"s of all other vertices are evaluated numerically.
    Beware of rounding errors...
  *)
  BEGIN
    <* ASSERT side # 0 *>
    WITH
      plf = pl.f,
      pl0 = FLOAT(plf[0], LONGREAL),
      pl1 = FLOAT(plf[1], LONGREAL),
      pl2 = FLOAT(plf[2], LONGREAL),
      pl3 = FLOAT(plf[3], LONGREAL)
    DO
      WHILE vls # NIL DO
        WITH v = vls.v, vc = v.wcs.c DO
          IF NOT v.mark THEN
            WITH
              t0 = pl0 * FLOAT(vc[0], LONGREAL),
              t1 = pl1 * FLOAT(vc[1], LONGREAL),
              t2 = pl2 * FLOAT(vc[2], LONGREAL),
              t3 = pl3 * FLOAT(vc[3], LONGREAL),

              dot = t0 + t1 + t2 + t3,
              err = 2.0d-7 * (ABS(t0) + ABS(t1) + ABS(t2) + ABS(t3))
            DO
              IF ABS(dot) <= err THEN
                v.val := 0.0d0 
              ELSIF side = +1 THEN 
                v.val := dot
              ELSE
                v.val := -dot
              END
            END
          END;
        END;
        vls := vls.rest
      END
    END;
    FOR i := 0 TO LAST(null) DO
      WITH v = null[i] DO 
        IF v # NIL THEN v.mark := TRUE; v.val := 0.0d0 END
      END
    END
  END EvalVertices;
  
PROCEDURE AuxCutScene(
    VAR sc: T;
    negF, zerF, posF: BOOLEAN;
    VAR neg, zer, pos: T;
  ) =
  (*
    Cuts the elements of "sc" by a plane, and prepends the negative,
    zero, and positive parts onto "neg", "zer", and "bS", respectively.
    
    The negative parts are returned only if "negF = TRUE";
    otherwise those parts are discarded, and the "neg" 
    scene is not changed.  In the same way, "zerF" and "posF"
    control the disposal of the zero and positive parts.

    This procedure assumes that the "val" fields of all vertices have
    been evaluated, relative to the desired cutting plane.  Note that
    this assumption applies to all vertices in "sc.V", to the
    endpoints of all edges in "sc.e", and to the corners of all
    triangles in "sc.t".
    
    The procedure also assumes that if "ln" is a "Line"
    record reachable through "sc" which has "ln.mark=TRUE" 
    then the "ln.cut" field is a vertex at the intersection
    of the line and the current cutting plane.  Similarly,
    if "pl" is a reachable "Plane" record with "pl.mark=TRUE",
    then "pl.cut" is the line where "pl" meets the cutting plane.
    
    Beware that the variable "sc" is cleared, and the nodes of the 
    input lists may be reused in the output lists.
    
    The variables "sc", "neg", "zer", and "pos" need not be distinct
    from each other. 
  *)
  VAR tmp: T;
  BEGIN
    (* We must clear "sc" right away, in case it is aliased with *)
    (* any of the output variables "neg", "zer", or "pos". *)
    tmp := sc; sc := T{NIL,NIL,NIL};

    (* Cut elements, in order of dimension: *)
    CutVertices  (tmp.v, negF, zerF, posF, neg.v, zer.v, pos.v);
    IF Debug THEN
      DebugSlices(stderr, "CutVertices:  ", neg, zer, pos)
    END;
    CutEdges     (tmp.e, negF, zerF, posF, zer.v, neg.e, zer.e, pos.e);
    IF Debug THEN
      DebugSlices(stderr, "CutEdges:     ", neg, zer, pos)
    END;
    CutTriangles (tmp.t, negF, zerF, posF, neg.e, zer.e, pos.e, neg.t, zer.t, pos.t);
    IF Debug THEN
      DebugSlices(stderr, "CutTriangles: ", neg, zer, pos)
    END;
  END AuxCutScene;
  
PROCEDURE DebugSlices(
    wr: Wr.T;
    msg: TEXT;
    neg, zer, pos: T
  ) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    Wr.PutText(wr, msg); 
    Wr.PutText(wr, "\n"); 
    Wr.PutText(wr, "    neg = "); PrintCounts(wr, neg); Wr.PutText(wr, "\n"); 
    Wr.PutText(wr, "    zer = "); PrintCounts(wr, zer); Wr.PutText(wr, "\n"); 
    Wr.PutText(wr, "    pos = "); PrintCounts(wr, pos); Wr.PutText(wr, "\n"); 
  END DebugSlices;

PROCEDURE CutVertices(
    VAR vin: VertexList;
    negF, zerF, posF: BOOLEAN;           (* Which parts to return *)
    VAR vneg, vzer, vpos: VertexList; (* (I/O) classified vertices *)
  ) =
  (*
    Prepends to "vneg", "vzer", and "vpos" the vertices in "vin"
    that lie on the positive, null, and negative part of 
    the plane "pl", respectively. 

    The negative parts are returned only if "negF = TRUE";
    otherwise those parts are discarded, and the "vneg" 
    list is not changed.  In the same way, "zerF" and "posF"
    control the disposal of the zero and positive parts.

    Assumes that the "val" fields of all vertices have been 
    evaluated, relative to the desired cutting plane.

    Beware that the variable "vin" is cleared, and the nodes of the 
    input list may be reused in the output lists.
    
    The variables "vin", "vneg", "vzer", and "vpos" need not be distinct
    from each other. 
  *)
  VAR tmp, rst: VertexList;
  BEGIN
    (* Beware that "vin" may be aliased with "vneg", "vzer", or "vpos". *)
    (* Copy "vin" to a local variable "tmp", and clear it out: *)
    tmp := vin; vin := NIL;
    WHILE tmp # NIL DO
      rst := tmp.rest; tmp.rest := NIL;
      WITH 
        v = tmp.v, 
        val = v.val
      DO
        IF val < 0.0d0 THEN
          IF negF THEN tmp.rest := vneg; vneg := tmp END
        ELSIF val > 0.0d0 THEN
          IF posF THEN tmp.rest := vpos; vpos := tmp END
        ELSE
          IF zerF THEN tmp.rest := vzer; vzer := tmp END
        END
      END;
      tmp := rst
    END
  END CutVertices;

PROCEDURE CutEdges(
    VAR ein: EdgeList;
    negF, zerF, posF: BOOLEAN;
    VAR vzer: VertexList;
    VAR eneg, ezer, epos: EdgeList;
  ) =
  (*
    Prepends to "eneg", "ezer", and "epos" the pieces of edges in "ein"
    that lie on the positive, null, and negative part of 
    the current cutting plane, respectively.  Also prepends to "vzer" 
    all 'new' vertices, which are the points where the 
    edges in "ein" cross the plane.

    The negative parts are returned only if "negF = TRUE";
    otherwise those parts are discarded, and the "eneg"
    list is not changed.  In the same way, "zerF" and "posF"
    control the disposal of the zero and positive parts.

    As a side effect, sets the "part" fields of all edges
    in "ein".  Any edge "e" that crosses the plane is split, and the
    two pieces stored in the "part" fields.  Any edge "e" that does
    not cross the plane gets both "part" fields set to "e" itself.
    In any case, the parts are numbered and oriented so that 
    "e.part[i].vert[i] = e.vert[i]" for "i=0..1".

    Beware that the variable "ein" is cleared, and the nodes of the 
    input list may be reused in the output lists.
    
    The variables "ein", "eneg", "ezer", and "epos" need not be distinct
    from each other. 
  *)
  VAR tmp, rst: EdgeList;
  BEGIN
    (* Beware that "ein" may be aliased with "eneg", "ezer", or "epos". *)
    (* Copy "ein" to a local variable "tmp", and clear it out: *)
    tmp := ein; ein := NIL;
    
    WHILE tmp # NIL DO
      rst := tmp.rest; tmp.rest := NIL;
      WITH 
        e = tmp.e,
        s0 = e.vert[0].val,
        s1 = e.vert[1].val
      DO
        e.part[0] := e; e.part[1] := e;
        (* Now return the two parts: *)
        IF s0 = 0.0d0 AND s1 = 0.0d0 THEN
          IF zerF THEN tmp.rest := ezer; ezer := tmp END
        ELSIF s0 <= 0.0d0 AND s1 <= 0.0d0 THEN
          IF negF THEN tmp.rest := eneg; eneg := tmp END
        ELSIF s0 >= 0.0d0 AND s1 >= 0.0d0 THEN
          IF posF THEN tmp.rest := epos; epos := tmp END
        ELSE
          (* Split edge: *)
          VAR ne, pe: REF Edge;
              zv: REF Vertex;
          BEGIN
            SplitEdge(e, ne, zv, pe, negF, posF);
            BEGIN
              IF negF THEN
                tmp.e := ne; tmp.rest := eneg; eneg := tmp; tmp := NIL
              END;
              IF zerF THEN
                WITH vnew = NEW(VertexList) DO
                  vnew.v := zv; vnew.rest := vzer; vzer := vnew
                END
              END;
              IF posF THEN
                IF tmp = NIL THEN tmp := NEW(EdgeList) END;
                tmp.e := pe; tmp.rest := epos; epos := tmp; tmp := NIL
              END;
            END
          END
        END
      END;
      tmp := rst
    END;
  END CutEdges;  
    
PROCEDURE SplitEdge(
    e: REF Edge;
    VAR ne: REF Edge;
    VAR zv: REF Vertex;
    VAR pe: REF Edge;
    negF, posF: BOOLEAN;
  ) =
  (*
    Breaks edge "e" in two new edges by the current plane "pl".
    Stores the two half-edges in "e.part[0..1]", 
    and returns them in "ne" (negative) and "pe" (positive).  Also returns
    the cut vertex "zv". 
    
    The negative part is created only if "negF=TRUE", otherwise "ne"
    is set to "NIL.  In the same way, "posF" determines the 
    creation of the the positive part.

    Assumes that the endpoint values rel. to "pl" have been computed
    and have strictly opposite signs. 
  *)
  VAR make0, make1: BOOLEAN;
  BEGIN
    WITH
      v0 = e.vert[0], s0 = v0.val,
      v1 = e.vert[1], s1 = v1.val,
      ln = e.line^
    DO
      <* ASSERT s0*s1 < 0.0d0 *>
      IF NOT ln.mark THEN
        (* Compute cut vertex for this line: *)
        WITH
          w0 = ABS(s0),
          w1 = ABS(s1),
          t = 1.0d0/(w0 + w1),
          cutpt = H3.Point{c := R4.Mix(w1*t, v0.wcs.c, w0*t, v1.wcs.c)},
          cut = MakeVertex(e.org, cutpt)
        DO
          cut.val := 0.0d0;
          ln.cut := cut;
        END;
        ln.mark := TRUE
      END;
      zv := ln.cut; <* ASSERT zv # NIL *>
      IF s0 < 0.0d0 THEN 
        make0 := negF; make1 := posF
      ELSE
        make0 := posF; make1 := negF
      END;
      IF make0 THEN 
        e.part[0] := MakeEdge(e.org, v0, zv, e.line)
      ELSE
        e.part[0] := NIL
      END;
      IF make1 THEN 
        e.part[1] := MakeEdge(e.org, zv, v1, e.line)
      ELSE
        e.part[1] := NIL
      END;
      IF s0 < 0.0d0 THEN 
        ne := e.part[0]; pe := e.part[1]
      ELSE
        ne := e.part[1]; pe := e.part[0]
      END
    END
  END SplitEdge;

PROCEDURE CutTriangles(
    VAR tin: TriangleList;
    negF, zerF, posF: BOOLEAN;
    VAR eneg, ezer, epos: EdgeList;
    VAR tneg, tzer, tpos: TriangleList;
  ) = 
  (*
    Prepends to "tneg", "tzer", and "tpos" the pieces of triangles
    in "tin" that lie on the positive, null, and negative part of 
    the current cutting plane, respectively.
    
    Assumes that the "val" fields of the vertices of all triangles in
    "tin" have been evaluated relative to the current plane, and that
    the edges have already been cut and their "part" fields set
    appropriately.

    Also prepends to "eneg", "ezer" and "epos" all 'new' edges,
    which are the segments where the triangles cross the plane,
    and any extra edges needed to triangulate non-triangular
    pieces.

    The negative parts are returned only if "negF = TRUE";
    otherwise those parts are discarded, and the "tneg"
    list is not changed.  In the same way, "zerF" and "posF"
    control the disposal of the zero and positive parts.

    Beware that the variable "tin" is cleared, and the nodes of the 
    input list may be reused in the output lists.
    
    The variables "tin", "tneg", "tzer", and "tpos" need not be distinct
    from each other. 
  *)
  VAR tmp, rst: TriangleList;
  BEGIN
    (* Beware that "tin" may be aliased with "tneg", "tzer", or "tpos". *)
    (* Copy "tin" to a local variable "tmp", and clear it out: *)
    tmp := tin; tin := NIL;
    
    WHILE tmp # NIL DO
      rst := tmp.rest; tmp.rest := NIL;
      WITH 
        t = tmp.t,
        e0 = t.edge[0], d0 = t.edir[0], s0 = e0.vert[d0].val,
        e1 = t.edge[1], d1 = t.edir[1], s1 = e1.vert[d1].val,
        e2 = t.edge[2], d2 = t.edir[2], s2 = e2.vert[d2].val
      DO
        IF s0 = 0.0d0 AND s1 = 0.0d0 AND s2 = 0.0d0 THEN
          IF zerF THEN tmp.rest := tzer; tzer := tmp END
        ELSIF s0 <= 0.0d0 AND s1 <= 0.0d0 AND s2 <= 0.0d0 THEN
          IF negF THEN tmp.rest := tneg; tneg := tmp END
        ELSIF s0 >= 0.0d0 AND s1 >= 0.0d0 AND s2 >= 0.0d0 THEN
          IF posF THEN tmp.rest := tpos; tpos := tmp END
        ELSE
          (* Split triangle: *)
          VAR nt0, nt1, pt0, pt1: REF Triangle;
              ne, ze, pe: REF Edge;
              iv: [0..2];
          BEGIN
            IF s0*s1 <= 0.0d0 AND s0*s2 < 0.0d0 THEN
              iv := 0;
            ELSIF s1*s2 <= 0.0d0 AND s1*s0 < 0.0d0 THEN
              iv := 1
            ELSIF s2*s0 <= 0.0d0 AND s2*s1 < 0.0d0 THEN
              iv := 2
            ELSE
              <* ASSERT FALSE *>
            END;
            SplitTriangle(t, iv, nt0, ne, nt1, ze, pt0, pe, pt1, negF, posF);
            IF negF THEN
              tmp.t := nt0; tmp.rest := tneg; tneg := tmp; tmp := NIL;
              IF nt1 # NIL THEN
                <* ASSERT ne # NIL *>
                WITH enew = NEW(EdgeList) DO
                  enew.e := ne; enew.rest := eneg; eneg := enew
                END;
                tmp := NEW(TriangleList);
                tmp.t := nt1; tmp.rest := tneg; tneg := tmp; tmp := NIL
              END;
            END;
            IF zerF THEN
              <* ASSERT ze # NIL *>
              WITH enew = NEW(EdgeList) DO
                enew.e := ze; enew.rest := ezer; ezer := enew
              END
            END;
            IF posF THEN
              IF tmp = NIL THEN tmp := NEW(TriangleList) END;
              tmp.t := pt0; tmp.rest := tpos; tpos := tmp; tmp := NIL;
              IF pt1 # NIL THEN
                <* ASSERT pe # NIL *>
                WITH enew = NEW(EdgeList) DO
                  enew.e := pe; enew.rest := epos; epos := enew
                END;
                tmp := NEW(TriangleList);
                tmp.t := pt1; tmp.rest := tpos; tpos := tmp; tmp := NIL
              END;
            END;
          END
        END
      END;
      tmp := rst
    END
  END CutTriangles;    

PROCEDURE SplitTriangle(
    t: REF Triangle;
    iv: [0..2];
    VAR nt0: REF Triangle;
    VAR ne: REF Edge;
    VAR nt1: REF Triangle;
    VAR ze: REF Edge;
    VAR pt0: REF Triangle;
    VAR pe: REF Edge;
    VAR pt1: REF Triangle;
    negF, posF: BOOLEAN;
  ) =
  (* 
    Cuts a triangle "t" by the current plane "pl".
    Returns the negative part of "t" in 
    "nt0", "ne", "nt1", the zero part in "ze", and the positive
    part in "pt0", "pe", "pt1".  The edge "ne" and the triangle 
    "nt1" are both NIL or both non-NIL; ditto for "pe" and "pt1".
    
    The negative parts are created only if "negF=TRUE"; otherwise
    "nt0", "ne", and "nt1" are set to NIL.  In the same way, "posF"
    determines the creation of the the positive parts "pt0", "pe", "pt1".

    The procedure assumes the vertex values have been computed, and
    the edges that cross "pl" have been split.
    
    Also assumes that the sign of vertex "iv" is different from
    that of "iv+1", and strictly opposite to that of "iv+2".
  *)
  VAR xv: REF Vertex;   (* Cut vertex on edge "av--bv" *)
      yv: REF Vertex;   (* Cut vertex on edge "av--cv" *)
      ke: REF Edge;     (* Diagonal edge "bv--yv", or NIL *)
      at: REF Triangle; (* Triangle next to vertex "av" *)
      ct: REF Triangle; (* Triangle next to vertex "cv" *)
      zt: REF Triangle; (* Triangle next to edge "ze", or NIL *)
      makea: BOOLEAN;   (* TRUE to create "at" *)
      makec: BOOLEAN;   (* TRUE to create "ct" and "zt" *)
  BEGIN
    WITH
      a = iv, 
      b = (iv+1) MOD 3, 
      c = (iv+2) MOD 3,
      ae = t.edge[a], ad = t.edir[a], av = ae.vert[ad], as = av.val,
      be = t.edge[b], bd = t.edir[b], bv = be.vert[bd], bs = bv.val,
      ce = t.edge[c], cd = t.edir[c], cv = ce.vert[cd], cs = cv.val
    DO
      <* ASSERT as*bs <= 0.0d0 AND as*cs < 0.0d0 *>
      <* ASSERT ce.part[0] # ce.part[1] *>
      <* ASSERT be.part[0] = be.part[1] *>
      <* ASSERT be.part[0] = be *>
      
      <* ASSERT negF OR posF *>
      IF as <= 0.0d0 THEN
        makea := negF; makec := posF
      ELSE
        makea := posF; makec := negF
      END;

      (* Cutting edge on cutting plane: *)
      IF bs = 0.0d0 THEN
        xv := bv
      ELSIF makea THEN
        xv := ae.part[ad].vert[1-ad]
      ELSE
        xv := ae.part[1-ad].vert[ad]
      END;
      IF makea THEN
        yv := ce.part[1-cd].vert[cd]
      ELSE
        yv := ce.part[cd].vert[1-cd];
      END;
      <* ASSERT xv.val = 0.0d0 *>
      <* ASSERT yv.val = 0.0d0 *>
      (* Create cut line, if needed: *)
      WITH pl = t.plane^ DO
        IF NOT pl.mark THEN
          WITH 
            cutln = NEW(REF Line, ln := H3.LineFromTwoPoints(xv.wcs, yv.wcs))
          DO
            cutln.cut := NIL;
            cutln.mark := TRUE;
            pl.cut := cutln
          END;
          pl.mark := TRUE
        END;
        ze := MakeEdge(t.org, xv, yv, pl.cut);
      END;
      
      (* Triangle next to vertex "av": *)
      IF makea THEN
        at := MakeTriangle(
          t.org, 
          ae.part[ad], ad,
          ze, 0,
          ce.part[1-cd], cd,
          plane := t.plane,
          lights := t.lights,
          rootv := a
        )
      ELSE
        at := NIL
      END;

      IF makec THEN
        IF bs = 0.0d0 THEN
          (* Only two triangles: *)
          <* ASSERT ae.part[0] = ae.part[1] *>
          ke := NIL; 
          (* Triangle next to vertex "cv": *)
          ct := MakeTriangle(
            t.org,
            ze, 1,
            be, bd,
            ce.part[cd], cd,
            plane := t.plane,
            lights := t.lights,
            rootv := a
          );
          zt := NIL;
        ELSE
          (* Three triangles: *)
          <* ASSERT ae.part[0] # ae.part[1] *>
          ke := MakeEdge(t.org, bv, yv);
          (* Triangle next to vertex "cv": *)
          ct := MakeTriangle(
            t.org,
            ke, 1,
            be, bd,
            ce.part[cd], cd,
            plane := t.plane,
            lights := t.lights,
            rootv := a
          );
          zt := MakeTriangle(
            t.org,
            ae.part[1-ad], ad,
            ke, 0,
            ze, 1,
            plane := t.plane,
            lights := t.lights,
            rootv := a
          )
        END
      ELSE
        ke := NIL;
        ct := NIL;
        zt := NIL
      END;
      (* Return in appropriate variables: *)
      IF as <= 0.0d0 THEN
        nt0 := at; ne := NIL; nt1 := NIL;
        pt0 := ct; pe := ke; pt1 := zt;
      ELSE
        nt0 := ct; ne := ke; nt1 := zt;
        pt0 := at; pe := NIL; pt1 := NIL;
      END
    END
  END SplitTriangle;

PROCEDURE CopyScene(src: T;  VAR dst: T) =
  BEGIN
    CopyVertexList  (src.v, dst.v);
    CopyEdgeList    (src.e, dst.e);
    CopyTriangleList(src.t, dst.t);
  END CopyScene;

PROCEDURE CopyVertexList(src: VertexList; VAR dst: VertexList) =
  (* Make a top-level copy of "src", prepends it to "dst": *)
  BEGIN
    WHILE src # NIL DO
      WITH nt = NEW(VertexList) DO
        nt.v := src.v; nt.rest := dst; dst := nt
      END;
      src := src.rest
    END;
  END CopyVertexList;

PROCEDURE CopyEdgeList(src: EdgeList; VAR dst: EdgeList) =
  (* Make a top-level copy of "src", prepends it to "dst": *)
  BEGIN
    WHILE src # NIL DO
      WITH nt = NEW(EdgeList) DO
        nt.e := src.e; nt.rest := dst; dst := nt
      END;
      src := src.rest
    END;
  END CopyEdgeList;
  
PROCEDURE CopyTriangleList(src: TriangleList; VAR dst: TriangleList) =
  (* Make a top-level copy of "src", prepends it to "dst": *)
  BEGIN
    WHILE src # NIL DO
      WITH nt = NEW(TriangleList) DO
        nt.t := src.t; nt.rest := dst; dst := nt
      END;
      src := src.rest
    END;
  END CopyTriangleList;

PROCEDURE MoveScene(VAR src: T;  VAR dst: T) =
  BEGIN
    MoveVertexList  (src.v, dst.v);
    MoveEdgeList    (src.e, dst.e);
    MoveTriangleList(src.t, dst.t);
  END MoveScene;

PROCEDURE MoveVertexList(VAR src: VertexList; VAR dst: VertexList) =
  (* Destructively prepend "src" onto "dst": *)
  BEGIN
    WHILE src # NIL DO
      VAR r := src.rest; BEGIN src.rest := dst; dst := src; src := r END;
    END;
  END MoveVertexList;

PROCEDURE MoveEdgeList(VAR src: EdgeList; VAR dst: EdgeList) =
  (* Destructively prepend "src" onto "dst": *)
  BEGIN
    WHILE src # NIL DO
      VAR r := src.rest; BEGIN src.rest := dst; dst := src; src := r END;
    END;
  END MoveEdgeList;

PROCEDURE MoveTriangleList(VAR src: TriangleList; VAR dst: TriangleList) =
  (* Destructively prepend "src" onto "dst": *)
  BEGIN
    WHILE src # NIL DO
      VAR r := src.rest; BEGIN src.rest := dst; dst := src; src := r END;
    END;
  END MoveTriangleList;

PROCEDURE CountVertices(lst: VertexList): CARDINAL =
  VAR n: CARDINAL := 0;
  BEGIN
    WHILE lst # NIL DO INC(n); lst := lst.rest END;
    RETURN n
  END CountVertices;
  
PROCEDURE CountEdges(lst: EdgeList): CARDINAL =
  VAR n: CARDINAL := 0;
  BEGIN
    WHILE lst # NIL DO INC(n); lst := lst.rest END;
    RETURN n
  END CountEdges;
  
PROCEDURE CountTriangles(lst: TriangleList): CARDINAL =
  VAR n: CARDINAL := 0;
  BEGIN
    WHILE lst # NIL DO INC(n); lst := lst.rest END;
    RETURN n
  END CountTriangles;
  
PROCEDURE PrintCounts(wr: Wr.T; sc: T) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    Wr.PutText(wr, "  { v = ");
    Wr.PutText(wr, Fmt.Int(CountVertices(sc.v)));
    Wr.PutText(wr, "  e = ");
    Wr.PutText(wr, Fmt.Int(CountEdges(sc.e)));
    Wr.PutText(wr, "  t = ");
    Wr.PutText(wr, Fmt.Int(CountTriangles(sc.t)));
    Wr.PutText(wr, " }");
  END PrintCounts;
  
PROCEDURE PrintNum(wr: Wr.T; r: REFANY; pad: CARDINAL := 0) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  VAR x: TEXT;
  BEGIN
    TYPECASE r OF
    | NULL            => x := "NIL";
    | REF Vertex(v)   => x := Fmt.Int(v.num) & "v";
    | REF Edge(e)     => x := Fmt.Int(e.num) & "e";
    | REF Triangle(t) => x := Fmt.Int(t.num) & "t";
    ELSE                 x := "???";
    END;
    Wr.PutText(wr, Fmt.Pad(x, pad))
  END PrintNum;
  
PROCEDURE PrintBoolean(wr: Wr.T; b: BOOLEAN) =
  CONST Boole = ARRAY BOOLEAN OF CHAR{'F', 'T'};
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    Wr.PutChar(wr, Boole[b])
  END PrintBoolean;
  
PROCEDURE PrintReal(wr: Wr.T; x: REAL) =
  CONST SigDigits = 6;
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    Wr.PutText(wr, Fmt.Real(x, prec := SigDigits, style := Fmt.Style.Auto))
  END PrintReal;
  
PROCEDURE PrintPoint(wr: Wr.T; READONLY p: H3.Point) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    FOR i := 0 TO 3 DO
      IF i = 0 THEN Wr.PutText(wr, "[") ELSE Wr.PutText(wr, ", ") END;
      PrintReal(wr, p.c[i])
    END;
    Wr.PutText(wr, "]");
  END PrintPoint;
  
PROCEDURE PrintVertex(wr: Wr.T; v: REF Vertex) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    PrintNum(wr, v, 5);
    Wr.PutText(wr, "  org = ");
    PrintNum(wr, v.org);
    Wr.PutText(wr, "  mark = ");
    PrintBoolean(wr, v.mark);
    Wr.PutText(wr, "  wcs = ");
    PrintPoint(wr, v.wcs);
    Wr.PutText(wr, "  val = ");
    Wr.PutText(wr, Fmt.Pad(Fmt.LongReal(v.val, Fmt.Style.Sci, 6), 14));
  END PrintVertex;
  
PROCEDURE PrintEdge(wr: Wr.T; e: REF Edge) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    PrintNum(wr, e, 5);
    Wr.PutText(wr, "  org = ");
    PrintNum(wr, e.org);
    Wr.PutText(wr, "  mark = ");
    PrintBoolean(wr, e.mark);
    Wr.PutText(wr, "  vert = {");
    FOR i := 0 TO 1 DO
      IF i > 0 THEN Wr.PutText(wr, ", ") END;
      PrintNum(wr, e.vert[i], 5);
    END;
    Wr.PutText(wr, "}");
    Wr.PutText(wr, "  part = {");
    FOR i := 0 TO 1 DO
      IF i > 0 THEN Wr.PutText(wr, ", ") END;
      PrintNum(wr, e.part[i], 5);
    END;
    Wr.PutText(wr, "}");
  END PrintEdge;

PROCEDURE PrintTriangle(wr: Wr.T; t: REF Triangle) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    PrintNum(wr, t, 5);
    Wr.PutText(wr, "  org = ");
    PrintNum(wr, t.org);
    Wr.PutText(wr, "  mark = ");
    PrintBoolean(wr, t.mark);
    Wr.PutText(wr, "  edge = {");
    FOR i := 0 TO 2 DO
      IF i > 0 THEN Wr.PutText(wr, ", ") END;
      PrintNum(wr, t.edge[i], 5);
      Wr.PutChar(wr, ':');
      Wr.PutText(wr, Fmt.Int(t.edir[i]))
    END;
    Wr.PutText(wr, "}");
    Wr.PutText(wr, "  light = ");
    FOR i := 0 TO MaxLights-1 DO 
      PrintBoolean(wr, t.lights[i])
    END;
  END PrintTriangle;

PROCEDURE PrintScene(wr: Wr.T; sc: T) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    Wr.PutText(wr, "-- vertices --\n");
    VAR lst: VertexList := sc.v;
    BEGIN
      WHILE lst # NIL DO 
        PrintVertex(wr, lst.v);
        Wr.PutText(wr, "\n");
        lst := lst.rest;
      END;
    END;
    
    Wr.PutText(wr, "-- edges --\n");
    VAR lst: EdgeList := sc.e;
    BEGIN
      WHILE lst # NIL DO 
        PrintEdge(wr, lst.e);
        Wr.PutText(wr, "\n");
        lst := lst.rest;
      END;
    END;
    
    Wr.PutText(wr, "-- triangles --\n");
    VAR lst: TriangleList := sc.t;
    BEGIN
      WHILE lst # NIL DO 
        PrintTriangle(wr, lst.t);
        Wr.PutText(wr, "\n");
        lst := lst.rest;
      END;
    END;
  END PrintScene;

PROCEDURE PlotScene(
    filename: TEXT; 
    sc: T;
    xmin, xmax: LONGREAL;
    ymin, ymax: LONGREAL;
    READONLY w2i: H3.PMap;
  ) =
  BEGIN
    WITH
      ps = NEW(PSPlot.EPSFile).open(filename & ".eps")
    DO
      ps.setScale(PSPlot.Axis.X, xmin, xmax);
      ps.setScale(PSPlot.Axis.Y, ymin, ymax);
      ps.setLineColor(PSPlot.Invisible);
      ps.setFillColor(PSPlot.Gray(0.75));

      VAR tl: TriangleList := sc.t;
      BEGIN
        WHILE tl # NIL DO 
          WITH f = tl.t DO
            PlotTriangle(ps, f^, w2i)
          END;
          tl := tl.rest
        END;
      END;

      ps.setLineColor(PSPlot.Black);

      VAR el: EdgeList;
      BEGIN
        el := sc.e;
        WHILE el # NIL DO 
          WITH e = el.e DO
            PlotEdge(ps, e^, w2i)
          END;
          el := el.rest
        END;
      END;

      (* Finish off: *)
      ps.close()
    END
  END PlotScene;
  
PROCEDURE PlotTriangle(
    ps: PSPlot.File;
    READONLY f: Triangle;
    READONLY w2i: H3.PMap;
  ) =
  BEGIN
    WITH
      av = f.edge[0].vert[f.edir[0]],
      bv = f.edge[1].vert[f.edir[1]],
      cv = f.edge[2].vert[f.edir[2]],
      
      avi = H3.MapPoint(av.wcs, w2i),
      bvi = H3.MapPoint(bv.wcs, w2i),
      cvi = H3.MapPoint(cv.wcs, w2i),

      ax = FLOAT(avi.c[1]/avi.c[0], LONGREAL),
      ay = FLOAT(avi.c[2]/avi.c[0], LONGREAL),
      bx = FLOAT(bvi.c[1]/bvi.c[0], LONGREAL),
      by = FLOAT(bvi.c[2]/bvi.c[0], LONGREAL),
      cx = FLOAT(cvi.c[1]/cvi.c[0], LONGREAL),
      cy = FLOAT(cvi.c[2]/cvi.c[0], LONGREAL)
    DO
      ps.triangle(ax, ay, bx, by, cx, cy);
    END
  END PlotTriangle;

PROCEDURE PlotEdge(
    ps: PSPlot.File;
    READONLY e: Edge;
    READONLY w2i: H3.PMap;
  ) =
  BEGIN
    WITH
      av = e.vert[0],
      bv = e.vert[1],
      
      avi = H3.MapPoint(av.wcs, w2i),
      bvi = H3.MapPoint(bv.wcs, w2i),

      ax = FLOAT(avi.c[1]/avi.c[0], LONGREAL),
      ay = FLOAT(avi.c[2]/avi.c[0], LONGREAL),
      bx = FLOAT(bvi.c[1]/bvi.c[0], LONGREAL),
      by = FLOAT(bvi.c[2]/bvi.c[0], LONGREAL)
    DO
      ps.segment(ax, ay, bx, by);
    END
  END PlotEdge;
  
BEGIN
END TrArScene.