/*===========================================================================
  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_fastfillholes
 *
 */
/*---------------------------------------------------------------------------
  Includes
  --------------------------------------------------------------------------*/
/**
   @name  flames_fastfillholes()  
   @short background fitting and subtraction, foreach x,order locate fibre 
   start-end
   "quick and dirty" no y shift correction flavour
   @author G. Mulas  -  ITAL_FLAMES Consortium. Ported to CPL by A. Modigliani


   @param allflatsin, 
   @param ordpos,  
   @param backname, 
   @param bxdegree, 
   @param bydegree, 
   @param bkgswitch, 
   @param bkgswitch2,
   @param badwinxsize, 
   @param badwinysize, 
   @param badfracthres,
   @param badtotthres, 
   @param kappa2,  
   @param maxbackiters, 
   @param maxdiscardfract,
   @param fracslicesthres,
   @param OUTPUTI)

   @return
   @doc

   DRS Functions called:          
   ordselect()
   standard()

   Pseudocode:                                                             
   for ord_start< ord <ord_end                                           
   if (no adjacent order overlapping)                                
   standard extract ord                                           
   endif                                                             
   endfor                                                               


   @note
 */


#ifdef HAVE_CONFIG_H
#  include <config.h>
#endif

#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <math.h>
#include <flames_midas_def.h>
#include <flames_scatter.h>
#include <flames_def_drs_par.h>
#include <flames_uves.h>
#include <flames_freeback.h>
#include <flames_fastfillholes.h>
#include <flames_fastlocatefibre.h>
#include <flames_getordpos.h>
#include <flames_computeback.h>
#include <flames_newmatrix.h>
#include <flames_readback.h>


flames_err 
fastfillholes(allflats *allflatsin, 
              orderpos *ordpos,  
              char *backname, 
              int bxdegree, 
              int bydegree, 
              scatterswitch bkgswitch, 
              scatterswitch2 bkgswitch2,
              int badwinxsize, 
              int badwinysize, 
              double badfracthres,
              int badtotthres, 
              double kappa2,  
              int32_t maxbackiters, 
              double maxdiscardfract,
              double fracslicesthres,
              int *OUTPUTI)
{
    char output[CATREC_LEN+1];
    flames_frame *framebuffer=0;
    frame_data **backframe;

    double x=0;
    double ordercentre=0;
    double totslices=0;
    double fracslices=0;
    int32_t ix=0;
    int32_t iy=0;
    int32_t iorder=0;
    int32_t iframe=0;
    int32_t lfibre=0;
    int32_t ifibre=0;
    int32_t ifibreixtotlimit=0;
    int32_t goodpixels=0;
    flames_err status=0;
    int32_t *newfibres=0;
    int32_t newnumfibres=0;
    frame_data pixelvalue=0;

    int32_t iyixend=0;
    int32_t iyixindex=0;
    int32_t iorderifibreindex=0;
    int32_t iorderifibreixindex=0;
    int32_t ifibreixindex=0;
    int32_t *lvecbuf1=0;


    frame_data *fdvecbuf1=0;
    frame_data *fdvecbuf2=0;
    frame_mask *fmvecbuf1=0;
    frame_mask *fmvecbuf2=0;

    int actvals=0;
    char drs_verbosity[10];
    int mid_stat=0;

    memset(drs_verbosity, 0, 10);
    if ((mid_stat=SCKGETC(DRS_VERBOSITY, 1, 3, &actvals, drs_verbosity))
                    != 0) {
        /* the keyword seems undefined, protest... */
        return MAREMMA;
    }

    memset(output, 0, CATREC_LEN+1);

    /* since it happens to have almost null fibre values on unmasked bad
     columns, try to avoid them */

    SCTPUT("Initialising internal structures...");

    /* compute the halfibrewidth in pixel units */
    ifibreixtotlimit = (allflatsin->maxfibres*allflatsin->subcols)-1;
    /* begin looping over orders at first */
    for (iorder=0; iorder <= ((ordpos->lastorder)-(ordpos->firstorder));
                    iorder++) {
        double order = (double) (iorder+(ordpos->firstorder));
        /* loop over x... */
        /* since we are at it, completely clean up normfactors and normsigmas */
        fdvecbuf1 = allflatsin->normfactors[iorder][0];
        fdvecbuf2 = allflatsin->normsigmas[iorder][0];
        for (ix=0; ix<=ifibreixtotlimit; ix++) {
            fdvecbuf1[ix] = 1;
            fdvecbuf2[ix] = 0;
        }
        for (ix=0; ix<=(allflatsin->subcols-1); ix++) {
            /* convert the ix pixel coordinate to the x world coordinate */
            x = allflatsin->substartx+(allflatsin->substepx)*((double) ix);
            /* find the unshifted central position and slope of this order at
	 this x */
            /* bail out if the function call return an error status */
            if ((status = get_ordpos(ordpos, order, x, &ordercentre))!=NOERR) {
                return(status);
            }
            for (iframe=0; iframe<=allflatsin->nflats-1;iframe++) {
                /* loop over lit fibres in this frame */
                for (lfibre=0; lfibre<=((allflatsin->flatdata)[iframe]).numfibres-1;
                                lfibre++) {
                    ifibre = (((allflatsin->flatdata)[iframe]).fibres)[lfibre];
                    /* find this fibre centre and boundaries */
                    if ((status=fastlocatefibre(allflatsin, ordpos, ordercentre,
                                    iorder, ifibre, ix)) != NOERR) {
                        return status;
                    }
                }
            }
        }
    }

    /* split the order loop in two parts, so that we can fit and subtract the
     background from fibre FF frames */

    SCTPUT("Background fitting and subtraction starting...");

    /* background fitting and subtraction section starts here */
    framebuffer = calloc(1, sizeof(flames_frame));
    framebuffer->subrows = allflatsin->subrows;
    framebuffer->subcols = allflatsin->subcols;
    framebuffer->substepx = allflatsin->substepx;
    framebuffer->substepy = allflatsin->substepy;
    framebuffer->substartx = allflatsin->substartx;
    framebuffer->substarty = allflatsin->substarty;
    framebuffer->nflats = 0;
    framebuffer->maxfibres = allflatsin->maxfibres;
    framebuffer->fibremask = cvector(0, framebuffer->maxfibres-1);
    iyixend = (allflatsin->subrows*allflatsin->subcols)-1;
    /* allocate and initialise the frame which will contain the
     estimated background */
    backframe = fdmatrix(0, framebuffer->subrows-1, 0, framebuffer->subcols-1);
    memset(backframe[0], 0,
           framebuffer->subrows*framebuffer->subcols*sizeof(frame_data));
    for (iframe=0; iframe<=allflatsin->nflats-1;iframe++) {
        /* build the fibremask for the framebuffer */
        for (ifibre=0; ifibre<=framebuffer->maxfibres-1; ifibre++)
            framebuffer->fibremask[ifibre]=FALSE;
        for (lfibre=0; lfibre<=allflatsin->flatdata[iframe].numfibres-1; lfibre++)
            framebuffer->fibremask[allflatsin->flatdata[iframe].fibres[lfibre]]=TRUE;
        /* assign pointers to the frame data, sigma and badpixel masks */
        framebuffer->frame_array = allflatsin->flatdata[iframe].data;
        framebuffer->frame_sigma = allflatsin->flatdata[iframe].sigma;
        framebuffer->badpixel = allflatsin->flatdata[iframe].badpixel;
        /* read in the background table */
        if ((status=readback(&(framebuffer->back), backname, bxdegree, bydegree))
                        != NOERR) {
            /* something went wrong while reading the background table */
            SCTPUT("Error while reading the background table");
            return flames_midas_fail();
        }
        if (framebuffer->back.Window_Number > 0) {
            if (scatter(framebuffer, ordpos, bkgswitch, bkgswitch2, badwinxsize,
                            badwinysize, badfracthres, badtotthres, kappa2,
                            maxbackiters, maxdiscardfract, OUTPUTI)) {
                SCTPUT("Error executing the scatter function");
                return flames_midas_fail();
            }
            /* compute the estimated background */
            if (computeback(framebuffer, backframe) != NOERR) {
                SCTPUT("Error computing fitted background");
                return flames_midas_fail();
            }
            if (freeback(&(framebuffer->back)) != NOERR) {
                SCTPUT("Error freeing framebuffer->back");
                return flames_midas_fail();
            }
            /* subtract the estimated background from the data frame */
            /* the error of the estimated background is assumed to be negligible */
            SCTPUT("Subtracting fitted background from fibre FF frame\n");
            fdvecbuf1 = framebuffer->frame_array[0];
            fdvecbuf2 = backframe[0];
            for (iyixindex=0; iyixindex<=iyixend; iyixindex++) {
                fdvecbuf1[iyixindex] -= fdvecbuf2[iyixindex];
            }
            /* 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 */
            fmvecbuf1 = framebuffer->badpixel[0];
            fdvecbuf1 = framebuffer->frame_array[0];
            fdvecbuf2 = framebuffer->frame_sigma[0];
            for (iyixindex=0; iyixindex<=iyixend; iyixindex++) {
                if (fmvecbuf1[iyixindex]==0 && (pixelvalue=fdvecbuf1[iyixindex])<0) {
                    if ((pixelvalue*pixelvalue)<4*(fdvecbuf2[iyixindex])) {
                        fdvecbuf1[iyixindex]=0;
                    }
                    else {
                        fmvecbuf1[iyixindex]=1;
                    }
                }
            }
        }
        else {
            SCTPUT("Error: no regions available for background estimation");
            return flames_midas_fail();
        }
    }
    /* free the background frame */
    free_fdmatrix(backframe, 0, framebuffer->subrows-1, 0,
                  framebuffer->subcols-1);
    /* free the framebuffer */
    free_cvector(framebuffer->fibremask, 0, framebuffer->maxfibres-1);
    free(framebuffer);

    SCTPUT("Background fitting and subtraction done");

    /* begin looping over orders at first */
    for (iorder=0; iorder <= ((ordpos->lastorder)-(ordpos->firstorder));
                    iorder++) {
        if ( strcmp(drs_verbosity,"LOW") == 0 ){
        } else {
            sprintf(output, "First cleaning step, order %d", iorder+1);
            SCTPUT(output);
        }
        lvecbuf1 = allflatsin->lowfibrebounds[iorder][0];
        int32_t* lvecbuf2 = allflatsin->highfibrebounds[iorder][0];
        fmvecbuf1 = allflatsin->goodfibres[iorder][0];
        /* loop over FF frames */
        for (iframe=0; iframe<=allflatsin->nflats-1;iframe++) {
            fdvecbuf1 = allflatsin->flatdata[iframe].data[0];
            fdvecbuf2 = allflatsin->flatdata[iframe].sigma[0];
            fmvecbuf2 = allflatsin->flatdata[iframe].badpixel[0];
            /* loop over lit fibres in this frame */
            for (lfibre=0; lfibre<=((allflatsin->flatdata)[iframe]).numfibres-1;
                            lfibre++) {
                ifibre = (((allflatsin->flatdata)[iframe]).fibres)[lfibre];
                iorderifibreindex = (iorder*allflatsin->maxfibres)+ifibre;
                /* loop over ix... */
                for (ix=0; ix<=(allflatsin->subcols-1); ix++){
                    ifibreixindex = (ifibre*allflatsin->subcols)+ix;
                    /* initialise data structures */
                    if (fmvecbuf1[ifibreixindex]==GOODSLICE) {
                        iorderifibreixindex = (iorderifibreindex*allflatsin->subcols)+ix;
                        goodpixels = 0;
                        /* loop over pixels belonging to the fibre */
                        for (iy=lvecbuf1[ifibreixindex];
                                        iy<=lvecbuf2[ifibreixindex];
                                        iy++){
                            iyixindex = (iy*allflatsin->subcols)+ix;
                            if (fmvecbuf2[iyixindex] == 0) {
                                /* this pixel is good, add to accumulators */
                                goodpixels++;
                            }
                        }
                        /* check goodpixels coverage of the fibre */
                        if ((double)goodpixels*allflatsin->substepy /
                                        (2*allflatsin->halfibrewidth) < allflatsin->minfibrefrac) {
                            /* the coverage of this fibre at this ix is too low, mark it
		 bad and disregard it altogether */
                            fmvecbuf1[iorderifibreixindex] = BADSLICE;
                            for (iy=lvecbuf1[ifibreixindex];
                                            iy<=lvecbuf2[ifibreixindex];
                                            iy++) {
                                iyixindex = (iy*allflatsin->subcols)+ix;
                                fmvecbuf2[iyixindex] = 1;
                            }
                        }
                    }
                }
            }
        }
    }

    /* ok, now let's check whether there is some fibre which is nominally
     lit but which is actually broken or otherwise unlit */
    fmvecbuf1 = allflatsin->goodfibres[0][0];
    for (ifibre=0; ifibre<=(allflatsin->maxfibres-1); ifibre++) {
        if (allflatsin->fibremask[ifibre]==TRUE) {
            totslices = fracslices = 0;
            for (iorder=0; iorder<=(ordpos->lastorder-ordpos->firstorder);
                            iorder++) {
                iorderifibreindex = (iorder*allflatsin->maxfibres)+ifibre;
                for (ix=0; ix<=(allflatsin->subcols-1); ix++) {
                    iorderifibreixindex = (iorderifibreindex*allflatsin->subcols)+ix;
                    totslices++;
                    if (fmvecbuf1[iorderifibreixindex] == GOODSLICE) {
                        /* there is at least one good slice for this fibre, keep it */
                        fracslices++;
                    }
                }
            }
            fracslices /= totslices;
            if (fracslices < fracslicesthres) {
                /* we were unable to find any good slice for this fibre, discard it */
                allflatsin->fibremask[ifibre] = FALSE;
                ordpos->fibremask[ifibre] = FALSE;
                sprintf(output, "fracslices=%f Warning: ignoring fibre %d which has \
poor coverage", fracslices, ifibre+1);
                SCTPUT(output);
            }
        }
    }

    /* reset the fibres and numfibres terms within allflatsin, they may have
     been changed */

    for (iframe=0; iframe<=(allflatsin->nflats-1); iframe++) {
        newnumfibres = 0;
        newfibres = lvector(0,allflatsin->maxfibres-1);
        for (lfibre=0; lfibre<=allflatsin->flatdata[iframe].numfibres-1;
                        lfibre++) {
            ifibre = allflatsin->flatdata[iframe].fibres[lfibre];
            if (allflatsin->fibremask[ifibre]==TRUE) {
                /* this fibre is ok, add it */
                newfibres[newnumfibres] = ifibre;
                newnumfibres++;
            }
        }
        /* ok, free the old fibres and set it to the new vector */
        free_lvector(allflatsin->flatdata[iframe].fibres,
                     0, allflatsin->maxfibres-1);
        allflatsin->flatdata[iframe].fibres = newfibres;
        allflatsin->flatdata[iframe].numfibres = newnumfibres;
    }


    /* free dynamically allocated temporary variables here */

    return NOERR;

}
/**@}*/