/* @(#)lnhough.c	17.1.1.1 (ES0-DMD) 01/25/02 17:53:59 */
/*===========================================================================
  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>

main ()
{

    char           catal[81], transf[81], line[81];
    char           colx[11], colw[11], mode[11];
    int            npix_hg[3], ndim, kunit;
    double         start_hg[3], step_hg[3];
    int            null, status, npar, iav, actvals;
    int            inpi[2], inter, ntrace, scan[2];
    int            linx, catwav;
    int            nbcol, nbrow, nsort, allcol, allrow;
    int            nbcolcat, nbrowcat, nsortcat, allcolcat, allrowcat;
    int            tidlin, tidcat, rowlin, rowcat;
    double         x, wave;
    float          unity = 1., rpar[5], avdisp, wcenter, delta;
    float          range, xnorm, xc, findmax(), max;
    int            col1, col2, col3;

    SCSPRO ("hough");

    SCKRDI ("INPUTI", 4, 1, &actvals, &ndim, &kunit, &null);
    SCKRDI ("INPUTI", 1, 3, &actvals, npix_hg, &kunit, &null);
    SCKRDD ("INPUTD", 1, 3, &actvals, start_hg, &kunit, &null);
    SCKRDD ("INPUTD", 4, 3, &actvals, step_hg,  &kunit, &null);
    SCKRDR ("INPUTR", 1, 4, &actvals, rpar,  &kunit, &null);
    wcenter = rpar[0], avdisp = rpar[1], xc = rpar[2], range = rpar[3];

    SCKGETC ("IN_A",  1, 60, &actvals, line);
    SCKGETC ("IN_B",  1, 60, &actvals, catal);
    SCKGETC ("OUT_A", 1, 60, &actvals, transf);
    SCKGETC ("OUT_B", 1, 10, &actvals, mode);

    SCKGETC ("INPUTC", 1,  10, &actvals, colx); 
    SCKGETC ("INPUTC", 10, 10, &actvals, colw); 

    TCTOPN(line,   F_I_MODE, &tidlin);
    TCIGET(tidlin, &nbcol, &nbrow, &nsort, &allcol, &allrow);
    TCCSER(tidlin, colx, &linx);

    if (catal[0] != '@') {
       TCTOPN(catal, F_I_MODE, &tidcat);
       TCIGET(tidcat, &nbcolcat, &nbrowcat, &nsortcat, &allcolcat, &allrowcat);
       TCCSER(tidcat, colw, &catwav);
    }
    else
       TCCSER(tidlin, colw, &catwav);


    create_hough(transf, npix_hg, start_hg, step_hg, ndim);

    for (rowlin=1; rowlin<=nbrow; rowlin++) {
        TCERDD(tidlin, rowlin, linx,   &x, &null);

    if (catal[0] != '@') {
    for (rowcat=1; rowcat<=nbrowcat; rowcat++) {    
        TCERDD(tidcat, rowcat, catwav, &wave, &null); 
        increment_hough(x,wave,unity,mode,avdisp,range);
      }}
    else {
       TCERDD(tidlin, rowlin, catwav, &wave, &null); 
       increment_hough(x,wave,unity,mode,avdisp,range);
    }

    }

    max = findmax(&col1,&col2,&col3);
    printf ("Found maximum %f in %d %d %d\n",max,col1,col2,col3);
    SCKWRI("OUTPUTI",&col1,3,1,&kunit);
    SCKWRI("OUTPUTI",&col2,4,1,&kunit);
    SCKWRI("OUTPUTI",&col3,5,1,&kunit);
    SCKWRR("OUTPUTR",&max,1,1,&kunit);

    close_hough();
    if (catal[0] != '@') TCTCLO(tidcat);
    TCTCLO(tidlin), SCSEPI();

  }