/* @(#)nechough.c	17.1.1.1 (ES0-DMD) 01/25/02 17:51:29 */
/*===========================================================================
  Copyright (C) 1995 European Southern Observatory (ESO)
 
  This program is free software; you can redistribute it and/or 
  modify it under the terms of the GNU General Public License as 
  published by the Free Software Foundation; either version 2 of 
  the License, or (at your option) any later version.
 
  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
 
  You should have received a copy of the GNU General Public 
  License along with this program; if not, write to the Free 
  Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
  MA 02139, USA.
 
  Corresponding concerning ESO-MIDAS should be addressed as follows:
	Internet e-mail: midas@eso.org
	Postal address: European Southern Observatory
			Data Management Division 
			Karl-Schwarzschild-Strasse 2
			D 85748 Garching bei Muenchen 
			GERMANY
===========================================================================*/

/* Program : hough.c                                              */
/* Author  : P. Ballester   -     ESO Garching                    */
/* Date    : 20.11.90      Version 1.0                            */
/*           22.03.91      Version 2.0                            */
/*           27.03.91      Version 3.0                            */
/*                                                                */
/* Purpose :                                                      */
/*                                                                */
/* a) Apply Hough Transform on the central region of an image     */
/*                                                                */
/* Input:                                                         */  
/*             - name of input image            : IN_A            */
/*             - name of output image           : IN_B            */
/*             - Inter-trace width              : INPUTI(1)       */
/*             - Number of perpendicular traces : INPUTI(2)       */
/*             - Number of columns of Hough Tr. : INPUTI(3)       */
/*             - Number of rows of Hough Tr.    : INPUTI(4)       */
/*             - Lower scan limit               : INPUTI(5)       */
/*             - Upper scan limit               : INPUTI(6)       */
/*             - Start of Hough Transform       : INPUTD(1),(2)   */
/*             - Step of Hough Transform        : INPUTD(3),(4)   */
/*                                                                */
/* Limits:                                                        */
/*                                                                */
/*  a) This algorithm must process only positive values.          */
/*  b) The row 0 of the transformed space is not computed         */
/*                                                                */
/* Algorithm:
 *
 *   a) Only a certain number of columns of the image are processed.
 *   Given the distance between two columns and the number of columns
 *   to process (INPUTI(1) and (2)), the positions of the different
 *   columns are given by the function icol(ncol,trace,ntrace,inter).
 *
 *   b) The image is first filtered and a constant value (background)
 *   is subtracted. The filtering consists of a median estimate over a
 *   kernel 3*5 pixels. The constant value corresponds to the minimum
 *   of all median estimates. The result is stored in a buffer of size
 *   ntrace * nrow (number of columns processed, number of rows of the
 *   original image). The first and last row (1 and n) of the filtered 
 *   image are copies of the adjacent rows (2 and n-1).
 *
 *   c) The Hough Transform is performed, in which the incremented value
 *   corresponds to the value of the pixel in the processed image. For
 *   more details on the Hough transform, see P. Ballester, 1991. Proceedings
 *   of the 3rd Data Analysis Workshop.
 */

#include <math.h>
#include <midas_def.h>
#include <osparms.h>
#include <time.h>

#include <hough.h>
#define icol(ncol,trace,ntrace,inter)  ncol/2.-0.5+(trace-((ntrace+1.)/2.))*inter

int  step_prgs=20, next_prgs=20;

main ()
{

    char           frame[TEXT_LEN+1], transf[TEXT_LEN+1], line[TEXT_LEN];
    char           ident[TEXT_LEN+1] , cunit[16*3 + 1];
    float          *pntra, *pntrb;
    int            imnoa, imnob, naxis, npix[2], npix_hg[2];
    double         start[2] , step[2], start_hg[2], step_hg[2];
    float          high_thres;
    int            null, status, npar, iav, actvals, kunit;
    int            inpi[2], inter, ntrace, scan[2];
    float          ctime;
    unsigned int   nobyt;
    float          *buffer;

    SCSPRO ("hough");

    if (SCKRDI ("INPUTI", 1, 2, &actvals, inpi, &kunit, &null))
         SCTPUT ("Error while reading keyword INPUTI");

    inter      = inpi[0];
    ntrace     = inpi[1];

    SCKRDI ("INPUTI", 3, 2, &actvals, npix_hg,  &kunit, &null);
    SCKRDI ("INPUTI", 5, 2, &actvals, &scan[0],  &kunit, &null);

    /* Because arrays start at 0 */
    scan[0]  -= 1;
    scan[1]  -= 1;

    SCKRDR ("INPUTR", 1, 1, &actvals, &high_thres, &kunit, &null);

    SCKRDD ("INPUTD", 1, 2, &actvals, start_hg, &kunit, &null);
    SCKRDD ("INPUTD", 3, 2, &actvals, step_hg,  &kunit, &null);

    if (SCKGETC ("IN_A", 1, 60, &actvals, frame))
         SCTPUT ("Error while reading keyword IN_A");

    if (SCKGETC ("IN_B", 1, 60, &actvals, transf))
         SCTPUT ("Error while reading keyword IN_C");
 
    if (SCFOPN(frame, D_R4_FORMAT, 1, F_IMA_TYPE, &imnoa))
                       SCTPUT ("Error while opening input frame");

    SCDRDI(imnoa, "NPIX",  1, 2,  &actvals, npix,  &kunit, &null);
    SCDRDD(imnoa, "START", 1, 2,  &actvals, start, &kunit, &null);
    SCDRDD(imnoa, "STEP",  1, 2,  &actvals, step,  &kunit, &null);

    strcpy(cunit,"Slope           Ordin. InterceptCell Value      ");
    strcpy(ident,"Hough transform image");

    if (SCIPUT (transf, D_R4_FORMAT, F_O_MODE, F_IMA_TYPE,
                2, npix_hg, start_hg, step_hg, ident, cunit, (char **)&pntrb,
                &imnob))
         SCTPUT ("Error while opening output frame");

    /* Allocate a buffer to store the filtered columns being
     * used to compute the Hough Transform */

    nobyt = ntrace*npix[1]*sizeof(float); /* Number of bytes to be read */
    buffer = (float *) osmmget(nobyt);

    /* The input image is filtered and a constant background is subtracted 
    This step is not used. See below.
    SCTPUT("Filtering image.");
    prepare_image(imnoa, npix, inter, ntrace, buffer, scan); */


    /* SCTPUT("Performing Hough Transform."); */

    correct_image(imnoa, npix, inter, ntrace, buffer, scan);
    hough_transform (buffer, pntrb, npix, npix_hg, start_hg, 
                            step_hg, inter, ntrace, high_thres, scan);


    SCFCLO(imnob);
    SCFCLO(imnoa);

    osmmfree((char *)buffer);       /*  Bug on big frames  */


    SCSEPI();

  }


prepare_image(imnoa,  npix, inter, ntrace, buffer, scan)

    int imnoa, npix[2], scan[2];
    int inter, ntrace;
    float    *buffer;

{
    int             trace, rowread=(-1);
    int             nobyt, actsize, first_pass=0;
    float           *previous,*line,*next, *tmpbuf;
    float           histo[15], median, minimum;
    int             count, pos, row, col, i, j, posmax, nh;
    char            text[TEXT_LEN];
    extern int step_prgs, next_prgs;

    /* Allocate three buffers corresponding each to a row of pixels 
       to read the image row by row and process it */
    nobyt = npix[0]*sizeof(float); /* Number of bytes to be read */

    previous = (float *) osmmget(nobyt);  /* Row n-1 */
    line     = (float *) osmmget(nobyt);  /* Row n   */
    next     = (float *) osmmget(nobyt);  /* Row n+1 */

    /* Initialization: Read the two first rows of the image and store
       in previous and line */

    pos = ipos(1, scan[0], npix[0]); /* First pixel */
    SCFGET(imnoa, pos, npix[0], &actsize, (char *)previous);
    pos += npix[0];
    SCFGET(imnoa, pos, npix[0], &actsize, (char *)line);
    next_prgs = step_prgs;

    for (row=(scan[0]+2); row<=scan[1]; row++) {

       display_progress(row, npix[1]);

       pos = ipos(1, row, npix[0]);
       SCFGET(imnoa, pos, npix[0], &actsize, (char *)next);

       for (trace=1; trace<=ntrace; trace++) {
          col = icol(npix[0],trace,ntrace,inter);

          /* Store neighbouring values in a buffer histo to 
           * estimate the median */

          nh = 0;
          if ((col-2) >= 0) {
            histo[nh++] = previous[col-2];
            histo[nh++] = line[col-2];
            histo[nh++] = next[col-2];
	  }

          if ((col-1) >= 0) {
             histo[nh++] = previous[col-1];
             histo[nh++] = line[col-1];
             histo[nh++] = next[col-1];
	   }

          histo[nh++] = previous[col];
          histo[nh++] = line[col];
          histo[nh++] = next[col];

          if ((col+1) < npix[0]) {
             histo[nh++] = previous[col+1];
             histo[nh++] = line[col+1];
             histo[nh++] = next[col+1];
	   }

          if ((col+2) < npix[0]) {
             histo[nh++] = previous[col+2];
             histo[nh++] = line[col+2];
             histo[nh++] = next[col+2];
	   }

          /* Finds the median value out of 15=2*7+1 values stored in histo */
          sort(nh, histo);
          median = histo[(nh-1)/2];

          pos = ipos((trace-1), (row-1), ntrace);
          buffer[pos] = median;
          if (first_pass == 0) {
               minimum = median;
               first_pass = 1;
	     }
          if (median < minimum) minimum = median;
        } /* Matches for (trace =1... */

        /* Circular permutation of the buffers */

        tmpbuf   = previous;
        previous = line;
        line     = next;
        next     = tmpbuf;

     } /* matches for (row=2... */

      /* Copies the first and last row of buffer */
      for (trace=0; trace<ntrace; trace++) {
        buffer[ipos(trace, scan[0], ntrace)] = buffer[ipos(trace, (scan[0]+1), ntrace)];
        buffer[ipos(trace, scan[1], ntrace)] = buffer[ipos(trace, (scan[1]-1), ntrace)];
      }

      posmax = ntrace * npix[1];
      for (pos=0; pos<posmax; pos++)
                    buffer[pos] -= minimum;

      sprintf(text,"Subtracted constant value %f from the frame",minimum);
      SCTPUT(text);

      osmmfree((char *)previous);
      osmmfree((char *)line);
      osmmfree((char *)next);
}


correct_image(imnoa,  npix, inter, ntrace, buffer, scan)

    int imnoa, npix[2], scan[2];
    int inter, ntrace;
    float    *buffer;

{
    int             trace, rowread=(-1), rowmin, rowmax;
    int        nobyt, actsize, first_pass=0;
    float           *line;
    float           minimum, median;
    int        count, pos, row, col, posmax;
    char            text[TEXT_LEN];
    extern int step_prgs, next_prgs;

    /* Allocate a buffer  corresponding to a row of pixels 
       to read the image row by row and process it */
    nobyt = npix[0]*sizeof(float); /* Number of bytes to be read */

    line     = (float *) osmmget(nobyt);  /* Row n   */


    /* The minimum is estimated in the inner half of the scanned interval */

    rowmin = 3.*scan[0]/4. + scan[1]/4.;
    rowmax = scan[0]/4.    + 3.*scan[1]/4.;

    for (row=scan[0]; row<=scan[1]; row++) {

       pos = ipos(1, row, npix[0]);
       SCFGET(imnoa, pos, npix[0], &actsize, (char *)line);

       for (trace=1; trace<=ntrace; trace++) {

          col = icol(npix[0],trace,ntrace,inter);
          median = line[col];
          pos = ipos((trace-1), row, ntrace);
          buffer[pos] = median;

 
          if (row>rowmin && row<rowmax) {
            if (first_pass == 0) {
               minimum = median;
               first_pass = 1;
	     }

          if (median < minimum )  minimum = median;
	  }  

        } /* Matches for (trace =1... */
     } /* matches for (row=... */

      posmax = ntrace * npix[1];
      for (pos=0; pos<posmax; pos++)
                    buffer[pos] -= minimum;

      sprintf(text,"Subtracted constant value %f from the frame",minimum);
      SCTPUT(text);
      osmmfree((char *)line);
}


hough_transform (buffer, pntrb, npix, npix_hg, start_hg, 
                              step_hg, inter, ntrace, high_thres, scan)

    float    *buffer, *pntrb;
    int      npix[2], npix_hg[2];
    double   start_hg[2], step_hg[2], high_thres;
    int      inter, ntrace, scan[2];

{

    int            nrow_img, row, col, row_hg, col_hg;
    int       pos, actsize, rowstart;
    double         theta, slope, origin;
    float          increment;
    int            trace;
    extern int next_prgs, step_prgs;

    nrow_img = npix[1];
    step_prgs = 10;
    next_prgs = step_prgs;


/* Computes the local slope aint  the orders */

        for (row=scan[0]; row<=scan[1]; row++) {

           display_progress(row, npix[1]);

           for (trace=1; trace<=ntrace; trace++) {
                col = icol(npix[0],trace,ntrace,inter);
                pos = ipos((trace-1), row, ntrace);
                increment = buffer[pos];

                if (increment < high_thres) {

                for (col_hg=0; col_hg<npix_hg[0]; col_hg++) {  

                    slope  = start_hg[0] + step_hg[0] * col_hg;
                    origin = (row+1) - (col+1) * slope;
                    row_hg = (int) ((origin - start_hg[1])/step_hg[1] + 0.5);

                    if (row_hg > 0 && row_hg < npix_hg[1]) /* Add = */
                        pntrb[ipos(col_hg,row_hg,npix_hg[0])] += increment;

	  }
}
}
}

}
    


display_progress(row, nrow)

  int row, nrow;

{
  float   prgs;
  extern  int     next_prgs;
  extern  int     step_prgs;

  char       date[28], text[TEXT_LEN];
  struct tm  date_struct;

  prgs = 100. * (float) row / (float) nrow;
  if (prgs > next_prgs) {
     if (oshdate(date,&date_struct)) date[0] = '\0';
     sprintf(text,"%s  %d %% performed...", date, next_prgs);
     next_prgs += step_prgs;
     SCTPUT(text);
   }
}

/*  Sorting algorithm:      Heapsort method                      */
/*  From: Numerical Recipes, Cambridge University Press, p. 226  */

sort(n,ra)
int n;
float ra[];
{
        int l,j,ir,i;
        float rra;

        l=(n >> 1)+1;
        ir=n;
        for (;;) {
                if (l > 1)
                        rra=ra[--l];
                else {
                        rra=ra[ir];
                        ra[ir]=ra[1];
                        if (--ir == 1) {
                                ra[1]=rra;
                                return;
                        }
                }
                i=l;
                j=l << 1;
                while (j <= ir) {
                        if (j < ir && ra[j] < ra[j+1]) ++j;
                        if (rra < ra[j]) {
                                ra[i]=ra[j];
                                j += (i=j);
                        }
                        else j=ir+1;
                }
                ra[i]=rra;
        }
}