PROCEDURE AuxVisible(
    VAR ins, bds: TrArScene.T;
    READONLY obs: HLR3.Point;
    invF: BOOLEAN;
    VAR vis, inv: TriangleList;
  ) = 
  (*
    Given a closed scene with opaque triangles, prepends onto "vis"
    and "inv" the parts of triangles that are visible and hidden from
    the point "obs" by other parts of the same scene.
    
    The invisible parts are returned only if "invF" is TRUE.

    Assumes that the scene has been clipped to some ray cone "R"
    with apex "obs", and is being presented in two sections: the
    parts "bds" on the walls of "R", and the parts "ins" strictly
    inside "R".  
    
    Note that "bds" is closed by itself, but the closure of "ins" 
    may include elements of "bds".

    The procedure will reuse the list nodes of "ins" and "bds".
  *)
  
  VAR 
    coins: Random.T := Random.New();
    nCalls: CARDINAL := 0;  (* Calls to AuxAuxVisible *)

  PROCEDURE AuxAuxVisible(
      VAR inx: TrArScene.T;    (* Elements inside "R". *)    
      VAR bdx: TrArScene.T;    (* Elements on the boundary of "R". *)
      VAR bg: TriangleList;  (* Background triangle, or NIL. *)
      depth: CARDINAL;       (* Recursion depth. *)
      dir: [0..2];           (* Preferred split direction. *)
    ) = 
    (*
      Same as "AuxVisible(Union(inx, bg), bdx), obs, vis)".
      The parameter "bg", if not NIL, is a list containing a
      single triangle that lies entirely inside
      the region "R", covers its cross-section completely, 
      and lies behind all elements of "inx".
    *)

    VAR debugCall: BOOLEAN := Debug AND depth >= DebugMinDepth AND depth <= DebugMaxDepth;
        ivis: CARDINAL;  (* Debug *)
        call: CARDINAL := nCalls;  (* Debug *)

    PROCEDURE SplitByTriangle() =
      (* 
        Split scene by plane of first triangle in "inx.t",
        assumed to be a cross-section of "R".
        Discards what is behind, and calls AuxAuxVisible 
        on the rest.
      *)
      VAR bdn, bdp: TrArScene.T := TrArScene.Empty;
          inn, inp: TrArScene.T := TrArScene.Empty;
          pl: HLR3.Plane;
          valid: BOOLEAN;
      BEGIN
        IF debugCall THEN Wr.PutText(stderr, "-- triangle split --\n") END;
        <* ASSERT inx.e = NIL *>
        REPEAT
          Wr.PutChar(stderr, 'T');
          IF inx.t = NIL AND bg = NIL THEN
            RETURN
          ELSIF inx.t = NIL THEN
            IF bg # NIL THEN bg.rest := vis; vis := bg; bg := NIL END;
          ELSIF inx.t.rest = NIL AND bg = NIL THEN
            inx.t.rest := vis; vis := inx.t; inx.t := NIL;
          ELSE
            bg := inx.t; inx.t := bg.rest; bg.rest := NIL;
            WITH tsplit = bg.t DO
              pl := TrianglePlane(tsplit^);
              valid := PlaneValid(pl);
              IF valid THEN
                IF HLR3.Side(obs, pl) = -1 THEN pl.f := R4.Neg(pl.f) END;
                (* Evaluate all vertices relative to splitting plane: *)
                TrArScene.EvalVertices(bdx.v, pl);
                TrArScene.EvalVertices(inx.v, pl);
                (* Force splitting triangle to lie on splitting plane: *)
                FOR j := 0 TO 2 DO
                  FOR i := 0 TO 1 DO tsplit.edge[j].vert[i].val := 0.0d0 END
                END;
              ELSE
                (* Discard triangle, try again: *)
                IF invF THEN bg.rest := inv; inv := bg; bg := NIL END
              END;
            END
          END
        UNTIL valid;
        (* Now cut scene: *)
        (* Note that we must cut "bdx" before cutting "inx", because *)
        (*   the triangles of "inx" may have edges in "bdx". *)
        TrArScene.CutScene(bdx, invF, TRUE, TRUE, bdn, bdp, bdp);
        TrArScene.CutScene(inx, invF, TRUE, TRUE, inn, inp, inp);
        (* Things on negative side are invisible: *)
        IF invF THEN
          TrArScene.MoveTriangleList(inn.t, inv); 
          inn.e := NIL; 
          inn.v := NIL;
        ELSE
          <* ASSERT inn = TrArScene.Empty *>
        END;
        <* ASSERT bdn.t = NIL *>
        bdn.e := NIL;
        bdn.v := NIL;
        <* ASSERT bdp.t = NIL *>
        (* Note that the cut may have created new edges inside "R". *)
        (* Recurse on parts that are left: *)
        AuxAuxVisible(inp, bdp, bg, depth + 1, dir)
      END SplitByTriangle;
      
    PROCEDURE PlaneValid(READONLY pl: HLR3.Plane): BOOLEAN =
      BEGIN
        RETURN pl.f # R4.Zero
      END PlaneValid;

    PROCEDURE SplitByEdge(): BOOLEAN =
      (* 
        Splits "R" by plane through observer,
        and applies AuxAuxVisible to each half.
        Returns FALSE if was unable to split;
        in that case "inx.e" will be NIL.
      *)
      CONST MaxTrials = 5;
      VAR bdn, bdz, bdp: TrArScene.T := TrArScene.Empty;
          inn, inz, inp: TrArScene.T := TrArScene.Empty;
          pl: HLR3.Plane;
          valid: BOOLEAN;
          VAR trial: CARDINAL := 0;
      BEGIN
        IF debugCall THEN Wr.PutText(stderr, "-- edge split --\n") END;
        REPEAT 
          INC(trial);
          IF inx.e = NIL THEN RETURN FALSE END;
          (* Put "bg" back into "inx": *)
          IF bg # NIL THEN
            bg.rest := inx.t; inx.t := bg; bg := NIL
          END;
          Wr.PutChar(stderr, '-');
          WITH 
            vsplit = PickVertexPair(inx.e, trial, MaxTrials)
          DO
            pl := ObsPlane(obs, vsplit[0]^, vsplit[1]^);
            IF debugCall THEN 
              Wr.PutText(stderr, "plane = ");
              PrintPoint(stderr, HLR3.Point{c := pl.f});
              Wr.PutText(stderr, "\n");
            END;
            valid := PlaneValid(pl);
            IF valid THEN
              (* Evaluate all vertices relative to splitting plane: *)
              TrArScene.EvalVertices(bdx.v, pl);
              TrArScene.EvalVertices(inx.v, pl);
              (* Force splitting vertices to lie on splitting plane: *)
              FOR i := 0 TO 1 DO vsplit[i].val := 0.0d0 END;
            ELSE
              (* Remove any edges or triangles incident to both these two vertices: *)
              RemoveBadElements(vsplit, inx, inv);
            END
          END;
        UNTIL valid AND GoodEnoughCut(inx.e, trial, MaxTrials);
        (* Now cut scene: *)
        (* Note that we must cut "bdx" before cutting "inx", because *)
        (*   the triangles of "inx" may have edges in "bdx". *)
        TrArScene.CutScene(bdx, TRUE, TRUE, TRUE, bdn, bdz, bdp);
        TrArScene.CutScene(inx, TRUE, TRUE, TRUE, inn, inz, inp);
        IF debugCall THEN
          Wr.PutText(stderr, "\n");
          Wr.PutText(stderr, "split result:\n");
          Wr.PutText(stderr, "  inn = "); PrintCounts(stderr, inn); Wr.PutText(stderr, "\n");
          Wr.PutText(stderr, "  inz = "); PrintCounts(stderr, inz); Wr.PutText(stderr, "\n");
          Wr.PutText(stderr, "  inp = "); PrintCounts(stderr, inp); Wr.PutText(stderr, "\n");
        END;
        (* Triangles on boundary are invisible: *)
        IF invF THEN
          TrArScene.MoveTriangleList(inz.t, inv)
        ELSE
          inz.t := NIL;
        END;
        (* Edges and vertices on cutting plane are on the boundary of each half: *)
        <* ASSERT bdz.t = NIL *>
        <* ASSERT inz.t = NIL *>
        TrArScene.CopyScene(bdz, bdn); TrArScene.MoveScene(bdz, bdp); 
        TrArScene.CopyScene(inz, bdn); TrArScene.MoveScene(inz, bdp); 
        (* Recurse: *)
        <* ASSERT bg = NIL *>
        AuxAuxVisible(inn, bdn, bg, depth + 1, (dir + 1) MOD 3);
        AuxAuxVisible(inp, bdp, bg, depth + 1, (dir + 1) MOD 3);
        RETURN TRUE
      END SplitByEdge;    
      
    PROCEDURE PickVertexPair(
        els: EdgeList; 
        trial: [1..LAST(CARDINAL)]; 
        expTrials: [1..LAST(CARDINAL)];
      ): TrArScene.VertexPair =
      (*
        Odd trials are edges of "els".
        Even trials are random pairs of endpoints.
        "expTrials" is expected number of trials.
      *)
      BEGIN
        IF trial MOD 2 = 1 THEN
          RETURN PickRandomEdge(els, (trial + 1) DIV 2, (expTrials + 1) DIV 2)
        ELSE
          RETURN PickRandomVertexPair(els, trial DIV 2, expTrials DIV 2)
        END
      END PickVertexPair;
      
    PROCEDURE PickRandomEdge(
        els: EdgeList; 
        trial: [1..LAST(CARDINAL)]; 
        expTrials: [1..LAST(CARDINAL)];
      ): TrArScene.VertexPair =
      (*
        Stratified random sample of "els".
      *)
      VAR ne: CARDINAL;
          t: EdgeList;
          vp: TrArScene.VertexPair;
      BEGIN
        (* Count edges: *)
        t := els; ne := 0;
        WHILE t # NIL DO INC(ne); t := t.rest END;
        <* ASSERT ne > 0 *>

        (* Select range: *)
        WITH
          ie = PickRandomIndex(ne, trial, expTrials)
        DO
          t := els;
          FOR i := 0 TO ie-1 DO t := t.rest END;
          FOR iv := 0 TO 1 DO vp[iv] := t.e.vert[iv] END
        END;
        RETURN vp
      END PickRandomEdge;

    PROCEDURE PickRandomVertexPair(
        els: EdgeList; 
        trial: [1..LAST(CARDINAL)]; 
        expTrials: [1..LAST(CARDINAL)];
      ): TrArScene.VertexPair =
      (*
        Stratified random sample of the endpoint pairs.
      *)
      VAR ne: CARDINAL;
          t: EdgeList;
          vp: TrArScene.VertexPair;
      BEGIN
        (* Count edges: *)
        t := els; ne := 0;
        WHILE t # NIL DO INC(ne); t := t.rest END;
        <* ASSERT ne > 0 *>

        WITH
          nv = 2 * ne,              (* Number of endpoints *)
          np = ne * (ne - 1) DIV 2, (* Number of edge pairs *)
          iefuv = PickRandomIndex(4*np, trial, expTrials),
          ief = UnpackPairIndex(iefuv DIV 4),
          iuv = ARRAY [0..1] OF [0..1]{iefuv MOD 2, iefuv DIV 2 MOD 2}
        DO
          <* ASSERT ief[1] < ief[0] AND ief[0] < ne *>
          FOR i := 0 TO 1 DO
            WITH
              eskip = ief[i],
              iend = iuv[i]
            DO
              t := els;
              FOR k := 1 TO eskip DO t := t.rest END;
              vp[i] := t.e.vert[iend]
            END
          END
        END;
        RETURN vp
      END PickRandomVertexPair;
      
    PROCEDURE UnpackPairIndex(ip: CARDINAL): ARRAY [0..1] OF CARDINAL =
      VAR iu, iv: CARDINAL; 
      BEGIN
        iu := 0;
        WHILE iu * (iu + 1) DIV 2 <= ip DO INC(iu) END;
        iv := ip - iu * (iu - 1) DIV 2;
        <* ASSERT iu * (iu - 1) DIV 2 + iv = ip *>
        RETURN ARRAY [0..1] OF CARDINAL{iu, iv}
      END UnpackPairIndex;
      
    PROCEDURE PickRandomIndex(
        n: CARDINAL;
        trial: [1..LAST(CARDINAL)];
        expTrials: [1..LAST(CARDINAL)]
      ): CARDINAL =
      (*
        Picks a random index in [0..n-1].  
        Tries to keep the first "expTrials" distinct.
      *)
      BEGIN
        WITH
          step = MAX(1, (n + expTrials - 1) DIV expTrials),
          lo = ((trial - 1) * step) MOD n,
          hi = lo + step - 1
        DO
          RETURN Random.Subrange(coins, lo, hi) MOD n
        END
      END PickRandomIndex;
      
    PROCEDURE GoodEnoughCut(
        els: EdgeList; 
        trial: [1..LAST(CARDINAL)]; 
        expTrials: [1..LAST(CARDINAL)];
      ): BOOLEAN =
      (*
        TRUE if the current "val" fields define a reasonably
        balanced cut of the edges in "els".
      *)
      <* FATAL Wr.Failure, Thread.Alerted *>
      BEGIN
        (* If there are at most two edges, the cut is good enough: *)
        IF els = NIL
        OR els.rest = NIL
        OR els.rest.rest = NIL
        THEN 
          RETURN TRUE
        ELSE
          (* Count edges on each side of cut: *)
          VAR ne, nn, np: CARDINAL := 0;
              t: EdgeList := els;
          BEGIN
            WHILE t # NIL DO 
              WITH
                e = t.e,
                s0 = e.vert[0].val,
                s1 = e.vert[1].val
              DO
                INC(ne);
                IF s0*s1 < 0.0d0 THEN
                  (* Edge will be split: *)
                  INC(nn); INC(np)
                ELSE
                  (* Edge will not be split: *)
                  IF s0 < 0.0d0 OR s1 < 0.0d0 THEN
                    INC(nn)
                  ELSIF s0 > 0.0d0 OR s1 > 0.0d0 THEN
                    INC(np)
                  ELSE
                    (* Edge lies on cutting plane *)
                  END
                END
              END;
              t := t.rest
            END;
            IF Debug THEN
              Wr.PutChar(stderr, '[');
              Wr.PutText(stderr, Fmt.Int(ne)); 
              Wr.PutChar(stderr, '>');
              Wr.PutText(stderr, Fmt.Int(nn)); 
              Wr.PutChar(stderr, ':');
              Wr.PutText(stderr, Fmt.Int(np));
              Wr.PutChar(stderr, ']');
            END;
            
            <* ASSERT nn + np <= 2*ne - 2 *>

            (* If we are over budget, accept anything: *)
            IF trial > expTrials THEN
              RETURN TRUE
            ELSE
              (* Decide whether cut is balanced: *)
              WITH 
                frac = FLOAT(trial)/FLOAT(expTrials),
                lo = FLOOR(FLOAT(ne) * 0.5 * (1.0 - frac)),
                hi = CEILING(FLOAT(ne) * 0.5 * (1.0 + frac))
              DO
                RETURN MIN(nn, np) >= lo AND MAX(nn,np) <= hi 
              END
            END
          END
        END
      END GoodEnoughCut;
      
    PROCEDURE RemoveBadElements(
        READONLY vp: TrArScene.VertexPair;
        VAR sc: TrArScene.T;
        VAR inv: TriangleList;
      ) =
      (*
        Given two vertices "vp[0]" and "vp[1]" that are collinear with
        "obs", removes from the scene "sc" any edges or triangles that are 
        incident to both "vp[0]" and "vp[1]", since those elements are 
        beeing viewed edge-on.
      *)
      BEGIN
        (* Discard bad edges: *)
        VAR ls, lsn: TrArScene.EdgeList := sc.e;
            lsa: TrArScene.EdgeList := NIL; 
        BEGIN
          WHILE ls # NIL DO
            lsn := ls.rest;
            WITH e = ls.e, u = e.vert[0], v = e.vert[1] DO
              IF (u = vp[0] OR u = vp[1]) AND (v = vp[0] OR v = vp[1]) THEN
                IF lsa = NIL THEN sc.e := ls.rest ELSE lsa.rest := ls.rest END;
                ls.rest := NIL
              ELSE
                lsa := ls
              END
            END;
            ls := lsn
          END
        END;
      
        (* Discard bad triangles: *)
        VAR ls: TrArScene.TriangleList := sc.t;
            lsa, lsn: TrArScene.TriangleList := NIL;
            badvs: CARDINAL;
        BEGIN
          WHILE ls # NIL DO
            lsn := ls.rest;
            WITH t = ls.t DO
              badvs := 0;
              FOR i := 0 TO 2 DO
                WITH 
                  j = (i+1) MOD 3, 
                  u = t.edge[i].vert[1-t.edir[i]], 
                  v = t.edge[j].vert[t.edir[j]]
                DO 
                  IF u = vp[0] OR u = vp[1] OR v = vp[0] OR v = vp[1] THEN
                    INC(badvs)
                  END
                END
              END;
              IF badvs >= 2 THEN
                IF lsa = NIL THEN sc.t := ls.rest ELSE lsa.rest := ls.rest END;
                IF invF THEN ls.rest := inv; inv := ls ELSE ls.rest := NIL END;
              ELSE
                lsa := ls
              END
            END;
            ls := lsn
          END
        END;

      END RemoveBadElements;

    <* FATAL Wr.Failure, Thread.Alerted *>
    BEGIN
      INC(nCalls);
      Wr.PutChar(stderr, '(');
      IF debugCall THEN
        ivis := CountTriangles(vis);
        Wr.PutText(stderr, "\n");
        Wr.PutText(stderr, "Entering AuxAuxVisible(" & Fmt.Int(call) & "):\n");
        Wr.PutText(stderr, "  bdx = "); PrintCounts(stderr, bdx); Wr.PutText(stderr, "\n");
        Wr.PutText(stderr, "  inx = "); PrintCounts(stderr, inx); Wr.PutText(stderr, "\n");
        Wr.PutText(stderr, "  bg = " & Fmt.Int(CountTriangles(bg)) & " triangles.\n");
        IF depth <= DebugMinDepth + 3 THEN
          DebugPlotScene("inx", call, inx, obs)
        END;
      END;
      <* ASSERT bdx.t = NIL *>
      IF NOT SplitByEdge() THEN SplitByTriangle() END;
      IF debugCall THEN
        WITH ovis = CountTriangles(vis) DO
          Wr.PutText(stderr, "\n");
          Wr.PutText(stderr, "Exiting AuxAuxVisible: " & Fmt.Int(ovis) & " triangles.\n");
          Wr.PutText(stderr, "\n");
          <* ASSERT ovis >= ivis *>
        END
      END;
      Wr.PutChar(stderr, ')');
    END AuxAuxVisible;

  VAR bgs: TriangleList;
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    Wr.PutText(stderr, "\n");
    (* Any boundary triangles are invisible: *)
    IF invF THEN
      TrArScene.MoveTriangleList(bds.t, inv)
    ELSE
      bds.t := NIL;
    END;
    bgs := NIL;
    AuxAuxVisible(ins, bds, bgs, 0, 0);
    Wr.PutText(stderr, "\n");
    Wr.PutText(stderr, Fmt.Int(nCalls) & " calls\n");
  END AuxVisible;
  
PROCEDURE TrianglePlane(
    READONLY t: TrArScene.Triangle
  ): HLR3.Plane =
  BEGIN
    WITH
      a = t.edge[0].vert[t.edir[0]].wcs,
      b = t.edge[1].vert[t.edir[1]].wcs,
      c = t.edge[2].vert[t.edir[2]].wcs
    DO
      RETURN HLR3.Join(a, b, c)
    END
  END TrianglePlane;

PROCEDURE ObsPlane(
    READONLY obs: HLR3.Point; 
    READONLY a, b: TrArScene.Vertex
  ): HLR3.Plane =
  BEGIN
    RETURN HLR3.Join(obs, a.wcs, b.wcs)
  END ObsPlane;