MODULE PSPlotMesh;

IMPORT PSPlot, LongReal, LR3, LR3x3;

(*DEBUG*) IMPORT Wr, Thread, Fmt;
(*DEBUG*) FROM Stdio IMPORT stderr;

CONST
  Huge = LongReal.MaxFinite;

PROCEDURE Square(
    ps: PSPlot.File;              (* Postscript file. *)
    READONLY g2i: LR3x3.T;        (* Grid-to-image transformation. *)
    READONLY i2g: LR3x3.T;        (* Image-to-grid transformation. *)
    xLo, xHi, yLo, yHi: LONG;     (* Clipping rectangle, in image coordinates. *)
    color: TileColorProc;         (* Maps tile indices to colors. *)
    defined: TileDefProc := NIL;  (* Predicate "tile is defined"; NIL is always "yes". *)
    outlineWidth: REAL := 0.0;    (* If positive, outlines the painted area. *)
  ) =
  BEGIN
    WITH
      (* Windowing lines in grid coordinates: *)
      side = ARRAY [0..3] OF LR3.T{
        LR3x3.MapCol(g2i, LR3.T{-xLo, +1.0d0, 00.0d0}), 
        LR3x3.MapCol(g2i, LR3.T{+xHi, -1.0d0, 00.0d0}), 
        LR3x3.MapCol(g2i, LR3.T{-yLo, 00.0d0, +1.0d0}), 
        LR3x3.MapCol(g2i, LR3.T{+yHi, 00.0d0, -1.0d0}) 
      },
      outline = (outlineWidth > 0.0)
    DO

      PROCEDURE FindXRange(VAR xMin, xMax: LONG) =
        (*
          Computes the min and max grid "x" for the image rectangle *)
        BEGIN
          xMin := +Huge; xMax := -Huge;
          WITH
            corner = ARRAY [0..3] OF LR3.T{
              LR3x3.MapRow(LR3.T{1.0d0, xLo, yLo}, i2g),
              LR3x3.MapRow(LR3.T{1.0d0, xHi, yLo}, i2g),
              LR3x3.MapRow(LR3.T{1.0d0, xLo, yHi}, i2g),
              LR3x3.MapRow(LR3.T{1.0d0, xHi, yHi}, i2g)
            }
          DO
            FOR i := 0 TO 3 DO
              WITH x = corner[i][1]/corner[i][0] DO
                xMin := MIN(xMin, x);
                xMax := MAX(xMax, x);
              END
            END;
          END;
        END FindXRange;

      PROCEDURE FindYRange(x: LONG; VAR yMin, yMax: LONG) =
        (*
          Computes the min and max "y" along the vertical line
          at abscissa "x" *)
        BEGIN
          yMin := -Huge;
          yMax := +Huge;
          FOR i := 0 TO 3 DO
            WITH 
              sW = side[i][0], sX = side[i][1], sY = side[i][2]
            DO
              IF sY # 0.0d0 THEN
                WITH y = -(sW + sX*x)/sY DO 
                  IF sY > 0.0d0 THEN 
                    yMin := MAX(yMin, y) 
                  ELSE 
                    yMax := MIN(yMax, y)
                  END
                END
              END
            END
          END;
        END FindYRange;

      PROCEDURE InWindow(x, y: INT): BOOL =
        (*
          TRUE if tile "(x,y)" has center inside the plotting window. *)
        BEGIN
          WITH
            c = LR3.T{1.0d0, FLOAT(x, LONG) + 0.5d0, FLOAT(y, LONG) + 0.5d0}
          DO
            FOR i := 0 TO 3 DO 
              IF LR3.Dot(c, side[i]) < 0.0d0 THEN RETURN FALSE END
            END
          END;
          RETURN TRUE
        END InWindow;
      
      PROCEDURE Visible(x, y: INT): BOOL =
        (*
          TRUE if tile "(x,y)" is defined and is to be plotted. *)
        BEGIN
          IF defined # NIL AND NOT defined(x,y) THEN
            RETURN FALSE
          ELSE
            RETURN InWindow(x,y)
          END
        END Visible;
        
      PROCEDURE PlotTile(x, y: INT) =
        (* 
          Assumes "(x+.5,y+.5)" is inside plotting rectangle. *)
        BEGIN
          IF defined = NIL OR defined(x,y) THEN
            IF NOT InWindow(x, y) THEN
              Debug(x, y, yMin, yMax)
            END;
            ps.setLineColor(PSPlot.Invisible);
            ps.setFillColor(color(x, y));
            WITH
              gx = FLOAT(x,LONG),
              gy = FLOAT(y,LONG),
              a = LR3x3.MapRow(LR3.T{1.0d0, gx,       gy      }, g2i),
              b = LR3x3.MapRow(LR3.T{1.0d0, gx+1.0d0, gy      }, g2i),
              c = LR3x3.MapRow(LR3.T{1.0d0, gx      , gy+1.0d0}, g2i),
              d = LR3x3.MapRow(LR3.T{1.0d0, gx+1.0d0, gy+1.0d0}, g2i),
              
              ra = a[1]/a[0], sa = a[2]/a[0],
              rb = b[1]/b[0], sb = b[2]/b[0],
              rc = c[1]/c[0], sc = c[2]/c[0],
              rd = d[1]/d[0], sd = d[2]/d[0]
            DO
              ps.triangle(ra, sa, rb, sb, rd, sd);
              ps.triangle(ra, sa, rd, sd, rc, sc);
              IF outline THEN
                IF NOT Visible(x-1,y) THEN 
                  ps.setLineColor(PSPlot.Black);
                  ps.segment(ra, sa, rc, sc)
                END;
                IF NOT Visible(x,y+1) THEN 
                  ps.setLineColor(PSPlot.Black);
                  ps.segment(rc, sc, rd, sd)
                END;
                IF NOT Visible(x+1,y) THEN 
                  ps.setLineColor(PSPlot.Black);
                  ps.segment(rd, sd, rb, sb)
                END;
                IF NOT Visible(x,y-1) THEN 
                  ps.setLineColor(PSPlot.Black);
                  ps.segment(rb, sb, ra, sa)
                END;
              END
            END
          END
        END PlotTile;

      VAR xMin, xMax, yMin, yMax: LONG;

      PROCEDURE Debug(x, y: INTEGER; yMin, yMax: LONG) =
        <* FATAL Wr.Failure, Thread.Alerted *>
        BEGIN
          Wr.PutText(stderr, "x = " & Fmt.Int(x) & "\n");
          Wr.PutText(stderr, "y = " & Fmt.Int(y) & "\n");
          Wr.PutText(stderr, "yMin = " & Fmt.LongReal(yMin) & "\n");
          Wr.PutText(stderr, "yMax = " & Fmt.LongReal(yMax) & "\n");
          FOR i := 0 TO 3 DO 
            Wr.PutText(stderr, "side = [" & Fmt.LongReal(side[i][0]));
            Wr.PutText(stderr, ", " & Fmt.LongReal(side[i][1]));
            Wr.PutText(stderr, ", " & Fmt.LongReal(side[i][2]));
            Wr.PutText(stderr, "]\n");
          END;
        END Debug;
        
      BEGIN
        IF outline THEN 
          ps.setLineWidth(outlineWidth)
        ELSE
          ps.setLineWidth(0.0)
        END;
        FindXRange(xMin, xMax);
        FOR x := CEILING(xMin - 0.5d0) TO FLOOR(xMax - 0.5d0) DO
          FindYRange(FLOAT(x, LONG) + 0.5d0, yMin, yMax);
          FOR y := CEILING(yMin - 0.5d0) TO FLOOR(yMax - 0.5d0) DO
            PlotTile(x, y)
          END
        END;
      END
    END
  END Square;

BEGIN
END PSPlotMesh.