MODULE PZShape;
(* Last edited on 1999-08-06 22:17:06 by hcgl *)

IMPORT PZLR3Chain, PZLRChain, PZFourier, LR3;

FROM Math IMPORT atan2, sqrt;
FROM PZTypes IMPORT LONG, NAT;

PROCEDURE Signal(READONLY c: PZLR3Chain.T): REF PZLRChain.T =
  BEGIN
    WITH
      N = NUMBER(c),
      rs = NEW(REF PZLRChain.T, N), s = rs^,
      L = Length(c)
    DO
      
      PROCEDURE Convert(i, j: NAT; f: LONG) =
        (*
          Converts the curve "c[i..j]" to signal "s[i+1..j-1]", 
          assuming "s[i]" and "s[j]" have been defined.
        *)
        BEGIN
          IF i < j-1 THEN
            <* ASSERT (i-j) MOD 2 = 0 *>
            WITH 
              k = (i + j) DIV 2,
              h = - f * Angle(c[i], c[k], c[j])
            DO
              s[k] := (s[i] + s[j] + h)/2.0d0;
              Convert(i,k,f/2.0d0); Convert(k,j,f/2.0d0)
            END
          END
        END Convert;
      
      BEGIN
        s[0] := 0.0d0;
        s[N-1] := 0.0d0;
        Convert(0,N-1,L/2.0d0)
      END;
      RETURN rs
    END
  END Signal;
  
PROCEDURE Length(READONLY c: PZLR3Chain.T): LONG =
  (* Assumes the steps are all equal except for rounding error *)
  (* and turning bias: *)
  VAR s: LONG := 0.0d0;
  BEGIN
    FOR i := 1 TO LAST(c) DO
      WITH d2 = LR3.DistSqr(c[i], c[i-1]) DO s := s + d2 END;
    END;
    RETURN sqrt(s*FLOAT(NUMBER(c)-1, LONG))
  END Length;
  
PROCEDURE Angle(READONLY a, b, c: LR3.T): LONG =
  (*
    Computes the signed external angle between the vectors "b-a" and "c-b".
    Assumes the three points lie on the plane "Z=0".
  *)
  BEGIN
    <* ASSERT a[2] =0.0d0 *>
    <* ASSERT b[2] =0.0d0 *>
    <* ASSERT c[2] =0.0d0 *>
    WITH
      u = LR3.Sub(b,a),
      v = LR3.Sub(c,b),
      s = u[0]*v[1] - u[1]*v[0],
      c = u[0]*v[0] + u[1]*v[1]
    DO
      RETURN atan2(s,c)
    END
  END Angle;

PROCEDURE ComputeSpectrum(READONLY p: PZLRChain.T): REF Spectrum =
  BEGIN
    WITH 
      N = PowerOfTwoNotLessThan(NUMBER(p)-1)
    DO
      IF NUMBER(p) # N THEN
        <* ASSERT NUMBER(p) = N+1 *>
        <* ASSERT p[0] = p[LAST(p)] *>
      END;
      WITH
        q = NEW(REF PZLRChain.T, N)^,
        ra = NEW(REF Spectrum, N), a = ra^
      DO
        q := SUBARRAY(p, 0, N);
        PZFourier.FFT(q,a);
        RETURN ra
      END
    END
  END ComputeSpectrum;

PROCEDURE ComputePowerSpectrum(READONLY a: Spectrum): REF PowerSpectrum  =
  VAR w: LONG;
  BEGIN
    WITH
      N = NUMBER(a),
      fMax = N DIV 2,
      rpwr = NEW(REF PZLRChain.T, fMax+1), pwr = rpwr^
    DO
      <* ASSERT N MOD 2 = 0 *>
      FOR f := 0 TO fMax DO
        WITH t = a[f] DO w := t*t END;
        IF f > 0 AND f < fMax THEN
          WITH t = a[N-f] DO w := w + t*t END
        END;
        pwr[f] := w;
      END;
      RETURN rpwr;
    END
  END ComputePowerSpectrum;


PROCEDURE PowerOfTwoNotLessThan(n: NAT): NAT = 
  (* Smallest power of two not less than "n" *)
  VAR m: NAT := 1;
  BEGIN
    <* ASSERT n > 0 *>
    WHILE m < n DO m := m+m END;
    RETURN m
  END PowerOfTwoNotLessThan;

BEGIN
END PZShape.
(*
  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.
*)