/*                                                                              *
 *   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: 2012-03-02 16:23:31 $
 * $Revision: 1.53 $
 * $Name: not supported by cvs2svn $
 *
 */

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

/*----------------------------------------------------------------------------*/
/**
 * @defgroup uves_backsub   Substep: Background Subtraction
 *
 *  Subtract the background from an echelle image using method spline (recommended),
 *  polynomial or smooth.
 */
/*----------------------------------------------------------------------------*/


#include <uves_backsub.h>

#include <uves_parameters.h>
#include <uves_pfits.h>
#include <uves_dump.h>
#include <uves_utils.h>
#include <uves_utils_wrappers.h>
#include <uves_utils_cpl.h>
#include <uves_error.h>
#include <uves_msg.h>
#include <uves.h>

#include <cpl.h>
#include <string.h>
#include <stdbool.h>
#include <float.h>
/*-----------------------------------------------------------------------------
                            Functions prototypes
 -----------------------------------------------------------------------------*/
static int first_order(const polynomial *order_locations, int nx);
static int last_order (const polynomial *order_locations, int nx, int ny);
static cpl_error_code lower_to_average(cpl_image *image, int RADX, int RADY);
static double sample_background(const cpl_image *image, int x0, double y_0,
                int radius_x, int radius_y, int nx, int ny,
                background_measure_method BM_METHOD);
static cpl_error_code subtract_background(cpl_image *image, cpl_image *background_im, 
                      const polynomial *background_pol);

/*-----------------------------------------------------------------------------
                                Defines
 -----------------------------------------------------------------------------*/

/* This is sort of ugly, because we fine tune parameters depending on
   wavelength and also different for masterflat/science exposures.
   A 'perfect' background subtraction algorithm should not need to
   know about its context.
*/ 

#define BACKSUB_FLAT_SMOOTHX_BLUE (25.0/4096)
#define BACKSUB_FLAT_SMOOTHX_RED  (50.0/4096)
#define BACKSUB_FLAT_SMOOTHY_BLUE (100.0/2048)
#define BACKSUB_FLAT_SMOOTHY_RED  (300.0/2048)

#define BACKSUB_SCI_SMOOTHX_BLUE  (300.0/4096)
#define BACKSUB_SCI_SMOOTHX_RED   (300.0/4096)
#define BACKSUB_SCI_SMOOTHY_BLUE  (200.0/2048)
#define BACKSUB_SCI_SMOOTHY_RED   (500.0/2048)

#define BACKSUB_SMOOTHY_WLEN 859.9

/**@{*/

/*-----------------------------------------------------------------------------
                            Implementation
 -----------------------------------------------------------------------------*/

/*----------------------------------------------------------------------------*/
/**
   @brief    Define parameters used for spline background subtraction
   @return   Newly allocated parameter list for this substep.

   The parameters defined are mmethod, npoints, radiusy, sdegree, smoothx, smoothy.
   See source code or 'esorex --man-page' for a description of each parameter.
*/
/*----------------------------------------------------------------------------*/

cpl_parameterlist *
uves_backsub_define_parameters(void)
{

    char *full_name = NULL;
    cpl_parameterlist *parameters = NULL;
    cpl_parameter *p = NULL;

    parameters = cpl_parameterlist_new();
    
    //
    const char* name = "mmethod";
    full_name = uves_sprintf("%s.%s", UVES_BACKSUB_ID, name);
    
    uves_parameter_new_enum(p, full_name,
                   CPL_TYPE_STRING,
                   "Background measuring method. If equal to 'median' "
                   "the background is sampled using the median of a subwindow. "
                   "If 'minimum', the subwindow minimum value is used. "
                   "If 'no', no background subtraction is done.",
                   UVES_BACKSUB_ID,
                   "median",                        /* Default */
                   3,                               /* Number of options */
                   "median", "minimum", "no");      /* List of options */
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
    cpl_parameterlist_append(parameters, p);
    cpl_free(full_name);

    //
    name = "npoints";
    full_name = uves_sprintf("%s.%s", UVES_BACKSUB_ID, name);
    uves_parameter_new_range(p, full_name,
                 CPL_TYPE_INT,
                 "This is the number of columns in interorder space "
                 "used to sample the background.",
                 UVES_BACKSUB_ID,
                 82, 0, INT_MAX);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
    cpl_parameterlist_append(parameters, p);
    cpl_free(full_name);
    
    //
    name = "radiusy";
    full_name = uves_sprintf("%s.%s", UVES_BACKSUB_ID, name);
    uves_parameter_new_range(p, full_name,
                CPL_TYPE_INT,
                "The height (in pixels) of the background sampling "
                "window is (2*radiusy + 1). "
                "This parameter is not corrected for binning.",
                UVES_BACKSUB_ID,
                2, 0, INT_MAX);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
    cpl_parameterlist_append(parameters, p);
    cpl_free(full_name);
    
    //
    name = "sdegree";
    full_name = uves_sprintf("%s.%s", UVES_BACKSUB_ID, name);
    uves_parameter_new_range(p, full_name,
                 CPL_TYPE_INT,
                 "Degree of interpolating splines. Currently "
                 "only degree = 1 is supported",
                 UVES_BACKSUB_ID,
                 1, 0, INT_MAX);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
    cpl_parameterlist_append(parameters, p);
    cpl_free(full_name);

    //
    name = "smoothx";
    full_name = uves_sprintf("%s.%s", UVES_BACKSUB_ID, name);
    uves_parameter_new_range(p, full_name,
                 CPL_TYPE_DOUBLE,
                 "If spline interpolation is used to measure the background, "
                 "the x-radius of the post-smoothing window is "
                 "(smoothx * image_width). Here, 'image_width' is the image "
                 "width after binning. If negative, the default values are used: "
                 make_str(BACKSUB_FLAT_SMOOTHX_BLUE) " for blue flat-field frames, "
                 make_str(BACKSUB_FLAT_SMOOTHX_RED) " for red flat-field frames, "
                 make_str(BACKSUB_SCI_SMOOTHX_BLUE) " for blue science frames and "
                 make_str(BACKSUB_SCI_SMOOTHX_RED) " for red science frames.",
                 UVES_BACKSUB_ID,
                 -1.0, -DBL_MAX, DBL_MAX);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
    cpl_parameterlist_append(parameters, p);
    cpl_free(full_name);
    
    //
    name = "smoothy";
    full_name = uves_sprintf("%s.%s", UVES_BACKSUB_ID, name);
    uves_parameter_new_range(p, full_name,
                 CPL_TYPE_DOUBLE,
                 "If spline interpolation is used to measure the "
                 "background, the y-radius of the post-smoothing "
                 "window is (smoothy * image_height). Here, "
                 "'image_height' is the image height after binning. "
                 "If negative, the default values are used: "
                 make_str(BACKSUB_FLAT_SMOOTHY_BLUE) " for blue flat-field frames, "
                 make_str(BACKSUB_FLAT_SMOOTHY_RED) " for red flat-field frames, "
                 make_str(BACKSUB_SCI_SMOOTHY_BLUE) " for blue science frames and "
                 make_str(BACKSUB_SCI_SMOOTHY_RED) " for red science frames.",
                 UVES_BACKSUB_ID,
                 -1.0, -DBL_MAX, DBL_MAX);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
    cpl_parameterlist_append(parameters, p);
    cpl_free(full_name);

    if (cpl_error_get_code() != CPL_ERROR_NONE)
    {
        cpl_msg_error(__func__, "Creation of spline background subtraction "
              "parameters failed: '%s'", cpl_error_get_where());
        cpl_parameterlist_delete(parameters);
        return NULL;
    }
    else
    {
        return parameters;
    }
}

/*----------------------------------------------------------------------------*/
/**
  @brief    Read background measuring method from parameter list
  @param    parameters             The parameter list
  @param    context                Context of parameter (or NULL)
  @param    subcontext             Subcontext of parameter
                                           
  @return   The measure method as read from the parameter context.subcontext.mmethod
  
**/
/*----------------------------------------------------------------------------*/
background_measure_method
uves_get_bm_method(const cpl_parameterlist *parameters, const char *context, 
           const char *subcontext)
{
    const char *bm = "";
    background_measure_method result = 0;

    check( uves_get_parameter(parameters, context, subcontext, "mmethod", CPL_TYPE_STRING, &bm),
       "Could not read parameter");
    
    if      (strcmp(bm, "median" ) == 0) result = BM_MEDIAN;
    else if (strcmp(bm, "minimum") == 0) result = BM_MINIMUM;
    else if (strcmp(bm, "no"     ) == 0) result = BM_NO;
    else
    {
        /* Impossible */ assure(false, CPL_ERROR_ILLEGAL_INPUT, 
                    "No such background measuring method: '%s'", bm);
    }
    
  cleanup:
    return result;
}
    
/*----------------------------------------------------------------------------*/
/**
   @brief   Subtract background of an echelle image
   @param   image          Subtract background of this image
   @param   raw_header     The raw image header.
   @param   ordertable     The order table, used to get minimum and maximum order number
   @param   order_locations   Polynomial describing the order locations.
   @param   parameters      The parameters used for background subtraction.
                Add parameters by calling @c uves_propagate_parameters_step() from
                the @c recipe_create() function
   @param   context         Use @em parameters belonging to this context
   @param   chip            CCD chip
   @param   flat_field      Use back-sub default parameters best suited 
                            for flat-field (true) or science frames (false).
   @param   background (out) The background image which was subtracted
   @return  CPL_ERROR_NONE iff OK.

   This function estimates and subtracts the background of an echelle image
   using spline interpolation. The input @em image is sampled at half-integer order locations
   using the x step size (image_width / @em NPOINTS). The median or minimum of a
   subwindow of size (2 * @em radius_x + 1) * (2 * @em radius_y + 1) is used as the background
   sample value.

   After the spline interpolation, the background is filtered using an average filter
   with radius (@em smooth_x, @em smooth_y).

   After subtraction, the input image is thresholded to be everywhere non-negative

   Currently, only splines of degree 1 and smooth factor 0
   is supported (i.e. linear interpolation).

   If the mmethod parameter is 'no', no background subtraction is done.

*/
/*----------------------------------------------------------------------------*/

cpl_error_code
uves_backsub_spline(cpl_image *image, const uves_propertylist *raw_header,
            const cpl_table *ordertable, const polynomial *order_locations,
            const cpl_parameterlist *parameters, const char *context, 
            enum uves_chip chip,
            bool flat_field,
            cpl_image **background)
{
    /* Recipe parameters */
    background_measure_method BM_METHOD;
    int npoints;
    int radius_y;
    int bin_x=1;
    int bin_y=1;

    int sdegree;
    double SMOOTHX;
    double SMOOTHY;
    
    /* Local variables */
    int nx, ny;
    int x, y;
    int stepx;
    int radius_x;
    int smooth_x, smooth_y;        /* Window radius in pixels */
    
    passure( image != NULL, " ");
    passure( raw_header != NULL, " ");
    passure( ordertable != NULL, " ");
    passure( order_locations != NULL, " ");
    passure( parameters != NULL, " ");
    passure( context != NULL, " ");
    passure( uves_polynomial_get_dimension(order_locations) == 2, 
         "%d", uves_polynomial_get_dimension(order_locations));
    passure( background != NULL, " ");

    /* Get recipe parameters */
    check( BM_METHOD = uves_get_bm_method(parameters, context, UVES_BACKSUB_ID),
       "Error getting background measuring method");
    
    check( uves_get_parameter(parameters, context, UVES_BACKSUB_ID,
                  "npoints", CPL_TYPE_INT   , &npoints) , "Could not read parameter");
    check( uves_get_parameter(parameters, context, UVES_BACKSUB_ID,
                  "radiusy", CPL_TYPE_INT   , &radius_y), "Could not read parameter");

    check(bin_x=uves_pfits_get_binx(raw_header),"error getting %s",UVES_BINX);
    check(bin_y=uves_pfits_get_biny(raw_header),"error getting %s",UVES_BINY);

    radius_y = uves_round_double((double)radius_y/bin_y);
 
    check( uves_get_parameter(parameters, context, UVES_BACKSUB_ID,
                  "sdegree", CPL_TYPE_INT   , &sdegree) , "Could not read parameter");
    check( uves_get_parameter(parameters, context, UVES_BACKSUB_ID,
                  "smoothx", CPL_TYPE_DOUBLE, &SMOOTHX) , "Could not read parameter");
    check( uves_get_parameter(parameters, context, UVES_BACKSUB_ID,
                  "smoothy", CPL_TYPE_DOUBLE, &SMOOTHY) , "Could not read parameter");

   
    /* Get other parameters */
    nx = cpl_image_get_size_x(image);
    ny = cpl_image_get_size_y(image);


    if (BM_METHOD == BM_NO)
    {
        uves_msg("Skipping background subtraction");

        /* Calculate a zero-background */
        check( *background = cpl_image_new(nx, ny, CPL_TYPE_DOUBLE),
           "Error allocating image");
    }
    else {
    /* If negative, set default values for smoothx, smoothy */
    if (SMOOTHX < 0)
        {
        if (chip == UVES_CHIP_BLUE)
            {
            SMOOTHX = (flat_field) ? 
                BACKSUB_FLAT_SMOOTHX_BLUE : BACKSUB_SCI_SMOOTHX_BLUE;
            }
        else
            {
            SMOOTHX = (flat_field) ? 
                BACKSUB_FLAT_SMOOTHX_RED : BACKSUB_SCI_SMOOTHX_RED;
            }
        }
    if (SMOOTHY < 0)
        {
        double wlen;
        
        /* Read wavelength from raw header */
        
        check( wlen = uves_pfits_get_gratwlen(raw_header, chip),
               "Error reading central wavelength");
        
        /* The criterion is not if the chip is BLUE/RED,
           but whether the wlen is < 860A */
        if (wlen < BACKSUB_SMOOTHY_WLEN)
            {
            SMOOTHY = (flat_field) ? 
                BACKSUB_FLAT_SMOOTHY_BLUE : BACKSUB_SCI_SMOOTHY_BLUE;
            }
        else
            {
            SMOOTHY = (flat_field) ? 
                BACKSUB_FLAT_SMOOTHY_RED : BACKSUB_SCI_SMOOTHY_RED;
            }
        }
    
    assure( 0 < SMOOTHX, CPL_ERROR_ILLEGAL_INPUT, "Illegal smoothx factor: %e", SMOOTHX);
    assure( 0 < SMOOTHY, CPL_ERROR_ILLEGAL_INPUT, "Illegal smoothy factor: %e", SMOOTHY);
    
    smooth_x = uves_round_double(SMOOTHX * nx - 0.5);
    smooth_y = uves_round_double(SMOOTHY * ny - 0.5);
    
    assure( 0 < npoints, CPL_ERROR_ILLEGAL_INPUT, 
        "Illegal number of sample points: %d", npoints);
    stepx = nx / npoints;
    assure( 0 < stepx, CPL_ERROR_ILLEGAL_INPUT, "Illegal step size: %d", stepx);
    radius_x = stepx/2;
    assure( 0 < radius_x, CPL_ERROR_ILLEGAL_INPUT, "Illegal x sample radius: %d", radius_x);
    assure( 0 < radius_y, CPL_ERROR_ILLEGAL_INPUT, "Illegal y sample radius: %d", radius_y);
    assure( 0 < smooth_x, CPL_ERROR_ILLEGAL_INPUT, "Illegal x sample smooth: %d", smooth_x);
    assure( 0 < smooth_y, CPL_ERROR_ILLEGAL_INPUT, "Illegal y sample smooth: %d", smooth_y);
    assure( sdegree == 1, CPL_ERROR_UNSUPPORTED_MODE, 
        "Spline degree must be 1. It is %d", sdegree);
    
    uves_msg("Sample window (pixels): radx, rady = %d, %d", radius_x, radius_y);
    
    check( *background = cpl_image_new(nx, ny, CPL_TYPE_DOUBLE), 
           "Error allocating background image");
    
    /* Process */
    
    for (x = stepx; x <= nx; x += stepx) {
        int order, minorder, maxorder;
        /* Find min. and max. order where background positions are inside image  */
            
        minorder = cpl_table_get_column_min(ordertable, "Order");
            
        /* If outside image, move to inside image */
        while (uves_round_double(
               uves_polynomial_evaluate_2d(order_locations, x + radius_x, minorder - 0.5)
               ) - radius_y < 1 ||
           uves_round_double(
               uves_polynomial_evaluate_2d(order_locations, x - radius_x, minorder - 0.5))
           - radius_y < 1  )
        {
            int sign;

            for (sign = -1; sign <= 1; sign += 2)
            {
                assure( 
                uves_polynomial_evaluate_2d(order_locations,
                                x + sign*radius_x, minorder+1 - 0.5) >
                uves_polynomial_evaluate_2d(order_locations,
                                x + sign*radius_x, minorder   - 0.5),
                CPL_ERROR_ILLEGAL_INPUT,
                "Order polynomial is not well-formed: "
                "p(%d, %f) = %e; p(%d, %f) = %e",
                x + sign*radius_x, minorder+1 - 0.5, uves_polynomial_evaluate_2d(
                    order_locations, x + sign*radius_x, minorder+1 - 0.5
                    ),
                x + sign*radius_x, minorder   - 0.5, uves_polynomial_evaluate_2d(
                    order_locations, x + sign*radius_x, minorder   - 0.5)
                );
            }

            minorder += 1;
        }
            
        maxorder = cpl_table_get_column_max(ordertable, "Order");
        
        /* If outside image, move to inside image */
        while (uves_round_double( 
               uves_polynomial_evaluate_2d(order_locations, x + radius_x, maxorder + 0.5)
               ) + radius_y > ny ||
           uves_round_double( 
               uves_polynomial_evaluate_2d(order_locations, x - radius_x, maxorder + 0.5)
               ) + radius_y > ny  ) {
        int sign;
        for (sign = -1; sign <= 1; sign += 2)
            {
            assure( 
                uves_polynomial_evaluate_2d(
                order_locations, x + sign*radius_x, maxorder-1 - 0.5) <
                uves_polynomial_evaluate_2d(order_locations, 
                            x + sign*radius_x, maxorder   - 0.5), 
                CPL_ERROR_ILLEGAL_INPUT,
                "Order polynomial is not well-formed: "
                "p(%d, %f) = %e; p(%d, %f) = %e",
                x + sign*radius_x, maxorder-1 - 0.5, uves_polynomial_evaluate_2d(
                order_locations, x + sign*radius_x, maxorder-1 - 0.5),
                x + sign*radius_x, maxorder   - 0.5, uves_polynomial_evaluate_2d(
                order_locations, x + sign*radius_x, maxorder   - 0.5)
                );
            }
                
        maxorder -= 1;
        }
        
            /* Move to min. order inside image */
            while (uves_round_double(uves_polynomial_evaluate_2d(
                                         order_locations, x + radius_x, minorder - 1.5)
                       ) - radius_y >= 1 &&
                   uves_round_double(uves_polynomial_evaluate_2d(
                                         order_locations, x - radius_x, minorder - 1.5)
                       ) - radius_y >= 1  )
                {
                    int sign;
                    for (sign = -1; sign <= 1; sign += 2)
                        {
                            assure( 
                                uves_polynomial_evaluate_2d(
                                    order_locations, x + sign*radius_x, minorder-1 - 1.5) <
                                uves_polynomial_evaluate_2d(
                                    order_locations, x + sign*radius_x, minorder   - 1.5), 
                                CPL_ERROR_ILLEGAL_INPUT,
                                "Order polynomial is not well-formed: "
                                "p(%d, %f) = %e ; p(%d, %f) = %e",
                                x + sign*radius_x, minorder-1 - 1.5, 
                                uves_polynomial_evaluate_2d(
                                    order_locations, x + sign*radius_x, minorder-1 - 1.5),
                                x + sign*radius_x, minorder   - 1.5,
                                uves_polynomial_evaluate_2d(
                                    order_locations, x + sign*radius_x, minorder   - 1.5));
                        }
                    
                    minorder -= 1;
                }
            
            /* Move to max. order inside image */
            while (uves_round_double( uves_polynomial_evaluate_2d(
                                          order_locations, x + radius_x, maxorder + 1.5)
                       ) + radius_y <= ny &&
                   uves_round_double( uves_polynomial_evaluate_2d(
                                          order_locations, x - radius_x, maxorder + 1.5)
                       ) + radius_y <= ny  ) {
                int sign;
                for (sign = -1; sign <= 1; sign += 2)
                    {
                        assure( 
                            uves_polynomial_evaluate_2d(
                                order_locations, x + sign*radius_x, maxorder+1 + 1.5)
                            >
                            uves_polynomial_evaluate_2d(
                                order_locations, x + sign*radius_x, maxorder   + 1.5),
                            CPL_ERROR_ILLEGAL_INPUT,
                            "Order polynomial is not well-formed: "
                            "p(%d, %f) = %e ; p(%d, %f) = %e",
                            x + sign*radius_x, maxorder+1 + 1.5,
                            uves_polynomial_evaluate_2d(
                                order_locations, x + sign*radius_x, maxorder+1 + 1.5),
                            x + sign*radius_x, maxorder   + 1.5,
                            uves_polynomial_evaluate_2d(
                                order_locations, x + sign*radius_x, maxorder   + 1.5));
                    }
                
                maxorder += 1;
            }
        
        uves_msg_debug("(x, order) = (%d, %f - %f)  ", x, minorder-.5, maxorder+.5);
        
        for (order = minorder; order <= maxorder; order++) {
        int ylo, yhi;
        double backlo, backhi;
            
        /* Sample background above and below order using the median of a window
         * with size (2*radius_x + 1) * (2*radius_y + 1)
         */
            
        ylo = uves_round_double( 
            uves_polynomial_evaluate_2d(order_locations, x, order - 0.5) );
        yhi = uves_round_double(
            uves_polynomial_evaluate_2d(order_locations, x, order + 0.5) );

        /* Fail cleanly if input polynomial is corrupted */
        assure( yhi > ylo, CPL_ERROR_ILLEGAL_INPUT,
            "Order polynomial is not well-formed: "
            "p(%d, %f) = %d ; p(%d, %f) = %d",
            x, order - 0.5, ylo,
            x, order + 0.5, yhi);
            
            
        check( backlo = 
               sample_background(
               image, x, ylo, radius_x, radius_y, nx, ny, BM_METHOD),
               "Error sampling background level");
            
        check( backhi = sample_background(
               image, x, yhi, radius_x, radius_y, nx, ny, BM_METHOD),
               "Error sampling background level");
                   
        uves_msg_debug("Background sample at (x, y, order) = (%d, %d, %f) = %f",
                   x, ylo, order-0.5, backlo);
        uves_msg_debug("Background sample at (x, y, order) = (%d, %d, %f) = %f",
                   x, yhi, order+0.5, backhi);
    
        /* Extrapolate (linearly, or constant if MIDAS) if first order */
        if (order == minorder) {
            for (y = 1; y <= ylo; y++) {
            double back = backlo + (backhi - backlo)*(y - ylo)/(yhi - ylo);
            cpl_image_set(*background, x, y, back);

                        cpl_image_set(*background, x, y, back);
            }
        }
            
        /* Make a linear interpolation (1-degree, no-smooth spline) from ylo to yhi */
        for (y = ylo; y <= yhi; y++) {
            double back;
            back = backlo + (backhi - backlo) * (y - ylo) / (yhi - ylo);
            /* We know that yhi > ylo */
            cpl_image_set(*background, x, y, back);
        }
            
        /* Extrapolate (linearly, or constant if MIDAS) if last order */
        if (order == maxorder) {
            for (y = yhi; y <= ny; y++) {
            double back;
            back = backlo + (backhi - backlo) * (y - ylo) / (yhi - ylo);

                        cpl_image_set(*background, x, y, back);
                    }
        }
        }
    }/* For column...  */

    /* Now interpolate between columns */
    for (y = 1; y <= ny; y++) {
        int col;
        for (col = stepx; col+stepx <= nx; col += stepx) {
        int pis_rejected; /* Not used, all pixels read are good; they've just been set */
            
        double backlo, backhi;
            
        /* Read this and next column */
        backlo = cpl_image_get(*background, col      , y, &pis_rejected);
        backhi = cpl_image_get(*background, col+stepx, y, &pis_rejected);
        
        /* Extrapolate (linear) before first column */
        if (col == stepx)
            for (x = 1; x <= col; x++)
            {
                double back = backlo + (backhi - backlo) * (x - col) / stepx;
                cpl_image_set(*background, x, y, back);
            }
            
        /* Interpolate between columns */
        for (x = col; x <= col + stepx; x++)
            {
            double back = backlo + (backhi - backlo) * (x - col) / stepx;
            cpl_image_set(*background, x, y, back);
            }

        /* Extrapolate (linear) after last column */
        if (col+stepx+stepx > nx)
            for (x = col; x <= nx; x++)
            {
                double back = backlo + (backhi - backlo) * (x - col) / stepx;
                cpl_image_set(*background, x, y, back);
            }
        }
    }

    /* All pixels in background image have been set.
     * Smooth background. 
     */

    uves_msg("Smoothing window (pixels): smox, smoy = %d, %d", smooth_x, smooth_y);
    check( uves_filter_image_average(*background, smooth_x, smooth_y), 
           "Error applying average filter to background image");

    uves_msg("Subtracting background image");

    check( subtract_background(image, *background, NULL),
           "Error subtracting background image");


   } /* BM_METHOD was not 'no' */
 

  cleanup:
    return cpl_error_get_code();
}

/*----------------------------------------------------------------------------*/
/**
   @brief   Subtract background of an echelle image
   @param   image   Subtract background of this image
   @param   orders  Table describing the orders locations
                    using the columns @em 'Slope' and @em 'Intersept'. Table
            must be sorted according to @em 'Intersept' or an error is set.
            May be NULL.
   @param   order_locations   Polynomial describing the order locations. May be NULL.

   @param   NPOINTS For efficiency reasons, the image is sampled on a grid with
                    mask size (@em image_width / NPOINTS , @em image_height / NPOINTS)
   @param   radius_y Vertical sample radius (in pixels)

   @param   BM_METHOD Background measure method. If equal to BM_NO, no subtraction
                      is done. Also see @c sample_background().

   @param   DEGX    Degree (in x) of polynomial
   @param   DEGY    Degree (in y) of polynomial
   @param   KAPPA   Value of kappa used in kappa-sigma clipping. If negative, no
                    clipping is done.

   @return  CPL_ERROR_NONE iff OK.

   This function estimates and subtracts the background of an echelle image by fitting
   a low degree polynomial to a subset of the image pixels and continuously rejecting
   points (such as the signal itself) that have large positive residuals (one-sided
   kappa-sigma clipping).

   At least one of the parameters @em orders or @em order_locations must be NULL. 

   - If both are NULL, the input @em image is sampled on a regular grid with mask size
   (@em image_width / @em NPOINTS , @em image_height / @em NPOINTS).
   This is used to subtract the background of an echelle image before the order 
   locations are known.

   - If @em orders is non-NULL, the background is sampled 
   (at locations separated by @em image_width / NPOINTS)
   between the order lines defined in the table
   (defined by the columns @em 'Slope' and @em 'Intersept').

   - If the polynomial @em order_locations is non-NULL, the background is sampled (at locations
   separated by @em image_width / NPOINTS) between the order lines defined by the polynomial.

*/
/*----------------------------------------------------------------------------*/
cpl_error_code
uves_backsub_poly(cpl_image *image,
          const cpl_table *orders, const polynomial *order_locations, 
          background_measure_method BM_METHOD,
          int NPOINTS,
          int radius_y,
          int DEGX, 
          int DEGY,
          double KAPPA)
{
    cpl_table  *t          = NULL;
    polynomial *background = NULL;
    int nx, ny;
    int stepx, stepy;                   /* Step size */
    int radius_x;                       /* Sample window x-radius */
    double mse, rmse;                   /* mse, rms of fit */
    cpl_size total_clipped = 0;
    
    if (BM_METHOD == BM_NO)
    {
        uves_msg("Skipping background subtraction");
    }
    else
    {
        passure( image != NULL, " ");
        passure( orders == NULL || order_locations == NULL, " ");
        
        nx = cpl_image_get_size_x(image);
        ny = cpl_image_get_size_y(image);
        
        assure( NPOINTS < nx, CPL_ERROR_ILLEGAL_INPUT,
            "Number of sample columns (%d) larger than image width (%d pixels)", 
            NPOINTS, nx);
        
        stepx = nx/NPOINTS;
        stepy = ny/NPOINTS;

        radius_x = stepx/2;
    
        /* First sample background */
        if (orders != NULL)
        {
            /* Using the order table */

            int x, ordersrow, row;
        
            /* Check input */
            passure( cpl_table_has_column(orders, "Slope"), " ");
            passure( cpl_table_has_column(orders, "Intersept"), " ");

            passure( cpl_table_get_column_type(orders, "Slope") == CPL_TYPE_DOUBLE,
                 "%s", 
                 uves_tostring_cpl_type(cpl_table_get_column_type(orders, "Slope")));
        
            passure( cpl_table_get_column_type(orders, "Intersept") == CPL_TYPE_DOUBLE,
                 "%s",
                 uves_tostring_cpl_type(cpl_table_get_column_type(orders, "Slope")));
        
            /* This check is computationally cheap because 
               there are never very many order lines */
            passure( uves_table_is_sorted_double(orders, "Intersept", false), " ");
        
            /* Need at least two lines to identify inter-order region */
            assure ( cpl_table_get_nrow(orders) >= 2, CPL_ERROR_ILLEGAL_INPUT, 
                 "Only %" CPL_SIZE_FORMAT " line(s) in order table", cpl_table_get_nrow(orders));
        
            t = cpl_table_new( (nx/stepx + 1)*(cpl_table_get_nrow(orders) + 1) );
            cpl_table_new_column(t, "X", CPL_TYPE_INT);
            cpl_table_new_column(t, "Y", CPL_TYPE_INT);
            cpl_table_new_column(t, "Z", CPL_TYPE_DOUBLE);
        
            row = 0;
            for (ordersrow = -1; ordersrow < cpl_table_get_nrow(orders); ordersrow++)
            {
                double slope, intersept;
            
                /* Sample positions between this and the next orderline */
            
                /* Lowest and highest orders are special cases */
                if (ordersrow == -1)
                {
                    slope     = cpl_table_get_double(
                    orders, "Slope"    , 0, NULL);

                    /* Interorder space below lowest order line is at: 
                       intersept0 - (intersept1-intersept0)/2 */
                    intersept =    
                    0.5*cpl_table_get_double(orders, "Intersept", 0, NULL) -
                    0.5*cpl_table_get_double(orders, "Intersept", 1, NULL) ;
                }
                else if (ordersrow == cpl_table_get_nrow(orders) - 1)
                {
                    slope     = cpl_table_get_double(
                    orders, "Slope"    , ordersrow, NULL);
                    
                    /* Interorder space above highest order line is at:
                       intersept(N) + (intersept(N)-intersept(N-1))/2 */
                    intersept =    
                    0.5*cpl_table_get_double(
                        orders, "Intersept", ordersrow, NULL) -
                    0.5*cpl_table_get_double(
                        orders, "Intersept", ordersrow-1, NULL) ;
                }
                else   /* The most common case */
                {
                    slope = 
                    (cpl_table_get_double(
                        orders, "Slope", ordersrow  , NULL) +
                     cpl_table_get_double(
                         orders, "Slope", ordersrow+1, NULL) ) / 2;
                    
                    intersept      = 
                    (cpl_table_get_double(
                        orders, "Intersept", ordersrow  , NULL) +
                     cpl_table_get_double(
                         orders, "Intersept", ordersrow+1, NULL) ) / 2;
                }
            
                /* Sample the interorder space */
                for (x = 1 + stepx/2; x <= nx; x += stepx)
                {
                    int y = uves_round_double(intersept + slope * x);
                
                    if (1 <= y && y <= ny)
                    {
                        double z;
                    
                        check( z = sample_background(
                               image, 
                               x, y,
                               radius_x, radius_y,
                               nx, ny,
                               BM_METHOD),
                           "Error sampling background "
                           "(x, y) = (%d, %d)", x, y);

                        cpl_table_set_int   (t, "X" , row, x);
                        cpl_table_set_int   (t, "Y" , row, y);
                        cpl_table_set_double(t, "Z" , row, z);
                        row++;
                    }
                }
            } /* for ordersrow... */
        
            cpl_table_set_size(t, row);

        }/* if  orders != NULL */
        
        else if (order_locations != NULL)
        {
            /* Sample background using the polynomial */

            int x, minorder, maxorder, order;
            int row;        /* Pointing to row in temporary table */
        
            /* Check input */
            assure( uves_polynomial_get_dimension(order_locations) == 2, 
                CPL_ERROR_ILLEGAL_INPUT,
                "Order location polynomial must be 2d. It is %d!", 
                uves_polynomial_get_dimension(order_locations));
            
            check(( minorder = first_order(order_locations, nx),
                maxorder = last_order(order_locations, nx, ny)),
               "Error getting min. and max. order numbers");

            t = cpl_table_new( (nx/stepx + 1) * (maxorder-minorder+1));
            cpl_table_new_column(t, "X", CPL_TYPE_INT);
            cpl_table_new_column(t, "Y", CPL_TYPE_INT);
            cpl_table_new_column(t, "Z", CPL_TYPE_DOUBLE);
        
            row = 0;
            for (order = minorder; order <= maxorder; order++) {
            /* Sample the interorder space from (minorder+0.5) to (maxorder+0.5) */
            for (x = 1+stepx/2; x <= nx; x += stepx) {
                int y = uves_round_double(
                uves_polynomial_evaluate_2d(order_locations, x, order + 0.5));
                
                if (1 <= y && y <= ny) {
                double z;
                
                check( z = sample_background(image, 
                                 x, y,
                                 radius_x, radius_y,
                                 nx, ny,
                                 BM_METHOD),
                       "Error sampling background (x, order) = (%d, %d+0.5)",
                       x, order);
                
                cpl_table_set_int   (t, "X" , row, x);
                cpl_table_set_int   (t, "Y" , row, y);
                cpl_table_set_double(t, "Z" , row, z);
                row++;
                }
            }
            }
            
            cpl_table_set_size(t, row);
        }
        else
        { 
            /* Grid sampling (order positions unknown) */
            int x, y, row;
        
            t = cpl_table_new((nx/stepx + 1) * (ny/stepy + 1));
            cpl_table_new_column(t, "X" , CPL_TYPE_INT);
            cpl_table_new_column(t, "Y" , CPL_TYPE_INT);
            cpl_table_new_column(t, "Z" , CPL_TYPE_DOUBLE);
        
            row = 0;
            for (y = 1 + stepy/2; y <= ny; y += stepy) 
            {
                for (x = 1+stepx/2; x <= nx; x += stepx) 
                {
                    double z;
                
                    check( z = sample_background(image, 
                                 x, y,
                                 radius_x, radius_y,
                                 nx, ny,
                                 BM_METHOD),
                       "Error sampling background (x, y) = (%d, %d)", x, y);
                
                    cpl_table_set_int   (t, "X" , row, x);
                    cpl_table_set_int   (t, "Y" , row, y);
                    cpl_table_set_double(t, "Z" , row, z);
                    row++;
                }
            }
            cpl_table_set_size(t, row);
        }
        
        /* Sampling done. Fit poly. */

        total_clipped = 0;
        {
        cpl_size n_clipped;
        cpl_size deg_xy=(DEGX + 1)*(DEGY + 1);
        do {
            assure( cpl_table_get_nrow(t) > (DEGX + 1)*(DEGY + 1), 
                CPL_ERROR_ILLEGAL_OUTPUT,
                "Too few sample points available (%" CPL_SIZE_FORMAT " point(s)) to make the fit "
                "(more than %" CPL_SIZE_FORMAT " points needed). "
                "Increase number of sample points or increase kappa",
                cpl_table_get_nrow(t),  deg_xy);
        
            /* Fit, calculate Zfit */
            uves_polynomial_delete(&background);
            check( background = uves_polynomial_regression_2d(
                   t, "X", "Y", "Z", NULL,
                   DEGX, DEGY, "Zfit", NULL, NULL, &mse,
                   NULL, NULL, -1, -1),
               "Error fitting polynomial");
        
            /* Residual := Z - Zfit */
            cpl_table_duplicate_column(t, "Residual", t, "Z");
            cpl_table_subtract_columns(t, "Residual", "Zfit");
        
            /* Compute residuals w.r.t. median of Z 
               (i.e. subtract median(residual) from all residuals),
               then get stdev based on this new mean/median value.
               This is to make kappa sigma clipping more robust */

            cpl_table_subtract_scalar(t, "Residual", 
                          cpl_table_get_column_median(t, "Residual"));
            rmse = cpl_table_get_column_stdev(t, "Residual");

            /* One-sided kappa-sigma clipping */
            if (KAPPA > 0)
            {
                check( n_clipped = uves_select_table_rows(
                       t,  "Residual", CPL_GREATER_THAN, KAPPA * rmse),
                   "Error selecting rows");
            }
            else
            {
                n_clipped = 0;
            }
            
            total_clipped += n_clipped;
        
            uves_msg_debug("RMS = %f. %" CPL_SIZE_FORMAT " of %" CPL_SIZE_FORMAT " points rejected in kappa-sigma clipping",
                   rmse, n_clipped, cpl_table_get_nrow(t));
            
            cpl_table_erase_selected(t);

            if (n_clipped > 0)
            {
                cpl_table_erase_column(t, "Zfit");
                cpl_table_erase_column(t, "Residual");
            }
        
        } while (n_clipped > 0);
        }

        /* Try to do some quality checking of the background subtraction.
           The number of rejected points (the signal) is often around 10-20 %  */
        {
        double percentage = 
            100.0 * ( (double)total_clipped ) / (total_clipped + cpl_table_get_nrow(t));
        
        if (KAPPA > 0) {
            uves_msg("%" CPL_SIZE_FORMAT " of %" CPL_SIZE_FORMAT " points (%.2f %%) were rejected in "
                 "kappa-sigma clipping. RMS = %.2f ADU", 
                 total_clipped,
                 cpl_table_get_nrow(t) + total_clipped,
                 percentage,
                 sqrt(mse));
        }
        
        /* For grid sampling: */
        if (orders == NULL && order_locations == NULL) 
            {
            if (total_clipped == 0)
                {
                uves_msg_warning("No points rejected during background "
                         "estimation. Background subtraction is "
                         "uncertain. Try to decrease KAPPA "
                         "(current value is %f)", KAPPA);
                }
            if (percentage > 40)
                {
                uves_msg_warning("%f %% of the sample points were "
                         "rejected during "
                         "background estimation", percentage);
                }
            }
        }
        
        check( subtract_background(image, NULL, background),
           "Error subtracting background polynomial");
    } /* BM_METHOD wasn't 'no' */
    
  cleanup:
    uves_free_table(&t);
    uves_polynomial_delete(&background);
    
    return cpl_error_get_code();
}

/*----------------------------------------------------------------------------*/
/**
   @brief   Subtract background of an echelle image
   @param   image   Subtract background of this image
   @param   RADX   x-radius of filter window
   @param   RADY   y-radius of filter window
   @param   ITER   number of iterations
   @return  CPL_ERROR_NONE iff OK.

   This function imitates the MIDAS command BACKGROUND/SMOOTH: The background
   of the image is estimated by repeatedly lowering the image fluxes to the
   local average of a box of size (2*RADX+1)*(2*RADY+1). The background is then
   subtracted from the original image. 
*/
/*----------------------------------------------------------------------------*/
/* Recipe parameter creation code for this function
/ * Backsmoothx, Backsmoothy * /
    uves_parameter_new_range(p, uves_orderpos.preproc.backsmoothx,
                CPL_TYPE_INT,
                "Radius of window used for average filtering in the "
                "background subtraction (mode=smooth) step",
                uves_orderpos.preproc,
                5, 0, INT_MAX);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "backsmoothx");
    cpl_parameterlist_append(recipe->parameters, p);
    
    uves_parameter_new_range(p, uves_orderpos.preproc.backsmoothy,
                CPL_TYPE_INT,
                "Radius of window used for average filtering in the "
                "background subtraction (mode=smooth) step",
                uves_orderpos.preproc,
                30, 0, INT_MAX);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "backsmoothy");
    cpl_parameterlist_append(recipe->parameters, p);

/ * Backsmoothiter * /
    uves_parameter_new_range(p, uves_orderpos.preproc.backsmoothiter,
                CPL_TYPE_INT,
                "Number of iterations when estimating the background "
                "(mode=smooth)",
                uves_orderpos.preproc,
                10, 1, INT_MAX);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "backsmoothiter");
    cpl_parameterlist_append(recipe->parameters, p);
*/
cpl_error_code
uves_backsub_smooth(cpl_image *image, int RADX, int RADY, int ITER)
{
    cpl_image  *background  = NULL;
    int i;
    
    assure( RADX >= 0 && RADY >= 0, CPL_ERROR_ILLEGAL_INPUT,
        "Negative radius ((%d)x(%d))", RADX, RADY);
    assure( ITER >= 1, CPL_ERROR_ILLEGAL_INPUT, 
        "Non-positive number of iterations (%d)", ITER);
    
    /* First estimate background */
    background = cpl_image_duplicate(image);
    
    for (i = 0; i < ITER; i++) {
      //uves_msg_debug("i=%d,%d ...",i, ITER);
    uves_msg("i = %d", i);
    check( lower_to_average(background,
                RADX, RADY), "Error smoothing image");
    }
    
    /* Then subtract background */
    check( cpl_image_subtract(image, background), "Could not subtract background image");
    
  cleanup:
    uves_free_image(&background);

    return cpl_error_get_code();
}

/*----------------------------------------------------------------------------*/
/**
   @brief   Sample background level
   @param   image   The image to look at.
   @param   x0      x-center of window.
   @param   y0      y-center of window.
   @param   radius_x Sample window half width
   @param   radius_y Sample window half height
   @param   nx       Image width
   @param   ny       Image height
   @param   BM_METHOD Background measure method. Supported values are BM_MEDIAN or BM_MINIMUM.
                      This parameter defines whether to use the median or minimum of the
              specified window.                                             
   @return  The median/minimum flux of the specified window, or undefined on error.

   The window center is at (@em x0, @em y0), where @em y0 is determined
   from the polynomial. The window slope follows the polynomial.
   
*/
/*----------------------------------------------------------------------------*/

static double
sample_background(const cpl_image *image, int x0, double y_0,
          int radius_x, int radius_y, int nx, int ny,
          background_measure_method BM_METHOD)
{
    double result = 0;
    /* Use a table to calculate the median. Invalid rows are ignored */
    cpl_table *temp = NULL;
    bool found_good = false;
    int row;
    int x, y;

    check( 
    (temp = cpl_table_new( (2*radius_x + 1) * (2*radius_y + 1) ),
     row = 0,
     cpl_table_new_column(temp, "Flux", CPL_TYPE_DOUBLE)),
    "Error allocating table");

    for(y = y_0 - radius_y; y <= y_0 + radius_y; y++)
    {
        for (x = x0 - radius_x; x <= x0 + radius_x; x++)
        {
            if (1 <= x && x <= nx &&
            1 <= y && y <= ny)
            {
                int pis_rejected;
                double flux = cpl_image_get(image, x, y, &pis_rejected);
                if( !pis_rejected )
                {
                    cpl_table_set(temp, "Flux", row, flux);
                    found_good = true;
                }
                else
                {
                    cpl_table_set_invalid(temp, "Flux", row);
                }
            }
            else
            {
                cpl_table_set_invalid(temp, "Flux", row);
            }
            
            row++;
        }
    }

    assure( found_good, CPL_ERROR_ILLEGAL_INPUT, "No valid pixels in sample window");

    if (BM_METHOD == BM_MEDIAN)
    {
        result = cpl_table_get_column_median(temp, "Flux");
    }
    else if (BM_METHOD == BM_MINIMUM)
    {
        result = cpl_table_get_column_min(temp, "Flux");
    }
    else
    {
        assure( false, CPL_ERROR_UNSUPPORTED_MODE,
            "Unsupported background sample method: %d", BM_METHOD);
    }

  cleanup:
    uves_free_table(&temp);
    return result;
}

/*----------------------------------------------------------------------------*/
/**
   @brief   Get minimum order number
   @param   order_locations   Polynomial defining the order positions
   @param   nx                Image width
   @return  The minimum order number

   This function returns the minimum (first) order partially inside the image.
*/
/*----------------------------------------------------------------------------*/
static int
first_order(const polynomial *order_locations, int nx)
{
    int result;
    
    result = 0;
    while (uves_polynomial_evaluate_2d(order_locations, 1 , result + 0.5) < 1 ||
       uves_polynomial_evaluate_2d(order_locations, nx, result + 0.5) < 1 )
    {
        result++;
    }

    while (uves_polynomial_evaluate_2d(order_locations, 1 , result - 0.5) >= 1 ||
       uves_polynomial_evaluate_2d(order_locations, nx, result - 0.5) >= 1 )
    {
        result -= 1;
        
        /* Fail cleanly even if 'order_locations' is corrupted */
        assure( result > -100000, 
            CPL_ERROR_CONTINUE,
            "Invalid polynomial: p(x=1, order=%d) = %f  p(x=%d, order=%d) = %f",
            result, uves_polynomial_evaluate_2d(order_locations, 1.0, result),
            nx, result, uves_polynomial_evaluate_2d(order_locations, nx, result));
    }
    
  cleanup:
    return result;
}


/*----------------------------------------------------------------------------*/
/**
   @brief   Get maximum order number
   @param   order_locations   Polynomial defining the order positions
   @param   nx                Image width
   @param   ny                Image height
   @return  The maximum order number

   This function returns the maximum (last) order partially inside the image.
*/
/*----------------------------------------------------------------------------*/
static int
last_order(const polynomial *order_locations, int nx, int ny)
{
    int result;
    
    result = 0;
    while (uves_polynomial_evaluate_2d(order_locations, 1 , result - 0.5) > ny ||
       uves_polynomial_evaluate_2d(order_locations, nx, result - 0.5) > ny )
    {
        result--;
    }

    while (uves_polynomial_evaluate_2d(order_locations, 1 , result + 1.5) <= ny ||
       uves_polynomial_evaluate_2d(order_locations, nx, result + 1.5) <= ny )
    {
        result += 1;
        
        /* Fail cleanly even if 'order_locations' is corrupted */
        assure( result < 100000, 
            CPL_ERROR_CONTINUE,
            "Invalid polynomial: p(x=1, order=%d) = %f  p(x=%d, order=%d) = %f",
            result, uves_polynomial_evaluate_2d(order_locations, 1.0, result),
            nx, result, uves_polynomial_evaluate_2d(order_locations, nx, result));
    }
    
  cleanup:
    return result;
}

/*----------------------------------------------------------------------------*/
/**
   @brief   Smooth an image
   @param   image            Image to smooth
   @param   RADX             X-radius of filter window
   @param   RADY             Y-radius of filter window
   @return  CPL_ERROR_NONE iff OK.

   Each pixel in the image is lowered to the average of the local neighbourhood. If
   the pixel is below average, it is left unchanged. The size of the window used for calculating
   the average is (2*RADX + 1, 2*RADY + 1).
*/
/*----------------------------------------------------------------------------*/
static cpl_error_code
lower_to_average(cpl_image *image, int RADX, int RADY)
{
    cpl_image  *average = NULL;
    double *image_data = NULL;
    double *average_data = NULL;
    int nx, ny;
    int x, y;
    
    passure( image != NULL, "Null image");
    nx = cpl_image_get_size_x(image);
    ny = cpl_image_get_size_y(image);
    
    /* Create smoothed image */
    uves_msg("Filtering...");
    check( average    = cpl_image_duplicate(image), "Error copying image");
    check( uves_filter_image_average(average, RADX, RADY), "Error applying average filter");
    uves_msg("done");
    
    image_data   = cpl_image_get_data(image);
    average_data = cpl_image_get_data(average);
    uves_msg("Lowering...");
    for (y = 0; y < ny; y++)
    {
        for (x = 0; x < nx; x++)
        {
            if (image_data[x + y*nx] > average_data[x + y*nx]) 
            {
                image_data[x + y*nx] = average_data[x + y*nx];
            }
        }
    }
    uves_msg("done");
    
  cleanup:
    uves_free_image(&average);
    
    return cpl_error_get_code();
}

/*----------------------------------------------------------------------------*/
/**
   @brief   Subtract the previously defined background 
   @param   image            Image to be background subtracted
   @param   background_im    Background image, may be NULL. If non-NULL,
                             this is updated to contain the flux
                 values that were actually subtracted.
   @param   background_pol   Background polynomial, may be NULL
   @return  CPL_ERROR_NONE iff OK.

   Exactly one of @em background_im and @em background_pol must be non-NULL.
*/
/*----------------------------------------------------------------------------*/
    
static cpl_error_code
subtract_background(cpl_image *image, cpl_image *background_im, 
            const polynomial *background_pol)
{
    int nx, ny;
    int x, y;

    double *image_data;
    double *background_data = NULL;

    passure(image != NULL, " ");
    /* Exactly one of 'background_im' and 'background_pol' must be non-NULL */
    passure((background_im == NULL) != (background_pol == NULL), " ");

    /* For efficiency, don't call cpl_image_get() */
    /* The following check is too strict. It can be avoided to solve PIPE-4893
    assure(cpl_image_count_rejected(image) == 0, 
       CPL_ERROR_UNSUPPORTED_MODE, "Input image contains bad pixels");
       */
    assure(cpl_image_get_type(image) == CPL_TYPE_DOUBLE,
       CPL_ERROR_UNSUPPORTED_MODE, 
       "Input image is of type %s. double expected", 
       uves_tostring_cpl_type(cpl_image_get_type(image)));

    if (background_im != NULL)
    {
        assure(cpl_image_count_rejected(background_im) == 0, 
           CPL_ERROR_UNSUPPORTED_MODE, "Background image contains bad pixels");
        assure(cpl_image_get_type(background_im) == CPL_TYPE_DOUBLE, 
           CPL_ERROR_UNSUPPORTED_MODE, 
           "Background image is of type %s. double expected", 
           uves_tostring_cpl_type(cpl_image_get_type(background_im)));
    }

    image_data = cpl_image_get_data_double(image);
    if (background_im != NULL)
    {
        background_data = cpl_image_get_data_double(background_im);
    }

    nx = cpl_image_get_size_x(image);
    ny = cpl_image_get_size_y(image);

    for (y = 1; y <= ny; y++)
    {
        for (x = 1; x <= nx; x++)
        {
            double back;
            double flux, new_flux;
            
            if (background_im != NULL)
            {
                /* Slow:  back = cpl_image_get(background_im, x, y, &pis_rejected); */
                back = background_data[(x-1) + (y-1) * nx]; 
            }
            else
            {
                /* Evaluate at (x,y) */
                back = uves_polynomial_evaluate_2d(background_pol, 
                                   x,
                                   y);
            }
            
            /* Slow: flux = cpl_image_get(image, x, y, &pis_rejected);  */
            flux = image_data[(x-1) + (y-1) * nx];
            
/* Exclude these sanity checks for backwards compatibility */
#if 0
            /* Make sure the estimated background is between zero and flux-value */
            if (back < 0)
            {
                back = 0.0;
            }
            if (back > flux)
            {
                back = flux;
            }
            
            /* Then subtract the background.
             * Pixel flux may be negative. Make sure the result is non-negative.
             */
                    new_flux = uves_max_double(0, flux - back);
#else
            new_flux = flux-back;            
#endif
            
            /* Slow: cpl_image_set(image, x, y, new_flux); */
            image_data[(x-1) + (y-1) * nx] = new_flux; 
            
            if (background_im != NULL)
            {
                /* Slow: cpl_image_set(background_im, x, y, flux - new_flux); */
                background_data[(x-1) + (y-1) * nx] = flux - new_flux;
            }
        }
    }/* for each pixel... */
    
  cleanup:
    return cpl_error_get_code();
}
/**@}*/