MODULE PZCurvHistogram EXPORTS Main;

(* Generates histograms of numeric curvature values *)
(* Last edited on 1999-11-12 23:23:03 by hcgl *)

(*
  Reads a bunch of segments of numeric curvature chains, and computes
  their histograms.

  IMPORTANT: These segments must not overlap corners (blurred
  by the appropriate blurring length), nor any smmooth (non-fracture) 
  parts of the contours.
*)

IMPORT ParseParams, Process, FileRd, FileWr, Wr;
IMPORT OSError, Thread, Stdio, Fmt;
IMPORT PZProc, PZLRChain, PZSegment, PZHistogram;

FROM Math IMPORT sqrt, pow;
FROM Stdio IMPORT stderr;
FROM PZTypes IMPORT LONG, NAT;

<* FATAL Wr.Failure, Thread.Alerted, OSError.E *>

TYPE
  Options = RECORD
      input: TEXT;           (* Name of segment file, sans extension. *)
      maxSegs: NAT;          (* Number of segments to use. *)
      chainDir: TEXT;        (* Directory where chains can be found. *)
      chainPrefix: TEXT;     (* Chain file name prefix *)
      band: NAT;             (* Nominal band width (lambda) for file names. *)
      extension: TEXT;       (* Extension for chains (".fcv" usually) *)
      maxCurv: LONG;         (* Maximum curvature in histogram. *)
      (* Parameters for "squeezing" of curvature values: *)
      sigmaCurv: LONG;       (* Nominal standard dev. or curvature for band 1. *)
      fractalDim: LONG;      (* Nominal fractal dimension of curves. *)
      output: TEXT;          (* Output histogram file name, sans extension. *)
    END;

TYPE
  Segments = PZSegment.List;

PROCEDURE Main() =
  CONST NBins = 25; (* On each side of 0 *)
  BEGIN
    WITH
      o = GetOptions(),
      segData = PZSegment.Read(FileRd.Open(o.input & ".seg")),
      seg = SUBARRAY(segData.s^, 0, MIN(NUMBER(segData.s^), o.maxSegs)),
      curves = PZSegment.GetCurves(seg)^,
      chData = PZLRChain.ReadAll(o.chainPrefix, o.band, o.extension, 
        sel := curves, dir := o.chainDir, headerOnly := FALSE
      ),
      lambda = MAX(0.5d0, FLOAT(o.band, LONG)),
      plainHist = PZHistogram.New(-o.maxCurv, +o.maxCurv, 2*NBins + 1)^,
      squeezeHist = PZHistogram.New(-1.2d0, +1.2d0, 2*NBins + 1)^
    DO
      ComputeHistograms(chData.chain^, seg,
        lambda, o.maxCurv, o.sigmaCurv, o.fractalDim,
        plainHist, squeezeHist
      );
      OutputHistogram(o.output & "-plain", plainHist);
      OutputHistogram(o.output & "-squeezed", squeezeHist);
    END
  END Main;

PROCEDURE GetOptions(): Options =
  VAR o: Options;
  BEGIN
    WITH
      pp = NEW(ParseParams.T).init(stderr)
    DO
      TRY
        pp.getKeyword("-input");
        o.input := pp.getNext();
        
        IF pp.keywordPresent("-maxSegs") THEN
          o.maxSegs := pp.getNextInt(1,LAST(NAT))
        ELSE
          o.maxSegs := LAST(NAT)
        END;
        
        IF pp.keywordPresent("-chainDir") THEN
          o.chainDir := pp.getNext()
        ELSE
          o.chainDir := "."
        END;
        
        pp.getKeyword("-chainPrefix");
        o.chainPrefix := pp.getNext();
        
        pp.getKeyword("-band");
        o.band := pp.getNextInt();
          
        IF pp.keywordPresent("-extension") THEN 
          o.extension:= pp.getNext();
        ELSE 
          o.extension := ".fcv";
        END;
        
        pp.getKeyword("-maxCurv");
        o.maxCurv := pp.getNextLongReal(1.0d-10, 1.0d+10);
          
        pp.getKeyword("-sigmaCurv");
        o.sigmaCurv := pp.getNextLongReal(1.0d-10, 1.0d+10);
          
        pp.getKeyword("-fractalDim");
        o.fractalDim := pp.getNextLongReal(1.0d0, 2.0d0);
          
        pp.getKeyword("-output");
        o.output := pp.getNext();

        pp.finish();
      EXCEPT
      | ParseParams.Error =>
          Wr.PutText(stderr, "Usage: PZCurvHistogram \\\n");
          Wr.PutText(stderr, "  -input NAME \\\n");
          Wr.PutText(stderr, "  [ -chainDir DIR ] -chainPrefix NAME \\\n");
          Wr.PutText(stderr, "  -band NUMBER \\\n");
          Wr.PutText(stderr, "  [ -extension EXT ] \\\n");
          Wr.PutText(stderr, "  -maxCurv NUMBER \\\n");
          Wr.PutText(stderr, "  -sigmaCurv NUMBER -fractalDim NUMBER \\\n");
          Wr.PutText(stderr, "  -output NAME\n");
        Process.Exit(1);
      END;
    END;
    RETURN o
  END GetOptions;

PROCEDURE ComputeHistograms(
    READONLY cvd: ARRAY OF PZLRChain.ReadData;
    READONLY seg: Segments;   (* List of segments to digest. *)
    lambda: LONG;
    maxCurv: LONG;
    sigmaCurv: LONG;
    fractalDim: LONG;
    (*OUT*)
    VAR plainHist: PZHistogram.T;
    VAR squeezeHist: PZHistogram.T;
  ) =
  VAR sxx: LONG := 0.0d0;  (* Sum of squared samples *)
      ns: NAT := 0;        (* Number of samples *)
  BEGIN
    WITH
      sigma = sigmaCurv/pow(lambda, fractalDim)
    DO
      Wr.PutText(stderr, ": input segments = " & Fmt.Int(NUMBER(seg)) & "\n");
      Wr.PutText(stderr, ": lambda = " & FLR(lambda,7,2) & "\n");
      Wr.PutText(stderr, ": sigmaCurv = " & FLR(sigmaCurv,7,2) & "\n");
      Wr.PutText(stderr, ": fractalDim = " & FLR(fractalDim,7,2) & "\n");
      Wr.PutText(stderr, ": max plain curvature = " & FLR(maxCurv,14,6) & "\n");
      Wr.PutText(stderr, ": expected deviation = " & FLR(sigma,14,6) & "\n");
      FOR i := 0 TO LAST(seg) DO
        WITH
          si = seg[i],
          c = cvd[si.cvx].c^
        DO
          FOR i := 0 TO si.ns-1 DO
            WITH
              x = c[(si.ini + i) MOD si.tot],
              y = PZProc.Squeeze(x, sigma)
            DO
              sxx := sxx + x*x;
              INC(ns);
              PZHistogram.Add(plainHist, x);
              PZHistogram.Add(squeezeHist, y)
            END
          END
        END
      END;
      WITH sdev = sqrt(sxx/FLOAT(ns,LONG)) DO
        Wr.PutText(stderr, ": observed deviation =  " & FLR(sdev,14,6) & "\n");
        Wr.PutText(stderr, ": total samples = " & Fmt.Int(ns) & "\n");
      END;
    END
  END ComputeHistograms;

PROCEDURE OutputHistogram(output: TEXT; READONLY h: PZHistogram.T) =
  BEGIN
    WITH
      fileName = output & ".plt"
    DO
      Wr.PutText(stderr, "writing " & fileName & "\n");
      WITH
        wr = FileWr.Open(fileName)
      DO
        PZHistogram.Output(wr, h);
        Wr.Close(wr)
      END
    END
  END OutputHistogram;

PROCEDURE FLR(x: LONG; w, p: NAT): TEXT =
  BEGIN
    RETURN Fmt.Pad(Fmt.LongReal(x, Fmt.Style.Fix, prec := p), w)
  END FLR;

BEGIN
  Main()
END PZCurvHistogram.

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