INTERFACE PixelWeightTable;
(* Last edited on 1998-10-25 23:12:34 by stolfi *)

(* Routines to create pixel weight tables of various kinds. *)

IMPORT Wr, ParseParams;

TYPE
  LONG = LONGREAL;
  BOOL = BOOLEAN;
  
  T = ARRAY OF ARRAY OF LONG;  
    (* 
      A table of pixel weights.

      A weight table of any kind has "NY" rows and "NX" columns, where
      "NX" and "NY are odd integers. The center pixel is thus
      "w[HY,HX]" where "HX=(NX-1)/2" and "HY=(NY-1)/2". When "NX=NY"
      and "HX=HY", we will denote these parameters by "N" and "H",
      respectively

      The weight of each pixel is usually defined in terms of a system
      of pixel coordinates "(x,y)" whose origin is at the center of
      the midle pixel, with "x" horizontal and "y" vertical. Thus, the
      pixel centers have coordinates in the range
      "[-HX..+HX]×[-HY..+HY]"; and the weight for the pixel with
      center at coordinates "(x,y)" is "w[HY-y,HX+x]".
   *)

CONST
  MaxWidth = 1023;

PROCEDURE NormalizeSum(VAR w: T);
  (*
    Normalizes the table entries so that they add to 1. *)

PROCEDURE NormalizeSumOfSquares(VAR w: T);
  (* 
    Normalizes "w" so that the sum of the squares is 1. *)

PROCEDURE NormalizeSumOfAbs(VAR w: T);
  (* 
    Normalizes "w" so that the sum of the absolute values
    of all elements is 1. *)

PROCEDURE NormalizeSignedSum(VAR w: T);
  (* 
    Normalizes "w" so that the sum of the positive elements is at
    most 1, the sum of negative elements is at least -1, and at least
    one of the two limits is tight. *)

PROCEDURE Balance(VAR w: T; READONLY u: T);
  (*
    Adds to "w" a multiple of "u" such that the resulting table has
    zero sum.  Requires "u" to have non-zero sum. *)

(*
  SPECIFIC TABLES
  
  In the following procedures, the dimensions "NX,NY" of the resulting
  table are always odd integers. If "uniDim" is false, the table
  dimensions are determined implicitly from the other parameters. If
  "uniDim" is true, the table will be reduced to a single row, equal
  to the middle row of the full table. *)

PROCEDURE Flat(R: LONG; uniDim: BOOL := FALSE): REF T;
  (*
    Circular uniform ("pillbox") weight function: "w[HY-y,HX+x]" is 1
    if pixel centered at "(x,y)" is entirely inside the circle of
    radius "R", is 0 if the pixel is entirely outside. 
    Pixels that straddle the circle's boundary will be antialiased.
    The full table will have "NX = NY = 2*CEILING(R) + 1". *)

PROCEDURE Tent(R: LONG; uniDim: BOOL := FALSE): REF T;
  (*
    Tent function, the convolution of two "Flat" distributions with
    radius "R".  In other words, "w[H-y,H+x]" is the area of overlap
    of two circles of radius "R", one centered at the origin, the
    other centered at "(x,y)". The full table will have 
    "N = 2*CEILING(2*R)+1". *)

PROCEDURE Gauss(R: LONG; cutoff: LONG; uniDim: BOOL := FALSE): REF T;
  (*
    Truncated Gaussian weight distribution with nominal radius "R",
    "w[H-y,H+x] = (exp(-r^2/R^2) - exp(-H^2/R^2))/(Pi*R^2)" for "r < H", 0
    for "r > H", where "r^2 = x^2 + y^2", and
    "H" is such that the total mass of the original Gaussian that lies
    outside the circle of radius "H" is equal to "cutoff". *)
    
PROCEDURE Binomial(RX, RY: CARDINAL; cutoff: LONG := 0.0d0): REF T;
  (* 
    A truncated 2D binomial weight distribution, "w[HY-y,HX+x] =
    (B(RX,x) - B(RX,HX+1))*(B(RY,Y) - (RY,HY+1))" for "|x| < HX" and
    "|y| < HY", where "B(R,k) = choose(2*R,R+k)/4^R", "HX" is such
    that the sum of "B(RX,k)" for "|k| > HX" is "cutoff/2" or less,
    and analogously for "HY". Thus the total mass of the original
    (non-truncated) 2D binomial distribution that lies outside the
    domain of "w" is approximately "cutoff" (for small "cutoff").
    
    Note that the distribution tends to a gaussian for large "RX = RY".  
    The default "cutoff=0" implies the table will have
    "NX=2*RX+1" and "NY=2*RY+1". *)

PROCEDURE Power(e: LONG; cutoff: LONG; uniDim: BOOL := FALSE): REF T;
  (*
    An inverse-power weight function, "w[H-y,H+x] = P(r) - P(H)" for
    "r < H", 0 for "r > H", where "r = sqrt(x^2 + y^2)", 
    "P(r) = 1/(1 + r^2)^(e/2)", and "H" is such that "P(r) = cutoff".
    Note that the integral of the non-truncated distribution "P(r)" 
    diverges when "e" is 2 or less. *)

PROCEDURE FromParams(
    pp: ParseParams.T
  ): REF T RAISES {ParseParams.Error};
  (*
    Parses a weight table description from the command line arguments.
    The arguments are assumed to follow the syntax defined in the
    "ParamsHelp" string below.  
    
    The "cutoff" parameter is mandatory for infinte distributions
    like "gauss" and "power". For other distributions, it defauts
    to 0.0. 
    
    For the meaning of the arguments, see the procedures "Gauss",
    "Flat", ""Power", etc above.  *)
    
CONST
  ParamsHelp = 
    "    [ 2d | 1d ] \\\n" &
    "    { flat RADIUS | \\\n" &
    "      tent RADIUS | \\\n" &
    "      power EXPT cutoff EPS \\\n" &
    "      gauss RADIUS cutoff EPS | \\\n" &
    "      binomial NX [ NY ] [ cutoff EPS ] | \\\n" &
    "    } \n";

PROCEDURE Print(wr: Wr.T; READONLY w: T);
  (*
    Prints "w" formatted to the writer "wr". *)

END PixelWeightTable.