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

/* @(#)splinadd.c	17.1.1.1  (ESO)  01/25/02  17:54:01 */
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
.IDENT        splinadd.c
.MODULE       main program -- splinadd.exe
.LANGUAGE     C
.AUTHOR       Cristian Levin - ESO La Silla
.PURPOSE      This program makes the search line process over
              an image given initial "guess" lines in a row.
.KEYWORDS     line searching.
.VERSION 1.0  1-Sep-1991   Implementation
.ENVIRONMENT  UNIX
------------------------------------------------------------*/

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

#include <midas_def.h>

#include <gl_defs.h>
#include <spec_defs.h>

#define GRAVITY		0
#define GAUSSIAN	1
#define MAXIMUM		2

/* Limits for the line table */
#define COLSLIN		15
#define ROWSLIN		800

/* max no. of rows in LINEADD.tbl */
#define MAXADD		100

/* Descriptors of image */
int	Npix[2];
float	Start[2], Step[2];

/* Search parameters */
int	Method;		/* GRAVITY, GAUSSIAN, MAXIMUM */
int	Window, Ybin;
int	Ystart;

char 	Image[81], LineTable[81], LineAddTable[81];
int 	IdIma, IdLineTab, IdAddTab;

int 	ColLineTab[3];	/* :X, :Y, :PEAK */
int 	LineRows;	/* no. of rows in LINE.tbl */
int 	LineCurrRow;	/* current updating row of LINE.tbl */
float	*LineX, *LineY;	/* array of :X, :Y of LINE.tbl */

int 	Width;			/* width of the window */
double 	*Xgaus, *Ygaus, *A;	/* Aux. variables for gaussian fitting */

int	GaussErr = FALSE;

void 	init_midas();
void 	end_midas();
void 	free_data();
void 	read_parameters();
void 	read_add_table();
void 	read_line_table();
void 	update_table();
void 	find_lines();

int *	ivector();
float *	fvector();
double *dvector();

main()
{
    int  nulval; /* garbage */
    int n_add, i, done = 1, nrows = 2, slice;
    float *x, *intens;
    float x_add[MAXADD], x_found[MAXADD], int_found[MAXADD];
    int inside_win[MAXADD];
    char msg[80];

    init_midas();
    read_parameters();
    read_line_table();
    read_add_table( x_add, &n_add );

    x      = fvector(0, Npix[0] - 1); /* X-coordinates */
    intens = fvector(0, Npix[0]*Npix[1] - 1); /* Total intensity values */

    for ( i = 0; i < Npix[0]; i++ )
	x[i]  = Start[0] + Step[0] * i;

    SCFGET( IdIma, 1, Npix[0]*Npix[1], &nulval, (char *)intens );

    slice = Npix[1]/5;

    for ( i = 0; i < n_add; i++ ) /* initial guess: x-coords of LINEADD */
	x_found[i] = x_add[i];

    find_lines(x, intens+(Ystart-1)*Npix[0], x_found, int_found, 
	       inside_win, n_add, Ystart);
    update_table(x_found, int_found, inside_win, n_add, Ystart);

    /* add lines from 'Ystart'+1 to Npix[1] */
    for ( i = Ystart+1; i <= Npix[1]; i++, nrows++ ) { 
        find_lines(x, intens+(i-1)*Npix[0], x_found, int_found, 
		   inside_win, n_add, i);
        update_table(x_found, int_found, inside_win, n_add, i);

        if ( nrows == slice * done && done != 5 ) {
	    sprintf(msg, "%3d%% done...", done * 20);
	    SCTPUT(msg);
	    done++;
        }
    }
    for ( i = 0; i < n_add; i++ )
	x_found[i] = x_add[i];

    /* add lines from 'Ystart'-1 to 1 */
    for ( i = Ystart-1; i >= 1; i--, nrows++ ) {
        find_lines(x, intens+(i-1)*Npix[0], x_found, int_found, 
	           inside_win, n_add, i);
        update_table(x_found, int_found, inside_win, n_add, i);

        if ( nrows == slice * done && done != 5 ) {
	    sprintf(msg, "%3d%% done...", done * 20);
	    SCTPUT(msg);
	    done++;
        }
    }

    sprintf( msg, "100%% done..." );
    SCTPUT(msg);

    free_fvector(x,      0, Npix[0] - 1);
    free_fvector(intens, 0, Npix[0]*Npix[1] - 1);
    free_data();
    end_midas();
}

void find_lines( x, intens, x_found, intens_found, inside_win, n_add, row )
float *x; 		/* IN    : x-coords of current row pixels */
float *intens; 		/* OUT   : intensity of current row pixels */
float *x_found; 	/* IN/OUT: x-coords of lines found */
float *intens_found;	/* OUT   : intensity of lines found */
int *inside_win;	/* OUT   : TRUE if xfound[i] is not in the borders */
int n_add;		/* IN    : no. of lines to add */
int row;		/* IN    : current row */
{
    int i, j, k, l, ipix, factor, imax, av_rows;
    float max_intens, sum, *aux_intens, *line;
    float a, b, shift, aleft, aright;
    float median();

    aux_intens = fvector(1, Width);
    line = fvector(0, Npix[0] - 1);

    /* smooth by averaging of 2*Ybin+1 rows */
    if ( Ybin == 0 || row - Ybin < 1 || row + Ybin > Npix[1] )
	for ( i = 0; i < Npix[0]; i++ )
	    line[i] = intens[i];
    else {
	av_rows = 2*Ybin + 1;
	for ( i = 0; i < Npix[0]; i++ ) {
	    sum = 0;
	    for ( j = -Ybin; j <= Ybin; j++ )
		sum += intens[i+Npix[0]*j];
	    line[i] = sum / av_rows;
	}
    }

    for ( i = 0; i < n_add; i++ ) { 
	inside_win[i] = TRUE;
        ipix = (x_found[i] - Start[0]) / Step[0];
        if ( ipix + Window >= Npix[0] || ipix - Window < 0 ) {
	    inside_win[i] = FALSE;
	    continue;
	}
	max_intens = line[ipix];
        imax = ipix;
	for ( j = ipix - Window; j <= ipix + Window; j++ )
	    if ( line[j] > max_intens ) {
		max_intens = line[j];
		imax = j;
	    }
	switch ( Method ) {
	      case GAUSSIAN:
    		A[1] = line[ipix]; 
		A[2] = x[ipix]; 
		A[3] = Step[0];
		for ( k = 1, j = ipix-Window; j <= ipix+Window; j++, k++ ) {
		    Xgaus[k] = x[j];
		    Ygaus[k] = line[j];
		}
		fit_gauss( Xgaus, Ygaus, Width, A );
		if ( GaussErr ) {
		    GaussErr = FALSE;
		    continue;
		}
		intens_found[i] = A[1];
		x_found[i] = A[2];
		break;
	      case GRAVITY:
		if ( imax == 0 || imax == Npix[0]-1 ) {
	    	    inside_win[i] = FALSE;
		    continue;
		}
		aleft  = line[imax-1];
		aright = line[imax+1];
		factor = 1;
		if ( aleft >= aright ) {
		    aleft  = intens[imax+1];
		    aright = intens[imax-1];
		    factor = -1;
		}
		a = line[imax] - aleft;
		b = aright - aleft;
		shift = (a+b == 0.0) ? 0.0 : Step[0] * b / (a+b);
		x_found[i] =  x[imax] + factor * shift;
		intens_found[i] = line[imax];
		break;
	      case MAXIMUM:
		intens_found[i] = line[imax];
		x_found[i] = x[imax];
		break;
	}
    }
    free_fvector(aux_intens, 1, Width);
    free_fvector(line, 0, Npix[0] - 1);
}

void init_midas()
{
    SCSPRO( "SPADDL" );
}

void end_midas()
{
    SCSEPI();
}

void free_data()
{
    if ( LineRows != 0 ) {
    	free_fvector( LineX, 0, LineRows-1 );
    	free_fvector( LineY, 0, LineRows-1 );
    }
}

void read_parameters()
{
    int unit;		/* useless */
    int actval, nulval;	/* useless */
    char msg[80], meth[21];
    float inputr[3];

    SCKGETC( "P1", 1, 80, &actval, Image );
    SCKGETC( "P4", 1, 20, &actval, meth );
    SCKGETC( "P5", 1, 80, &actval, LineAddTable );
    SCKGETC( "P6", 1, 80, &actval, LineTable );
    SCKRDI( "INPUTI", 1, 1, &actval, &Ystart, &unit, &nulval );
    SCKRDR( "INPUTR", 1, 2, &actval, inputr, &unit, &nulval );
    Window = inputr[0];
    Ybin   = inputr[1];

    Method = GRAVITY;			/* default */
    if ( !strncmp(meth, "GAU", 3) || !strncmp(meth, "gau", 3) )
	Method = GAUSSIAN;
    else if ( !strncmp(meth, "MAX", 3) || !strncmp(meth, "max", 3) )
	Method = MAXIMUM;

    if ( SCFOPN( Image, D_R4_FORMAT, 0, F_IMA_TYPE, &IdIma ) != 0 ) {
	sprintf( msg, "Frame %s invalid...", Image );
	SCTPUT( msg );
	end_midas();
    }

    SCDRDI( IdIma, "NPIX", 1, 2, &actval, Npix, &unit, &nulval );
    SCDRDR( IdIma, "START", 1, 2, &actval, Start, &unit, &nulval );
    SCDRDR( IdIma, "STEP", 1, 2, &actval, Step, &unit, &nulval );

    Width  = 2 * Window + 1;
    Xgaus  = dvector( 1, Width );
    Ygaus  = dvector( 1, Width );
    A      = dvector( 1, 3 );
    A[3]   = Step[0]; /* initial guess for the deviation in gaussian fitting */
}

void read_line_table()
{
    int i, col;
    int  nulval, sortcol, aw, ar, ncols; 

    if ( !file_exists( LineTable, ".tbl" ) ) {
	LineCurrRow = LineRows = 0;
    	TCTINI( LineTable, F_TRANS, F_O_MODE, COLSLIN, ROWSLIN, &IdLineTab );
    	TCCINI( IdLineTab, D_R4_FORMAT, 1, "F10.2", "PIXEL", "X", &ColLineTab[0] );
    	TCCINI( IdLineTab, D_R4_FORMAT, 1, "F10.2", "PIXEL", "Y", &ColLineTab[1] );
    	TCCINI( IdLineTab, D_R4_FORMAT, 1, "E12.3", "DN", "PEAK", &ColLineTab[2] );
    	TCCINI( IdLineTab, D_C_FORMAT,  4, "A1", " ", "ERASED", &col );
	return;
    }

    TCTOPN( LineTable, F_IO_MODE, &IdLineTab );
    TCIGET( IdLineTab, &ncols, &LineRows, &sortcol, &aw, &ar );
    TCCSER( IdLineTab, ":X", &ColLineTab[0] );
    TCCSER( IdLineTab, ":Y", &ColLineTab[1] );
    TCCSER( IdLineTab, ":PEAK", &ColLineTab[2] );

    LineCurrRow = LineRows;
    LineX = fvector( 0, LineRows-1 );
    LineY = fvector( 0, LineRows-1 );

    for ( i = 0; i < LineRows; i++ ) {
	TCERDR( IdLineTab, i+1, ColLineTab[0], &LineX[i], &nulval );
	TCERDR( IdLineTab, i+1, ColLineTab[1], &LineY[i], &nulval );
    }
}

void read_add_table( x_add, n_add)
float *x_add;
int *n_add;
{
    int  nulval, sortcol, aw, ar, ncols; 
    int i, col;
    char str[MAXLINE];

    if ( TCTOPN( LineAddTable, F_I_MODE, &IdAddTab ) != 0 ) {
	sprintf( str, "Table %s couldn't be opened.", LineAddTable );
	SCTPUT( str );
	end_midas();
    }
    TCIGET( IdAddTab, &ncols, n_add, &sortcol, &aw, &ar );
    TCCSER( IdAddTab, ":X", &col );

    for ( i = 0; i < *n_add; i++ )
	TCERDR( IdAddTab, i+1, col, &x_add[i], &nulval );

    TCTCLO( IdAddTab );
}

void update_table( x_found, intens_found, inside_win, n_add, image_row )
float *x_found; 	/* x-coords of lines found */
float *intens_found; 	/* intensity of lines found */
int *inside_win;	/* TRUE if xfound[i] is not in the borders */
int n_add; 		/* no. of lines to add */
int image_row;		/* current image row */
{
    float value[3];
    int i, j;

    value[1] = Start[1] + Step[1] * (image_row - 1);

    for ( i = 0; i < n_add; i++ ) {
	if ( !inside_win[i] ) 
	    continue;
	/* check if the spectral line already exists */
        for ( j = 0; j < LineRows; j++ )
	    if ( LineY[j] == value[1] && (x_found[i] <= LineX[j] + Window) &&
		 (x_found[i] >= LineX[j] - Window) )
		break;
	if ( j < LineRows )
	    continue;

    	value[0] = x_found[i];
    	value[2] = intens_found[i];
	LineCurrRow++;
    	TCRWRR( IdLineTab, LineCurrRow, 3, ColLineTab, value );
    }
}

/* for management of math local routine errors */

/*
void nrerror( msg ) 
char *msg;
{
    GaussErr = TRUE;
}
*/