/* @(#)flames_tracing.c    */
/*===========================================================================
  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
===========================================================================*/
/* Program  : flames_tracing.c                                             */
/* Author   : I. Porceddu  -  ITAL_FLAMES Consortium                       */
/* Date     :                                                              */
/*                                                                         */
/* Purpose  :                                                              */
/* Follow orders. Create and fill in table ORDER                           */
/*                                                                         */
/* Input    :                                                              */
/*      IN_A      Name of Flat-Field frame                                 */
/*      IN_B      Name of input table                                      */
/*      OUT_A     Name of order table                                      */
/*      INPUTI(1) Step (pixels)                                            */
/*      INPUTI(2) Lower scan limit                                         */
/*      INPUTI(3) Upper scan limit                                         */
/*      INPUTR(1) Hot threshold                                            */
/*                                                                         */
/*                                                               */
/* Output:                                                       */
/*    tracing_orders table                                               */
/*                                                                         */
/* DRS Functions called:                                                   */
/* Missing                                                                 */
/*                                                                         */
/* Pseudocode:                                                             */
/* Missing                                                                 */ 
/*                                                                         */ 
/* Version  :                                                              */
/* Last modification date: 2002/08/05                                      */
/* Who     When        Why                Where                            */
/* AMo     02-08-05   Add header         header                            */
/*-------------------------------------------------------------------------*/
#ifdef HAVE_CONFIG_H
#  include <config.h>
#endif
#include <flames_tracing.h>
#include <flames.h>

#include <flames_midas_def.h>
#include <flames_def_drs_par.h>
#include <flames_uves.h>
#include <flames_newmatrix.h>
#include <flames_midas_tblsys.h>
#include <flames_midas_tblerr.h>
#include <flames_midas_tbldef.h>
#include <flames_midas_macrogen.h>
#include <flames_midas_atype.h>
#include <uves_msg.h>
#include <math.h>
#include <string.h>
#include <stdio.h>
#include <stdlib.h>


//jmlarsen: If we define this ipos macro correctly, then my gaussian fitting code
//          may start working!
//#define ipos(col,row,siz)   row * siz + col  /* Pointer value */
#define ipos(col,row,siz)  ((row)*(siz) + (col))
#define nint(f)             (int) (f + 0.5)  /* Nearest integer */


static int 
find_center (const int x, 
             const double slop,
             const double inter,
             const double low_th,
             const double high_th,
             const int direct,
             double * ycent);

static flames_err 
follow (const int order, 
        double slope_loc, 
        double intercept_loc, 
        const double ord_thres, 
        const double hot_thres, 
        const int step, 
        int xprev);

static flames_err 
store (const int x, const double y, const int pos);

static int 
where_start (const double slope_loc, 
             const double interc, 
             int * ordsta_loc, 
             int * ordend_loc, 
             int * ordlen_loc);

static double 
estim_thresh (const int order, 
              const int nbord, 
              const double hot_th, 
              const int xcent, 
              double * minith);
/*
static flames_err 
estim_thres2 (const double slope, const double interc);
 */

static flames_err 
update_out (const int tid, 
            const int order,
            const int fibre,
            const int ordcol,
            const int fibrecol,
            const int xcol,
            const int ycol,
            const int ordfibcol,
            const int stepix);


static float *pntra=0;
static int nrow;
static int ncol;
static int stkmin;
static int stkmax;
static int tbrow;
static int scan[2]={0,0};
static int maxorder;
static int midorder;
static float *slope;
static float *intercept;
//jmlarsen: not used, float *fwhm;
static float *userthres;
static int *ordsta;
static int *ordend;
static int *ordernum;
static int *fibrenum;
static int *selected;
static float *xstack;
static float *ystack;

static double 
Center  (const int col_loc, 
         int row_loc,
         const double low_th, 
         const double high_th, 
         int* ylow, 
         int* yupp)
{

    double pixel=0;
    double ycent=0;
    double ysum=0;
    double sum=0;

    //extern int nrow;
    //extern int ncol;
    //extern float *pntra;

    /* Moves to the lower limit of the order */
    while ((row_loc > 0) && (pntra[ipos(col_loc, row_loc, ncol)] >= low_th))  row_loc--;

    /* Computes center of gravity between lower and upper limit */
    sum  = 0.;
    ysum = 0.;
    *ylow = row_loc++;
    while ((row_loc < nrow) && ((pixel = pntra[ipos(col_loc, row_loc, ncol)]) >= low_th)) {
        if (pixel < high_th) {
            ysum += row_loc*(pixel-low_th);
            sum  += (pixel-low_th);
        }
        *yupp = row_loc++;
    }

    ycent = ysum/sum;

    *ylow += 1;   /* Convert from array index to pixel location */
    *yupp += 1;


    /* jmlarsen: Using the center of gravity is an inferior method.
     Fit a Gaussian. */
    if (*yupp - *ylow < 4) {
        *yupp += 2;
        *ylow += 2;
    }


    if ((*yupp - *ylow) > 4) {
        int N = *yupp - *ylow + 1;
        //const bool horizontal = true;
        cpl_vector *vx = cpl_vector_new(N);
        cpl_vector *vy = cpl_vector_new(N);
        double y;
        double sigma;
        double area;
        double offset;
        int i;

        for (i = 0; i < N; i++) {
            cpl_vector_set(vx, i, i + (*ylow-1));
            cpl_vector_set(vy, i, pntra[ipos(col_loc, i + (*ylow-1), ncol)]);

            /*          fprintf(stderr, "(%d, %d)   %f\n",
                  col_loc,
                  i + (*ylow-1),
                  pntra[ipos(col_loc, i + (*ylow-1), ncol)]);
             */
        }

        cpl_vector_fit_gaussian(vx, NULL,
                                vy, NULL,
                                CPL_FIT_ALL,
                                &y,
                                &sigma,
                                &area,
                                &offset,
                                NULL, NULL, NULL);

        cpl_vector_delete(vx);
        cpl_vector_delete(vy);
#if FLAMES_DEBUG
        fprintf(stderr, "Before, after, error = %f %f %d\n",
                ycent, y, cpl_error_get_code());
#endif

        if (cpl_error_get_code() != CPL_ERROR_NONE) {
            cpl_error_reset();
        }
        else if ((y-ycent)*(y-ycent) < 1*1) {
            /* Fit succeeded and new position is within one
             pixel from previous */
            ycent = y;
        }
    }

    return(ycent);
}


static int
ReadInput(char name[])
{
    int Tid=0;
    int slope_col=0;
    int inter_col=0;
    //  int fwhm_col=0;
    int order_col=0;
    int thres_col=0;
    int fibre_col=0;
    int nbcol=0;
    int nbrow=0;
    int null=0;
    int nb_order=0;

    int row;



    TCTOPN (name, F_I_MODE, &Tid);

    TCIGET (Tid, &nbcol, &nbrow);

    TCCSER (Tid, "SLOPE", &slope_col);
    TCCSER (Tid, "ORIG",  &inter_col);
    //  TCCSER (Tid, "FWHM",  &fwhm_col);
    TCCSER (Tid, "THRES", &thres_col);
    TCCSER (Tid, "NEWORD", &order_col);
    TCCSER (Tid, "FIBRE", &fibre_col);
    //  nb_order = nbrow;

    //jmlarsen:
    //Changed reading of table to be robust to unselected rows in middle of table

    nb_order = 0;
    for (row=1; row<=nbrow; row++) {
        TCSGET (Tid, row, &selected[nb_order+1]);
        if (selected[nb_order+1]) {
            nb_order += 1;
            TCERDR (Tid, row, slope_col, &slope[nb_order],     &null);
            TCERDR (Tid, row, inter_col, &intercept[nb_order], &null);
            //      TCERDR (Tid, row, fwhm_col,  &fwhm[nb_order],      &null);
            TCERDR (Tid, row, thres_col, &userthres[nb_order], &null);
            TCERDI (Tid, row, order_col, &ordernum[nb_order],  &null);
            if (null == 1) selected[nb_order] = FALSE;
            TCERDI (Tid, row, fibre_col, &fibrenum[nb_order],  &null);
            if (null == 1) selected[nb_order] = FALSE;
        }
    }


    TCTCLO(Tid);
    return(nb_order);
} 


static int 
where_start (const double slope_loc, 
             const double interc, 
             int * ordsta_loc, 
             int * ordend_loc, 
             int * ordlen_loc)
{

    int xsta=0;
    int ysta=0;
    int xend=0;
    int yend=0;
    int xcent=0;
    /* extern  int   nbrow, nbcol; */
    //extern int scan[2];

    if (slope_loc < -DEPSILON) {
        yend  = nint(interc);
        ysta  = nint(slope_loc*ncol+interc);
        xsta  = (yend > scan[0])    ?  1   : nint((scan[0]-interc)/slope_loc);
        xend  = (ysta < scan[1])    ? ncol : nint((scan[1]-interc)/slope_loc);
    }
    else if (slope_loc > DEPSILON) {
        ysta  = nint(interc);
        yend  = nint(slope_loc*ncol+interc);
        xsta  = (ysta > scan[0])    ?  1   : nint((scan[0]-interc)/slope_loc);
        xend  = (yend < scan[1])    ? ncol : nint((scan[1]-interc)/slope_loc);
    }
    else {
        ysta  = yend = nint(interc);
        if (ysta > scan[0] && yend < scan[1]) {
            xsta = 1;
            xend = ncol;
        }
        else {
            xsta = xend = scan[0]-1;
        }
    }
    xcent = (xsta + xend)/2;
    *ordlen_loc = xend - xsta;
    *ordsta_loc = xsta;
    *ordend_loc = xend;
    return(xcent);

}


static double 
estim_thresh (const int order, 
              const int nbord, 
              const double hot_th, 
              const int xcent, 
              double * minith)
{

    //extern int ncol;
    //extern float *pntra;
    //extern float *slope;
    //extern float *intercept;
    //extern int scan[];

    int row=0;

    int ylow=0;
    int yupp=0;
    int ycent=0;
    int ymax=0;
    int ymin=0;
    float mini=0;
    float maxi=0;
    double threshold=0;

    ycent = nint (slope[order] * xcent + intercept[order] - 1);
    ymax  = scan[1] - 1;
    ymin  = scan[0] - 1;

    if (order < nbord) {
        yupp  = nint (slope[order+1] * xcent + intercept[order+1])-1;
        yupp  = (ycent + yupp)/2;
    }
    else
        yupp  = ymax;   /* order == nbord */

    if (yupp > ymax) yupp = ymax;

    if (order > 1) {
        ylow  = nint (slope[order-1] * xcent + intercept[order-1])-1;
        ylow  = (ylow + ycent)/2;
    }
    else
        ylow  = ymin;  /* order == 1 */

    if (ylow < ymin)  ylow = ymin;

    /* Init */

    mini = maxi = pntra[ipos(xcent,ylow,ncol)];

    /* Find the minimum and maximum value */

    for (row=ylow; row<yupp; row++) {
        int pos = ipos(xcent,row,ncol);
        if (pntra[pos] > maxi) maxi   = pntra[pos];
        if (pntra[pos] < mini) mini = pntra[pos];
    }

    /* some tuning is in order here: to trace a larger portion
     of each fibre, the threshold needs to be low, but setting 
     it too low may cause the wrong fibre to be selected; a 
     decent value must be found by trial and error */

    threshold = (maxi - mini)*0.50 + mini;
    *minith   = (threshold - mini)*0.2 + mini;
    //jmlarsen: As in UVES, this minimum threshold is based on the center column
    //and too high
    *minith   = mini;


    return(threshold);

}


static flames_err 
follow (const int order, 
        double slope_loc, 
        double intercept_loc, 
        const double ord_thres, 
        const double hot_thres, 
        const int step, 
        int xprev)
{

    int xnext=0;
    int xini=0;
    int eof=0;
    int direct=1;

    double ycent=0;

    double ynext=0;
    double yini=0;
    double slope_ini=0;
    double inter_ini=0;
    double updateweight=.02;
    /* extern int  ncol; */

    eof = find_center (xprev, slope_loc, intercept_loc, ord_thres, hot_thres,
                    direct, &ycent);

    store(xprev, ycent, midorder);

    xini = xprev;
    yini = ycent;
    slope_ini = slope_loc;
    inter_ini = intercept_loc;

    for (direct = -1; direct<=1; direct += 2) {

        xprev = xini;
        double yprev = yini;
        int pos   = midorder;
        eof   = 0;      /* Init of EndOfFrame flag */
        slope_loc     = slope_ini;
        intercept_loc = inter_ini;

        while (!eof) {
            xnext = xprev + direct*step;
            eof = find_center (xnext, slope_loc, intercept_loc, ord_thres, hot_thres,
                            direct, &ynext);
            if (!eof) {
                pos += direct;
                store (xnext, ynext, pos);
                /* in the case of UVES, the orders are almost straight, with
       a tiny curvature; therefore, we adjust slope and intercept 
       slowly, to reduce the possibility to jump fibre altogether */
                slope_loc = (1-updateweight)*slope_loc+
                                updateweight*(ynext - yprev)/(xnext - xprev);
                intercept_loc = (1-updateweight)*intercept_loc+
                                updateweight*(xnext*yprev - xprev*ynext)/(xnext - xprev);
                /* we just adjust the intercept so that the straight line
       approximation of the order is forced to go through the 
       previously found point
                 */
                /* intercept_loc = yprev-slope_loc*xprev; */
                xprev = xnext;
                yprev = ynext;
            }
        }
    }
    return(0);
}

static int 
find_center (const int x, 
             const double slop,
             const double inter,
             const double low_th,
             const double high_th,
             const int direct,
             double * ycent)
{

    float ymid=0;



    int ylow=0;
    int yupp=0;

    /*  double  Center(); */

    //extern int ncol;
    //extern int scan[2];
    //extern float *pntra;

    ymid = slop*x+inter;

    if (1<=x && x<=ncol && scan[0]<=ymid && ymid<=scan[1]) {
        int col = x - 1;
        int row = nint(ymid) - 1;

        /* Find row of maximum signal connex to the center */
        int pass=1;
        int pos = 0;
        while(pass != 0) {
            pos  = ipos(col, row, ncol);
            pass = 0;
            if (pntra[ipos(col, (row+1), ncol)] > pntra[pos])  {pass=1; row++;}
            if (pntra[ipos(col, (row-1), ncol)] > pntra[pos])  {pass=1; row--;}
        }

        if (pntra[pos] > (float) low_th) {

            *ycent = Center(col,row,low_th,high_th,&ylow,&yupp)+1.;

            if (ylow > scan[0] && yupp < scan[1])
                return(0); /* Centering correctly done */
            else return(1);
        }

        else return(1); /* No more signal above the threshold */
    }

    else  return(1); /* Reached edge of frame */

}



static flames_err 
store (const int x, const double y, const int pos)
{
    //extern int stkmin;
    //extern int stkmax;
    //extern float *xstack;
    //extern float *ystack;

    xstack[pos] = x;
    ystack[pos] = y;

    if (pos < stkmin) stkmin = pos;
    if (pos > stkmax) stkmax = pos;
    return(0);
}


static flames_err 
update_out (const int tid, 
            const int order,
            const int fibre,
            const int ordcol,
            const int fibrecol,
            const int xcol,
            const int ycol,
            const int ordfibcol,
            const int stepix)

{
    //extern int tbrow;
    //extern int stkmin;
    //extern int stkmax;
    //extern float *xstack;
    //extern float *ystack;

    int stkpntr=0;
    int xempty=1;
    char orderfibre[20];

    memset(orderfibre, '\0', 20);
    sprintf(orderfibre,"%d,%d",order,fibre);

    while(xempty < xstack[stkmin]) {
        TCEWRI (tid, tbrow, fibrecol, &fibre);
        TCEWRI  (tid, tbrow,   ordcol, &order);
        TCEWRC (tid, tbrow, ordfibcol, orderfibre);
        TCEWRI  (tid, tbrow++, xcol,   &xempty);
        xempty += stepix;
    }

    for (stkpntr=stkmin; stkpntr<=stkmax; stkpntr++) {

        if(ystack[stkpntr] != 0) {
            TCEWRI (tid, tbrow, ordcol, &order);
            TCEWRI (tid, tbrow, fibrecol, &fibre);
            TCEWRC (tid, tbrow, ordfibcol, orderfibre);
            TCEWRR (tid, tbrow, xcol, &xstack[stkpntr]);
            TCEWRR (tid, tbrow, ycol, &ystack[stkpntr]);
            tbrow++;}
        else
            tbrow++;
    }

    xempty = xstack[stkmax] + stepix;
    while (xempty < ncol) {
        TCEWRI  (tid, tbrow,   ordcol, &order);
        TCEWRI (tid, tbrow, fibrecol, &fibre);
        TCEWRC (tid, tbrow, ordfibcol, orderfibre);
        TCEWRI  (tid, tbrow++, xcol,   &xempty);
        xempty += stepix;
    }

    TCEWRI  (tid, tbrow,   ordcol, &order);
    TCEWRI (tid, tbrow, fibrecol, &fibre);
    TCEWRC (tid, tbrow, ordfibcol, orderfibre);
    TCEWRI  (tid, tbrow++, xcol,   &ncol);

    stkmin = stkmax = midorder;
    return(0);
}

/*
static flames_err 
estim_thres2 (const double slope_loc, const double interc)

{

  int xsta=0;
  int ysta=0;
  int xend=0;
  int yend=0;
  int xcent=0;
  //extern int nrow;
  //extern int ncol;
  //extern int scan[2];

  ysta  = nint(interc);
  xsta  = (ysta > 1)    ?  1   : nint((scan[0]-interc)/slope_loc);

  yend  = nint(slope_loc*ncol+interc);
  xend  = (yend < nrow) ? ncol : nint((scan[1]-interc)/slope_loc);

  xcent = (xsta + xend)/2;
  return(xcent);

}
 */

/* Returns: table handle */
static
int flames_create_ordertable(const char *OTAB)
{
    char ordtab[TEXT_LEN+1];

    int unit=0;
    int null=0;
    int ordcol=0;
    int xcol=0;
    int ycol=0;
    int yfitcol=0;
    int deltacol=0;
    int fibrecol=0;
    int ordfibcol=0;
    int tid=0;
    int actvals=0;
    int maxrows=0;
    int maxcols=0;

    memset(ordtab, '\0', TEXT_LEN+1);

    //  SCSPRO("create_ordertable");

    SCKGETC(OTAB, 1, 60, &actvals, ordtab); /* User defined */

    SCKRDI(&MAXROWS, 1, 1, &actvals, &maxrows, &unit, &null);
    SCKRDI(&MAXCOLS, 1, 1, &actvals, &maxcols, &unit, &null);

    TCTINI (ordtab, F_O_MODE, maxrows, &tid);
    TCCINI (tid, D_R4_FORMAT, 1, "I6", "  ", "ORDER", &ordcol);
    TCCINI (tid, D_R4_FORMAT, 1, "I6", "  ", "X",     &xcol);
    TCCINI (tid, D_R4_FORMAT, 1, "I6", "  ", "Y",     &ycol);

    TCCINI(tid, D_R4_FORMAT, 1, "f7.2", " ", "YFIT", &yfitcol);
    TCCINI(tid, D_R4_FORMAT, 1, "f7.2", " ", "RESIDUAL", &deltacol);

    TCCINI (tid, D_I4_FORMAT, 1, "I6", "  ", "FIBRE", &fibrecol);

    TCCINI (tid, D_C_FORMAT, 21, "A", "  ", "ORDERFIB", &ordfibcol);


    //  TCTCLO(tid);
    //
    //  SCSEPI();
    //  return(0);
    return tid;
}


int flames_tracing(const char *IN_A,
                   const char *IN_B,
                   const char *OTAB,
                   const int *INPUTI,
                   const float *INPUTR,
                   const int *MAX_ORDER)
{
    char frame[TEXT_LEN+1];
    char inptab[TEXT_LEN+1];
    char ordtab[TEXT_LEN+1];
    //char ident[TEXT_LEN+1];
    //char cunit[16*3 + 1];
    char text[TEXT_LEN+1];


    int order=0;
    int nb_order=0;
    int stepix=0;
    int stkdif=0;
    int stklimit=0;
    int counter=0;
    int ordlen=0;

    //int imnoa=0;
    //int naxis=0;
    //int npix[2]={0,0};
    int null=0;
    int kunit=0;
    int actvals=0;
    int tid=0;
    int ordcol=0;
    int fibrecol=0;
    int xcol=0;
    int ycol=0;
    int xcent=0;
    int i=0;
    int ordfibcol=0;
    int yfitcol=0;
    int residcol=0;

    //extern int nrow;
    //extern int ncol;
    //extern float *pntra;
    //extern float *slope;
    //extern float *intercept;
    //extern int *ordsta;
    //extern int *ordend;
    //extern int scan[];
    //extern int *selected;

    float hot_thres=0;
    float thres[2]={0,0};
    double ord_thres=0;
    double mini_thres=0;

    //double start[2]={0,0};
    //double step[2]={0,0};

    char drs_verbosity[10];
    int status=0;

    //jmlarsen: Reset static memory so that the function
    //          is reentrant
    pntra=0;
    nrow=0;
    ncol=0;
    stkmin=0;
    stkmax=0;
    tbrow=1;
    //  scan[2]={0,0};
    scan[0]=0;
    scan[1]=0;
    maxorder=0;
    midorder=0;
    slope=0;
    intercept=0;
    //jmlarsen: not used, fwhm=0;
    userthres=0;
    ordsta=0;
    ordend=0;
    ordernum=0;
    fibrenum=0;
    selected=0;
    xstack=0;
    ystack=0;




    memset(frame, '\0', TEXT_LEN+1);
    memset(inptab, '\0', TEXT_LEN+1);
    memset(ordtab, '\0', TEXT_LEN+1);
    //memset(ident, '\0', TEXT_LEN+1);
    memset(text, '\0', TEXT_LEN+1);
    memset(text, '\0', 16*3+1);

    SCSPRO("flames_tracing");




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

    SCKGETC(IN_A, 1, 60,  &actvals, frame); /* middummi.bdf */
    SCKGETC(IN_B, 1, 60,  &actvals, inptab); /* middummr.tbl */
    SCKGETC(OTAB, 1, 60, &actvals, ordtab); /* User defined */

    if ( strcmp(drs_verbosity,"LOW") == 0 ){
    } else {
        char output[200];
        sprintf(output, "Frame: %s \t Inptab: %s \t Ordtab: %s \n",frame,
                        inptab,ordtab);
        SCTPUT(output);
    }
    SCKRDI(INPUTI, 1, 1, &actvals, &stepix,     &kunit, &null);
    SCKRDI(INPUTI, 2, 2, &actvals, &scan[0],    &kunit, &null);
    SCKRDR(INPUTR, 1, 1, &actvals, &thres[0],   &kunit, &null);

    hot_thres = thres[0];

    SCKRDI(MAX_ORDER, 1, 1, &actvals, &maxorder, &kunit, &null);
    if (maxorder<2) maxorder=2;
    midorder = maxorder/2;
    stkmin = stkmax = midorder;

    slope = calloc(maxorder, sizeof(float));
    intercept = calloc(maxorder, sizeof(float));
    //  fwhm =  calloc(maxorder, sizeof(float));
    userthres = calloc(maxorder, sizeof(float));
    xstack = calloc(maxorder, sizeof(float));
    ystack = calloc(maxorder, sizeof(float));
    ordsta = calloc(maxorder, sizeof(int));
    ordend = calloc(maxorder, sizeof(int));
    ordernum = calloc(maxorder, sizeof(int));
    fibrenum = calloc(maxorder, sizeof(int));
    selected = calloc(maxorder, sizeof(int));

    for (i=0; i<=maxorder-1; i++) {
        slope[i] = 0;
        intercept[i] = 0;
        //    fwhm[i] = 0;
        userthres[i] = 0;
        ordsta[i] = 0;
        ordend[i] = 0;
        ordernum[i] = 0;
        fibrenum[i] = 0;
        selected[i] = 0;
        xstack[i] = 0;
        ystack[i] = 0;
    }

    if (scan[1]<scan[0]) {
        nrow = scan[1];
        scan[1] = scan[0];
        scan[0] = nrow;
    }

    /* Fixed on  23/05/2000 to prevent going out of the frames boundaries */
    scan[0] = scan[0] + 2;
    scan[1] = scan[1] - 2;

    //jmlarsen: The following code loads an image buffer from the file 'frame'
    // to the array pntra
    //This is the only place where SCIGET is used, so don't bother to implement
    //this function
    //use CPL directly to do the same
    //not used: naxis, step, ident, cunit, imnoa

    /*  Old code:

      strcpy(ident, " ");
      strcpy(cunit, " ");
      SCIGET (frame, D_R4_FORMAT, F_I_MODE, F_IMA_TYPE, 2,
      &naxis, npix, start, step, ident, cunit, (char **)&pntra,
      &imnoa);

      nrow  =  npix[1];
      ncol  =  npix[0];

     */
    /* New code: */
    cpl_image *im=NULL;
    {
        int plane = 0;
        int extension = 0;

        im = cpl_image_load(frame, CPL_TYPE_FLOAT, plane, extension);
        pntra = cpl_image_get_data_float(im);
        nrow = cpl_image_get_size_y(im);
        ncol = cpl_image_get_size_x(im);
        //cpl_image_unwrap(im);
    }



    /* Opening table OTAB for output */
    //jmlarsen:
    //  TCTOPN(ordtab,F_IO_MODE,&tid);
    tid = flames_create_ordertable(ordtab);


    TCCSER (tid,"ORDER", &ordcol);
    TCCSER (tid, "FIBRE", &fibrecol);
    TCCSER (tid, "X", &xcol);
    TCCSER (tid, "Y", &ycol);
    TCCSER (tid, "ORDERFIB", &ordfibcol);
    TCCSER (tid, "YFIT", &yfitcol);
    TCCSER (tid, "RESIDUAL", &residcol);

    nb_order = ReadInput(inptab);

    // printf("Number of orders : %d\n",nb_order);
    uves_msg_debug("Number of orders : %d\n",nb_order);

    for (order=1; order<=nb_order; order++) {

        if (selected[order] == TRUE) {

            /* Defines the geometrical central position of the order/fibre */
            xcent     = where_start (slope[order], intercept[order],
                            &ordsta[order], &ordend[order], &ordlen);

            if (xcent>=1 && xcent<=ncol) {

                /* Defines the absolute threshold of the order */
                if (userthres[order]>0.001)  {
                    ord_thres  = userthres[order];
                    mini_thres = userthres[order];
                }
                else
                    ord_thres = estim_thresh (order, nb_order, hot_thres,
                                    xcent, &mini_thres);

                /* printf ("xcent, ordlen, thres, mini: %d %d %f %f\n",xcent,
       ordlen,ord_thres,mini_thres); */

                /* The threshold must be low enough to find at least */
                /* three quarters of the positions aint  each order  */

                stklimit = 9*ordlen/stepix/10;
                stkdif   = stklimit - 1;       /* Initial value */
                counter  = 0;                  /* Number of loops (maxi 10) */

                while (stkdif < stklimit && counter < 10) {

                    /* Follows the order to estimate order center and */
                    /* background center                              */
                    follow     (order, slope[order], intercept[order], ord_thres,
                                    hot_thres, stepix, xcent);

                    stkdif = stkmax - stkmin + 1;
                    counter++;

                    if (stkdif < stklimit)      ord_thres = 0.75 * ord_thres;
                    if (ord_thres < mini_thres) ord_thres = mini_thres;
                }

                /* Display information line */
                if ( strcmp(drs_verbosity,"LOW") == 0 ){
                } else {
                    sprintf (text,
                                    "Order: %3d  Fibre %3d  Threshold: %2f  Numb. of values: %d",
                                    ordernum[order], fibrenum[order], ord_thres, stkdif);
                    SCTPUT  (text);
                }
            }
            else {
                /* Display information line */
                if ( strcmp(drs_verbosity,"LOW") == 0 ){
                } else {

                    sprintf (text,"Order %3d,  Fibre %3d: out of boundaries",
                                    ordernum[order], fibrenum[order]);
                    SCTPUT  (text);
                }
            }

            /* Writes data in the output table */
            update_out(tid, ordernum[order], fibrenum[order], ordcol, fibrecol,
                       xcol, ycol, ordfibcol, stepix);

        }

    }

    SCDWRI(tid,"ORDSTA",&ordsta[1],1,nb_order,&kunit);
    SCDWRI(tid,"ORDEND",&ordend[1],1,nb_order,&kunit);

    TCTCLO(tid);
    cpl_image_delete(im);
    free(slope);
    free(intercept);
    //  free(fwhm);
    free(userthres);
    free(xstack);
    free(ystack);
    free(ordsta);
    free(ordend);
    free(fibrenum);
    free(ordernum);
    free(selected);

    return SCSEPI();

}