/*===========================================================================
  Copyright (C) 2001 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
  ===========================================================================*/
/*-------------------------------------------------------------------------*/
/**
 * @defgroup flames_backfit   Substep: fit the background 
 *
 */
/*-------------------------------------------------------------------------*/

/*---------------------------------------------------------------------------
  Includes
  --------------------------------------------------------------------------*/
#ifdef HAVE_CONFIG_H
#  include <config.h>
#endif


#include <ctype.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
/* MIDAS include files */
#include <flames_midas_def.h>
#include <flames_computeback.h>
/* FLAMES-UVES include files */ 
#include <flames_readframe.h>
#include <flames_scatter.h>
#include <flames_freeordpos.h>
#include <flames_freeframe.h>
#include <flames_writeback.h>
#include <flames_uves.h>
#include <flames_backfit.h>
#include <flames_newmatrix.h>
#include <flames_readordpos.h>
#include <flames_readback.h>

/**@{*/

/*---------------------------------------------------------------------------
  Implementation
  ---------------------------------------------------------------------------*/
/**
   @name  flames_backfit()  
   @short   Fit background 
   @author G. Mulas  -  ITAL_FLAMES Consortium. Ported to CPL by A. Modigliani

   @param MAXDISCARDFRACT  maximum fraction of windows/pixels which we are 
                           willing to discard by kappa-sigma clipping in 
                           each iteration of the background fitting loop 
   @param MAXBACKITERS     maximum number of iterations which we are willing 
                           to perform in correlation
   @param FRAMENAME        name of frame
   @param BACKTABLE        background table
   @param ORDERTABLE       order table
   @param BACKFRAME        background frame
   @param BKGPOL           array tohold X,Y polynomial degree
   @param BKGFITMETHOD     Background fitting method
   @param BKGBADSCAN       Background table bad pixel frame scanning switch:
                           possible values are: none, fraction, absolute
   @param BKGBADWIN        Background table bad pixel frame scanning window 
                           size. Defaults are: 50,50
   @param BKGBADMAXFRAC    Background table bad pixel frame scanning threshold 
                           fraction
   @param BKGBADMAXTOT     Background table bad pixel frame scanning threshold 
                           number
   @param SIGMA            Kappa value used in kappa-sigma clipping (it should 
                           be called kappa!)
   @param OUTPUTI

   @return initialized frame

   DRS Functions called:                                                   

   Pseudocode:   -init all relevant structures
                 -read all relevant inputs 
                 -call scatter function to compute inter-order background

                 -to make the thing more robust: scan the frame for any 
                  negative pixel values; if any are found, compare them to the
                  standard deviation: if the negative value is compatible with 
                  zero, set that pixel to zero, otherwise mark that pixel as 
                  bad==> BEWARE HERE IN SOME WAY MAY ALTERATE BACKGROUND!!


   @note
 */

/* C functions include files */ 


int 
flames_backfit(
		const double *MAXDISCARDFRACT,
		const int *MAXBACKITERS,
		const char *FRAMENAME,
		const char *BACKTABLE,
		const char *ORDERTABLE,
		const char *BACKFRAME,
		const int *BKGPOL,
		const char *BKGFITMETHOD,
		const char *BKGBADSCAN,
		const int *BKGBADWIN,
		const double *BKGBADMAXFRAC,
		const int *BKGBADMAXTOT,
		const double *SIGMA,
		int *OUTPUTI)

{ 

	char output[200];
	int i=0;
	int unit=0;
	int actvals=0;
	int bxdegree=0, bydegree=0;
	double kappa=0;
	double kappa2=0;
	frame_data **backframe;
	flames_err status=0;

	flames_frame *ScienceFrame=0;

	orderpos *Order=0;

	char backfile[CATREC_LEN+2];
	char infile[CATREC_LEN+2];
	char orderfile[CATREC_LEN+2];
	char backframefile[CATREC_LEN+2];
	int nval=0;
	int null=0;
	scatterswitch bkgswitch=USEALL;
	scatterswitch2 bkgswitch2=NOBADSCAN;
	int badwinxsize=0;
	int badwinysize=0;
	double badfracthres=0;
	int badtotthres=0;
	char keytype=0;
	int noelem=0;
	int bytelem=0;
	int maxbackiters=0;
	double maxdiscardfract=0;
	char bkgfitmethod[CATREC_LEN+1];
	int32_t negativepixels=0;






	memset(output, 0, 200);
	memset(bkgfitmethod, 0, CATREC_LEN+1);

	memset(backfile, 0, CATREC_LEN+2);
	memset(infile, 0, CATREC_LEN+2);
	memset(orderfile, 0, CATREC_LEN+2);
	memset(backframefile, 0, CATREC_LEN+2);

	bkgswitch = USEALL;

	SCSPRO("flames_backfit"); /* Get into MIDAS Env. */

	/* Read once and for all, here at the beginning, what MIDAS keywords we
     need */

	/* initialise MAXDISCARDFRACT from keyword */
	if (SCKRDD(MAXDISCARDFRACT, 1, 1, &actvals, &maxdiscardfract, &unit,
			&null) != 0) {
		/* problems reading MAXDISCARDFRACT */
		SCTPUT("Error reading the MAXDISCARDFRACT keyword");
		return flames_midas_fail();
	}

	/* initialise MAXBACKITERS from keyword */
	if (SCKRDI(MAXBACKITERS, 1, 1, &actvals, &maxbackiters, &unit,
			&null) != 0) {
		/* problems reading MAXBACKITERS */
		SCTPUT("Error reading the MAXBACKITERS keyword");
		return flames_midas_fail();
	}

	/* get input frame name from MIDAS env.*/
	SCKGETC(FRAMENAME,1,CATREC_LEN+1,&nval,infile);

	/* read input table name */
	SCKGETC(BACKTABLE,1,160,&nval,backfile);

	/* read order table name */
	SCKGETC(ORDERTABLE,1,160,&nval,orderfile);

	/* read output fitted background frame name */
	SCKGETC(BACKFRAME,1,160,&nval,backframefile);

	/* initialise background fitting scalars */
	if (SCKRDI(BKGPOL, 1, 1, &actvals, &bxdegree, &unit, &null)
			!= 0) {
		/* problems reading xdegree */
		SCTPUT("Error reading the x degree of the background polynomial");
		return flames_midas_fail();
	}
	if (SCKRDI(BKGPOL, 2, 1, &actvals, &bydegree, &unit, &null)
			!= 0) {
		/* problems reading xdegree */
		SCTPUT("Error reading the y degree of the background polynomial");
		return flames_midas_fail();
	}


	/* before trying to read it, make sure that BKGFITMETHOD exists and it is
     of the appropriate type */
	if (SCKFND_string(BKGFITMETHOD, &keytype, &noelem, &bytelem) != 0) {
		/* SCKFND failed, give up! */
		SCTPUT("Internal MIDAS error in mainopt: SCKFND failed");
		return flames_midas_fail();
	}
	switch(keytype) {
	case 'C':
		/* it exists and it is a character keyword, go ahead and read it */
		if (SCKGETC(BKGFITMETHOD, 1, CATREC_LEN, &nval, bkgfitmethod)
				!= 0) {
			/* problems reading BKGFITMETHOD */
			SCTPUT("Warning: error reading the BKGFITMETHOD keyword, falling back to \
default");
		}
		else {
			/* is the string long enough to be unambiguous? */
			if (nval<2) {
				/* no, fall back to the default */
				SCTPUT("Warning: BKGFITMETHOD is ambiguous, falling back to default");
			}
			else {
				/* convert bkgfitmethod to upper case, to ease the subsequent check */
				for (i=0; i<=(nval-1); i++) bkgfitmethod[i] = toupper(bkgfitmethod[i]);
				/* compare as many letters as we have with the expected values */
				if (strncmp("ALL", bkgfitmethod, (size_t) nval) == 0)
					bkgswitch = USEALL;
				else if (strncmp("MEDIAN", bkgfitmethod, (size_t) nval)
						== 0) bkgswitch = USEMEDIAN;
				else if (strncmp("MINIMUM", bkgfitmethod, (size_t) nval)
						== 0) bkgswitch = USEMINIMUM;
				else if (strncmp("AVERAGE", bkgfitmethod, (size_t) nval)
						== 0) bkgswitch = USEAVERAGE;
				else {
					SCTPUT("Warning: unsupported BKGFITMETHOD value, falling back to \
default");
				}
			}
		}
		break;
	case ' ':
		/* the keyword does not exist at all */
		SCTPUT("Warning: BKGFITMETHOD is undefined, falling back to default");
		break;
	default:
		/* the keyword is of the wrong type */
		SCTPUT("Warning: BKGFITMETHOD is not a string, falling back to default");
		break;
	}

	/* before trying to read it, make sure that BKGFITMETHOD exists and it is
     of the appropriate type */
	if (SCKFND_string(BKGBADSCAN, &keytype, &noelem, &bytelem) != 0) {
		/* SCKFND failed, give up! */
		SCTPUT("Internal MIDAS error in mainopt: SCKFND failed");
		return flames_midas_fail();
	}
	switch(keytype) {
	case 'C':
		/* it exists and it is a character keyword, go ahead and read it */
		if (SCKGETC(BKGBADSCAN, 1, CATREC_LEN, &nval, bkgfitmethod)
				!= 0) {
			/* problems reading BKGFITMETHOD */
			SCTPUT("Warning: error reading the BKGBADSCAN keyword, falling back to \
default");
		}
		else {
			/* is the string long enough to be unambiguous? */
			if (nval<1) {
				/* no, fall back to the default */
				SCTPUT("Warning: BKGBADSCAN is ambiguous, falling back to default");
			}
			else {
				/* convert bkgfitmethod to upper case, to ease the subsequent check */
				for (i=0; i<=(nval-1); i++) bkgfitmethod[i] = toupper(bkgfitmethod[i]);
				/* compare as many letters as we have with the expected values */
				if (strncmp("NONE", bkgfitmethod, (size_t) nval) == 0)
					bkgswitch2 = NOBADSCAN;
				else if (strncmp("FRACTION", bkgfitmethod, (size_t) nval)
						== 0) bkgswitch2 = FRACBADSCAN;
				else if (strncmp("ABSOLUTE", bkgfitmethod, (size_t) nval)
						== 0) bkgswitch2 = ABSBADSCAN;
				else {
					SCTPUT("Warning: unsupported BKGBADSCAN value, falling back to \
default");
				}
			}
		}
		break;
	case ' ':
		/* the keyword does not exist at all */
		SCTPUT("Warning: BKGBADSCAN is undefined, falling back to default");
		break;
	default:
		/* the keyword is of the wrong type */
		SCTPUT("Warning: BKGBADSCAN is not a string, falling back to default");
		break;
	}

	/* if neighborhood bad pixel scanning was requested, read the other
     keywords needed */
	if (bkgswitch2 == FRACBADSCAN) {
		if ((SCKRDI(BKGBADWIN, 1, 1, &actvals, &badwinxsize, &unit, &null) != 0)
				|| (SCKRDI(BKGBADWIN, 2, 1, &actvals, &badwinysize, &unit, &null)
						!= 0)) {
			SCTPUT("Error reading the BKGBADWIN keyword");
			return flames_midas_fail();
		}
		if (SCKRDD(BKGBADMAXFRAC, 1, 1, &actvals, &badfracthres, &unit,
				&null) != 0) {
			SCTPUT("Error reading the BKGBADMAXFRAC keyword");
			return flames_midas_fail();
		}
		/* check the values read for consistence */
		if ((badwinxsize < 0) || (badwinysize < 0)) {
			SCTPUT("Warning: BKGBADWIN values must be non negative, disabling \
BKGBADSCAN");
			bkgswitch2 = NOBADSCAN;
		}
		if (badfracthres < 0) {
			SCTPUT("Warning: BKGBADMAXFRAC value must be non negative, disabling \
BKGBADSCAN");
			bkgswitch2 = NOBADSCAN;
		}
	}
	else if (bkgswitch2 == ABSBADSCAN) {
		if ((SCKRDI(BKGBADWIN, 1, 1, &actvals, &badwinxsize, &unit, &null) != 0)
				|| (SCKRDI(BKGBADWIN, 2, 1, &actvals, &badwinysize, &unit, &null)
						!= 0)) {
			SCTPUT("Error reading the BKGBADWIN keyword");
			return flames_midas_fail();
		}
		if (SCKRDI(BKGBADMAXTOT, 1, 1, &actvals, &badtotthres, &unit,
				&null) != 0) {
			SCTPUT("Error reading the BKGBADMAXTOT keyword");
			return flames_midas_fail();
		}
		/* check the values read for consistence */
		if ((badwinxsize < 0) || (badwinysize < 0)) {
			SCTPUT("Warning: BKGBADWIN values must be non negative, disabling \
BKGBADSCAN");
			bkgswitch2 = NOBADSCAN;
		}
		if (badtotthres < 0) {
			SCTPUT("Warning: BKGBADMAXTOT value must be non negative, disabling \
BKGBADSCAN");
			bkgswitch2 = NOBADSCAN;
		}
	}

	/* read the kappa factor to be used later in kappa-sigma clipping */
	if ((status=SCKRDD(SIGMA, 1, 1, &actvals, &kappa, &unit, &null))!=0) {
		/* something went wrong while reading the kappa-sigma factor */
		sprintf(output, "Error %d while reading SIGMA keyword", status);
		SCTPUT(output);
		return flames_midas_fail();
	}
	/* compute once and for all the square of kappa, as we will be using that */
	kappa2 = kappa*kappa;


	/* Link the MIDAS names of the frames to the physical ones */

	if(!(ScienceFrame = calloc(1, sizeof(flames_frame)))) {
		SCTPUT("Allocation error during ScienceFrame memory allocation");
		return flames_midas_fail();
	}

	/* Read the data of the frames and put them in the right C structures */

	SCTPUT("Background fitting started\n");

	/* Get the Midas name of the frame to be read */

	/* let's read the Science Frame */

	sprintf(output, "I'm reading the frame %s", infile);
	SCTPUT(output);

	if (readframe(ScienceFrame, infile) != NOERR) {
		SCTPUT("Error while reading the frame");
		return flames_midas_fail();
	}

	/* Read the table data and then put them in the C structures */
	sprintf(output,"Reading the order/fibre table...");
	SCTPUT(output);
	if(!(Order = calloc(1, sizeof(orderpos)))) {
		SCTPUT("Allocation error during the allocation of Order structure");
		return flames_midas_fail();
	}
	/* use the readordpos function to read the descriptors */
	if (readordpos(orderfile, Order) != NOERR) {
		SCTPUT("Error while reading the order table");
		return flames_midas_fail();
	}
	/* does the order table match the Science frame chip? */
	if(Order->chipchoice != ScienceFrame->chipchoice) {
		/* no, it doesn't */
		SCTPUT("Error: chip mismatch between Science frame and order table");
		return flames_midas_fail();
	}

	/* initialise firstorder and lastorder in ScienceFrame from Order */
	ScienceFrame->firstorder = Order->firstorder;
	ScienceFrame->lastorder = Order->lastorder;
	ScienceFrame->tab_io_oshift = Order->tab_io_oshift;

	sprintf(output, "I'm reading the background table %s", backfile);
	SCTPUT(output); /* Display text w/o storing it into MIDAS logfile */

	/* read in the background table */
	if ((status=readback(&(ScienceFrame->back), backfile, bxdegree, bydegree))
			!= NOERR) {
		/* something went wrong while reading the background table */
		sprintf(output, "Error %d while reading the background table", status);
		SCTPUT(output);
		return flames_midas_fail();
	}


	/* Calculate the fit model of the background */

	SCTPUT("Start the background fitting procedure");

	if (ScienceFrame->back.Window_Number > 0) {
		if (scatter(ScienceFrame, Order, bkgswitch, bkgswitch2, badwinxsize,
				badwinysize, badfracthres, badtotthres, kappa2,
				(int32_t) maxbackiters, maxdiscardfract, OUTPUTI)) {
			SCTPUT("Error executing the scatter function");
			return flames_midas_fail();
		}
		/* allocate and initialise the frame which will contain the
       estimated background */
		backframe = fdmatrix(0, ScienceFrame->subrows-1, 0,
				ScienceFrame->subcols-1);
		memset(&backframe[0][0], 0,
				ScienceFrame->subrows*ScienceFrame->subcols*sizeof(frame_data));
		/* compute the estimated background */
		if (computeback(ScienceFrame, backframe) != NOERR) {
			SCTPUT("Error computing fitted background");
			return flames_midas_fail();
		}
		/* some more black magic to make the thing more robust: scan the frame
       for any negative pixel values; if any are found, compare them to the
       standard deviation: if the negative value is compatible with zero, 
       set that pixel to zero, otherwise mark that pixel as bad */
		negativepixels=0;
		frame_data* fdvecbuf1 = ScienceFrame->frame_array[0];
		frame_data* fdvecbuf2 = ScienceFrame->frame_sigma[0];
		frame_mask* fmvecbuf1 = ScienceFrame->badpixel[0];
		frame_data* fdvecbuf3 = backframe[0];
		frame_data pixelvalue=0;
		for (i=0; i<=((ScienceFrame->subrows*ScienceFrame->subcols)-1); i++) {
			if (fmvecbuf1[i]==0 && (pixelvalue=(fdvecbuf1[i]-fdvecbuf3[i]))<0) {
				if ((pixelvalue*pixelvalue)>fdvecbuf2[i]) {
					negativepixels++;
				}
			}
		}
		if (negativepixels!=0) {
			sprintf(output, "Warning: %d pixels result lower than fitted background",
					negativepixels);
			SCTPUT(output);
			SCTPUT("either they are unmasked bad pixels or some contamination is");
			SCTPUT("skewing background determination");
		}
	}
	else {
		SCTPUT("Error: no regions available for background estimation");
		return flames_midas_fail();
	}

	sprintf(output, "Writing fitted background frame to %s", backframefile);
	SCTPUT(output);
	if (writeback(ScienceFrame, backframefile, backframe)!=NOERR) {
		SCTPUT("Error writing background frame to disk");
		return flames_midas_fail();
	}


	/* free the memory for internal ScienceFrame members */
	if (freeframe(ScienceFrame) != NOERR) {
		SCTPUT("Warning: error freeing frame structure memory");
	}
	/* free the ScienceFrame structure itself */
	free(ScienceFrame);

	/* free the memory for internal Order members */
	if (freeordpos(Order) != NOERR) {
		SCTPUT("Warning: error freeing Order structure memory");
	}
	/* free the Order structure itself */
	free(Order);

	/* free the fitted background frame */
	free_fdmatrix(backframe, 0, ScienceFrame->subrows-1, 0,
			ScienceFrame->subcols-1);


	SCTPUT("\n*** Background fitting complete ***\n");

	return SCSEPI();

}
/**@}*/