/* Generate list of straight segments in symbolic chains */
/* Last edited on 2015-01-20 16:46:31 by stolfilocal */

/* 
  Looks for contour segments that have low average curvature
  in the current scale, and a specified minimum length.
  
  These segments are broadened to include the adjacent blurred
  corners.  Remember to use a negative "blurFactor" when mapping these
  segments to other scales.
*/
#include <stdio.h>

#include <argparser.h>

#include <pz_symbol_chain.h>
#include <pz_symbol.h>
#include <pz_segment.h>
#include <pz_types.h>


CONST
  ExpectedSegsPerChain == 1;

TYPE

  Options = struct ??? {
    char *chainDir;         /* Directory where to find chain data */
    char *chainPrefix;      /* Invariant chain file name prefix. */
    unsigned nChains;           /* Number of chains. */
    char *output;           /* Output segment file name. */
    unsigned band;              /* Filtering "band" (wavelength parameter). */
    double step;             /* Sampling step size. */
    double minLength;        /* Min length of straight segments. */
    double minCurvature;     /* Min curvature to consider pairing. */
    double blurFactor;       /* Corner broadening factor. */
    unsigned printSegs; /* Max segments to print. */
  } ???;

TYPE
  VTable == double_vec_t;


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

    { /* with */
      ??? o = pz_get_options(int argc, char **argv);
      ??? lambda = ((double)o.band);
      ??? sel = SelectAll(o.nChains)^;
      ??? chAllData = pz_symbol_chain_read_all(o.chainPrefix, o.band, ".cvc", sel;
        header_only := FALSE, dir := o.chainDir
      );
      ??? dch = pz_symbol_make_decode_table(chAllData.sigmaMin)^;
      ??? ech = pz_symbol_make_error_var_table(chAllData.sigmaMin)^;
      seg = ComputeStraightSegments(
        chAllData.chData^;
        dch := dch;
        ech := ech;
        minCurvature := o.minCurvature;
        minLength := EnsurePositive(o.minLength - 2.0e0 * o.blurFactor * lambda);
        step := o.step;
        /* Extend the result to cover the blurred corners too: */
        extraLength := 4.0e0 * o.blurFactor * lambda
      )^;
      /* do */
      affirm(chAllData.sigmaMin == chAllData.sigmaMax );

      { /* with */
        ??? wr = open_write(o.output & ".seg", TRUE);
        /* do */
        pz_segment_write(wr, SegComment(o), seg);
      };
      PrintSegments(chAllData.chData^, seg, o.band, o.printSegs);;
    };
  } /* Main */
  

double  EnsurePositive(
    double segLength 
  )
  {
    if (( segLength < 0.0e0 )){
      fprintf(stderr, "warning: requested segment length too small for this scale\n");
      fflush(stderr);
      Process.Exit(1);;
    };
    return max(0.0e0, segLength)
  } /* EnsurePositive */


bool_vec_t *SelectAll(
    unsigned n 
  ) 
  {

    { /* with */
      ??? rb = bool_vec_new(n), b = rb^;
      /* do */
      for (i = 0; i <= n-1 ; i++){ b[i] = TRUE ;};
      return rb;
    }
  } /* SelectAll */


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
        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, "-nChains");
        o.nChains = argparser_get_next_int(pp);

        argparser_get_keyword(pp, "-band");
        o.band = argparser_get_next_int(pp);

        argparser_get_keyword(pp, "-step");
        o.step = argparser_get_next_double(pp, 0.0e0);

        argparser_get_keyword(pp, "-output");
        o.output = argparser_get_next(pp);

        argparser_get_keyword(pp, "-minLength");
        o.minLength = argparser_get_next_double(pp, 0.0e0);

        argparser_get_keyword(pp, "-blurFactor");
        o.blurFactor = argparser_get_next_double(pp, 0.0e0);

        argparser_get_keyword(pp, "-minCurvature");
        o.minCurvature = argparser_get_next_double(pp, 0.0e0);

        if (( argparser_keyword_present(pp, "-printSegs") )){
          o.printSegs = argparser_get_next_int(pp, 1)
        }else{
          o.printSegs = 1000;
        };

        argparser_finish(pp);
      EXCEPT
      | ParseParams.Error ==>
          fprintf(stderr, "Usage: pz_get_straight_segs \\\n");
          fprintf(stderr, "  [ -chainDir DIR ] -chainPrefix NAME \\\n");
          fprintf(stderr, "  -nChains NUMBER \\\n");
          fprintf(stderr, "  -band NUMBER -step NUMBER \\\n");
          fprintf(stderr, "  -output NAME \\\n");
          fprintf(stderr, "  -minLength NUMBER  -blurFactor NUMBER \\\n");
          fprintf(stderr, "  -minCurvature NUMBER \\\n");
          fprintf(stderr, "  [ -printSegs NUMBER ]\n");
        Process.Exit(1);;
      };;
    };
    return o
  } /* GetOptions */


PROCEDURE ComputeStraightSegments(
    ARRAY *OF pz_symbol_chain_read_data chain,
    VTable *dch,              /* Curvature decoding table. */
    VTable *ech,              /* Quantization noise variances. */
    double minCurvature,      /* Min curvature to consider */
    double minLength,         /* Min length of straight segment (
    already shrunk 
  ) */
    double step,              /* Sampling step */
    extraLength: double;       /* How much to extend the result */
  ): REF pz_segment_vec_t ==
  VAR
    s: REF pz_segment_vec_t;
    nSegs : unsigned = 0;
  {

    { /* with */
      ??? nChains = (chain.ne);
      /* do */
      s = NEW(REF pz_segment_vec_t, nChains * ExpectedSegsPerChain);
      for (k = 0; k <= nChains-1 ; k++){

        { /* with */
          ??? c = chain[k].c^;
          ??? minSteps = FLOOR(minLength / step);
          ??? halfStretch = FLOOR(0.5e0 * extraLength / step );
          /* do */
          ProcessChain(
            c, k,
            s, nSegs,
            dch = dch,
            ech = ech,
            minSteps = minSteps,
            minCurvature = minCurvature,
            halfStretch = halfStretch
          );
        };;
      };

      { /* with */
        ??? st = pz_segment_vec_new(nSegs);
        /* do */
        st^ =  SUBARRAY(s^, 0, nSegs);
        return st;;
      };;
    };
  } /* ComputeStraightSegments */


void ProcessChain(
    pz_symbol_chain_t *c,
    unsigned k,
    REF *pz_segment_vec_t s,
    unsigned *nSegs,
    VTable *dch        /* Curvature decoding table */
    VTable *ech,       /* Quantization noise variances. */
    unsigned minSteps,              /* Minimum number of steps in segment. */
    double minCurvature,         /* Min curvature to consider */
    unsigned halfStretch, /* Number of steps to add at each end */
  )
  /*
    Consider the set of indices "ic" of the chain "c".
    This set is a priodic sequence with period "mc==(c.ne)".

    We compute a boolean vector "ok", such that "ok[ic]" tells whether
    the the "minSteps+1" elements beginning with element "r" would be
    an acceptable segment. Specifically, "ok[ic] == TRUE" iff the root
    mean of the squared curvatures "c[ic+k]^2" is less than
    "minCurvature".  All averages are taken for "k" ranging in
    [0..minSteps]".

    The segments are then found by identifying the maximal substrings
    of consecutive "TRUE"s in the "ok" vector.  Finally, "halfStretch"
    steps are added at each end of every segment found.
  */
  VAR c2Sum: double;
  {

    { /* with */
      ??? mc = (c.ne);
      ??? nOld = nSegs + 0;
      /* Properties of a minimal-size segment: */
      ??? minSamples = ((double)minSteps + 1);
      ??? minTotValueSqr = minSamples * minCurvature*minCurvature;
      /* Work vectors: */
      ??? c2s = double_vec_new(mc)^;  /* Accumulated "c[ia]^2" */
      ??? ok = bool_vec_new(mc)^;
      /* do */
      fprintf(stderr, "chain: " & Fmt.Int(k) & "[" & Fmt.Int(mc) & "]");
      fprintf(stderr, " minSteps: " & Fmt.Int(minSteps));
      /*
        Accumulates squared values,
        so that "c2s[r] == SUM { c[k]^2 : k IN [0..r] }".
      */
      c2Sum = 0.0e0;
      for (r = 0; r <= mc-1 ; r++){

        { /* with */
          ??? cv = dch[c[r]], e2 = ech[c[r]];
          /* do */
          c2Sum = c2Sum + cv*cv + e2; c2s[r] = c2Sum;;
        };;
      };
      /*
        Mark in "ok" all positions "r" along the diagonal
        that are the beginnings of substrings, consisting OF
        "minSteps+1" consecutive samples, which have
        sufficiently small root mean square curvatures:
      */
      for (r = 0; r <= mc-1 ; r++){

        { /* with */
          ??? rLo = r - 1;
          ??? rmLo = rLo MOD mc;
          ??? qLo = ((double)rLo DIV mc);

          ??? rHi = r + minSteps;
          ??? rmHi = rHi MOD mc;
          ??? qHi = ((double)rHi DIV mc);

          ??? c2Tot = (c2s[rmHi] + qHi*c2Sum) - (c2s[rmLo] + qLo*c2Sum);
          /* do */
          ok[r] = c2Tot < minTotValueSqr;
        };;
      };
      BreakChain(
        c, k,
        ok,
        s, nSegs,
        dch = dch, ech = ech,
        minSteps = minSteps,
        minCurvature = minCurvature, 
        halfStretch = halfStretch
      );
      fprintf(stderr,"  nSegs: " & Fmt.Int(nSegs));
      fprintf(stderr, " (+%s)\n", Fmt.Int(nSegs - nOld) );;
    };
  } /* ProcessChain */


void BreakChain(
    pz_symbol_chain_t *c,
    unsigned k,
    bool_vec_t *ok,
    REF *pz_segment_vec_t s,
    unsigned *nSegs,
    VTable *dch        /* Curvature decoding table */
    VTable *ech,       /* Quantization noise variances. */
    unsigned minSteps,              /* Minimum number of steps in segment. */
    double minCurvature,         /* Min curvature to consider */
    unsigned halfStretch, /* Number of steps to add at each end */
  )
  VAR
    t: unsigned;
    count, tStart, tIni, tLen: unsigned;
    inx: unsigned;
  {

    { /* with */
      ??? mc = (c.ne);
      ??? maxSteps = mc-1;
      /* do */
      affirm(mc > 1 );
      /* Look for a TRUE element "ok[tStart]" right after a FALSE element: */
      t = 0;
      while ( ok[t MOD mc] == TRUE ) ANDAND ( t < mc ){ t++ ;};
      if (( t == mc )){
        /* All TRUE */
        tStart = 0
      }else{
        inx = 0;
        while ( ok[t MOD mc] == FALSE ) ANDAND ( inx < mc ){ t++; inx++ ;};
        if (( inx == mc )){ /* All FALSE */ return ;};
        tStart = t MOD mc;;
      };
      /* Collect runs of TRUEs": */
      t = tStart;
      REPEAT
        /* Here begins a new run */
        tIni = t;
        count = 0;
        REPEAT
          t++;
          count++;
        UNTIL ok[t MOD mc] == FALSE ) || ( count == mc;

        /* Compute end of segment, clipping to "maxSteps": */

        { /* with */
          ??? nSteps = count - 1 + minSteps;
          ??? excess = MAX(0, nSteps - maxSteps);
          /* do */
          tIni = (tIni + excess DIV 2) MOD mc;
          tLen = min(maxSteps, nSteps) + 1;;
        };

        /* Remove high-curvature points at either end: */
        TrimCrookedEnds(c, minCurvature*minCurvature, dch, ech, tIni, tLen);

        /* Extend segment in both directions, clipping to "maxSteps": */
        tIni = (tIni - halfStretch) MOD mc;
        tLen = tLen + 2*halfStretch;

        { /* with */
          ??? excess = MAX(0, tLen - 1 - maxSteps);
          /* do */
          tIni = (tIni + excess DIV 2) MOD mc;
          tLen = min(maxSteps + 1, tLen);;
        };

        /* Finally, store segment if long enough: */
        if (( tLen - 1 >= minSteps )){

          /* Make room for a new segment: */
          if (( nSegs >= (s^.ne) )){

            { /* with */
              ??? st = NEW(REF pz_segment_vec_t, nSegs * 2);
              /* do */
              SUBARRAY(st^, 0, (s^.ne)) = s^;
              s = st;;
            };;
          };

          /* Save this run as a new segment: */
          s[nSegs] = pz_segment_t{
            cvx = k,
            tot = mc,
            ini = tIni MOD mc,
            ns = tLen,
            rev = FALSE
          };
          nSegs++;
        };
        while ( ok[t MOD mc] == FALSE ){ t++ ;};
      UNTIL (t MOD mc) == tStart;
    };
  } /* BreakChain */


void TrimCrookedEnds(
    pz_symbol_chain_t *c,
    double minCurvSqr,           /* Min curvature to be non-straight, squared. */
    VTable *dch,       /* Curvature decoding table */
    VTable *ech,       /* Quantization noise variances. */
    VAR /*IO*/ tIni, tLen: unsigned;
  )
  VAR tFin := tIni + tLen - 1;
  {

    { /* with */
      ??? mc = (c.ne);
      /* do */
      while ( tIni < tFin ){

        { /* with */
          ??? ic = tIni MOD mc;
          ??? cvi = dch[c[ic]], e2i = ech[c[ic]];

          ??? fc = tFin MOD mc;
          ??? cvf = dch[c[fc]], e2f = ech[c[fc]];
          /* do */
          if (( cvi*cvi + e2i >= minCurvSqr )){
            tIni++
          }else if (( cvf*cvf + e2f >= minCurvSqr )){
            tFin--
          }else{
            tLen = tFin - tIni + 1;
            return;
          };
        };
      };
    }
  } /* TrimCrookedEnds */


void PrintSegments (
    ARRAY *OF pz_symbol_chain_read_data chain,
    pz_segment_vec_t *seg,
    unsigned band,
    unsigned maxSegs 
  )
  {
    if (( maxSegs == 0 )){ return ;};
    fprintf(stderr, "===== band == %s=====\n", Fmt.Int(band) );
    for (i = 0; i <= min(maxSegs, (seg.ne))-1 ; i++){

      { /* with */
        ??? s = seg[i];
        /* do */
        fprintf(stderr, "seg[%s]\n", Fmt.Int(i) );
        pz_segment_print(stderr, s);

        { /* with */
          ??? ch = chain[s.cvx].c;
          /* do */
          if (( ch != NULL )){
            affirm((ch^.ne) == s.tot );
            pz_symbol_chain_print_segment(stderr, ch^, s);
          };
        };
      };
    };
    fflush(stderr);
  } /* PrintSegments */


char *SegComment(
    Options *o 
  )
  {
    return
      "pz_get_straight_segs\n" &
      "    chainPrefix " & o.chainPrefix & "\n" &
      "    nChains == " & Fmt.Int(o.nChains) & "\n" &
      "    band == " & Fmt.Int(o.band) & "\n" &
      "    minCurvature == " & Fmt.LongReal(o.minCurvature) & "\n" &
      "    minLength == " & Fmt.LongReal(o.minLength) & "\n" &
      "    blurFactor == " & Fmt.LongReal(o.blurFactor) & "\n"
  } /* SegComment */

{

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