/*
  st-steps - finds horizontal intensity steps in a PGM image
  
  usage: st-steps PARMSFILE.parms [ < ] IMGFILE.pgm > ISTEPS.ist
  
  Format of ISTEPS.ist lines is:
  
    x y p1 p2
    
  where 
  
    x, y are the image coordinates of the intensity step,
    p1, p2 are the pixel values at the beginning and end of the step.
    
*/

#include <sys/time.h>
#include <sys/resource.h>
#include <stdio.h>
#include <string.h>
#include <pgm.h>
#include <st-types.h>
#include <st-parms.h>

/* INTERNAL PROTOTYPES */

int main(int argc, char* argv[]);
static void parse_options(
  int *argc, char **argv,
  char **parmsnameP, char **imgnameP
);
static FILE *open_rd(char *name);
static void writeistep(FILE *f, float y, IStep *ist);
int istep_value(int x_beg, int x_end, gray p_beg, gray p_end, int max_step_w);
float istep_middle(int r, int s, gray row[]);
void find_isteps(
    gray row[], int cols,
    Parms *p, int maxval, 
    IStep ist[], int *nstepsP
  );

/* ROUTINES */
 
int main(int argc, char **argv)
  { 
    char *imgname, *parmsname;
    FILE *imgfile;
    int cols, rows;
    gray maxval;
    int format;
    int nsteps;
    gray *row;
    IStep *ist;
    Parms *p;
    int r, i, j; 
    float y, y_top;

    FILE *istfile; /*DEBUG */
    gray *srow; /*DEBUG */

    parse_options(&argc, argv, &parmsname, &imgname);
    p = readparms(parmsname);
    imgfile = (imgname == NULL ? stdin : open_rd(imgname));
    pgm_readpgminit(imgfile, &cols, &rows, &maxval, &format);
    y_top = ((float)(rows-1))/2.0;
    row = (gray*)(malloc(cols*sizeof(gray)));
    if (row == NULL) pm_error("out of memory");

    istfile = fopen("debug.pgm", "w"); /*DEBUG */
    pgm_writepgminit(istfile, cols, rows, maxval, 0); /*DEBUG */
    srow = (gray*)(malloc(cols*sizeof(gray))); /*DEBUG */
    if (srow == NULL) pm_error("out of memory"); /*DEBUG */
    
    ist = (IStep*)(malloc(p->max_steps*sizeof(IStep)));
    if (ist == NULL) pm_error("out of memory");
    for (r = 0; r < rows; r++)
      { y = y_top - (float)r;
        pgm_readpgmrow(imgfile, row, cols, maxval, format);
        find_isteps(row, cols, p, maxval, ist, &nsteps);
        for (j = 0; j < cols; j++) srow[j] = maxval; /*DEBUG */
        for (i = 0; i < nsteps; i++)
          { writeistep(stdout, y, &(ist[i]));
            for (j = ist[i].x_beg; j < ist[i].x_end; j++) srow[j] = row[j]; /*DEBUG */
          }
        fflush(stdout);
        pgm_writepgmrow(istfile, srow, cols, maxval, 0); /*DEBUG */
      }
    fclose(imgfile);
    fclose(stdout);
    fclose(istfile); /*DEBUG */
    return(0);
  }
  
static void writeistep(FILE *f, float y, IStep *ist)
  {
    fprintf(f, "%6.2f %6.2f %5.3f %5.3f\n", 
      ist->x_mid, y, ist->p_beg, ist->p_end
    );
  }

int istep_value(int x_beg, int x_end, gray p_beg, gray p_end, int max_step_w)
  {
    return(abs(p_end - p_beg) * (max_step_w - (x_end - x_beg)));
  }

float istep_middle(int r, int s, gray row[])
  {
    return (((float)(r+s))/2.0);
  }
  
void find_isteps(
    gray row[], int cols,
    Parms *p, int maxval, 
    IStep ist[], int *nstepsP
  )
  {
    int ix;
    float x_left =  -((float)(cols-1))/2.0;
    float pscale = 1.0/(float)maxval;

    float p1, p2 = 0.0;
    IStep s;            /* last gap (if n_runs > 0), possibly not definitive; */
    int run_beg;        /* index of first pixel of current run */
    int run_length = 0; /* current width of current run */
    int n_runs = 0;     /* number of significant runs before current one */
    int n = 0;          /* number of complete gaps found so far  */

    n = 0;
    n_runs = 0;
    run_beg = 0;
    run_length = 1;
    for(ix = 1; ix <= cols; ix++)
      {
        /* 
          Invariant: 
            n_runs is the number of complete significant runs before this one;
            run_length is the current length of the current run;
            run_beg is the end of the last gap/beginning of current run;
            s.x_beg (if n_runs > 0) is the beginning of the last gap;
            s.p_beg is its intensity;
            s.x_end (if n_runs > 0) is the end of the last gap, same as run_beg;
            s.p_end is its intensity;
            ist[0..n-1] are the complete gaps found, sorted by decreasing val.
        */
	p1 = pscale * row[ix-1];
	if (ix < cols) p2 = pscale * row[ix];
	if ((ix < cols) && (abs(p2-p1) <= 2*p->delta))
	  { run_length++; }
	else
	  { /* run ended: */
            if (run_length >= p->min_run_w)
              { /* this was a significant run */
                if ((n_runs > 0) && (run_beg - s.x_beg < p->max_step_w))
                  { /* save previous gap, if good enough: */
                    int j;
                    s.x_mid = istep_middle(s.x_beg, s.x_end, row) + x_left;
                    s.val = istep_value(s.x_beg, s.x_end, s.p_beg, s.p_end, p->max_step_w);
                    for (j = n; (j > 0) && (ist[j-1].val < s.val); j--)
                      { if (j < p->max_steps) ist[j] = ist[j-1]; }
                    if (j < p->max_steps) { ist[j] = s; }
                    if (n < p->max_steps) n++;
                  }
                /* new gap started: */
                n_runs++;
                s.x_beg = ix-1;
                s.p_beg = p1;
              }
            /* new run started: */
            run_beg = ix;
            run_length = 1;
            /* gap provisionally ended here: */
            s.x_end = ix;
            s.p_end = p2;
          }
      }
    (*nstepsP) = n;
  }

static FILE*
open_rd(char *name)
  {
    FILE *f = fopen(name, "r");
    if (f == NULL) pm_error("couldn't open input file");
    return (f);
  }

static void 
parse_options(
    int *argc, char **argv,
    char **parmsnameP, 
    char **imgnameP
  )
  {
    int argn;
    char* usage = " PARMSFILE [ PGMFILE ]";
    pgm_init(argc, argv);
    argn = 1;
    (*imgnameP) = NULL;
    (*parmsnameP) = NULL;
    
    if (argn >= (*argc)) pm_usage(usage);
    (*parmsnameP) = argv[argn];
    ++argn;
    if (argn != (*argc))
      { (*imgnameP) = argv[argn]; argn++; }
    if (argn != (*argc)) pm_usage(usage);
  }