/*                                                                              *
 *   This file is part of the ESO UVES Pipeline                                 *
 *   Copyright (C) 2004,2005 European Southern Observatory                      *
 *                                                                              *
 *   This library 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, 51 Franklin St, Fifth Floor, Boston, MA  02111-1307  USA       *
 *                                                                              */

/*
 * $Author: amodigli $
 * $Date: 2013-10-09 12:11:14 $
 * $Revision: 1.34 $
 * $Name: not supported by cvs2svn $
 * $Log: not supported by cvs2svn $
 * Revision 1.33  2012/03/02 16:43:12  amodigli
 * fixed warning related to upgrade to CPL6
 *
 * Revision 1.32  2011/12/08 13:59:43  amodigli
 * Fox warnings with CPL6
 *
 * Revision 1.31  2011/03/25 07:44:19  amodigli
 * cleaned output
 *
 * Revision 1.30  2011/03/23 12:27:31  amodigli
 * changed QC key as user likes
 *
 * Revision 1.29  2011/03/23 10:08:31  amodigli
 * added QC to better characterize wave accuracy
 *
 * Revision 1.28  2010/09/24 09:32:10  amodigli
 * put back QFITS dependency to fix problem spot by NRI on FIBER mode (with MIDAS calibs) data
 *
 * Revision 1.26  2007/08/21 13:08:26  jmlarsen
 * Removed irplib_access module, largely deprecated by CPL-4
 *
 * Revision 1.25  2007/06/06 08:17:34  amodigli
 * replace tab with 4 spaces
 *
 * Revision 1.24  2007/05/02 13:20:01  jmlarsen
 * Take error bars into account in line searching if arclamp was flat-fielded
 *
 * Revision 1.23  2007/04/24 12:50:29  jmlarsen
 * Replaced cpl_propertylist -> uves_propertylist which is much faster
 *
 * Revision 1.22  2007/04/20 14:46:45  jmlarsen
 * Added commented out code
 *
 * Revision 1.21  2007/03/05 10:25:08  jmlarsen
 * Include slope in Gaussian fit
 *
 * Revision 1.20  2007/02/23 13:33:38  jmlarsen
 * Added code to test unweighted fitting
 *
 * Revision 1.19  2007/02/22 15:38:26  jmlarsen
 * Use linear background term in Gaussian fit
 *
 * Revision 1.18  2006/11/15 15:02:15  jmlarsen
 * Implemented const safe workarounds for CPL functions
 *
 * Revision 1.16  2006/11/15 14:04:08  jmlarsen
 * Removed non-const version of parameterlist_get_first/last/next which is already
 * in CPL, added const-safe wrapper, unwrapper and deallocator functions
 *
 * Revision 1.15  2006/11/06 15:19:42  jmlarsen
 * Removed unused include directives
 *
 * Revision 1.14  2006/08/18 13:51:01  jmlarsen
 * Moved one message from info to debug level
 *
 * Revision 1.13  2006/08/17 13:56:53  jmlarsen
 * Reduced max line length
 *
 * Revision 1.12  2006/08/17 09:18:47  jmlarsen
 * Removed CPL2 code
 *
 * Revision 1.11  2006/08/11 14:38:24  jmlarsen
 * Minor text output change
 *
 * Revision 1.10  2006/08/11 09:01:17  jmlarsen
 * Set unextracted bins to zero flux rather than marking them as bad
 *
 * Revision 1.9  2006/07/14 12:45:58  jmlarsen
 * Removed a few messages
 *
 * Revision 1.8  2006/07/03 13:29:30  jmlarsen
 * Adapted to new 1d-fitting function interface
 *
 * Revision 1.7  2006/06/13 12:05:11  jmlarsen
 * Shortened max line length
 *
 * Revision 1.6  2006/05/12 15:06:30  jmlarsen
 * Killed code for method = gravity
 *
 * Revision 1.5  2006/04/24 09:34:26  jmlarsen
 * Adapted to new interface of gaussian fitting routine
 *
 * Revision 1.4  2006/03/03 13:54:11  jmlarsen
 * Changed syntax of check macro
 *
 * Revision 1.3  2006/02/15 13:19:15  jmlarsen
 * Reduced source code max. line length
 *
 * Revision 1.2  2006/02/08 07:52:16  jmlarsen
 * Added function returning library version
 *
 * Revision 1.34  2006/01/12 15:41:14  jmlarsen
 * Moved gauss. fitting to irplib
 *
 * Revision 1.33  2005/12/20 08:11:44  jmlarsen
 * Added CVS  entry
 *
 */

/*----------------------------------------------------------------------------*/
/**
 * @addtogroup uves_wavecal
 */
/*----------------------------------------------------------------------------*/
/**@{*/

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

#include <uves_wavecal_search.h>
#include <uves_utils.h>
#include <uves_utils_wrappers.h>
#include <uves_utils_cpl.h>
#include <uves_pfits.h>
#include <uves_dump.h>
#include <uves_error.h>
#include <uves_msg.h>
#include <uves_qclog.h>

#include <cpl.h>
#include <float.h>

#define FIT_SLOPE 1
#define WEIGHTED_FIT 1   /* Define to zero to get unweighted fit of emmision line
                            (like MIDAS) */

static double
xcenter(const cpl_image *image, const cpl_image *noise, int xlo, int xhi, int row,
    centering_method CENTERING_METHOD, int bin_disp,
    double *sigma, double *intensity, double *dx0, double *slope, double *background);

static cpl_error_code
detect_lines(const cpl_image *spectrum, const cpl_image *noise, 
         const uves_propertylist *spectrum_header, 
             bool flat_fielded,
         int RANGE, double THRESHOLD, centering_method CENTERING_METHOD, 
             int bin_disp,
         const polynomial *order_locations, cpl_image *arcframe, 
         cpl_table *linetable, 
         int *ndetected, int *nrows);

/*----------------------------------------------------------------------------*/
/**
  @brief    Search for a given number of emission lines
  @param    spectrum          The spectrum to search.
  @param    noise             The noise (1 sigma) of the spectrum.
  @param    spectrum_header   Header of spectrum image. The (relative) order 
                              numbers are read from this header
  @param    flat_fielded      was the arclamp frame flat-fielded?
  @param    order_locations   The order position bi-polynomial (defines the
                              order locations as a function of @em x and
                  relative order number)
  @param    arcframe          The raw input image (before spectrum extraction). The
                              results of the line search will be drawn on this image.
  @param    RANGE             Width (in pixels) of search window is 2 * @em RANGE + 1
  @param    MINLINES          The minimum number of lines to search for
  @param    MAXLINES          The maximum number of lines to search for
  @param    CENTERING_METHOD  The method used to calculate the line position.
  @return   A newly allocated line table containing the detected lines, 
            or NULL on error

  This function tries to detect a number of emission lines in the already extracted @em
  spectrum. On success the number of detected lines is in the range from @em MINLINES
  to @em MAXLINES (both inclusive). This is achieved by calling @c detect_lines() .
  with varying detection threshold levels.

  The result line table contains the columns @em "Y" (the relative order number),
  @em "X" (the line positions), @em "Xwidth" (the detected line widths
  (stddev in pixels)) and @em "Peak" (the detected emission line intensities).

 */
/*----------------------------------------------------------------------------*/
cpl_table *
uves_wavecal_search(const cpl_image *spectrum, const cpl_image *noise,
                    const uves_propertylist *spectrum_header,
                    bool flat_fielded,
                    const polynomial *order_locations, cpl_image *arcframe,
                    int RANGE, int MINLINES, int MAXLINES,
                    centering_method CENTERING_METHOD,int bin_disp,
                    const int trace, const int window,  cpl_table* qclog)
{
    cpl_table *linetable = NULL;       /* Result */

    int nx, ny, norders;               /* Dimensions of raw image, number of orders */
    double threshold_low;              /* Threshold limits used for binary search */
    double threshold_high;
    double threshold = 0;             /* Current threshold */
    int lines_in_table;               /* Number of lines in line table */
    int lines_detected;               /* Number of lines actually found */
    bool max_thresh_found = false;    /* Is 'threshold_high' large enough? */
  

    passure( spectrum        != NULL, "Null input spectrum");
    passure( order_locations != NULL, "Null polynomial");
    passure( arcframe        != NULL, "Null raw image");

    if (flat_fielded) {
        assure( cpl_image_get_type(spectrum) == CPL_TYPE_DOUBLE,
                CPL_ERROR_TYPE_MISMATCH,
                "Spectrum image type is %s, must be double",
                uves_tostring_cpl_type(cpl_image_get_type(spectrum)));
    }
    
    check(( nx      = cpl_image_get_size_x(spectrum),
        norders = cpl_image_get_size_y(spectrum)), "Error reading input spectrum");
    check( ny      = cpl_image_get_size_y(arcframe), "Error reading input image");
    assure(nx == cpl_image_get_size_x(arcframe), CPL_ERROR_INCOMPATIBLE_INPUT, 
       "Spectrum and image widths are different (%d and %" CPL_SIZE_FORMAT ")",
       nx, cpl_image_get_size_x(arcframe));
    
    assure( MINLINES <= MAXLINES, CPL_ERROR_ILLEGAL_INPUT, 
        "minlines=%d maxlines=%d", MINLINES, MAXLINES );
    
    /* Initialize result line table */
    check(( linetable = cpl_table_new(MAXLINES),
        cpl_table_new_column(linetable, "X"     , CPL_TYPE_DOUBLE),
        cpl_table_new_column(linetable, "dX"    , CPL_TYPE_DOUBLE),
        cpl_table_new_column(linetable, "Xwidth", CPL_TYPE_DOUBLE),
        cpl_table_new_column(linetable, "Y"     , CPL_TYPE_INT),
        cpl_table_new_column(linetable, "Peak"  , CPL_TYPE_DOUBLE),
        cpl_table_new_column(linetable, "Background" , CPL_TYPE_DOUBLE),
        cpl_table_new_column(linetable, "Slope" , CPL_TYPE_DOUBLE)),
      "Could not create line table");
    cpl_table_set_column_unit(linetable,"X","pix");
    cpl_table_set_column_unit(linetable,"dX","pix");
    cpl_table_set_column_unit(linetable,"Xwidth","pix");
    cpl_table_set_column_unit(linetable,"Y","pix");
    cpl_table_set_column_unit(linetable,"Peak","ADU");
    cpl_table_set_column_unit(linetable,"Background","ADU");
    cpl_table_set_column_unit(linetable,"Slope","  ");
    uves_msg("Searching for emission lines");

    threshold_low  = 0.0;

    /* This start (guess) value is doubled until too few lines are detected */
    if (flat_fielded) {
        threshold_high = 10.0; /* dimensionless, number of stdevs above continuum */
    }
    else {
        threshold_high = cpl_image_get_mean(spectrum);

        assure( threshold_high > 0, CPL_ERROR_ILLEGAL_INPUT,
                "Spectrum median flux is %e. Must be positive",
                cpl_image_get_median(spectrum));
    }
    
    max_thresh_found = false;

    /* Detect lines and adjust threshold
       until MINLINES <= lines_detected <= MAXLINES */
    lines_detected = 0;


    char qc_key[40];
 
    int kk=0;
    while( (lines_detected < MINLINES || MAXLINES < lines_detected) && 
        fabs(threshold_low - threshold_high) > DBL_EPSILON )
    {
       kk++;
        threshold = (threshold_low + threshold_high)/2.0;

        check( detect_lines(spectrum, noise, spectrum_header,
                                flat_fielded,
                RANGE, threshold, CENTERING_METHOD,
                                bin_disp,
                order_locations,
                NULL,
                linetable,
                &lines_detected,
                &lines_in_table),
           "Could not search for emission lines");
        
        /* Update threshold */
        /* 'threshold_high' is doubled until the threshold is too high.
           Then a binary search is performed. */
        if (lines_detected < MINLINES)
        {
            max_thresh_found = true;
            threshold_high = threshold;
        }
        else if (MAXLINES < lines_detected) 
        {
            if (!max_thresh_found)
            {
                threshold_high *= 2;
            }
            else
            {                
                threshold_low = threshold;
            }
        }
        sprintf(qc_key,"QC TRACE%d WIN%d NLINDET%d",trace,window,kk);
        uves_msg_debug("ThAr lamp on trace %d window %d detected lines %d",
                 trace,window,lines_detected);
        ck0_nomsg(uves_qclog_add_int(qclog,qc_key,lines_detected,
                                        "ThAr lamp detected lines","%d"));




    } /* end while loop */

    sprintf(qc_key,"QC TRACE%d WIN%d NLINDET NITERS",trace,window);
    ck0_nomsg(uves_qclog_add_int(qclog,qc_key,kk,"Number of iterations","%d") );
    assure( MINLINES <= lines_detected && lines_detected <= MAXLINES, 
        CPL_ERROR_CONTINUE,
        "Could not detect between %d and %d lines. Try to increase search range",
        MINLINES, MAXLINES);
    
    /* Draw detections on input image  */
    check( detect_lines(spectrum, noise, spectrum_header,
                        flat_fielded,
            RANGE, threshold, CENTERING_METHOD,
                        bin_disp,
            order_locations,
            arcframe,
            linetable,
            &lines_detected,
            &lines_in_table),
       "Could not search for emission lines");
    
    /* Remove the last part of the line table (garbage) */
    check( cpl_table_set_size(linetable, lines_in_table), 
       "Could not resize line table");
    
    uves_sort_table_1(linetable, "X", false);
    
  cleanup:
#if 0 /* if flat-fielded */
    uves_free_image(&temp);
#endif
    if (cpl_error_get_code() != CPL_ERROR_NONE)
    {
        uves_free_table(&linetable);
    }
    else
    {
        /* Returned is... */
        passure( cpl_table_get_ncol(linetable) == 7, "%" CPL_SIZE_FORMAT "",
             cpl_table_get_ncol(linetable));
        passure( cpl_table_has_column(linetable, "X"     ), " ");
        passure( cpl_table_has_column(linetable, "dX"    ), " ");
        passure( cpl_table_has_column(linetable, "Xwidth"), " ");
        passure( cpl_table_has_column(linetable, "Y"     ), " ");
        passure( cpl_table_has_column(linetable, "Peak"  ), " ");
        passure( cpl_table_has_column(linetable, "Background" ), " ");
        passure( cpl_table_has_column(linetable, "Slope" ), " ");
        
    }
    return linetable;
}

/*----------------------------------------------------------------------------*/
/**
  @brief    Find emission lines above a certain threshold
  @param    spectrum          The spectrum to search.
  @param    noise             The noise (1 sigma) of the spectrum
  @param    flat_fielded      was the spectrum flat-fielded?
  @param    spectrum_header   Header of spectrum image. The (relative) order 
                              numbers are read from this header
  @param    RANGE             Width (in pixels) of search window is 2 * @em RANGE + 1
  @param    THRESHOLD         Detection threshold (relative to local background level)
  @param    CENTERING_METHOD  The method used to calculate the line position.
  @param    bin_disp          binning in dispersion direction
  @param    arcframe          The (raw) arc frame image. Iff non-NULL, the results
                              of the line search will be drawn on this image
                  using the provided @em order_locations.
  @param    order_locations   The order position bi-polynomial (defines the 
                              order locations as a function of @em x and
                  relative order number)
  @param    linetable         The line table where detected emission lines 
                              will be written
  @param    ndetected         (output) Number of lines detected
  @param    nrows             (output) Number of lines written to the line table
  @return   CPL_ERROR_NONE iff okay

  This function searches for emission lines at each row of the input image 
  (an extracted spectrum).

  The detected lines will be written to the columns @em "Y" (the relative
  order number), @em "X" (the line position), @em "Xwidth" (the detected line 
  width (stddev in pixels)) and @em "Peak" (the detected emission line intensity).

  A line is detected iff the flux at a position is more than @em THRESHOLD
  plus the local background level, which is defined as the median of a window
  of width (2 * @em RANGE + 1) centered on the current position. If flat_fielding
  was done @em THRESHOLD is the peak height in units of error bars.

  After the line is initially detected, the line position is calculated depending
  on the @em CENTERING_METHOD. See also @c xcenter() .

  The number of detected lines, @em ndetected, will be different from the number
  of rows written to the line table, @em nrows, only if the table is not large
  enough to contain all the detected lines.

  Finally, doublets (i.e. lines with positions within 2.0 pixels) are removed 
  from the line table. The maximum number of lines to search for is defined by
  the number of available rows in the line table.

  If the provided @em arcframe is not NULL, the results of the line search 
  will be marked on this image (used for debugging purposes).
  
 */
/*----------------------------------------------------------------------------*/
static cpl_error_code
detect_lines(const cpl_image *spectrum, const cpl_image *noise, 
         const uves_propertylist *spectrum_header, 
             bool flat_fielded,
         int RANGE, double THRESHOLD, centering_method CENTERING_METHOD, 
             int bin_disp,
         const polynomial *order_locations, cpl_image *arcframe, 
         cpl_table *linetable, 
         int *ndetected, int *nrows)
{
    int norders;      /* Number of orders */
    int minorder;     /* Relative order number of first row in spectrum image */
    int MAXLINES;     /* Number of rows in line table (max no. of 
             lines to search for) */
    int nx;           /* Width of spectrum (and raw image) */
    int x, order;     /* 'order' always counts from 1 */
    
    const double *spectrum_data;
    const double *noise_data;

    /* Check input */
    passure( spectrum        != NULL, " ");
    passure( noise           != NULL, " ");
    passure( spectrum_header != NULL, " ");
    nx      = cpl_image_get_size_x(spectrum);
    norders = cpl_image_get_size_y(spectrum);
    
    /* For efficiency reasons, get direct pointer to buffer,
       support only CPL_TYPE_DOUBLE */
    assure( cpl_image_get_type(spectrum) == CPL_TYPE_DOUBLE,
        CPL_ERROR_UNSUPPORTED_MODE,
        "Image type must be double. It is %s", 
        uves_tostring_cpl_type(cpl_image_get_type(spectrum)));

    spectrum_data = cpl_image_get_data_double_const(spectrum);
    noise_data    = cpl_image_get_data_double_const(noise);

    passure( RANGE > 0, "%d", RANGE);

    if (arcframe != NULL)
    {
        passure( order_locations != NULL, " ");
        passure( nx == cpl_image_get_size_x(arcframe), 
             "%d %" CPL_SIZE_FORMAT "", nx, cpl_image_get_size_x(arcframe));
    }
    
    passure( linetable != NULL, " ");
    MAXLINES = cpl_table_get_nrow(linetable);
    passure( cpl_table_get_ncol(linetable) == 7, "%" CPL_SIZE_FORMAT "",
         cpl_table_get_ncol(linetable));
    passure( cpl_table_has_column(linetable, "X"     ), " ");
    passure( cpl_table_has_column(linetable, "dX"    ), " ");
    passure( cpl_table_has_column(linetable, "Xwidth"), " ");
    passure( cpl_table_has_column(linetable, "Y"     ), " ");
    passure( cpl_table_has_column(linetable, "Peak"  ), " ");
    passure( cpl_table_has_column(linetable, "Background" ), " ");
    passure( cpl_table_has_column(linetable, "Slope" ), " ");
    
    assure( THRESHOLD > 0, CPL_ERROR_ILLEGAL_INPUT, "Illegal threshold: %e",
        THRESHOLD);

    check( minorder = uves_pfits_get_crval2(spectrum_header), 
       "Error reading order number of first row");

    *ndetected = 0;    /* Counts the number of lines detected so far. */
    *nrows = 0;        /* A pointer to the first unused row in the
                  line table */
    
    for (order = minorder; order < minorder + norders; order++) {
        int spectrum_row = order - minorder + 1;
        int ndetected_order = 0;
        for (x = 1; x <= nx; x++) {
        double flux, dflux;
        int peak_width = 0;
        int xlo, xhi;
        double local_median;
        
        /* Check if there is a peak and determine its position and width */
        // flux = cpl_image_get(spectrum, x, spectrum_row, &pis_rejected);
        flux  = spectrum_data[(x-1) + (spectrum_row - 1) * nx];
                dflux = noise_data   [(x-1) + (spectrum_row - 1) * nx];
        
        xlo = uves_max_int(x - RANGE, 1);
        xhi = uves_min_int(x + RANGE, nx);
        
        local_median = cpl_image_get_median_window(
            spectrum,
            uves_max_int(xlo, 1 ), spectrum_row,
            uves_min_int(xhi, nx), spectrum_row);
        
        while(x <= nx && 
                      (
                          (!flat_fielded && flux - local_median > THRESHOLD) 
                          ||
                          (flat_fielded && (flux - local_median) > THRESHOLD * dflux)
                          )
                    ) {
#if WANT_BIG_LOGFILE
            uves_msg_debug("threshold = %f\tx = %d\tflux = %f\tmedian = %f", 
                   THRESHOLD, x, flux, local_median);
#endif
            
            x += 1;
            peak_width += 1;
            
            if (x <= nx) {
            /* flux =
               cpl_image_get(spectrum, x,
               spectrum_row, &pis_rejected);
            */
            flux = spectrum_data[(x-1) + (spectrum_row - 1) * nx];
            xlo = uves_max_int(x - RANGE, 1);
            xhi = uves_min_int(x + RANGE, nx);
            local_median = cpl_image_get_median_window(
                spectrum,
                uves_max_int(xlo, 1 ), spectrum_row,
                uves_min_int(xhi, nx), spectrum_row);
            }
        }
        /* x is now first position that is below (median + threshold) */
        
        if (peak_width > 0) {
            double x_peak, dx = 0, sigma, slope, back;
            check( x_peak = xcenter(spectrum, noise,
                        uves_max_int(1, x - peak_width), 
                        /* First position above threshold */ 
                        uves_max_int(1, x - 1),          
                        /* Last  position above threshold */ 
                        spectrum_row,
                        CENTERING_METHOD,
                                            bin_disp,
                        &sigma,
                        &flux,
                        &dx,
                                            &slope,
                                            &back),
               "Could not locate peak center");
            
#if WANT_BIG_LOGFILE
            uves_msg_debug("(Order, x, flux) = (%d, %f, %f)", 
                   order, x_peak, flux);
#endif        
            /* Add line to line table, but only if less
               lines that MAXLINES have been detected */
            if (*nrows < MAXLINES) {
            check(( cpl_table_set_int   (linetable, "Y"     , *nrows, order),
                cpl_table_set_double(linetable, "X"     , *nrows, x_peak),
                cpl_table_set_double(linetable, "dX"    , *nrows, dx),
                cpl_table_set_double(linetable, "Xwidth", *nrows, sigma),
                cpl_table_set_double(linetable, "Peak"  , *nrows, flux),
                cpl_table_set_double(linetable, "Background" , *nrows, back),
                cpl_table_set_double(linetable, "Slope" , *nrows, slope)),
                  "Could not update line table row %d", *nrows);
            (*nrows)++;
            }
            
            ndetected_order++;
            (*ndetected)++;
            
            if (arcframe != NULL) {
            int x1;
            int pen = 0;  /* Value to write */
            int ny = cpl_image_get_size_y(arcframe);
            /* We already know 'nx' from the width of the spectrum image */
            
            for (x1  = uves_max_int(
                 1 , uves_round_double(
                     x_peak - peak_width - 0*RANGE/2.0)); 
                 x1 <= uves_min_int(
                 nx, uves_round_double(
                     x_peak + peak_width + 0*RANGE/2.0)); 
                 x1++) {
                check( cpl_image_set(
                       arcframe,
                       x1,
                       uves_min_int(
                       ny, 
                       uves_max_int(
                           1, 
                           (int) uves_polynomial_evaluate_2d(
                           order_locations, x1, order)     
                           )), 
                       pen),
                   "Error writing input image");
                check( cpl_image_set(
                       arcframe,
                       uves_min_int(
                       nx,
                       uves_max_int((int) x_peak, 1)),
                       uves_min_int(
                       ny, 
                       uves_max_int(
                           1,
                           (int) uves_polynomial_evaluate_2d(
                           order_locations, x1, order)
                           - 10)), 
                       pen),
                   "Error writing input image");
            }
            }
        } /* line found */
        }/* for x */
        if (arcframe != NULL) uves_msg_debug("Order #%d: %d lines detected", 
                         order, ndetected_order);
    }/* for order */
    
    /* Remove doublets */
    {
    int i;
    int doublets_removed = 0;
    for (i = 0; i+1 < *nrows; i++) {
        if (fabs(cpl_table_get_double(linetable, "X", i  , NULL) - 
             cpl_table_get_double(linetable, "X", i+1, NULL))  < 2.0) 
        {
            /* If a doublet was found, delete it. 
               Make sure the table stays the same size
               by adding two rows at the end. */

            check( cpl_table_erase_window(linetable, i, 2),
               "Error removing rows");
            *nrows -= 2;
            *ndetected -= 2;
            
            check( cpl_table_set_size(linetable, 
                          cpl_table_get_nrow(linetable) + 2),
               "Could not resize line table");
            
            doublets_removed++;
        }
    }
    if (doublets_removed > 0)
        {
        uves_msg_debug("%d doublet%s removed", 
                   doublets_removed, doublets_removed > 1 ? "s" : "");
        }
    }
    
            uves_msg("Range = %d pixels; threshold = %.2f %s; %d lines detected", 
         RANGE, THRESHOLD, flat_fielded ? "stdev" : "ADU", *ndetected);
    
  cleanup:    
    return cpl_error_get_code();
}

/*----------------------------------------------------------------------------*/
/**
  @brief    Refine the center position of an initially detected emission line
  @param    image             The input image
  @param    noise             The noise (1 sigma) associated with the input image
  @param    xlo               Left x-coordinate of search window
  @param    xhi               Right x-coordinate of search window
  @param    row               Image row to look at
  @param    CENTERING_METHOD  The method used to calculate the position
  @param    bin_disp          CCD binning, used to finetune size of fitting window
  @param    sigma             (output) Peak width
  @param    intensity         (output) Peak height (relative to 0)
  @param    dx0               (output) The statistical uncertainty of peak center
  @param    slope             (output) fitted background slope
  @param    background        (output) fitted background
  @return   The refined center location, or undefined on error

  This function locates the center of an emission peak in the image window
  (@em xlo, @em y)-(@em xhi, @em y) by calculating guess values for the center, 
  the stddev, the background level and the height of the peak, then calling @c
  uves_fit_gaussian_1d_image() with the "guess" values.

  The @em centering_method indicates whether to use the 'gravity' method or gaussian
  fitting and should always be set to @em CENTERING_GAUSSIAN .
**/
/*----------------------------------------------------------------------------*/
static double
xcenter(const cpl_image *image, const cpl_image *noise, int xlo, int xhi, int row,
    centering_method CENTERING_METHOD, int bin_disp,
    double *sigma, double *intensity, double *dx0, double *slope, double *background)
{
    double x0;           /* Result */
    cpl_matrix *covariance = NULL;
    const double *image_data;
    bool converged;
    int lo_r, hi_r;
    
    int nx = cpl_image_get_size_x(image);

    passure(cpl_image_get_type(image) == CPL_TYPE_DOUBLE, " ");

    image_data = cpl_image_get_data_double_const(image);

    /* Make sure fit window is 13-19 pixels
       (7-9 pixels for binned CCD) */
    lo_r = 6;
    hi_r = 8;
    if (bin_disp >= 2) 
        {
            lo_r = 4;
            hi_r = 5;
        }

    {
        int xm = (xlo+xhi)/2;

        xlo = uves_max_int(1, xm - lo_r);
        xhi = uves_min_int(nx, xm + lo_r);
    }

    /* Increase fit window (up to 19 pixels) */
    do {
        converged = true;
        if (1 < xlo && 0 <
            //cpl_image_get(image, xlo - 1, row, &pis_rejected) &&
            //cpl_image_get(image, xlo - 1, row, &pis_rejected) <
            //cpl_image_get(image, xlo    , row, &pis_rejected) )
            image_data[(xlo-1-1) + (row - 1) * nx] &&
            image_data[(xlo-1-1) + (row - 1) * nx] <
            image_data[(xlo  -1) + (row - 1) * nx] )
            {
                converged = false;
                xlo -= 1;
            }
        
        if (xhi < nx && 0 <
            //cpl_image_get(image, xhi + 1, row, &pis_rejected) && 
            //cpl_image_get(image, xhi + 1, row, &pis_rejected) <
            //cpl_image_get(image, xhi    , row, &pis_rejected) )
            image_data[(xhi+1-1) + (row - 1) * nx] &&
            image_data[(xhi+1-1) + (row - 1) * nx] <
            image_data[(xhi  -1) + (row - 1) * nx] )
            {
                converged = false;
                xhi += 1;
            }

        if ((xhi-xlo+1) >= hi_r) 
            {
                converged = true;
            }

    } while (!converged);

    /* Get precise location */
    if (CENTERING_METHOD == CENTERING_GAUSSIAN)
    {
#if WEIGHTED_FIT
        uves_fit_1d_image(image, noise, NULL,
#else /* Unweighted fit like MIDAS which gives larger sigma */
        uves_fit_1d_image(image, NULL, NULL,
#endif
                  true, false, false,
                  xlo, xhi, row,
                  &x0, sigma, intensity, background, slope,
#if WEIGHTED_FIT
                  NULL, NULL, &covariance,
#else
                  NULL, NULL, NULL,
#endif

#if FIT_SLOPE
                              uves_gauss_linear, uves_gauss_linear_derivative, 5);
#else
                              uves_gauss, uves_gauss_derivative, 4);
            *slope = 0;
#endif

        /* The fitting routine sometimes (i.e. regularly) fails
         * because of low statistics.
         * Recover from specific fitting errors */
        if (cpl_error_get_code() == CPL_ERROR_NONE)
        {
            /* Variance is guaranteed to be positive */
#if WEIGHTED_FIT
            *dx0 = sqrt(cpl_matrix_get(covariance, 0, 0));
#else
            *dx0 = *sigma / sqrt(*intensity);
#endif

#if WANT_BIG_LOGFILE
             uves_msg_debug("Gaussian fit succeeded at (x, row, N) = (%f, %d, %d)",
                    x0, row, xhi-xlo+1);
#endif
        }
        else if (cpl_error_get_code() == CPL_ERROR_CONTINUE)
        {
            /* Fitting failed */
            uves_error_reset();
#if WANT_BIG_LOGFILE
            uves_msg_debug("Gaussian fit failed at (x, row, N) ="
                   " (%f, %d, %d), using centroid", 
                   x0, row, xhi-xlo+1);
#endif
            *dx0 = *sigma / sqrt(*intensity);
        }
        else if (cpl_error_get_code() == CPL_ERROR_SINGULAR_MATRIX)
        {
            uves_error_reset();
            
            /* Fitting succeeded but covariance computation failed */
            uves_msg_debug("Covariance matrix computation failed");
            *dx0 = *sigma / sqrt(*intensity);
        }
        
        assure(cpl_error_get_code() == CPL_ERROR_NONE, cpl_error_get_code(),
           "Gaussian fitting failed");

#if WANT_BIG_LOGFILE
        uves_msg_debug("Fit   = (x0=%f, sigma=%f, norm=%f, backg=%f, N=%d)", 
               x0,
               *sigma,
               *intensity,
               background,
               xhi - xlo + 1);
#endif

        /* 'intensity' is the norm (area) of the gaussian fit.
           But we need to return the height above zero
           (for MIDAS compatibility). 
           height = f(x0) = background + norm/sqrt(2pi sigma^2)
        */

        *intensity = *background + (*intensity)/(sqrt(2*M_PI) * (*sigma));

    }
    else   /*  if (CENTERING_METHOD == CENTERING_GRAVITY) */
    {
        assure (false, CPL_ERROR_UNSUPPORTED_MODE,
            "Centering method (no. %d) is unsupported", 
            CENTERING_METHOD);
    }

  cleanup:
    uves_free_matrix(&covariance);
    return x0;
}

/**@}*/