(* Last edited on 1999-04-19 03:29:54 by stolfi *)
MODULE PGMMakeKernel EXPORTS Main;

(*
  Writes various convolution kernels. 
*)

IMPORT Wr, Fmt, Math, Thread, Process, ParseParams;
IMPORT PGM, PGMImageFilter, PixelWeightTable;
FROM Stdio IMPORT stdout, stderr;

TYPE
  LONG = LONGREAL;
  NAT = CARDINAL;
  BOOL = BOOLEAN;
  Weights = PixelWeightTable.T;
  
  Options = RECORD
      maxVal: NAT;           (* Maximum output value *)
      weights: REF Weights;  (* Pixel weight table. *)
      (* Trigonometric factor: *)
      wave: BOOL;            (* TRUE applies a trigonometric wave. *)
      sin: BOOL;             (* TRUE applies a "sin" wave, FALSE a "cos" wave. *)
      fx, fy: LONG;          (* Frequency vector for multiplying wave. *)
      (* Quadratic factor: *)
      quadratic: BOOL;       (* TRUE to apply a quadratic factor. *)
      p: LONG;               (* Coefficient of constant term. *)
      px, py: LONG;          (* Coefficients of linear terms. *)
      pxx, pxy, pyy: LONG;   (* Coefficients of quadratic terms. *)
      noCenter: BOOL;        (* TRUE to force center weight to be zero. *)
      signed: BOOL;          (* TRUE forces normalization as signed. *)
      printWeights: BOOL;    (* TRUE to print pixel weight table. *)
    END;
    
CONST
  TwoPi = 6.28318530717958647688d0;
    
PROCEDURE Main() =
  <* FATAL Wr.Failure, Thread.Alerted *>
  VAR maxAbs, scale: LONG;
  VAR bias: NAT;
  VAR DXMax, DYMax: NAT;
  BEGIN
    WITH 
      o = GetOptions(),
      w = o.weights^,
      NY = NUMBER(w),
      NX = NUMBER(w[0]),
      HX = NX DIV 2,
      HY = NY DIV 2,
      ot = NEW(REF ARRAY OF ARRAY OF NAT, NY, NX)^
    DO
      
      (* Apply modifier waves if requested: *)
      IF o.wave THEN
        ApplyWave(w, o.sin, o.fx, o.fy)
      END;
      IF o.quadratic THEN
        ApplyQuadratic(w, o.p, o.px, o.py, o.pxx, o.pxy, o.pyy)
      END;
      
      (* Remove center pixel if requeted: *)
      IF o.noCenter THEN w[HY,HX] := 0.0d0 END;
      
      (* Scale so that the maximum absolute value is 1: *)
      maxAbs := 0.0d0;
      FOR iy := 0 TO NY-1 DO 
        FOR ix := 0 TO NX-1 DO
          maxAbs := MAX(maxAbs, ABS(w[iy,ix]))
        END
      END;
      IF maxAbs # 0.0d0 THEN
        FOR iy := 0 TO NY-1 DO 
          FOR ix := 0 TO NX-1 DO
            w[iy,ix] := w[iy,ix]/maxAbs
          END
        END
      END;

      (* Print weights: *)
      IF o.printWeights THEN 
        PixelWeightTable.Print(stderr, w)
      END;
      
      (* Compute scaling parameters *)
      IF o.signed THEN
        bias := o.maxVal DIV 2; scale := FLOAT(o.maxVal DIV 2, LONG)
      ELSE
        bias := 0; scale := FLOAT(o.maxVal, LONG)
      END;
      
      (* Convert to integer image, and find true dimensions: *)
      DXMax := 0; DYMax := 0;
      FOR yp := 0 TO NY-1 DO
        FOR xp := 0 TO NX-1 DO
          ot[yp,xp] := bias + ROUND(w[yp,xp] * scale);
          IF ot[yp,xp] # bias THEN
            WITH dx = ABS(xp - HX) DO IF dx > DXMax THEN DXMax := dx END END;
            WITH dy = ABS(yp - HY) DO IF dy > DYMax THEN DYMax := dy END END;
          END
        END;
      END;
      
      (* Write image *)
      WITH
        MX = 2*DXMax + 1,
        MY = 2*DYMax + 1,
        otPGM = NEW(PGM.Writer).init(stdout, MX, MY, o.maxVal, binary := o.maxVal <= 255)
      DO
        Wr.PutText(stderr, "Actual image size = " & Fmt.Int(MX) & "×" & Fmt.Int(MY) & "\n");
        FOR yp := HY - DYMax TO HY + DYMax DO
          WITH otRow = SUBARRAY(ot[yp], HX - DXMax, MX) DO
            PGMImageFilter.WritePixels(otPGM, otRow,  mkWr := NIL, noPixel := LAST(NAT))
          END
        END;
        otPGM.finish()
      END
    END
  END Main;
  
PROCEDURE ApplyWave(
    VAR w: Weights;
    sin: BOOL;
    fx, fy: LONG; 
  ) =
  (*
    Multiplies the weights by a trigonometric wave of given
    frequency vector:
  *)
  VAR a: LONG;
  BEGIN
    <* ASSERT fx # 0.0d0 OR fy # 0.0d0 *>
    WITH
      NY = NUMBER(w),
      NX = NUMBER(w[0]),
      HX = NX DIV 2,
      HY = NY DIV 2
    DO
      FOR iy := 0 TO NY-1 DO 
        FOR ix := 0 TO NX-1 DO
          WITH
            x = FLOAT(ix - HX, LONG),
            y = FLOAT(iy - HY, LONG),
            r = x*fx + y*fy
          DO
            IF sin THEN 
              a := Math.sin(TwoPi*r)
            ELSE
              a := Math.cos(TwoPi*r)
            END;
            w[iy,ix] := a * w[iy,ix]
          END
        END
      END
    END
  END ApplyWave;
  
PROCEDURE ApplyQuadratic(
    VAR w: Weights;
    p, px, py, pxx, pxy, pyy: LONG; 
  ) =
  (*
    Multiplies the weights by a quadratic bivariate 
    polynomial with given coefficients.
  *)
  BEGIN
    WITH
      NY = NUMBER(w),
      NX = NUMBER(w[0]),
      HX = NX DIV 2,
      HY = NY DIV 2
    DO
      FOR iy := 0 TO NY-1 DO 
        FOR ix := 0 TO NX-1 DO
          WITH
            x = FLOAT(ix - HX, LONG),
            y = FLOAT(iy - HY, LONG),
            a = p + px*x + py*y + pxx * x*x + pxy*x*y + pyy*y*y
          DO
            w[iy,ix] := a * w[iy,ix]
          END
        END
      END
    END
  END ApplyQuadratic;
  
PROCEDURE GetOptions (): Options =
  <* FATAL Thread.Alerted, Wr.Failure *>
  CONST
    MaxMaxVal: NAT = 256*256-1;
  VAR o: Options;
  BEGIN
    WITH
      pp = NEW(ParseParams.T).init(stderr)
    DO
      TRY
        IF pp.keywordPresent("-maxVal") THEN
          o.maxVal := pp.getNextInt(1, MaxMaxVal)
        ELSE
          o.maxVal := MaxMaxVal
        END;
        
        pp.getKeyword("-type");
        o.weights := PixelWeightTable.FromParams(pp);
        
        o.signed := pp.keywordPresent("-signed");

        o.wave := FALSE;
        o.sin := FALSE;
        o.quadratic := FALSE;
        
        IF pp.keywordPresent("-modulate") THEN

          IF pp.testNext("sin") THEN
            o.wave := TRUE;
            o.sin := TRUE;
            o.fx := pp.getNextLongReal(); 
            o.fy := pp.getNextLongReal(); 
            o.signed := TRUE;
          ELSIF pp.testNext("cos") THEN
            o.wave := TRUE;
            o.sin := FALSE;
            o.fx := pp.getNextLongReal(); 
            o.fy := pp.getNextLongReal(); 
            o.signed := TRUE;
          ELSIF pp.testNext("quadratic") THEN
            o.quadratic := TRUE;
            o.p   := pp.getNextLongReal(); 
            o.px  := pp.getNextLongReal(); 
            o.py  := pp.getNextLongReal(); 
            o.pxx := pp.getNextLongReal(); 
            o.pxy := pp.getNextLongReal(); 
            o.pyy := pp.getNextLongReal(); 
            o.signed := o.signed OR (
              o.p < 0.0d0 OR
              o.px # 0.0d0 OR o.py # 0.0d0 OR 
              o.pxx < 0.0d0 OR o.pxy < 0.0d0 OR o.pyy < 0.0d0 
            );
          END
          
        END;
        
        o.noCenter := pp.keywordPresent("-noCenter");
        
        o.printWeights := pp.keywordPresent("-printWeights");
        
        pp.finish();                                       
      EXCEPT                                                            
      | ParseParams.Error =>                                              
          Wr.PutText(stderr, "Usage: PGMMakeKernel \\\n");
          Wr.PutText(stderr, "  [ -maxVal NUM ]");
          Wr.PutText(stderr, "  -type \\\n");
          Wr.PutText(stderr, PixelWeightTable.ParamsHelp & "\\\n");
          Wr.PutText(stderr, "  [ -modulate \\\n");
          Wr.PutText(stderr, "      { {sin|cos} FX FY } | \\\n");
          Wr.PutText(stderr, "      { quadratic P PX PY PXX PXY PYY } \\\n");
          Wr.PutText(stderr, "  ] \\\n");
          Wr.PutText(stderr, "  [ -noCenter ] \\\n");
          Wr.PutText(stderr, "  [ -signed ] \\\n");
          Wr.PutText(stderr, "  [ -printWeights ]\n");
          Process.Exit(1);
      END;
    END;
    RETURN o
  END GetOptions;

BEGIN 
  Main();
END PGMMakeKernel.