INTERFACE PixelWeightTable;
(* Last edited on 1999-04-19 17:36:20 by stolfi *)

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

IMPORT Wr, Rd, ParseParams, IA;

TYPE
  LONG = LONGREAL;
  BOOL = BOOLEAN;
  NAT = CARDINAL;
  Interval = IA.T;
  
  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 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. *)

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

PROCEDURE NormalizeMoment(VAR w: T; X, Y: LONG);
  (* 
    Scales the table entries so that its dot product with a ramp
    "P(x,y) = X*x + Y*y" is 1. The "(x,y)" coordinates are measured
    from the center of the central pixel. *)

PROCEDURE NormalizeWarp(VAR w: T; H, M, N: LONG);
  (*
    Scales the table entries so that its dot product with the quadric
    "H(x^2+y^2) + M(x^2-y^2) + N x y" is 1. The "(x,y)" coordinates
    are measured from the center of the central pixel. The quadric is 
    adjusted assuming a Gaussian pixel response function. *)

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

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

PROCEDURE NormalizeSignedSum(VAR w: T);
  (* 
    Scales "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 CoeffRange(READONLY w: T; pixel: Interval): Interval;
  (*
    Returns the range for the result of convolving the weights "w"
    with an image whose pixals range independently over the spacified
    range "pixRange". *)

(*
  SPECIFIC TABLES
  
  In the following procedures, the dimensions "NX,NY" of the resulting
  table are always odd integers. If "uniDim" is false, the height
  of the table is 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.
  
  The table entries are defined by a weight function "f(x,y)",
  that often depends on a characteristic radius "R".
  
  Note that "f" may be nonzero outside that radius. The actual width
  of the table is usually determined by the "cutoff" and "clip"
  parameters, when present. Typically, the "cutoff" level determines
  indirectly a ``significant radius'' "S", the maximum 
  distance from the center for which the weight function "f" is still
  significant. The actual window radius "H" is defined as
  "MIN(CEILING(S)-1, clip)". For non-isotropic functions like
  "Binomial", these criteria define separately an half-width "HX" and
  a half-height "HY"; otherwise it is assumed that "HX = HY = H".
  
  The weight function "f" gets modified, if necessary so that it is
  zero when "|x| >= HX+1" or "|y| >= HY+1". Since all the weights
  provided here are non-increasing along each ray, so the adjustment
  is performed by adding or subtracting the offset
  "MAX(ABS(f(HX+1,0)),ABS(f(0,HY+1)))". In any case, the returned
  table has "NX = 2*HX+1" columns and "NY = 2*HY+1" rows.

*)

CONST 
  DontClip = LAST(NAT);

PROCEDURE Flat(R: LONG; uniDim: BOOL := FALSE): REF T;
  (*
    Circular uniform ("pillbox") weight function: "f(x,y)" 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 significant radius is "S = R + 0.5". *)

PROCEDURE Tent(R: LONG; clip: NAT := DontClip; uniDim: BOOL := FALSE): REF T;
  (*
    Tent function, the convolution of two "Flat" distributions with
    radius "R".  In other words, "f(x,y)" is the area of overlap
    of two circles of radius "R", one centered at the origin, the
    other centered at "(x,y)". The significant radius is "S = 2*R". *)

PROCEDURE Gauss(
    R: LONG; 
    cutoff: LONG; 
    clip: NAT := DontClip; 
    uniDim: BOOL := FALSE;
  ): REF T;
  (*
    Gaussian weight distribution with nominal radius "R",
    "f(r) = exp(-r^2/R^2)/(Pi*R^2)" where where "r^2 = x^2 + y^2"
    The significant radius is "S" such that "f(S) = cutoff".
  *)
    
PROCEDURE Binomial(RX, RY: CARDINAL; cutoff: LONG := 0.0d0; clip: NAT := DontClip): REF T;
  (* 
    A truncated 2D binomial weight distribution with natural width "RX"
    and natural height "NY".
    
    Let "B(R,k) = choose(2*R,R+k)/4^R". The significant X-radius 
    "SX" is the smallest integer such that the sum of "B(RX,k)" for "|k| >= SX" is
    "cutoff/2" or less.  The significant Y-radius "SY" is defined in the same way.
    
    Note that the distribution tends to a gaussian for large "RX = RY".  
    The default "cutoff = 0", "clip = DontClip" implies the table will have
    "NX=2*RX+1" and "NY=2*RY+1". *)

PROCEDURE Power(
    e: LONG; 
    bias: LONG;
    cutoff: LONG; 
    clip: NAT := DontClip; 
    uniDim: BOOL := FALSE;
  ): REF T;
  (*
    An inverse-power weight function, approximately 
    "f(r) = 1/(bias^2 + r^2)^(e/2)" where "r = sqrt(x^2 + y^2)".
    The significant radius "S" is the radius for which "f(S) = cutoff".
    Note that the integral of the non-truncated distribution "f(r)" 
    diverges when "e" is 2 or less. *)

PROCEDURE FromFile(rd: Rd.T; lo, hi: LONG): REF T;
  (*
    Reads the weight table from the given image file, presently restricted to 
    PGM format.  The image pixels are converted linearly to float numbers
    between "lo" and "hi". *)

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. 
    
    In the "file" option, the "FILENAME" must include the extension
    (only ".pgm" implemented for now). Pixels will be mapped 
    linearly from their maximum range in the file to the range 
    "[lo_hi]" (default "[0_1]").
    
    For the meaning of the arguments, see the procedures "Gauss",
    "Flat", ""Power", etc above.  *)
    
CONST
  ParamsHelp = 
    "    [ 2d | 1d ] \\\n" &
    "    { flat RADIUS | \\\n" &
    "      tent RADIUS [ clip MAXRAD ] | \\\n" &
    "      power EXPT [ bias BS ] { cutoff EPS | clip MAXRAD } \\\n" &
    "      gauss RADIUS { cutoff EPS | clip MAXRAD }\\\n" &
    "      binomial NX [ NY ] [ cutoff EPS | clip MAXRAD ] | \\\n" &
    "      file FILENAME [ range LOW HIGH ] \\\n" &
    "    } \n";

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

END PixelWeightTable.