MODULE BugDragonInt EXPORTS Main;

IMPORT 
  Triang, TrArScene, TrArShade, PSPlot, R3, H3, R4, Process,
  ParseParams, Wr, Fmt, Stdio, Thread;
FROM Oct IMPORT 
  Rot, Sym, Lnext, Onext, Oprev;
FROM Triang IMPORT 
  Topology, Arc, OrgV, LeftF, Coords;
FROM TrArScene IMPORT MaxLights, TriangleList, Plane, CountTriangles;
FROM Stdio IMPORT stderr;

TYPE
  Options = RECORD
      inFile: TEXT;
      outFile: TEXT;
      normalize: BOOLEAN;  (* Normalize vertex coordinates *)
      drawEdges: BOOLEAN;  (* TRUE to draw edges *)
      mapEdges: REAL;  (* Thickness of original map edges *)
      triEdges: REAL;  (* Thickness of original triangulation edges *)
      silEdges: REAL;  (* Thickness of silhouette/intersection edges *)
      cutEdges: REAL;  (* Thickness of fragment edges *)
      obs: H3.Point;
      focus: H3.Point;
      w2i: H3.PMap;
      wireFrame: BOOLEAN;  (* Don't paint faces or remove hidden surfaces *)
      light: REF ARRAY OF LightSource;
    END;
    
  LightSource = RECORD
      pos: H3.Point;      (* Light position; [0,0,0,0] for ambient light *)
      rd2: LONGREAL;      (* Reference distance, squared. *)
      color: R3.T;        (* Color at reference distance. *)
      shadows: BOOLEAN;   (* TRUE if casts shadows *)
    END;
    
CONST
  XMin = -2.0; XMax = +2.0;
  YMin = -2.0; YMax = +2.0;
  ZMax = 1.0e10;

PROCEDURE Main() =
  VAR tri, opq: TriangleList;
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    WITH 
      o = GetOptions(),
      r = Triang.Read(o.inFile),
      t = r.t,
      coords = ComputeVertexCoordinates(t, r.coords^, o.normalize)^,
      bar = ComputeFaceBarycenters(t, coords)^,
      normal = ComputeFaceNormals(t, coords)^,
      thickness = ComputeEdgeThicknesses(
        t, coords, o.obs, 
        o.mapEdges, o.triEdges, o.silEdges
      )^
    DO
      (* Collect triangles: *)
      tri := CreateInitialTriangles(t, coords);
      Wr.PutText(stderr, Fmt.Int(CountTriangles(tri)) & " initial triangles.\n");
      
      (* Save a copy of the list to use as shadow casters: *)
      opq := NIL; TrArScene.CopyTriangleList(tri, opq);
      
      (* Clip against window: *)
      Wr.PutText(stderr, "clipping against window... ");
      TrArShade.ClipToWindow(tri, o.w2i, XMin, XMax, YMin, YMax, ZMax);
      Wr.PutText(stderr, Fmt.Int(CountTriangles(tri)) & " triangles.\n");
      
      (* Remove hidden surfaces: *)
      IF NOT o.wireFrame THEN
        Wr.PutText(stderr, "removing hidden surfaces... ");
        TrArShade.Visible(tri, opq, o.obs);
        Wr.PutText(stderr, Fmt.Int(CountTriangles(tri)) & " triangles.\n");
      
        (* Compute shadows: *)
        FOR i := 0 TO LAST(o.light^) DO
          IF o.light[i].shadows THEN
            Wr.PutText(stderr, "shadowing for light " & Fmt.Int(i) & "... ");
            TrArShade.Shadow(tri, opq, o.light[i].pos, i);
            Wr.PutText(stderr, Fmt.Int(CountTriangles(tri)) & " triangles.\n");
          END;
        END;
      END;
      
      (* Write result: *)
      WritePostscript(
        o.outFile, tri, 
        o.w2i, o.obs, 
        NOT o.wireFrame,
        o.light, 
        bar, normal, 
        o.drawEdges,
        thickness,
        o.silEdges,
        o.cutEdges
      );
    END
  END Main;
  
PROCEDURE ComputeVertexCoordinates(
    READONLY t: Topology;
    READONLY coords: Coords;
    normalize: BOOLEAN;
  ): REF Coords = 
  BEGIN
    WITH
      rc = NEW(REF Coords, t.NV),
      c = rc^
    DO
      FOR i := 0 TO LAST(coords) DO
        c[i] := coords[i];
      END;
      IF normalize THEN
        Triang.NormalizeCoords(t, c)
      END;
      RETURN rc
    END;
  END ComputeVertexCoordinates;

PROCEDURE ComputeFaceNormals(
    READONLY t: Topology;
    READONLY coords: Coords;
  ): REF ARRAY OF R3.T =
  BEGIN
    WITH
      rn = NEW(REF ARRAY OF R3.T, t.NF),
      n = rn^
    DO
      FOR i := 0 TO t.NF-1 DO
        WITH a = Rot(t.side[i]), f = LeftF(a) DO
          IF f.exists THEN
            n[i] := Triang.FaceNormal(a, coords)
          ELSE
            n[i] := R3.T{0.0, 0.0, 0.0}
          END
        END
      END;
      RETURN rn
    END;
  END ComputeFaceNormals;
  
PROCEDURE ComputeFaceBarycenters(
    READONLY t: Topology;
    READONLY coords: Coords;
  ): REF ARRAY OF H3.Point =
  BEGIN
    WITH
      rbar = NEW(REF ARRAY OF H3.Point, t.NF),
      bar = rbar^
    DO
      FOR i := 0 TO t.NF-1 DO
        WITH a = Rot(t.side[i]) DO
          IF LeftF(a).exists THEN
            WITH
              b = Lnext(a),
              c = Lnext(b)
            DO
              bar[i] := H3.FromCartesian(FaceBarycenter(a, b, c, coords))
            END
          ELSE
            bar[i] := H3.Point{c := R4.T{0.0, 0.0, 0.0, 0.0}}
          END
        END
      END;
      RETURN rbar
    END;
  END ComputeFaceBarycenters;
  
PROCEDURE FaceBarycenter(
    a, b, c: Arc; 
    READONLY coords: Coords;
  ): R3.T =
  BEGIN
    WITH
      av = coords[OrgV(a).num],
      bv = coords[OrgV(b).num],
      cv = coords[OrgV(c).num],
      bar = R3.Scale(1.0d0/3.0d0, R3.Add(R3.Add(av, bv), cv))
    DO
      RETURN bar
    END
  END FaceBarycenter;

PROCEDURE ComputeFacePlane(
    a, b, c: Arc;
    READONLY v: ARRAY OF REF TrArScene.Vertex;
  ): REF Plane =
  BEGIN
    WITH
      av = v[OrgV(a).num].wcs,
      bv = v[OrgV(b).num].wcs,
      cv = v[OrgV(c).num].wcs,
      r = NEW(REF Plane, pl := H3.PlaneFromThreePoints(av, bv, cv))
    DO
      RETURN r
    END
  END ComputeFacePlane;
  
PROCEDURE ComputeEdgeThicknesses(
    READONLY t: Topology;
    READONLY coords: Coords;
    READONLY obs: H3.Point;
    mapEdges: REAL;  (* Draw the map edges (those with nonzero radius) *)
    triEdges: REAL;  (* Draw the edges of the triangulation *)
    silEdges: REAL;  (* Draw the silhouette edges (those with "mark = TRUE" *)
  ): REF ARRAY OF REAL =
  BEGIN
    WITH
      rth = NEW(REF ARRAY OF REAL, t.NE),
      th = rth^
    DO
      FOR i := 0 TO t.NE-1 DO
        th[i] := ComputeEdgeThickness(
          t.edge[i], coords, obs, 
          mapEdges, triEdges, silEdges
        )
      END;
      RETURN rth
    END;
  END ComputeEdgeThicknesses;
  
PROCEDURE ComputeEdgeThickness(
    a: Arc;
    READONLY coords: Coords;
    READONLY obs: H3.Point;
    mapEdges: REAL;  (* Thickness of map edges (those with nonzero radius) *)
    triEdges: REAL;  (* Thickness of other triangulation edges *)
    silEdges: REAL;  (* Thickness of silhouette edges *)
  ): REAL =
  CONST Invisible = -1.0;
  VAR thickness: REAL := Invisible;
  BEGIN
    WITH e = NARROW(a.edge, Triang.Edge) DO
      IF e.exists THEN
        IF triEdges > thickness THEN
          thickness := triEdges
        END;
        IF mapEdges > thickness AND e.radius > 0.0 THEN
          thickness := mapEdges
        END;
        IF silEdges > thickness AND IsSilhouette(a, coords, obs) THEN
          thickness := silEdges
        END;
      END
    END;
    RETURN thickness 
  END ComputeEdgeThickness;
  
PROCEDURE IsSilhouette(
    a: Arc;
    READONLY coords: Coords;
    READONLY obs: H3.Point;
  ): BOOLEAN =
  (*
    TRUE iff "a.edge" is a silhouette edge 
  *)
  BEGIN
    WITH
      b = Onext(a),
      c = Oprev(a),

      u = H3.FromCartesian(coords[OrgV(a).num]),
      v = H3.FromCartesian(coords[OrgV(Sym(a)).num]),
      
      x = H3.FromCartesian(coords[OrgV(Sym(b)).num]),
      y = H3.FromCartesian(coords[OrgV(Sym(c)).num]),

      pl = H3.PlaneFromThreePoints(obs, u, v)
    DO
      RETURN R4.Dot(pl.f, x.c)*R4.Dot(pl.f, y.c) >= 0.0d0
    END
  END IsSilhouette;
  
PROCEDURE CreateInitialTriangles(
    READONLY t: Topology;
    READONLY coords: Coords;
  ): TriangleList =
  VAR tri: TriangleList := NIL;
  BEGIN
    WITH
      scv = NEW(REF ARRAY OF REF TrArScene.Vertex, t.NV)^,
      sca = NEW(REF ARRAY OF Arc, t.NE)^,
      sil = NEW(REF ARRAY OF BOOLEAN, t.NE)^,
      sce = NEW(REF ARRAY OF REF TrArScene.Edge, t.NE)^,
      sct = NEW(REF ARRAY OF REF TrArScene.Triangle, t.NF)^
    DO
      FOR i := 0 TO t.NV-1 DO
        WITH 
          tv = t.vertex[i],
          wcs = H3.FromCartesian(coords[i])
        DO
          scv[i] := TrArScene.MakeVertex(tv, wcs)
        END
      END;
      FOR i := 0 TO t.NE-1 DO
        WITH 
          ta = t.edge[i],
          u = scv[OrgV(ta).num],
          v = scv[OrgV(Sym(ta)).num]
        DO
          sce[i] := TrArScene.MakeEdge(ta.edge, u, v);
          sca[i] := ta;
        END;
      END;
      FOR i := 0 TO t.NF-1 DO
        WITH 
          tf = t.face[i],
          a = Rot(t.side[i]), ia = NARROW(a.edge, Triang.Edge).num,
          b = Lnext(a),       ib = NARROW(b.edge, Triang.Edge).num,
          c = Lnext(b),       ic = NARROW(c.edge, Triang.Edge).num
        DO
          sct[i] := TrArScene.MakeTriangle(
            tf,
            sce[ia], ComputeEDir(a, sca[ia]),
            sce[ib], ComputeEDir(b, sca[ib]),
            sce[ic], ComputeEDir(c, sca[ic]),
            plane := ComputeFacePlane(a, b, c, scv)
          );
          IF tf.exists THEN
            WITH tl = NEW(TriangleList) DO
              tl.t := sct[i]; tl.rest := tri; tri := tl
            END
          END
        END;
      END;
    END;
    RETURN tri
  END CreateInitialTriangles;
  
PROCEDURE ComputeEDir(a: Arc; r: Arc): TrArScene.EDir =
  (*
    Computes the direction of "a" relative to "r" (which must be
    either "a" or "Sym(a)". 
  *)
  BEGIN
    IF a = r THEN 
      RETURN 0
    ELSIF a = Sym(r) THEN
      RETURN 1
    ELSE
      <* ASSERT FALSE *>
    END
  END ComputeEDir;
  
PROCEDURE WritePostscript(
    filename: TEXT; 
    tri: TriangleList;
    READONLY w2i: H3.PMap;
    READONLY obs: H3.Point;
    paintFaces: BOOLEAN;
    light: REF ARRAY OF LightSource;
    READONLY bar: ARRAY OF H3.Point;
    READONLY normal: ARRAY OF R3.T;
    drawEdges: BOOLEAN;
    READONLY thickness: ARRAY OF REAL;
    <*UNUSED*> silEdges: REAL;
    cutEdges: REAL;
  ) =
  BEGIN
    WITH 
      ps = NEW(PSPlot.EPSFile).open(filename & ".eps")
    DO
      ps.setScale(PSPlot.Axis.X, LR(XMin), LR(XMax));
      ps.setScale(PSPlot.Axis.Y, LR(YMin), LR(YMax));
      
      IF paintFaces THEN
        VAR tl: TriangleList := tri;
        BEGIN
          WHILE tl # NIL DO 
            WITH
              f = tl.t,
              fn = NARROW(f.org, Triang.Face).num
            DO
              WritePostscriptTriangle(
                ps, f^, bar[fn], normal[fn], 
                w2i, obs, light^
              )
            END;
            tl := tl.rest
          END;
        END;
      END;
      
      IF drawEdges THEN
        VAR sc: TrArScene.T; 
            el: TrArScene.EdgeList;
            th: REAL;
        BEGIN
          (* Collect edges: *)
          sc.t := tri; sc.e := NIL; sc.v := NIL;
          TrArScene.CloseScene(sc);
          el := sc.e;
          WHILE el # NIL DO 
            WITH e = el.e DO
              IF e.org = NIL THEN
                th := cutEdges
              ELSE
                TYPECASE e.org OF
                | Triang.Edge(et) => th := thickness[et.num]
                ELSE th := cutEdges
                END              
              END;
              WritePostscriptEdge(ps, e^, th, w2i)
            END;
            el := el.rest
          END;
        END;
      END;
      
      (* Finish off: *)
      ps.close();
    END
  END WritePostscript;
  
PROCEDURE LR(x: REAL): LONGREAL = 
  BEGIN
    RETURN FLOAT(x, LONGREAL) (* Modula-3 is downright silly... *)
  END LR;
      
PROCEDURE WritePostscriptTriangle(
    ps: PSPlot.File;
    READONLY f: TrArScene.Triangle;
    READONLY bar: H3.Point;
    READONLY normal: R3.T;
    READONLY w2i: H3.PMap;
    READONLY obs: H3.Point;
    READONLY light: ARRAY OF LightSource;
  ) =
  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),
      
      fc = ComputeTriangleColor(f, bar, normal, obs, light)
    DO
      ps.setFillColor(fc);
      ps.triangle(ax, ay, bx, by, cx, cy);
    END
  END WritePostscriptTriangle;

PROCEDURE ComputeTriangleColor(
    READONLY f: TrArScene.Triangle;
    READONLY bar: H3.Point;
    READONLY normal: R3.T;
    READONLY obs: H3.Point;
    READONLY light: ARRAY OF LightSource;
  ): R3.T =
  VAR c: R3.T;
      atten: LONGREAL;
  BEGIN
    WITH
      vue = H3.Dir(bar, obs),
      vcos = R3.Dot(vue, normal)
    DO
      c := R3.Zero();
      FOR i := 0 TO LAST(light) DO
        IF f.lights[i] THEN
          WITH li = light[i] DO
            IF li.pos.c = R4.T{0.0, 0.0, 0.0, 0.0} THEN
              (* Ambient light *)
              c := R3.Add(c, li.color)
            ELSE
              WITH
                dir = H3.Dir(bar, li.pos),
                lcos = R3.Dot(dir, normal)
              DO
                IF lcos * vcos <= 0.0d0 THEN
                  (* Light on wrong side of triangle *)
                  atten := 0.0d0
                ELSIF li.pos.c[0] = 0.0 THEN
                  (* Light at infinity *)
                  atten := ABS(lcos)
                ELSE
                  (* Light at finite distance: *)
                  WITH ld2 = H3.DistSqr(bar, li.pos) DO
                    atten := ABS(lcos) * (li.rd2/ld2)
                  END
                END;
                c := R3.Add(c, R3.Scale(atten, li.color))
              END
            END
          END
        END
      END;
      RETURN c
    END
  END ComputeTriangleColor;
  
PROCEDURE WritePostscriptEdge(
    ps: PSPlot.File;
    READONLY e: TrArScene.Edge;
    thickness: REAL;
    READONLY w2i: H3.PMap;
  ) =
  BEGIN
    IF thickness >= 0.0 THEN
      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.setLineWidth(thickness);
        ps.segment(ax, ay, bx, by);
      END
    END
  END WritePostscriptEdge;
  
PROCEDURE GetOptions (): Options =
  <* FATAL Thread.Alerted, Wr.Failure *>
  VAR o: Options;
      up: H3.Point;
  BEGIN
    WITH
      pp = NEW(ParseParams.T).init(stderr)
    DO

      TRY
        pp.getKeyword("-inFile");                               
        o.inFile := pp.getNext();  
        o.outFile :=  o.inFile;

        pp.getKeyword("-obs");
        o.obs := ParsePointOption(pp);

        o.focus := H3.FromCartesian(R3.Zero());
        up := H3.Point{c := R4.T{0.0, 0.0, 0.0, 1.0}};
        o.w2i := H3.PerspMap(o.obs, o.focus, up);
        o.normalize := FALSE;

        o.mapEdges := -1.0;
        o.triEdges := -1.0;
        o.silEdges := -1.0;
        o.cutEdges := -1.0;

        o.drawEdges := 
          o.mapEdges >= 0.0 OR
          o.triEdges >= 0.0 OR
          o.silEdges >= 0.0 OR
          o.cutEdges >= 0.0;

        o.wireFrame := FALSE;
        o.light := DefaultUnshadedLights();

        pp.finish();                                       
      EXCEPT                                                            
      | ParseParams.Error =>                                              
          Wr.PutText(stderr, "Usage: BugDragonInt \\\n");
          Wr.PutText(stderr, "  -inFile <name> \\\n");
          Wr.PutText(stderr, "  -obs <w> <x> <y> <z>\\\n");
          Wr.PutText(stderr, "  ]\n");
          Process.Exit(1);
      END;
    END;
    RETURN o
  END GetOptions;
  
PROCEDURE ParsePointOption(pp: ParseParams.T): H3.Point RAISES {ParseParams.Error} =
  VAR p: H3.Point;
  BEGIN
    FOR i := 0 TO 3 DO
      p.c[i] := pp.getNextReal(-1.0e20, +1.0e20)
    END;
    RETURN p
  END ParsePointOption;
  
PROCEDURE DefaultUnshadedLights(): REF ARRAY OF LightSource =
  BEGIN
    WITH
      r = NEW(REF ARRAY OF LightSource, 1)
    DO
      WITH li = r[0] DO
        (* White ambient light: *)
        li.pos := H3.Point{c := R4.Zero()};
        li.color := R3.T{0.75, 0.75, 0.75};
        li.rd2 := 1.0d0;
        li.shadows := FALSE
      END;
      RETURN r
    END
  END DefaultUnshadedLights;
  
PROCEDURE DefaultShadedLights(): REF ARRAY OF LightSource =
  BEGIN
    WITH
      r = NEW(REF ARRAY OF LightSource, 2) 
    DO
      WITH li = r[0] DO
        (* Some ambient light... *)
        li.pos := H3.Point{c := R4.Zero()};
        li.color := R3.T{0.3, 0.3, 0.3};
        li.rd2 := 1.0d0;
        li.shadows := FALSE
      END;
      WITH li = r[1] DO
        (* ...and a light at infinity: *)
        li.pos := H3.Point{c := R4.T{0.0, 1.7, 2.3, 5.1}};
        li.color := R3.T{0.7, 0.7, 0.7};
        li.rd2 := 1.0d0;
        li.shadows := TRUE
      END;
     RETURN r
    END
  END DefaultShadedLights;

BEGIN 
  Main();
END BugDragonInt.