MODULE MakeTexture2D EXPORTS Main;

IMPORT Texture2DGray AS Tx;
IMPORT ParseParams, Process, Wr, Thread, Text, Fmt, FileWr, OSError;
IMPORT Random, Math;
FROM Stdio IMPORT stderr;

CONST
  Debug = FALSE;

TYPE
  LONG = LONGREAL;
  LONGS = ARRAY OF LONG;

  BOOL = BOOLEAN;
  NAT = CARDINAL;
  INT = INTEGER;

  PixelValue = LONG;
  PixelValues = ARRAY OF PixelValue;

TYPE
  Field = OBJECT
      NL: NAT;           (* Number of components (layers) *)
      pixMin: REF LONGS; (* The nominal min pixel value per layer *)
      pixMax: REF LONGS; (* The nominal max pixel value per layer *)
    METHODS
      print(wr: Wr.T);
        (*
          Prints a description of the field to "wr". *)

      eval(x, y: NAT; VAR pix: PixelValues);
        (*
          Computes the texture value "pix[i]" at "(x,y)" for each
          layer "i". *)
    END;

  Options = RECORD
      random: REAL;     (* Amount of random noise to add to input *)
      inName: TEXT;     (* Input texture file name, or "" if zero texture. *)
      outName: TEXT;    (* Output filename prefix. *)
      UX, VY: NAT;      (* Texture size *)
      VX: NAT;          (* Brick-row slant *)
      NX, NY: NAT;      (* Output replication factors *)
      field: Field;     (* The texture formula *)
    END;

CONST
  Pi = 3.141592653589793238d0;
  TwoPi = 2.0d0 * Pi;
  HalfPi = 0.5d0 * Pi;
  CompName = ARRAY OF TEXT{"A", "B", "C", "D", "E", "F", "G", "H"};

PROCEDURE Main() =
  <* FATAL Thread.Alerted, Wr.Failure, OSError.E *>
  BEGIN
    WITH
      o = GetOptions(),
      NL = o.field.NL,
      A = NEW(REF ARRAY OF Tx.T, o.field.NL)^,
      rnd = NEW(Random.Default).init(fixed := TRUE),
      zeroPixel = Tx.LongRealPixel{0.0d0}
    DO
      (* Allocate working textures: *)
      FOR i := 0 TO NL-1 DO
        IF Text.Equal(o.inName, "") THEN
          A[i] := Tx.New(o.UX, o.VX, o.VY, zeroPixel, zeroPixel);
          Tx.Zero(A[i]);
        ELSE
          WITH
            fName = o.inName & "-" & CompName[i] & ".pgm",
            a = Tx.LongRealPixel{0.5d0 * (o.field.pixMin[i] + o.field.pixMax[i])},
            s = Tx.LongRealPixel{0.5d0 * (o.field.pixMax[i] - o.field.pixMin[i])}
          DO
            A[i] := Tx.Read(fName, o.VX, a, s);
            <* ASSERT A[i].UX = o.UX *>
            <* ASSERT A[i].VX = o.VX *>
            <* ASSERT A[i].VY = o.VY *>
          END
        END;
        IF o.random # 0.0 THEN
          Tx.URandom(A[i], rnd,
            min := Tx.LongRealPixel{FLOAT(o.random, LONG) * o.field.pixMin[i]},
            max := Tx.LongRealPixel{FLOAT(o.random, LONG) * o.field.pixMax[i]},
            clear := FALSE
          )
        END;
      END;

      MapTexture(A, o.field);

      (* Write texture: *)
      WITH
        fName = o.outName,
        wr = FileWr.Open(fName & ".txt")
      DO
        o.field.print(wr);
        Wr.PutText(wr, "\n");
        Wr.Flush(wr);
        FOR i := 0 TO NL-1 DO
          WITH
            pixMin = o.field.pixMin[i],
            pixMax = o.field.pixMax[i]
          DO
            Wr.PutText(wr, CompName[i] & " range = [");
            Wr.PutText(wr, FLR(pixMin, 4) & " _ " & FLR(pixMax, 4) & "]");
            Wr.Flush(wr);
            Tx.Write(
              name := fName & "-" & CompName[i] & ".pgm",
              t := A[i],
              a := Tx.LongRealPixel{0.5d0 * (pixMin + pixMax)},
              s := Tx.LongRealPixel{0.5d0 * (pixMax - pixMin)},
              NX := o.NX, NY := o.NY
            );
          END
        END
      END;
    END
  END Main;

PROCEDURE MapTexture(READONLY A: ARRAY OF Tx.T; field: Field) =
  (*
    Modifies "A" as determiend by the given "field"..
  *)
  BEGIN
    <* ASSERT NUMBER(A) = field.NL *>
    IF field.NL = 0 THEN RETURN END;
    WITH
      UX = A[0].UX,
      VY = A[0].VY,
      NL = field.NL,
      pix = NEW(REF PixelValues, NL)^
    DO
      FOR i := 0 TO NL-1 DO
        <* ASSERT A[i].UX = UX *>
        <* ASSERT A[i].VY = VY *>
      END;
      FOR y := 0 TO VY-1 DO
        FOR x := 0 TO UX-1 DO
          FOR i := 0 TO NL-1 DO pix[i] := FLOAT(A[i].st[y,x], LONG) END;
          field.eval(x, y, pix);
          FOR i := 0 TO NL-1 DO A[i].st[y,x] := FLOAT(pix[i]) END
        END
      END
    END;
  END MapTexture;

PROCEDURE GetOptions (): Options =
  <* FATAL Thread.Alerted, Wr.Failure *>
  VAR o: Options;
  BEGIN
    WITH pp = NEW(ParseParams.T).init(stderr) DO
      TRY

        IF pp.keywordPresent("-random") THEN
          o.random := pp.getNextReal(0.0, 1000.0)
        ELSE
          o.random := 0.0
        END;

        IF pp.keywordPresent("-inName") THEN
          o.inName := pp.getNext();
        ELSE
          o.inName := "";
        END;

        pp.getKeyword("-outName");
        o.outName := pp.getNext();

        pp.getKeyword("-size");
        o.UX := pp.getNextInt(1, 4096);
        o.VY := pp.getNextInt(1, 4096);

        IF pp.keywordPresent("-offset") THEN
          o.VX := pp.getNextInt(0, o.UX-1)
        ELSE
          o.VX := 0
        END;
        
        IF pp.keywordPresent("-repeat") THEN
          o.NX := pp.getNextInt(1, 100);
          o.NY := pp.getNextInt(1, 100);
        ELSE
          o.NX := 1;
          o.NY := 1
        END;

        pp.getKeyword("-field");
        WITH sysName = pp.getNext() DO
          IF Text.Equal(sysName, "Wave") THEN
            WITH
              fX = pp.getNextInt(-10000, +10000),
              fY = pp.getNextInt(-10000, +10000)
            DO
              o.field := NEW(WaveField,
                NL := 2,
                pixMin := NewLongs(LONGS{-1.5d0, -1.5d0}),
                pixMax := NewLongs(LONGS{+1.5d0, +1.5d0}),
                fX := fX,
                fY := fY,
                wX := TwoPi * FLOAT(fX, LONG)/FLOAT(o.UX, LONG),
                wY := TwoPi * FLOAT(fY, LONG)/FLOAT(o.VY, LONG),
                dXdY := FLOAT(o.VX, LONG)/FLOAT(o.VY, LONG)
              )
            END
          ELSIF Text.Equal(sysName, "Swirl") THEN
            o.field := NEW(SwirlField,
              NL := 2,
              pixMin := NewLongs(LONGS{-1.5d0, -1.5d0}),
              pixMax := NewLongs(LONGS{+1.5d0, +1.5d0}),
              rX := FLOAT(o.UX-1, LONG)/2.0d0,
              rY := FLOAT(o.VY-1, LONG)/2.0d0
            )
          ELSIF Text.Equal(sysName, "Const") THEN
            WITH
              c = pp.getNextLongReal(-1000.0d0, 1000.0d0)
            DO
              o.field := NEW(ConstField,
                NL := 1,
                pixMin := NewLongs(LONGS{MIN(c, 0.0d0) - 1.5d0}),
                pixMax := NewLongs(LONGS{MAX(c, 0.0d0) + 1.5d0}),
                c := c
              )
            END
          ELSIF Text.Equal(sysName, "Squares") THEN
            WITH
              nLayers = pp.getNextInt(1, 10000),
              rMax = pp.getNextLongReal(0.5d0, 10000.0d0),
              theta = pp.getNextLongReal(-HalfPi, +HalfPi)
            DO
              o.field := NEW(SquaresField,
                NL := 1,
                pixMin := NewLongs(LONGS{- 1.5d0}),
                pixMax := NewLongs(LONGS{+ 1.5d0}),
                nLayers := nLayers,
                rMax := rMax,
                theta := theta,
                uX := Math.cos(theta),
                uY := Math.sin(theta),
                rX := FLOAT(o.UX-1, LONG)/2.0d0,
                rY := FLOAT(o.VY-1, LONG)/2.0d0
              )
            END
          ELSIF Text.Equal(sysName, "Rings") THEN
            WITH
              r1 = pp.getNextLongReal(0.0d0, 10000.0d0),
              r2 = pp.getNextLongReal(r1 + 1.0d0, 10000.0d0),
              half = pp.testNext("half"),
              rMax = MAX(r1, r2),
              rMin = MIN(r1, r2)
            DO
              o.field := NEW(RingsField,
                NL := 1,
                pixMin := NewLongs(LONGS{- 1.5d0}),
                pixMax := NewLongs(LONGS{+ 1.5d0}),
                rMax := rMax,
                rMin := rMin,
                rX := FLOAT(o.UX-1, LONG)/2.0d0,
                rY := FLOAT(o.VY-1, LONG)/2.0d0,
                half := half
              )
            END
          ELSIF Text.Equal(sysName, "Ripple") THEN
            WITH
              rMax = pp.getNextLongReal(0.5d0, 10000.0d0),
              alpha = pp.getNextLongReal(0.01d0, 3.0d0),
              ratio = 1.0d0 + alpha
            DO
              VAR rMin: LONG := rMax;
              BEGIN
                REPEAT rMin := rMin/ratio UNTIL alpha*rMin/ratio < 2.0d0;
                o.field := NEW(RippleField,
                  NL := 1,
                  pixMin := NewLongs(LONGS{- 1.5d0}),
                  pixMax := NewLongs(LONGS{+ 1.5d0}),
                  rMax := rMax,
                  rMin := rMin,
                  alpha := alpha,
                  coef := Pi / Math.log(1.0d0 + alpha),
                  rX := FLOAT(o.UX-1, LONG)/2.0d0,
                  rY := FLOAT(o.VY-1, LONG)/2.0d0
                )
              END
            END
          ELSIF Text.Equal(sysName, "FTest") THEN
            IF o.UX # 120 OR o.VY # 120 THEN 
              pp.error("the \"FTest\" image must be 120x120")
            END;
            o.field := NEW(FTestField,
              NL := 1,
              pixMin := NewLongs(LONGS{-1.5d0, -1.5d0}),
              pixMax := NewLongs(LONGS{+1.5d0, +1.5d0})
            )
          ELSE
            pp.error("invalid field name = " & sysName)
          END
        END;

        pp.finish();
      EXCEPT
      | ParseParams.Error =>
          Wr.PutText(stderr, "Usage: MakeTexture2D\\\n");
          Wr.PutText(stderr, "  [ -inName <name> ] [ -random <amt> ] \\\n");
          Wr.PutText(stderr, "  -outName <name> [ -repeat nx ny ] \\\n");
          Wr.PutText(stderr, "  -size <width> <height> [ -offset <offset> ] \\\n");
          Wr.PutText(stderr, "  -field \\\n");
          Wr.PutText(stderr, "    { Const <val> \\\n");
          Wr.PutText(stderr, "    | Wave <fx> <fy> \\\n");
          Wr.PutText(stderr, "    | Swirl \\\n");
          Wr.PutText(stderr, "    | Squares <nLayers> <rMax> <theta> \\\n");
          Wr.PutText(stderr, "    | Ripple <rMax> <alpha> \\\n");
          Wr.PutText(stderr, "    | Rings <rMax> <rMin> [ half ] \\\n");
          Wr.PutText(stderr, "    } \n");
          Process.Exit (1);
      END;
    END;
    RETURN o
  END GetOptions;

PROCEDURE FLR(x: LONG; prec: NAT := 6): TEXT =
  BEGIN
    RETURN Fmt.Pad(Fmt.LongReal(x, prec := prec, style := Fmt.Style.Fix), prec + 6)
  END FLR;

<*UNUSED*>
PROCEDURE FR(x: REAL; prec: NAT := 6): TEXT =
  BEGIN
    RETURN Fmt.Pad(Fmt.Real(x, prec := prec, style := Fmt.Style.Fix), prec + 6)
  END FR;

PROCEDURE NewLongs(READONLY a: LONGS): REF LONGS =
  (*
    Modula-3 is so silly....
  *)
  BEGIN
    WITH n = NEW(REF LONGS, NUMBER(a)) DO n^ := a; RETURN n END
  END NewLongs;

(* INTERESTING EQUATION FIELDS *)

TYPE
  WaveField = Field OBJECT
      fX, fY: NAT;   (* frequency *)
      wX, wY: LONG;  (* coefficients *)
      dXdY: LONG;    (* brick-row slant *)
    OVERRIDES
      print := WavePrint;
      eval := WaveEval
    END;

PROCEDURE WavePrint(s: WaveField; wr: Wr.T) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    Wr.PutText(wr, "Wave: sin/cos periodic wave\n");
    Wr.PutText(wr, "  A(x,y) = exp(TwoPi * I * (fX/UX * (x - dXdY*y) + fY/VY * y)\n");
    Wr.PutText(wr, "  fX =   " & Fmt.Pad(Fmt.Int(s.fX), 8) & "\n");
    Wr.PutText(wr, "  fY =   " & Fmt.Pad(Fmt.Int(s.fY), 8) & "\n");
    Wr.PutText(wr, "  dXdY = " & FLR(s.dXdY, 6) & "\n");
    Wr.Flush(wr)
  END WavePrint;

PROCEDURE WaveEval(
    s: WaveField;
    x, y: NAT;
    VAR pix: PixelValues
  ) =
  BEGIN
    WITH
      xx = FLOAT(x, LONG),
      yy = FLOAT(y, LONG),
      arg = s.wX * (xx - s.dXdY * yy) + s.wY * yy,
      c = Math.cos(arg),
      s = Math.sin(arg)
    DO
      pix[0] := c;
      pix[1] := s;
    END
  END WaveEval;

TYPE
  ConstField = Field OBJECT
      c: LONG;
    OVERRIDES
      print := ConstPrint;
      eval := ConstEval
    END;

PROCEDURE ConstPrint(s: ConstField; wr: Wr.T) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    Wr.PutText(wr, "Const: constant field\n");
    Wr.PutText(wr, "  A(x,y) = " & FLR(s.c, 6) & "\n");
    Wr.Flush(wr)
  END ConstPrint;

PROCEDURE ConstEval(
    s: ConstField;
    <*UNUSED*> x, y: NAT;
    VAR pix: PixelValues
  ) =
  BEGIN
    pix[0] := s.c;
  END ConstEval;

TYPE
  SwirlField = Field OBJECT
      rX, rY: LONG;
    OVERRIDES
      print := SwirlPrint;
      eval := SwirlEval
    END;

PROCEDURE SwirlPrint(s: SwirlField; wr: Wr.T) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    Wr.PutText(wr, "Swirl: swirling vector field\n");
    Wr.PutText(wr, "  A(x,y) = w((x-rX)/rX, (y-rY)/rY)*rot(dir(x-rX, y-rY))\n");
    Wr.PutText(wr, "  w(u,v) = (1 - u^2)^2 * (1 - v^2)^2\n");
    Wr.PutText(wr, "  dir(u,v) = (u,v)/sqrt(u^2+v^2)\n");
    Wr.PutText(wr, "  rot(u,v) = (-v,u)\n");
    Wr.PutText(wr, "  rX = " & FLR(s.rX, 6) & "\n");
    Wr.PutText(wr, "  rY = " & FLR(s.rY, 6) & "\n");
    Wr.Flush(wr)
  END SwirlPrint;

PROCEDURE SwirlEval(
    s: SwirlField;
    x, y: NAT;
    VAR pix: PixelValues
  ) =
  <* FATAL Thread.Alerted, Wr.Failure *>
  BEGIN
    WITH
      uX = (FLOAT(x, LONG) - s.rX),
      uY = (FLOAT(y, LONG) - s.rY),
      sX = uX/s.rX,
      sY = uY/s.rY,
      w = (1.0d0 - sX*sX)*(1.0d0 - sY*sY),
      L = Math.sqrt(uX*uX + uY*uY),
      m = w * w / MAX(1.0d-20, L)
    DO
      IF Debug THEN
        Wr.PutText(stderr, "x = " & FLR(uX, 3) & "\n");
        Wr.PutText(stderr, "y = " & FLR(uY, 3) & "\n");
        Wr.PutText(stderr, "L = " & FLR(L, 6) & "\n");
        Wr.PutText(stderr, "w = " & FLR(w, 6) & "\n");
        Wr.PutText(stderr, "m = " & FLR(m, 6) & "\n");
        Wr.PutText(stderr, "\n");
      END;
      pix[0] := - m * uY;
      pix[1] := + m * uX;
    END
  END SwirlEval;

TYPE
  SquaresField = Field OBJECT
      nLayers: NAT;
      rMax: LONG;
      theta: LONG;
      uX, uY: LONG; 
      rX, rY: LONG;
    OVERRIDES
      print := SquaresPrint;
      eval := SquaresEval
    END;

PROCEDURE SquaresPrint(s: SquaresField; wr: Wr.T) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    Wr.PutText(wr, "Squares: nested squares\n");
    Wr.PutText(wr, "  \"nLayers\" nested squares, radius \"rMax\", tilted by \"theta\"\n");
    Wr.PutText(wr, "  nLayers = " & Fmt.Pad(Fmt.Int(s.nLayers), 6) & "\n");
    Wr.PutText(wr, "  rMax =    " & FLR(s.rMax, 6) & "\n");
    Wr.PutText(wr, "  theta =   " & FLR(s.theta, 6) & "\n");
    Wr.Flush(wr)
  END SquaresPrint;

PROCEDURE SquaresEval(
    s: SquaresField;
    x, y: NAT;
    VAR pix: PixelValues
  ) =
  BEGIN
    WITH
      fX = FLOAT(x, LONG) - s.rX,
      fY = FLOAT(y, LONG) - s.rY,
      fU = s.uX*fX + s.uY*fY,
      fV = s.uY*fX - s.uX*fY,
      r = MAX(ABS(fU), ABS(fV)),
      relArea = 0.999999d0 - 0.999998d0 * (r*r)/(s.rMax*s.rMax),
      layer = 1 + FLOOR(FLOAT(s.nLayers, LONG)*relArea)
    DO
      IF r > s.rMax THEN 
        pix[0] := -1.0d0
      ELSE 
        pix[0] := 2.0d0 * FLOAT(layer MOD 2, LONG) - 1.0d0
      END;
    END
  END SquaresEval;

TYPE
  RingsField = Field OBJECT
      rMax, rMin: LONG;
      rX, rY: LONG;
      half: BOOL;
    OVERRIDES
      print := RingsPrint;
      eval := RingsEval
    END;

PROCEDURE RingsPrint(s: RingsField; wr: Wr.T) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    Wr.PutText(wr, "Rings: round rings\n");
    Wr.PutText(wr, "  outer radius \"rMax\", inner radius \"rMin\"\n");
    Wr.PutText(wr, "  rMax = " & FLR(s.rMax, 6) & "\n");
    Wr.PutText(wr, "  rMin = " & FLR(s.rMin, 6) & "\n");
    IF s.half THEN Wr.PutText(wr, "  (cut in half)\n") END;
    Wr.Flush(wr)
  END RingsPrint;

PROCEDURE RingsEval(
    s: RingsField;
    x, y: NAT;
    VAR pix: PixelValues
  ) =
  BEGIN
    WITH
      fX = FLOAT(x, LONG) - s.rX,
      fY = FLOAT(y, LONG) - s.rY,
      r2 = fX*fX + fY*fY
    DO
      IF r2 > s.rMax*s.rMax OR r2 < s.rMin*s.rMin 
      OR (s.half AND fX + fY < 0.0d0) THEN 
        pix[0] := -1.0d0 
      ELSE 
        pix[0] := +1.0d0
      END;
    END
  END RingsEval;

TYPE
  FTestField = Field OBJECT
    OVERRIDES
      print := FTestPrint;
      eval := FTestEval
    END;

PROCEDURE FTestPrint(<*UNUSED*> s: FTestField; wr: Wr.T) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    Wr.PutText(wr, "FTest: test pattern for local feature mappers.\n");
    Wr.Flush(wr)
  END FTestPrint;

PROCEDURE FTestEval(
    <*UNUSED*> s: FTestField;
    x, y: NAT;
    VAR pix: PixelValues
  ) =
  VAR c, i: NAT;
  BEGIN
    WITH
      row = y DIV 5,
      col = x DIV 5,
      case = 24 * row + col,
      u = x MOD 5 - 2,
      v = y MOD 5 - 2
    DO
      IF row >= 24 OR col >= 24 
      OR case >= 512 
      OR ABS(u) > 1 OR ABS(v) > 1 THEN 
        pix[0] := 0.0d0
      ELSE
        c := case;
        i := (3 * u + v) MOD 9;
        WHILE i > 0 DO c := c DIV 2; DEC(i) END;
        IF (c MOD 2) = 0 THEN 
          pix[0] := 0.0d0
        ELSE
          pix[0] := 1.0d0
        END
      END
    END
  END FTestEval;

TYPE
  RippleField = Field OBJECT
      rMax: LONG;
      rMin: LONG;
      alpha: LONG;
      coef: LONG;
      rX, rY: LONG;
    OVERRIDES
      print := RipplePrint;
      eval := RippleEval
    END;

PROCEDURE RipplePrint(s: RippleField; wr: Wr.T) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    Wr.PutText(wr, "Ripple: exponential ripples\n");
    Wr.PutText(wr, "  a circular ripple with exponential wavelength,\n");
    Wr.PutText(wr, "  radius \"rMax..rMin\", wavelength growth by \"alpha\"/cycle\n");
    Wr.PutText(wr, "  rMax =    " & FLR(s.rMax, 6) & "\n");
    Wr.PutText(wr, "  alpha =   " & FLR(s.alpha, 6) & "\n");
    Wr.Flush(wr)
  END RipplePrint;

PROCEDURE RippleEval(
    s: RippleField;
    x, y: NAT;
    VAR pix: PixelValues
  ) =
  BEGIN
    WITH
      fX = FLOAT(x, LONG) - s.rX,
      fY = FLOAT(y, LONG) - s.rY,
      r2 = fX*fX + fY*fY,
      m2 = s.rMax * s.rMax,
      u2 = s.rMin * s.rMin
    DO
      IF r2 >= m2 OR r2 <= u2 THEN 
        pix[0] := 0.0d0
      ELSE 
        pix[0] := Math.sin(s.coef * Math.log(r2/m2))
      END;
    END
  END RippleEval;


(*
  Equation Field Template:

TYPE
  ~Field = Field OBJECT
      ~: LONG;
      ~: LONG;
    OVERRIDES
      print := ~Print;
      eval := ~Eval
    END;

PROCEDURE ~Print(s: ~Field; wr: Wr.T) =
  <* FATAL Wr.Failure, Thread.Alerted *>
  BEGIN
    Wr.PutText(wr, "~:\n");
    Wr.PutText(wr, "  ~\n");
    Wr.PutText(wr, "  ~\n");
    Wr.PutText(wr, "  ~\n");
    Wr.PutText(wr, "  ~ = " & FLR(s.~, 6) & "\n");
    Wr.PutText(wr, "  ~ = " & FLR(s.~, 6) & "\n");
    Wr.Flush(wr)
  END ~Print;

PROCEDURE ~Eval(
    s: ~Field;
    x, y: NAT;
    VAR pix: PixelValues
  ) =
  BEGIN
    WITH
      ~
    DO
      pix[0] := ~;
      pix[1] := ~;
    END
  END ~Eval;

*)

BEGIN
  Main();
END MakeTexture2D.