INTERFACE PZLR3Chain;

(* Sequences of points of R^2 OR R^3, closed or open. *)
(* Last edited on 2001-11-15 17:28:15 by stolfi *)

IMPORT LR3, LR4x4, Rd, Wr;
IMPORT PZIntChain, PZLRChain, PZSegment, PZMatch, PZWindow;
FROM PZTypes IMPORT LONG, LONGS, NAT, BOOL, BOOLS;

TYPE
  T = ARRAY OF LR3.T;
    (*
      A sequence of points of R^3, often used to represent points on
      the plane (with the third component set to zero).  The elements
      may be also velocities or accelerations.  The chain may be
      interpreted as a polygonal or a spline, closed or open,
      depending on the context. *)

(* BASIC MANIPULATION *)

(*
  The "DoXXX" procedures below are analogous to the "XXX" 
  procedures, except that the result is returned in an array
  provided by the client (which must have the correct size). *)

PROCEDURE FromIntChain(READONLY c: PZIntChain.T): REF T;
  (*
    Converts "c" to floating point, adding the third coordinate "Z=0". *)

PROCEDURE ToIntChain(READONLY c: T): REF PZIntChain.T;
  (*
   Converts "c" to integer format, ignoring the "Z" coodinate
   and rounding "X" and "Y" to the nearest INT. *)
   
PROCEDURE FromLineSegment(READONLY p, q: LR3.T; n: NAT): REF T;
  (*
    Creates a chain with "n" points equally spaced
    along the segment "p--q". *)

PROCEDURE Trim(READONLY c: T; start, length: NAT): REF T;
PROCEDURE DoTrim(READONLY c: T; start, length: NAT; VAR x: T);
  (*
    Returns a new chain that is a copy of elements from "c[start]"
    to "c[start+length-1]".  Assumes the chain is periodic, i.e.
    "c[i+n]=c[i]" for all "i", where "n=NUMBER(c)". *)
    
PROCEDURE Reverse(VAR c: T); 
  (*
    Reverses the direction of traversal of the Chain "c". *)
  
PROCEDURE ExtractSegment(READONLY c: T; s: PZSegment.T): REF T;
PROCEDURE DoExtractSegment(READONLY c: T; s: PZSegment.T; VAR x: T);
  (*
    Extracts from "c" the part described by "s", 
    reversing it if "s.rev" is TRUE. *)

PROCEDURE GetWindow(READONLY c: T): PZWindow.T;
  (*
    The bounding box of chain "c", viewed as a polygonal chain. *)

PROCEDURE SampleBarycenter(READONLY c: T): LR3.T;
  (*
    The barycenter of all points of the chain. *)
    
PROCEDURE AreaBarycenter(READONLY c: T): LR3.T;
  (*
    Returns the barycenter "b" of the projection of "c" on the
    XY plane, seen as a polygonal *area*, with "b[2] = 0". *)

PROCEDURE Translate(READONLY c: T; d: LR3.T): REF T;
PROCEDURE DoTranslate(VAR c: T; d: LR3.T);
  (*
    Translates "c" by the vector "d". "Translate" makes a new copy and
    translates it, leaving "c" unchanged; "DoTranslate" modifies "c"
    itself. *)

PROCEDURE Map(READONLY c: T; READONLY m: LR4x4.T): REF T;
PROCEDURE DoMap(VAR c: T; READONLY m: LR4x4.T);
  (*
    Maps "c" by the given projective transformation.  "Map" makes a
    new copy and transforms it, leaving "c" unchanged; "DoMap"
    modifies "c" itself. *)
    
PROCEDURE AlignmentMatrix(READONLY a, b: T): LR4x4.T;
  (*
    Given two open chains "a" and "b", returns a matrix that moves the
    barycenter of "a" to the barycenter of "b", and rotates "a" so
    that its general direction is parallel to the general direction of "b". *)

(* CHAINS AS CLOSED POLYGONS *)

(* 
  The procedures in this section interpret the chain "p" as the
  vertices of a closed polygonal chain with "m = NUMBER(p)" sides.  It
  is assumed that each edge is traversed at constant speed, in such a
  way that the chain reaches vertex "p[i]" at time "t[i]", for "i IN
  [0..m-1]"; and the chain goes back to vertex "p[0]" at time "t[0] + tPeriod". *)
    
PROCEDURE LinearLengths(
    READONLY p: T;              (* "p[i]" is chain position at time "t[i]". *) 
    VAR s: PZLRChain.T;         (* "s[j]" is the length from "p[0] to "p[j] "j". *) 
  ): LONG;
  (*
    Computes the Euclidean length "s[j]" on the polygonal from
    vertex "p[0]" to vertex "p[j]", for "j" in "[0..LAST(p)]".
    Returns the total length of the polygonal as the result of the call. *)
    

PROCEDURE LinearUniformSample(
    READONLY t: PZLRChain.T;    (* Times of input samples. *) 
    tPeriod: LONG;              (* Time period *)
    READONLY p: T;              (* "p[i]" is chain position at time "t[i]". *) 
    tStart: LONG;               (* Time of first output sample *)
    VAR q: T;                   (* "q[j]" will be position sample number "j". *) 
  );
  (*
    Computes "n" samples along the chain, equally spaced in time, starting
    with time "tStart".  The samples are returned in "q[j]", "j IN [0..n-1]". *)

PROCEDURE LinearEquidistantSample(
    READONLY p: T;              (* "p[i]" is chain position at time "t[i]". *) 
    tStart: LONG;               (* Time of first output sample *)
    VAR q: T;                   (* "q[j]" will be position sample number "j". *) 
  );
  (*
    Same as "LinearUniformSample", assuming that the chain is traversed
    with constant velocity in unit time. *)

PROCEDURE LinearTimeSample(
    READONLY t: PZLRChain.T;    (* Times of input samples. *) 
    tPeriod: LONG;              (* Time period *)
    READONLY p: T;              (* Input chain positions at times "t[i]". *)
    READONLY r: PZLRChain.T;    (* Times of output samples. *)
    VAR q: T;                   (* Output chain positions at times "r[j]". *) 
  );
  (*
    Computes the positions "q[j]" at the specified times "tq[j]",
    "j IN [0..n-1]". *)
   
(* CHAINS AS CLOSED C_1 CUBIC SPLINES *)

(* 
  The procedures in this section interpret the chains "p" and "dp" as the
  positions and velocities of a closed spline chain with "m = NUMBER(p)" 
  cubic arcs, at the nodal times "t[0..m-1]".  It is assumed the chain 
  goes back to vertex "p[0]" and velocity "dp[0]" at time "t[0] + tPeriod". *)

 
PROCEDURE HermiteUniformSample(
    READONLY t: PZLRChain.T;    (* Times of input samples. *) 
    tPeriod: LONG;              (* Time period *)
    READONLY p: T;              (* "p[i]" is chain position at time "t[i]". *)
    READONLY dp: T;             (* "dp[i]" is velocity at time "t[i]". *)
    tStart: LONG;               (* Time of first output sample *)
    VAR q: T;                   (* "q[j]" will be position sample number "j". *) 
    VAR dq: T;                  (* "dq[j]" will be velocity at point "q[j]". *) 
  );
  (*
    Computes "n" samples "q[0..n-1]" along the chain, and their velocities "dq[0..n-1]",
    equally spaced in time, starting with time "tStart". *)

PROCEDURE HermiteEquidistantSample(
    READONLY p: T;              (* "p[i]" is chain position at time "t[i]". *) 
    tStart: LONG;               (* Time of first output sample *)
    VAR q: T;                   (* "q[j]" will be position sample number "j". *) 
    VAR dq: T;                  (* "dq[j]" will be velocity at point "q[j]". *) 
  );
  (*
    Same as "HermiteUniformSample", assuming that the chain is traversed
    with constant velocity in the total time "tPeriod". *)

PROCEDURE HermiteTimeSample(
    READONLY t: PZLRChain.T;    (* Times of input samples. *) 
    tPeriod: LONG;              (* Time period *)
    READONLY p: T;              (* "p[i]" is chain position at time "t[i]". *)
    READONLY dp: T;             (* "dp[i]" is velocity at time "t[i]". *)
    READONLY r: PZLRChain.T;    (* Times of output samples. *)
    VAR q: T;                   (* "q[j]" will be chain position at time "r[j]". *) 
    VAR dq: T;                  (* "dq[j]" will be velocity at time "r[j]". *) 
  );  
  (*
    Computes the positions "q[j]" and velocities "dq[j]" at the specified 
    times "r[j]", "j IN [0..n-1]". *)

PROCEDURE EstimateVelocities(
    READONLY t: PZLRChain.T;
    tPeriod: LONG;
    READONLY p: T;
    VAR dp: T;
    est: VelEstimator;
  );
  (*
    Estimates the velocity "dp[i]" at each point "p[i]",
    by applying the estimator "est" to three consecutive times. *)
    
TYPE
  VelEstimator = PROCEDURE (
      a: LONG; READONLY pa: LR3.T; 
      b: LONG; READONLY pb: LR3.T; 
      c: LONG; READONLY pc: LR3.T;
    ): LR3.T;

(* INPUT/OUTPUT *)

PROCEDURE AdjustUnit(givenUnit: LONG; READONLY c: T): LONG;
  (*
    Adjusts the "givenUnit" if needed to avoid overflow
    or excessive error when quantizing the values in "c".
    May print warnings to "stderr". *)

TYPE
  ReadData = RECORD 
      cmt: TEXT;      (* Comment text *)
      samples: NAT;   (* Number of samples in chain. *)
      c: REF T;       (* The samples *)
      unit: LONG;     (* Unit that was used when writing the data *)
    END;

PROCEDURE Read(rd: Rd.T; headerOnly: BOOL; centralize: BOOL): ReadData;
  (*
    Reads a chain (and its comments) from "rd". If "headerOnly=TRUE",
    reads only the parameters, and leaves the "c" field as NIL.
    If "centralize=TRUE", displaces the chain so that its "AreaBarycenter"
    lies at the origin. *)
 
PROCEDURE Write(wr: Wr.T; READONLY cmt: TEXT;  READONLY c: T; unit: LONG);  
  (*
    Writes chain "c" to "wr", prefixed with comments "cmt".
    The coordinates will be written as integer multiples of the given "unit". *)
    
TYPE
  ReadAllData = RECORD
      chData: REF ARRAY OF ReadData;
      unitMin: LONG;      (* Minimum quantization unit. *)
      unitMax: LONG;      (* Maximum quantization unit. *)
      cmt: TEXT;          (* Comment of first chain, if any. *)
    END;

PROCEDURE ReadAll(
    prefix: TEXT; 
    band: NAT;
    extension: TEXT;
    READONLY sel: BOOLS;
    headerOnly: BOOL;
    centralize: BOOL;
    dir: TEXT := ".";
  ): ReadAllData;
  (*
    Reads a set of numbered chains, returns a "ReadAllData" record
    "r". Chain number "i" is read if and only if "sel[i] = TRUE", and
    its "ReadData" is stored in "r.chain[i]"; otherwise "r.chain[i] "
    is set to a record with null "c" field.
    
    The data is read from file "DDD/KKKK/PPPBBBEEE", where
    "DDD" is the directory "dir", "KKKK" is the chain number "k" (4
    digits, zero padded), "PPP" is the given "prefix", "BBB" is the
    "band" number (3 digits, zero padded), and "EEE" is the given
    "extension".
    
    If "headerOnly" is TRUE, only the header of each file is read;
    the "c" fields are all set to NIL.
    
    The range of quantization "unit"s of the chains that were read is
    returned in "r.unitMin", "r.unitMax".  The comment of the first chain
    read, plus the call arguments, are stored in "r.cmt". *)

(* CHAIN MATCHING BY DYNAMIC PROGRAMMING *)

(*
  The procedures "FindBestMatch" and "RefineMatch" look for a monotone
  matching "mt" (a "PZMatch.T") between two chains that minimizes some
  measure of the distance between corresponding samples.
 
  See "PZMismatch" for the relevant concepts.  These procedures take
  the sample distance "d(a[r], b[s])" to be the Euclidean distance
  between the points "a[r] and "b[s]".
    
  The procedures return in "length" the total length "Len(mt)" of the
  matching (average of the lengths of the two chain segments spanned
  by "mt").  They also return in "lengthMatched" the total length of
  the matching steps whose mean chain distance "dist(S)" was strictly less than
  "maxDist".

  The procedures use internally a working array with "NUMBER(a)" rows
  and "NUMBER(b)" columns.  The client may optionally provide an array
  "rCost", at least this big, which will be automatically (re) allocated
  if omitted or too small. *)

PROCEDURE FindBestMatch(
    READONLY a: T;
    READONLY b: T;
    maxDistSqr: LONG;
    removeUnmatchedEnds: BOOL;
    VAR (*OUT*) mt: REF PZMatch.T;
    VAR (*OUT*) avgDist: LONG;
    VAR (*OUT*) length: LONG;
    VAR (*OUT*) matchedLength: LONG; 
    VAR (*WORK*) rCost: REF PZMatch.CostMatrix;
  ); 
  (*
    Computes (by dynamic programming) a complete monotone matching
    "mt" between two open chains "a" and "b" that minimizes the 
    integrated squared chain distance "IntgDistSqr(mt)".
    
    This procedure allows the chains to be sampled at arbitrary intervals,
    and takes the metric "h(a,i1,i2)" to be the Euclidean distance
    of the points "a[i1]" and "a[i2]".
    
    If "removeUnmatchedEnds" is true, the procedure removes any
    initial or final steps "s" whose distance "dist(s)" is equal to
    "maxDist".  Otherwise the matching "mt" will completely cover the
    two chains. *)
    
PROCEDURE RefineMatch(
    READONLY a: T;
    READONLY b: T;
    step: LONG;
    READONLY mtOld: PZMatch.T;
    maxAdjust: NAT;
    critDistSqr: LONG;
    maxDistSqr: LONG;
    VAR (*OUT*) mt: REF PZMatch.T;
    VAR (*OUT*) mismatch: LONG;
    VAR (*OUT*) length: LONG;
    VAR (*OUT*) matchedLength: LONG; 
    VAR (*OUT*) minAvgDist: LONGS;
    VAR (*WORK*) rCost: REF PZMatch.CostMatrix;
  ); 
  (*
    Finds (by bidirectional dynamic programming) a partial monotone
    matching "mt" between the polygonal chains "a" and "b", that is
    a refinement of an `input matching' "mtOld", and minimizes
    the mismatch "Mis(mt)".  
    
    See "PZMatch.Refine" for an explanation of the refinement concept,
    and the meaning of "mtOld" and "maxAdjust".  See "PZMismatch.i3"
    for the definition of the mismatch "Mis(mt)" and the meaning of the parameters
    "critDistSqr" and "maxDistSqr".
    
    This procedure assumes that the two chains are smooth chains that
    were sampled at uniform intervals of given size "step".
    Therefore, the metric "h(a,i1,i2)" is assumed to be
    "step*(i2-i1)".
    
    The procedure returns in "minAvgDist[K]", for "K" in
    "[0..na+nb-2]", the root mean square chain distance "dist(s)" of the
    optimum match among all matches with total step count "K".  (Note
    that, if the step count is fixed, the "critDist" term does not
    affect the optimality of the match.)  *)
  
(* FOURIER ANALYSIS AND SYNTHESIS *)
    
TYPE
  Spectrum = ARRAY OF LONGS;

PROCEDURE FourierTransform(READONLY c: T; VAR f: Spectrum);
  (*
    Computes the Fourier spectrum "f" of chain "c". 
    The length of the array "c" and of the rows of "f"
    must be the same power of two. *)
  
PROCEDURE InvFourierTransform(VAR f: Spectrum; VAR c : T);  
  (*
    Computes a chain "c" given its Fourier spectrum "f".
    Their lengths must be equal and a power of two.
    WARNING: destroys "f" in the process. *)
    
END PZLR3Chain.
(*
  Copyright © 2001 Universidade Estadual de Campinas (UNICAMP).
  Authors: Helena C. G. Leitão and Jorge Stolfi.
  
  This file can be freely distributed, used, and modified, provided
  that this copyright and authorship notice is preserved, and that any
  modified versions are clearly marked as such.
  
  This software has NO WARRANTY of correctness or applicability for
  any purpose. Neither the authors nor their employers chall be held
  responsible for any losses or damages that may result from its use.
*)