/* Generates histograms of numeric curvature values */
/* Last edited on 2015-01-20 16:44:34 by stolfilocal */

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

  IMPORTANT: For best results
IMPORT these segments must not overlap corners
  (even after the corners have been blurred by the appropriate
  blurring length), nor any smmooth (non-fracture) parts of the
  contours.
*/

#include <stdio.h>
#include <math.h>

#include <argparser.h>

#include <pz_proc.h>
#include <pz_double_chain.h>
#include <pz_segment.h>
#include <pz_histogram.h>
#include <pz_types.h>


TYPE

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

TYPE
  Segments == pz_segment_vec_t;


int main(int argc, char **argv 
  )
  CONST NBins == 25; /* On each side of 0 */
  {

    { /* with */
      ??? o = pz_get_options(int argc, char **argv);
      ??? segData = pz_segment_read(open_read(o.input & ".seg", TRUE));
      ??? seg = SUBARRAY(segData.s^, 0, MIN((segData.s^.ne), o.maxSegs));
      ??? chainUsed = pz_segment_chains_used(seg)^;
      ??? chAllData = pz_double_chain_read_all(o.chainPrefix, o.band, o.extension;
        sel := chainUsed, dir := o.chainDir, header_only := FALSE
      );
      ??? lambda = MAX(0.5e0, ((double)o.band));
      ??? plainHist = pz_histogram.New(-o.maxCurv, +o.maxCurv, 2*NBins + 1)^;
      ??? squeezeHist = pz_histogram.New(-1.2e0, +1.2e0, 2*NBins + 1)^;
      /* do */
      ComputeHistograms(chAllData.chData^, seg,
        lambda, o.maxCurv, o.sigmaCurv, o.fractalDim,
        plainHist, squeezeHist
      );
      OutputHistogram(o.output & "-plain", plainHist);
      OutputHistogram(o.output & "-squeezed", squeezeHist);;
    }
  } /* Main */


Options  pz_get_options(int argc, char **argv 
  )
  VAR o: Options;
  {
    argparser_t *pp = argparser_new(stderr, argc, argv);
    argparser_set_help(pp, PROG_NAME " version " PROG_VERS ", usage:\n" PROG_HELP);
    argparser_set_info(pp, PROG_INFO);
    argparser_process_help_info_options(pp);
    

    { /* with */
      /* do */
      TRY
        argparser_get_keyword(pp, "-input");
        o.input = argparser_get_next(pp);
        
        if (( argparser_keyword_present(pp, "-maxSegs") )){
          o.maxSegs = argparser_get_next_int(pp, 1,(unsigned.ne - 1))
        }else{
          o.maxSegs = (unsigned.ne - 1);
        };
        
        if (( argparser_keyword_present(pp, "-chainDir") )){
          o.chainDir = argparser_get_next(pp)
        }else{
          o.chainDir = ".";
        };
        
        argparser_get_keyword(pp, "-chainPrefix");
        o.chainPrefix = argparser_get_next(pp);
        
        argparser_get_keyword(pp, "-band");
        o.band = argparser_get_next_int(pp);
          
        if (( argparser_keyword_present(pp, "-extension") )){ 
          o.extension= argparser_get_next(pp);
        }else{ 
          o.extension = ".fcv";;
        };
        
        argparser_get_keyword(pp, "-maxCurv");
        o.maxCurv = argparser_get_next_double(pp, 1.0e-10, 1.0e+10);
          
        argparser_get_keyword(pp, "-sigmaCurv");
        o.sigmaCurv = argparser_get_next_double(pp, 1.0e-10, 1.0e+10);
          
        argparser_get_keyword(pp, "-fractalDim");
        o.fractalDim = argparser_get_next_double(pp, 1.0e0, 2.0e0);
          
        argparser_get_keyword(pp, "-output");
        o.output = argparser_get_next(pp);

        argparser_finish(pp);
      EXCEPT
      | ParseParams.Error ==>
          fprintf(stderr, "Usage: pz_curv_histogram \\\n");
          fprintf(stderr, "  -input NAME \\\n");
          fprintf(stderr, "  [ -chainDir DIR ] -chainPrefix NAME \\\n");
          fprintf(stderr, "  -band NUMBER \\\n");
          fprintf(stderr, "  [ -extension EXT ] \\\n");
          fprintf(stderr, "  -maxCurv NUMBER \\\n");
          fprintf(stderr, "  -sigmaCurv NUMBER -fractalDim NUMBER \\\n");
          fprintf(stderr, "  -output NAME\n");
        Process.Exit(1);;
      };;
    };
    return o
  } /* GetOptions */


void ComputeHistograms(
    ARRAY *OF pz_double_chain_read_data cvd,
    Segments *seg,   /* List of segments to digest. */
    double lambda,
    double maxCurv,
    double sigmaCurv,
    double fractalDim,
     /* (OUT) */
    pz_histogram.T *plainHist,
    pz_histogram.T *squeezeHist 
  )
  VAR sxx: double := 0.0e0;  /* Sum of squared samples */
      ns: unsigned = 0;        /* Number of samples */
  {

    { /* with */
      ??? sigma = sigmaCurv/pow(lambda, fractalDim);
      /* do */
      fprintf(stderr, ": input segments == %s\n", Fmt.Int((seg.ne)) );
      fprintf(stderr, ": lambda == %s\n", FLR(lambda,7,2) );
      fprintf(stderr, ": sigmaCurv == %s\n", FLR(sigmaCurv,7,2) );
      fprintf(stderr, ": fractalDim == %s\n", FLR(fractalDim,7,2) );
      fprintf(stderr, ": max plain curvature == %s\n", FLR(maxCurv,14,6) );
      fprintf(stderr, ": expected deviation == %s\n", FLR(sigma,14,6) );
      for (i = 0; i < (seg.ne ) ; i++){
        affirm(cvd[seg[i].cvx].c != NULL );

        { /* with */
          ??? si = seg[i];
          ??? c = cvd[si.cvx].c^;
          /* do */
          for (i = 0; i <= si.ns-1 ; i++){

            { /* with */
              ??? x = c[(si.ini + i) MOD si.tot];
              ??? y = pz_proc.Squeeze(x, sigma);
              /* do */
              sxx = sxx + x*x;
              ns++;
              pz_histogram.Add(plainHist, x);
              pz_histogram.Add(squeezeHist, y);
            };
          };
        };
      };

      { /* with */
        ??? sdev = sqrt(sxx/FLOAT(ns,double));
        /* do */
        fprintf(stderr, ": observed deviation ==  %s\n", FLR(sdev,14,6) );
        fprintf(stderr, ": total samples == %s\n", Fmt.Int(ns) );;
      };;
    }
  } /* ComputeHistograms */


void OutputHistogram(
    char *output, 
    pz_histogram.T *h 
  )
  {

    { /* with */
      ??? fileName = txtcat(output , ".plt");
      /* do */
      fprintf(stderr, "writing %s\n", fileName );

      { /* with */
        ??? wr = open_write(fileName, TRUE);
        /* do */
        pz_histogram.Output(wr, h);
        if (wr == stdout) { fflush(wr); } else { fclose(wr); }
      };
    }
  } /* OutputHistogram */


char *FLR(
    double x, 
    unsigned w,
    unsigned p 
  )
  {
    return Fmt.Pad(Fmt.LongReal(x, Fmt.Style.Fix, prec = p), w)
  } /* FLR */

{

/*
  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.
*/