MODULE PixelWeightTable;
(* Last edited on 1998-10-25 23:16:25 by stolfi *)

IMPORT LongReal;
FROM Math IMPORT sqrt, exp, log, acos;
IMPORT Text, Fmt, Wr, Thread, ParseParams;

CONST
  Tiny = LongReal.MinPos;
  Huge = LongReal.MaxFinite;
  Pi = 3.14159265358979324d0;

TYPE
  WeightFunction = PROCEDURE(x, y: LONG): LONG;
    (*
      Evaluates the weight function for the
      pixel whose center has coordinates "(x,y)". *)
      
PROCEDURE MakeSymmetricTable(
    func: WeightFunction; 
    NX, NY: CARDINAL; 
    uniDim: BOOL;
  ): REF T =
  (*
    
    If "uniDim" is false, builds table with "NX" columns and "NY" rows,
    whose entries are given by the x- and y-symmetric function "func"
    applied to the coordinates of the center of each pixel.
    
    The values of "func" are shifted by "bias" towards zero, and
    those less than "bias" in absolute value are set to zero.
    
    If "uniDim" is true, builds the middle row of the table.
  *)
  BEGIN
    IF uniDim THEN NY := 1 END;
    <* ASSERT NX MOD 2 = 1 *>
    <* ASSERT NY MOD 2 = 1 *>
    WITH
      HX = NX DIV 2,
      HY = NY DIV 2,
      wRef = NEW(REF T, NY, NX),
      w = wRef^
    DO
      FOR r := 0 TO MAX(HX,HY) DO
        FOR s := 0 TO MIN(MIN(HX,HY), r) DO
          WITH 
            x = FLOAT(r, LONG),
            y = FLOAT(s, LONG),
            e = func(x, y)
          DO
            IF r <= HY THEN 
              w[HY+r, HX+s] := e;
              w[HY-r, HX+s] := e;
              w[HY+r, HX-s] := e;
              w[HY-r, HX-s] := e
            END;
                       
            IF r <= HX THEN 
              w[HY+s, HX+r] := e;
              w[HY+s, HX-r] := e;
              w[HY-s, HX+r] := e;
              w[HY-s, HX-r] := e
            END;
          END
        END
      END;
      RETURN wRef
    END
  END MakeSymmetricTable;

PROCEDURE NormalizeSum(VAR w: T) =
  VAR sum: LONG := 0.0d0;
  BEGIN
    FOR yw := 0 TO LAST(w) DO
      WITH wy = w[yw] DO
        FOR xw := 0 TO LAST(wy) DO
          WITH e = wy[xw] DO sum := sum + e END
        END
      END
    END;
    FOR yw := 0 TO LAST(w) DO
      WITH wy = w[yw] DO
        FOR xw := 0 TO LAST(wy) DO
          WITH e = wy[xw] DO e := e/sum END
        END
      END
    END
  END NormalizeSum;

PROCEDURE NormalizeSumOfSquares(VAR w: T) =
  VAR sum: LONG := 0.0d0;
  BEGIN
    FOR y := 0 TO LAST(w) DO 
      WITH wy = w[y] DO
        FOR x := 0 TO LAST(wy) DO
          WITH wxy = wy[x] DO
            sum := sum + wxy*wxy
          END
        END
      END
    END;
    <* ASSERT sum > 0.0d0 *>
    WITH coef = 1.0d0/sqrt(sum) DO
      FOR y := 0 TO LAST(w) DO 
        WITH wy = w[y] DO
          FOR x := 0 TO LAST(wy) DO
            wy[x] := coef * wy[x]
          END
        END
      END
    END
  END NormalizeSumOfSquares;
  
PROCEDURE NormalizeSumOfAbs(VAR w: T) =
  VAR sum: LONG := 0.0d0;
  BEGIN
    FOR y := 0 TO LAST(w) DO 
      WITH wy = w[y] DO
        FOR x := 0 TO LAST(wy) DO
          WITH wxy = wy[x] DO
            sum := sum + ABS(wxy)
          END
        END
      END
    END;
    <* ASSERT sum > 0.0d0 *>
    WITH coef = 1.0d0/sum DO
      FOR y := 0 TO LAST(w) DO 
        WITH wy = w[y] DO
          FOR x := 0 TO LAST(wy) DO
            wy[x] := coef * wy[x]
          END
        END
      END
    END
  END NormalizeSumOfAbs;
  
PROCEDURE NormalizeSignedSum(VAR w: T) =
  VAR sumP, sumN: LONG := 0.0d0;
  BEGIN
    FOR y := 0 TO LAST(w) DO 
      WITH wy = w[y] DO
        FOR x := 0 TO LAST(wy) DO
          WITH wxy = wy[x] DO
            IF wxy >= 0.0d0 THEN sumP := sumP + wxy ELSE sumN := sumN + wxy END
          END
        END
      END
    END;
    <* ASSERT sumP > 0.0d0 OR sumN < 0.0d0 *>
    WITH coef = 1.0d0/MAX(sumP, -sumN) DO
      FOR y := 0 TO LAST(w) DO 
        WITH wy = w[y] DO
          FOR x := 0 TO LAST(wy) DO
            wy[x] := coef * wy[x]
          END
        END
      END
    END
  END NormalizeSignedSum;
  
PROCEDURE Balance(VAR w: T; READONLY g: T) =
  (*
    Adds to "w" a multiple of "g" such that the resulting table has
    zero sum.  This is the same as making "w/g" orthogonal to the
    constant image in the "g"-weighted sum metric. *)
  VAR wSum, gSum: LONG := 0.0d0;
  BEGIN
    FOR y := 0 TO LAST(w) DO 
      WITH wy = w[y], gy = g[y] DO
        FOR x := 0 TO LAST(wy) DO
          WITH wxy = wy[x], gxy = gy[x] DO
            wSum := wSum + wxy; gSum := gSum + gxy
          END
        END
      END
    END;
    <* ASSERT gSum > 0.0d0 *>
    WITH coef = wSum / gSum DO
      FOR y := 0 TO LAST(w) DO 
        WITH wy = w[y], gy = g[y] DO
          FOR x := 0 TO LAST(wy) DO
            WITH wxy = wy[x], gxy = gy[x] DO
              wxy := wxy - coef * gxy
            END
          END
        END
      END
    END;
  END Balance;
  
PROCEDURE Flat(R: LONG; uniDim: BOOL := FALSE): REF T =
  BEGIN
    WITH 
      H = CEILING(R),
      N = 2*H + 1,
      R2 = R*R
    DO

      PROCEDURE FlatEval(x, y: LONG): LONG =
        BEGIN
          WITH
            xu = MAX(0.0d0, ABS(x)-0.5d0),
            yu = MAX(0.0d0, ABS(y)-0.5d0),
            r2u = xu*xu+yu*yu,
            xv = ABS(x) + 0.5d0,
            yv = ABS(y) + 0.5d0,
            r2v = xv*xv+yv*yv
          DO
            IF r2u > R2 THEN
              RETURN 0.0d0
            ELSIF r2v < R2 THEN
              RETURN 1.0d0
            ELSE
              (* Rough antialising: linear interpolation... *)
              WITH s = (R2 - r2u)/(r2v - r2u), t = 1.0d0 - s DO
                IF s < 0.5d0 THEN 
                  RETURN 2.0d0 * s*s
                ELSE
                  RETURN 1.0d0 - 2.0d0 * t*t
                END
              END
            END
          END
        END FlatEval;

      BEGIN
        RETURN MakeSymmetricTable(FlatEval, N, N, uniDim)
      END
    END
  END Flat;

PROCEDURE Tent(R: LONG; uniDim: BOOL := FALSE): REF T =
  BEGIN
    WITH 
      D = 2.0d0 * R,
      D2 = D*D,
      H = CEILING(D),
      N = 2*H + 1,
      norm = D2
    DO

      PROCEDURE TentEval(x, y: LONG): LONG =
        BEGIN
          WITH
            d2 = x*x+y*y
          DO
            IF d2 >= D2 THEN
              RETURN 0.0d0
            ELSE
              WITH 
                r = sqrt(d2),
                cos = r/D,
                sin = sqrt(1.0d0 - cos*cos),
                T = cos * sin,
                S = acos(cos)
              DO
                RETURN norm * MAX(0.0d0, S-T)
              END
            END
          END
        END TentEval;

      BEGIN
        RETURN MakeSymmetricTable(TentEval, N, N, uniDim)
      END;
    END
  END Tent;

PROCEDURE Gauss(R: LONG; cutoff: LONG; uniDim: BOOL := FALSE): REF T =
  BEGIN
    <* ASSERT cutoff > 0.0d0 AND cutoff < 1.0d0 *>
    IF R = 0.0d0 THEN
      WITH wR = NEW(REF T, 1, 1), w = wR^ DO
        w[0,0] := 1.0d0;
        RETURN wR
      END
    ELSE
      (* 
        The total probability of the unmodified 2D Gaussian outside a
        circle of radius "H" is "exp(-(H/R)^2)".
      *)
      WITH 
        R2 = R*R,
        scale = 1.0d0/R2,
        norm = 1.0d0/(Pi*R2),
        H = CEILING(R*sqrt(-log(cutoff))),
        N = 2*H + 1,
        H2 = FLOAT(H*H, LONG),
        bias = exp(-scale*H2)
      DO

        PROCEDURE GaussEval(x, y: LONG): LONG =
          BEGIN
            WITH r2 = x*x+y*y DO
              IF r2 > H2 THEN 
                RETURN 0.0d0 
              ELSE 
                RETURN norm * (exp(-scale*r2) - bias)
              END
            END
          END GaussEval;

        BEGIN
          RETURN MakeSymmetricTable(GaussEval, N, N, uniDim)
        END;
      END
    END
  END Gauss;

PROCEDURE BinomialDistr(M: CARDINAL): REF ARRAY OF LONG =
  (* 
    Returns a vector "c" such that "c[k] = (1/2)^M * choose(M,k)"
  *)
  VAR p: LONG := 1.0d0; f: LONG := 1.0d0;
  BEGIN
    WITH
      R = (M-1) DIV 2,
      cR = NEW(REF ARRAY OF LONG, M+1), c = cR^
    DO
      FOR k := 1 TO M DO f := f/2.0d0 END;
      <* ASSERT f > 0.0d0 *>
      FOR k := 0 TO R DO 
        c[k] := p*f;
        c[M-k] := p*f;
        p := (p*FLOAT(M-k, LONG))/FLOAT(k+1, LONG)
      END;
      IF M MOD 2 = 0 THEN c[R+1] := p END;
      RETURN cR
    END
  END BinomialDistr;

PROCEDURE Binomial(RX, RY: CARDINAL; cutoff: LONG := 0.0d0): REF T =
  VAR HX, HY: CARDINAL;
      SX, SY: LONG;
      biasX, biasY: LONG;
  BEGIN
    WITH
      cX = BinomialDistr(2*RX)^,
      cY = BinomialDistr(2*RY)^
    DO
      (* Find true table dimensions "HX, HY": *)
      HX := RX; SX := 0.0d0; biasX := 0.0d0; 
      WHILE SX < 0.5d0 * cutoff DO 
        SX := SX + cX[RX+HX] + cX[RX-HX];
        biasX := MIN(cX[RX+HX], cX[RX-HX]);
        DEC(HX)
      END;
      HY := RY; SY := 0.0d0; biasY := 0.0d0; 
      WHILE SY < 0.5d0 * cutoff DO 
        SY := SY + cY[RY+HY] + cY[RY-HY];
        biasY := MIN(cY[RY+HY], cY[RY-HY]);
        DEC(HY)
      END;
            
      PROCEDURE BinomialEval(x, y: LONG): LONG =
        BEGIN
          WITH
            i = RX+ROUND(x),
            j = RY+ROUND(y)
          DO
            RETURN (cX[i]-biasX)*(cY[j]-biasY)
          END
        END BinomialEval;

      BEGIN
        RETURN MakeSymmetricTable(BinomialEval, 2*HX+1, 2*HY+1, (HY = 0))
      END;
    END
  END Binomial;

PROCEDURE Power(e: LONG; cutoff: LONG; uniDim: BOOL := FALSE): REF T =
  BEGIN
    <* ASSERT e > 0.0d0 *>
    <* ASSERT cutoff > 0.0d0 AND cutoff <= 1.0d0 *>
    WITH
      eh = e/2.0d0,
      tau = exp(-log(cutoff)/eh),
      H = CEILING(sqrt(tau - 1.0d0)),
      H2 = FLOAT(H*H, LONG),
      N = 2*H + 1,
      bias = exp(-log(1.0d0+H2)*eh)
    DO

      PROCEDURE PowerEval(x, y: LONG): LONG =
        BEGIN
          WITH r2 = x*x+y*y DO 
            IF r2 >= H2 THEN 
              RETURN 0.0d0
            ELSE
              RETURN exp(-log(1.0d0+r2)*eh) - bias
            END
          END
        END PowerEval;

      BEGIN
        RETURN MakeSymmetricTable(PowerEval, N, N, uniDim)
      END;
    END
  END Power;
  
PROCEDURE FromParams(pp: ParseParams.T): REF T RAISES {ParseParams.Error} =

  PROCEDURE ParseCutoff(required: BOOL): LONG RAISES {ParseParams.Error} =
    BEGIN
      IF pp.testNext("cutoff") THEN 
        RETURN pp.getNextLongReal(0.0d0, Huge)
      ELSE
        IF required THEN 
          pp.error("parameter \"-cutoff\" is required for this distribution");
          <* ASSERT FALSE *>
        ELSE
          RETURN 0.0d0
        END
      END
    END ParseCutoff;
    
  VAR uniDim: BOOL;
      maxRadius: CARDINAL;
  BEGIN
    IF pp.testNext("2d") THEN
      uniDim := FALSE
    ELSIF pp.testNext("1d") THEN
      uniDim := TRUE
    ELSE
      uniDim := FALSE
    END;
    
    IF uniDim THEN maxRadius := 10000 ELSE maxRadius := 100 END;
    
    WITH fam = pp.getNext() DO
      IF Text.Equal(fam, "flat") THEN
        WITH 
          R = pp.getNextLongReal(Tiny, FLOAT(maxRadius, LONG))
        DO
          RETURN Flat(R, uniDim)
        END
      ELSIF Text.Equal(fam, "tent") THEN
        WITH 
          R = pp.getNextLongReal(Tiny, FLOAT(maxRadius DIV 2, LONG))
        DO
          RETURN Tent(R, uniDim)
        END
      ELSIF Text.Equal(fam, "power") THEN
        WITH 
          e = pp.getNextLongReal(0.0d0, 10.0d0),
          cutoff = ParseCutoff(required := TRUE) 
        DO
          RETURN Power(e, cutoff, uniDim)
        END
      ELSIF Text.Equal(fam, "gauss") THEN
        WITH 
          R = pp.getNextLongReal(1.0d0, FLOAT(maxRadius, LONG)/2.0d0),
          cutoff = ParseCutoff(required := TRUE) 
        DO
          RETURN Gauss(R, cutoff, uniDim)
        END
      ELSIF Text.Equal(fam, "binomial") THEN
        VAR RX, RY: CARDINAL;
        BEGIN
          RX := pp.getNextInt(0, maxRadius);
          IF uniDim THEN RY := 0 ELSE RY := pp.getNextInt(0, maxRadius) END;
          WITH 
            cutoff = ParseCutoff(required := FALSE)
          DO
            RETURN Binomial(RX, RY, cutoff)
          END
        END
      ELSE
        pp.error("invalid distribution");
        RETURN NIL
      END;
    END
  END FromParams;
        
PROCEDURE Print(wr: Wr.T; READONLY w: T) =
  <* FATAL Thread.Alerted, Wr.Failure *>
  
  PROCEDURE FLR(x: LONG; w, p: CARDINAL): TEXT =
    BEGIN
      RETURN Fmt.Pad(Fmt.LongReal(x, Fmt.Style.Fix, p), w)
    END FLR;
  
  VAR sumw, sumww, sumwx, sumwy, sumwxx, sumwxy, sumwyy: LONG := 0.0d0;
      sumwP, sumwN: LONG := 0.0d0;
      allPos: BOOLEAN := TRUE;
  BEGIN
    WITH 
      NY = NUMBER(w),
      NX = NUMBER(w[0]),
      xc = FLOAT(NX-1, LONG)/2.0d0,
      yc = FLOAT(NY-1, LONG)/2.0d0
    DO
      Wr.PutText(wr, "NX = " & Fmt.Int(NX) & "\n");
      Wr.PutText(wr, "NY = " & Fmt.Int(NY) & "\n");
      FOR yw := NY-1 TO 0 BY -1 DO
        FOR xw := 0 TO NX-1 DO
          WITH 
            x = FLOAT(xw, LONG) - xc,
            y = FLOAT(yw, LONG) - yc,
            w = w[yw,xw]
          DO 
            Wr.PutText(wr, FLR(w, 8, 5));
            sumw := sumw + w;
            sumww := sumww + w*w;
            sumwx := sumwx + w*x;
            sumwy := sumwy + w*y;
            sumwxx := sumwxx + w*x*x;
            sumwxy := sumwxy + w*x*y;
            sumwyy := sumwyy + w*y*y;
            IF w >= 0.0d0 THEN
              sumwP := sumwP + w
            ELSE
              sumwN := sumwN + w;
              allPos := FALSE
            END
          END
        END;
        Wr.PutText(wr, "\n")
      END;
      Wr.PutText(wr, "sum w =    " & FLR(sumw,0,7) & "\n");
      Wr.PutText(wr, "sum w+ =   " & FLR(sumwP,0,7) & "\n");
      Wr.PutText(wr, "sum w- =   " & FLR(sumwN,0,7) & "\n");
      Wr.PutText(wr, "sum |w| =  " & FLR(sumwP - sumwN,0,7) & "\n");
      Wr.PutText(wr, "sum w^2 =  " & FLR(sumww,0,7) & "\n");
      IF sumw > 0.0d0 THEN
        WITH
          ax = sumwx/sumw,
          ay = sumwy/sumw
        DO
          Wr.PutText(wr, "center =   " & "(" & FLR(ax,0,2) & "," & FLR(ay,0,2) & ")" & "\n");
        END
      END;
      IF allPos THEN 
        WITH
          axx = sumwxx/sumw,
          axy = sumwxy/sumw,
          ayy = sumwyy/sumw,
          R = sqrt(axx + ayy)
        DO
          Wr.PutText(wr, "xx/xy/yy = " & FLR(axx,0,4) & "/" & FLR(axy,0,4) & "/" & FLR(ayy,0,4) & "\n");
          Wr.PutText(wr, "R =   " & FLR(R,0,4) & "\n");
        END
      END;
    END
  END Print;

BEGIN
END PixelWeightTable.