(* Last edited on 1999-03-14 07:07:42 by stolfi *)

INTERFACE QuadOps;

(* Second-order image operators based on a local quadratic approximation *)

IMPORT PixelWeightTable, PGMImageFilter, Random;

(*
  This interfaces provides some second-order image operators,
  including some rotational invariants.

  The operators can be expressed in terms of a weighted least-squares
  second-degree approximation P of the image in the neighborhood of a
  given pixel. Let "(x,y)" be Cartesian coordinates of a point
  relative to the center of the pixel in question. The approximation
  has the form

    P(x,y) = A (x^2+y^2-k)/4 + 
             B (x^2-y^2)/4 +
             C x y / 2 + 
             D x + E y + 
             F + 
             o(x^2+y^2)  (1)

  where k is some constant that makes the basis element (x^2 + y^2 - k)
  orthogonal to the constant element, in the discrete scalar
  product used in the least squares approximation. Note that the other
  pairs of basis elements are orthogonal in this metric, too.
  
  These coefficients are related to the value of P at (0,0)
  and its derivatives Px, Py, Pxx, Pxy, Pyy as follows:
  
    P   = F + A k        F = P - (Pxx + Pyy) k
    Px  = D              D = Px
    Py  = E              E = Py
    Pxx = (A + B)/2      A = Pxx + Pyy
    Pxy = C/2            B = Pxx - Pyy
    Pyy = (A - B)/2      C = 2 Pxy
  
*)

TYPE
  Coeffs = ARRAY [0..5] OF Interval; 
    (*
      Coefficients of the local quadratic approximation, in order
      F E D C B A.
    *)
TYPE
  Interval = PGMImageFilter.Interval;

  LONG = LONGREAL;
  NAT = CARDINAL;
  BOOL = BOOLEAN;

TYPE
  WeightTable = PixelWeightTable.T;
  WeightTables = ARRAY [0..5] OF REF WeightTable;
    (* 
      Pixel weights used to compute the coefficients 
       of the local quadratic approximation. *)
  
PROCEDURE BuildWeightTables(
    READONLY g: WeightTable;
    varScale: BOOL := FALSE;
    sumScale: BOOL := FALSE;
  ): WeightTables;
  (* 
    Returns pixel weight tables such that the coefficients A-F (in
    the "Coeffs" order, see above) are obtained by weighted sums of
    the pixels surrounding a given pixel. 
    
    The table "g" defines the size of the window and the 
    relative importance of each pixel in the neighborhood.
    
    If "varScale" and "sumScale" are false, the tables are scaled
    so that an image P that is a second-degree polynomial of the 
    form (1) will have "Coeffs" matching those of formula (1).
    
    If "varScale" is true, the tables are scaled so that the sum of the
    squares if 1. This means the variance of each coefficient will be
    1/12, when the input is white noise in "[0 _ 1]".
    
    If "sumScale" is true, the tables are scaled so that the sum of the
    absolute values of the terms is 1. This ensures that when the
    input pixels range in "[0 _ 1]" the range of the result will have
    width 1. Note that for the F coefficient (index 0) the range is
    centered at 1/2 whereas for the other coefficients it is centered
    at 0.
    
    Note that either normalization preserves the ratios D:E and B:C,
    but may not preserve the ratios between the groups {F}, {D,E},
    {B,C}, and {A}.  This fact must be taken into account when
    computing the Gaussian curvature, for example.
  *)
  
(* ROTATIONAL INVARIANTS *)

(* 
  If the image is rotated around the reference pixel by an angle
  "theta", the coefficients A and F remain unchanged; the gradient
  vector (Px,Py) = (D,E) rotates by "theta"; and the warp vector 
  (Pxx-Pyy, 2 Pxy) = (B,C) rotates by "2*theta".
  
  Therefore, from the coefficients A..F, we can derive the following
  rotational invarants:

    mv mean value     F                       P                             
                                                                        
    gr slope squared  D^2 + E^2               Px^2 + Py^2               
                                                                        
    lp Laplacian      A                       Pxx + Pyy                     
                                                                        
    wp warp           B^2 + C^2               (Pxx - Pyy)^2 + 4*Pxy^2     
                                                                        
    sd saddle dot     B (D^2-E^2) + 2 C D E   ...                             
    sc saddle cross   2 B D E - C (D^2 - E^2) ...                          

  The first four invariants are independent. The last two introduce
  only one degree of freedom; they are the dot and cross product of
  the warp vector (B,C) and the vector(D^2-E^2, 2DE), which is the
  complex square of the gradient vector (D, E). Unfortunately we
  cannot keep only one of them (we would miss a sign) nor reduce them
  to an angle (sometimes both are zero).

  From these invariants we can get also the Gaussian (intrinsic, product)
  curvature 
  
    cv = (A^2 - B^2 - C^2)/4 = Pxx Pyy - Pxy^2
    
  which is positive for bump/dent, negative for saddle-shape, and 0 
  for flat/cylinder. We can get also the the cylindricity indicator
  
    cy = A*sqrt(B^2 + C^2)/2 = (Pxx + Pyy)*sqrt((Pxx - Pyy)^2/4 + Pxy^2)
    
  which is positive for a valley, negative for a crest, and 0 for a
  saddle, bump, dent, or flat.  And, finally, the total quadratic
  component
  
    tq = (A^2 + B^2 + C^2)/4 = (Pxx^2 + Pyy^2)/4 + Pxy^2
    
  which is 0 for flat regions, positive for non-flat.
  Note that 
   
    cv^2 + cy^2 = tq^2
    
  so the ratios cv/tq and cy/tq are like sin and cos of an angle.
*)  

TYPE 
  Invariants = RECORD
      mv: Interval;   (* Mean value, F. *)
      gr: Interval;   (* Gradient length squared, D^2+E^2. *)
      lp: Interval;   (* Laplacian, A. *)
      wp: Interval;   (* Warp, B^2 + C^2. *)
      sd: Interval;   (* Saddle dot, B (D^2-E^2) + 2 C D E. *)
      sc: Interval;   (* Saddle cross, 2 B D E - C (D^2 - E^2). *)
      cv: Interval;   (* Gaussian curvature, (A^2 - B^2 - C^2)/4. *)
      cy: Interval;   (* Cylindricity, A*sqrt(B^2 + C^2)/2. *)
      tq: Interval;   (* Total quadratic component, (A^2 + B^2 + C^2)/4. *)
      tv: Interval;   (* Total variation, (A^2 + B^2 + C^2)/4 + (D^2 + E^2)/2. *)
      cx: Interval;   (* Cylindricity index, cy/tv. *)
    END;

PROCEDURE MeanValue(F: Interval): Interval;
PROCEDURE GradSqr(D, E: Interval): Interval;
PROCEDURE Laplacian(A: Interval): Interval;
PROCEDURE Warp(B, C: Interval): Interval;
PROCEDURE SaddleDot(B, C, D, E: Interval): Interval;
PROCEDURE SaddleCross(B, C, D, E: Interval): Interval;
PROCEDURE Curvature(A, B, C: Interval): Interval;
PROCEDURE Cylindricity(A, B, C: Interval): Interval;
PROCEDURE TotQuadratic(A, B, C: Interval): Interval;
PROCEDURE TotVariation(A, B, C, D, E: Interval): Interval;
PROCEDURE CylIndex(A, B, C, D, E: Interval): Interval;
  (*
    The rotational invariants, individually. *)

PROCEDURE InvariantsFromCoeffs(READONLY c: Coeffs): Invariants;
  (*
    Computes the rotational invariants from the given coefficients. *)
  
PROCEDURE CoeffRange(READONLY ww: WeightTables; pixelRange: Interval): Coeffs;
  (*
    Returns the coefficient ranges for the given tables, 
    assuming the input pixels are independently over the
    specified pixel range. *)
    
PROCEDURE CoeffSample(
    READONLY ww: WeightTables;
    range: Interval; 
    rnd: Random.T; 
    eps: LONG;
  ): Coeffs;
  (*
    Returns a coefficient vector for a random neighborhood.
    The pixels values will be random numbers uniformly distributed
    in the given "range", plus or minus "eps". *)
  
END QuadOps.