/* Generates a plot to estimate the fractal dimension of a curve. */
/* Last edited on 2015-01-20 16:46:16 by stolfilocal */

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

#include <argparser.h>
#include <r3.h>

#include <pz_r3_chain.h>
#include <pz_types.h>


TYPE
  MismatchOptions == pz_mismatch.Options;

  Options = struct ??? {
    char *inFile;      /* Input curve file (without ".flc") */
    char *outFile;     /* Output plot data file (without ".plt") */
    double minDist; /* Minimum beeline distance to consider */
  } ???;
  
CONST 
  Pi == 3.14159265358979323e0;
  Sqrt2 == 1.4142135623730950488e0;
  Log2 == 0.69314718055994530941e0;


int main(int argc, char **argv 
  )      
  { 

    { /* with */ 
      ??? o = pz_get_options(int argc, char **argv);
      ??? fc = pz_r3_chain_read(open_read(o.inFile & ".flc", TRUE), recenter := pz_ctr_NONE);
      ??? chain = fc.c^;
      /* do */ 
      WriteLengthScalingPlot(open_write(o.outFile & ".plt", TRUE), chain, o.minDist);;
    };
  } /* 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, "-inFile");
        o.inFile = argparser_get_next(pp);
        
        argparser_get_keyword(pp, "-outFile");
        o.outFile = argparser_get_next(pp);
        
        argparser_get_keyword(pp, "-minDist");
        o.minDist = argparser_get_next_double(pp, 0.1e0, 100000.0e0);
        
        argparser_finish(pp);                                       
      EXCEPT                                                            
      | ParseParams.Error ==>                                              
          fprintf(stderr, "Usage: pz_filter \\\n");
          fprintf(stderr, "  -inFile NAME \\\n");
          fprintf(stderr, "  -outFile NAME \\\n");
          fprintf(stderr, "  -minDist: NUM\n");
        Process.Exit(1);;
      };;
    };
    return o
  } /* GetOptions */

  

void WriteLengthScalingPlot(
    FILE *wr, 
    pz_r3_chain_t *p, 
    double minDist 
  )
  VAR j, k: int;
      /* These arrays are indexed by "ROUND(log_2(dij/minDist)/2)": */
      tSum: ARRAY [0..60] OF double;
      dSum: ARRAY [0..60] OF double;
      sSum: ARRAY [0..60] OF double;
      nSteps: ARRAY [0..60] OF unsigned;
  {

    { /* with */ 
      ??? n = (p.ne);
      ??? t = double_vec_new(n+1);
      /* do */ 
      fprintf(stderr, "number of samples == %s\n", Fmt.Int(n) );
      
      /* Store in "t[i]" the length of the chain from "p[0]" to "p[i]". */
      t[0] = 0.00e0; 
      for (i = 1; i <= n-1 ; i++){ 
        t[i] = t[i-1] + Dist(p[i-1], p[i]);
      };
      t[n] = t[n-1] + Dist(p[n-1], p[0]);
      fprintf(stderr, "Total length == %s pixels\n", Fmt.LongReal(t[n]) );
      
      /* Clear accumulators: */
      for (r = 0; r < (tSum.ne ) ; r++){ 
        tSum[r] = 0.0e0; dSum[r] = 0.0e0; sSum[r] = 0.0e0;
        nSteps[r] = 0;
      };
      
      /* Now generate a plot of "log(L)/log(D)": */
      k = 2;
      while ( k < n DIV 2 ){ 
        for (i = 0; i <= n - 1 BY k DIV 2 ; i++){ 
          j= i + k; 

          { /* with */ 
            /* Beeline distance (chord): */
            ??? dij = Dist(p[i], p[j MOD n]);
            
            /* Length along the chain: */
            ??? tij = t[j MOD n] + t[n]*((double)j DIV n) - t[i];
            
            /* Half of angle subtended by "tij" on a circle of length "t[n]": */
            ??? theta2 = M_PI * tij / t[n];
            
            /* Length of arc of circle with same angle and chord "dij": */
            ??? sij = dij * theta2 / ABS(sin(theta2));
            /* do */
            if (( tij > minDist/Sqrt2 ) ANDAND ( dij > minDist/Sqrt2 )){ 

              { /* with */
                ??? r = MAX(0, ROUND(2.0e0*log(dij/minDist)/Log2));
                /* do */
                tSum[r] = tSum[r] + tij;
                dSum[r] = dSum[r] + dij;
                sSum[r] = sSum[r] + sij;
                INC(nSteps[r]);
              };
            };
          };
        };
        k = max(k+1, ROUND(((double)k)*Sqrt2));
      };
      
      /* Now output the plot: */
      for (r = 0; r < (tSum.ne ) ; r++){ 
        if (( nSteps[r] >= 2 )){

          { /* with */
            ??? ns = ((double)nSteps[r]);
            ??? tAvg = tSum[r]/ns;
            ??? dAvg = dSum[r]/ns;
            ??? sAvg = sSum[r]/ns;
            ??? dDim = log(tSum[r])/log(dSum[r]);
            ??? sDim = log(tSum[r])/log(sSum[r]);
            /* do */
            fprintf(wr, FLR(tAvg, 7, 2));
            fprintf(wr, " ");
            fprintf(wr, FLR(dAvg, 7, 2)); 
            fprintf(wr, " ");
            fprintf(wr, FLR(sAvg, 7, 2)); 
            fprintf(wr, " ");
            fprintf(wr, FLR(dDim, 7, 4));
            fprintf(wr, " ");
            fprintf(wr, FLR(sDim, 7, 4));
            fprintf(wr, "\n");
            fflush(wr);;
          };
        };
      };
    }; 
  } /* WriteLengthScalingPlot */
  

char *FLR(
    double x, 
    unsigned w,
    unsigned p 
  )
  {
    return Fmt.Pad(Fmt.LongReal(x, style = 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.
*/