/*                                                                              *
 *   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 02110-1301 USA          *
 */
 
/*
 * $Author: amodigli $
 * $Date: 2011-11-02 15:43:23 $
 * $Revision: 1.62 $
 */

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

/*----------------------------------------------------------------------------*/
/**
 * @defgroup uves_mbias  Recipe: Master Bias
 *
 * This recipe calculates the master bias frame.
 * See man-page for details.
 */
/*----------------------------------------------------------------------------*/

/*-----------------------------------------------------------------------------
                                Includes
 -----------------------------------------------------------------------------*/
#include <uves_mbias_impl.h>

#include <uves_utils.h>
#include <uves_corrbadpix.h>
#include <uves_parameters.h>
#include <uves.h>
#include <uves_dfs.h>
#include <uves_pfits.h>
#include <uves_qclog.h>
#include <uves_recipe.h>
#include <uves_utils_wrappers.h>
#include <uves_error.h>
#include <irplib_mkmaster.h>
#include <uves_msg.h>

/* Library */
#include <cpl.h>
#include <float.h>
#include <ctype.h>
#include <string.h>
/*-----------------------------------------------------------------------------
                            Functions prototypes
 -----------------------------------------------------------------------------*/

static void uves_mbias_qclog(const cpl_imagelist* raw_imgs,
                             uves_propertylist **raw_headers,
                 enum uves_chip chip,
                 const cpl_image* mbia,
                 /* int sx_pix, Size of X bin in pix 
                 int sy_pix, Size of Y bin in pix */
                            cpl_table* qclog
    );

static void  
uves_mbias_qc_ron_raw(const cpl_image* rbia,
                      enum uves_chip chip,
                      const int x_cent_s,
                      const int x_cent_e,
                      const int y_cent_s,
                      const int y_cent_e,
                      cpl_table* qclog);

static int
uves_mbias_define_parameters(cpl_parameterlist *parameters);

/*-----------------------------------------------------------------------------
                            Recipe standard code
 -----------------------------------------------------------------------------*/
#define cpl_plugin_get_info uves_mbias_get_info
UVES_RECIPE_DEFINE(
    UVES_MBIAS_ID, UVES_MBIAS_DOM,
    /* Warning: if more parameters are added to this recipe, they 
       need to be propagated to uves_cal_mkmaster! */
    uves_mbias_define_parameters,
    "Jonas M. Larsen", "cpl@eso.org",
    "Creates the master bias frame", 
    "This recipe creates a master bias frame by computing the median of all input\n"
    "bias frames. All input frames must have same tag and size and must be either\n"
    "BIAS_BLUE or BIAS_RED.\n"
    "On blue input the recipe computes one master bias frame; on red input a \n"
    "master bias frame for each chip is produced. The average, standard deviation\n"
    "and median of the master bias image(s) are written to the FITS header(s)");


/*-----------------------------------------------------------------------------
                            Functions code
 -----------------------------------------------------------------------------*/
/**@{*/

/*----------------------------------------------------------------------------*/
/**
  @brief    Setup the recipe options    
  @param    parameters        the parameterlist to fill
  @param    recipe_id         name of calling recipe
  @return   0 if everything is ok
 */
/*----------------------------------------------------------------------------*/
int uves_mbias_define_parameters_body(cpl_parameterlist *parameters,
                       const char *recipe_id)
{

    /*****************
     *    General    *
     *****************/
    if (uves_define_global_parameters(parameters) != CPL_ERROR_NONE)
        {
            return -1;
        }

    /**************************************
     *  detector's trap correction        *
     **************************************/
 
    if (uves_corr_traps_define_parameters(parameters,recipe_id) 
        != CPL_ERROR_NONE)
        {
            return -1;
        }

    /**************************************
     *  Master stack generation           *
     **************************************/

    if (uves_master_stack_define_parameters(parameters,recipe_id) 
        != CPL_ERROR_NONE)
        {
            return -1;
        }

    return (cpl_error_get_code() != CPL_ERROR_NONE);
}



/*----------------------------------------------------------------------------*/
/**
  @brief    Setup the recipe options    
  @param    parameters        the parameterlist to fill
  @return   0 if everything is ok
 */
/*----------------------------------------------------------------------------*/
static int
uves_mbias_define_parameters(cpl_parameterlist *parameters)
{
    return uves_mbias_define_parameters_body(parameters, make_str(UVES_MBIAS_ID));
}





/*----------------------------------------------------------------------------*/
/**
  @brief    Process a single chip
  @param    raw_images   The input images
  @param    raw_headers  FITS headers of input images
  @param    mbias_header (output)
  @param    binx         x-binning
  @param    biny         y-binning
  @param    chip        CCD chip
  @return   The master bias image

  This function 
  - computes the master bias image by taking
  the pixel-by-pixel median of all input frames.
  - corrects bad pixels (the layout of bad pixels
  is hardcoded, valid for UVES chips).
 */
/*----------------------------------------------------------------------------*/
static cpl_image *
uves_mbias_process_chip(const cpl_imagelist *raw_images, 
                        uves_propertylist **raw_headers, 
                        uves_propertylist *mbias_header,
                        int binx, int biny,
                        enum uves_chip chip,
                        bool CLEAN_TRAPS,
                        const char* STACK_METHOD,
                        const double STACK_KLOW,
                        const double STACK_KHIGH,
                        const int STACK_NITER)
{
    cpl_image *master_bias        = NULL; /* Result */
    double exposure_time = 0;
    int badpixels_cleaned;
    int i;
    int nraw=0;
    bool red_ccd_is_new=false;
    cpl_vector* bias_levels=NULL;
    double bias_mean=0;

    uves_msg("Calculating master bias...");

    check_nomsg(red_ccd_is_new=uves_ccd_is_new(raw_headers[0]));
    /* Get the median at each pixel */
    if(strcmp(STACK_METHOD,"MEDIAN")==0) {
       uves_msg("method median");
       master_bias=irplib_mkmaster_median(raw_images,5.,5,1.e-5);
    } else {
       uves_msg("method mean");
       master_bias=irplib_mkmaster_mean(raw_images,5.,5,1.e-5,STACK_KLOW,STACK_KHIGH,STACK_NITER);
    }

    /* Set mbias exposure time to average of inputs */
    exposure_time = 0;
    nraw=cpl_imagelist_get_size(raw_images);
    for (i = 0; i < nraw; i++)
    {
        check( exposure_time += uves_pfits_get_exptime(raw_headers[i]), 
           "Error reading exposure time");
    }
    exposure_time /= nraw;

    uves_pfits_set_extname(mbias_header,"Master bias");
    
    check( uves_pfits_set_exptime(mbias_header, exposure_time),
       "Error setting master bias exposure time");
   if(CLEAN_TRAPS) {
      check( badpixels_cleaned = 
             uves_correct_badpix_all(master_bias, mbias_header, 
                                     chip, binx, biny, false,red_ccd_is_new),
             "Error replacing bad pixels");

      uves_msg("%d bad pixels replaced", badpixels_cleaned);
   }

  cleanup:
    if (cpl_error_get_code() != CPL_ERROR_NONE)
    {
        uves_free_image(&master_bias);
    }
    
    return master_bias;
}

/*----------------------------------------------------------------------------*/
/**
  @brief    Executor function
  @param    parameters  the parameters list
  @param    frames      the frames list
  @param    starttime   execution start
 */
/*----------------------------------------------------------------------------*/
static void
UVES_CONCAT2X(UVES_MBIAS_ID,exe)(cpl_frameset *frames, 
                   const cpl_parameterlist *parameters,
                   const char *starttime)
{
    uves_mbias_exe_body(frames, parameters, starttime, make_str(UVES_MBIAS_ID));
    return;
}

/*----------------------------------------------------------------------------*/
/**
  @brief    Get the command line options and execute the data reduction
  @param    parameters  the parameters list
  @param    frames      the frames list

  After computing the master bias frame, the pixel average, standard deviation
  and median values are also computed and written in appropriate keywords in the
  output image header.

 */
/*----------------------------------------------------------------------------*/
void
uves_mbias_exe_body(cpl_frameset *frames, 
            const cpl_parameterlist *parameters,
            const char *starttime,
            const char *recipe_id)
{
    /* Input */
    cpl_imagelist *raw_images[2]        = {NULL, NULL};
    uves_propertylist **raw_headers[2]   = {NULL, NULL};    /* Two arrays of pointers */

    cpl_table* qclog[2] = {NULL, NULL};

    /* Output */
    uves_propertylist *product_header[2] = {NULL, NULL};
    cpl_image *master_bias              = NULL;
    cpl_stats *mbias_stats              = NULL;
    
    /* Local variables */
    char *product_filename = NULL;
    bool blue;
    enum uves_chip chip;
    int binx, biny;
    const char* PROCESS_CHIP=NULL;
    bool CLEAN_TRAPS;
    int raw_index = 0;
    int i=0;

    const char* STACK_METHOD=NULL;
    double STACK_KLOW=0;
    double STACK_KHIGH=0;
    int STACK_NITER=0;

    /* Load and check raw bias images and headers, identify arm (blue/red) */
    /* On success, 'raw_headers' will be arrays with the same length as 'raw_images' */
    check( uves_load_raw_imagelist(frames, 
                   false,     /* FLAMES format? (no) */
                   UVES_BIAS(true), UVES_BIAS(false),
                   CPL_TYPE_DOUBLE,
                   raw_images, raw_headers, product_header, 
                   &blue), "Error loading raw frames");
    
    /* Get binning from first header (i.e. BLUE or REDL chip, first raw frame) */
    check( binx = uves_pfits_get_binx(raw_headers[0][0]), 
       "Could not get raw frame x-binning");
    check( biny = uves_pfits_get_biny(raw_headers[0][0]), 
       "Could not get raw frame y-binning");
    check( uves_get_parameter(parameters, NULL, "uves", "process_chip", CPL_TYPE_STRING, &PROCESS_CHIP),
               "Could not read parameter");
    uves_string_toupper((char*)PROCESS_CHIP);

    check( uves_get_parameter(parameters, NULL, recipe_id, "clean_traps", CPL_TYPE_BOOL, &CLEAN_TRAPS),
               "Could not read parameter");

    check( uves_get_parameter(parameters, NULL, recipe_id, "stack_method", CPL_TYPE_STRING, &STACK_METHOD),
               "Could not read parameter");
    uves_string_toupper((char*)STACK_METHOD);

    check( uves_get_parameter(parameters, NULL, recipe_id, "klow", CPL_TYPE_DOUBLE, &STACK_KLOW),
               "Could not read parameter");
    check( uves_get_parameter(parameters, NULL, recipe_id, "khigh", CPL_TYPE_DOUBLE, &STACK_KHIGH),
               "Could not read parameter");
    check( uves_get_parameter(parameters, NULL, recipe_id, "niter", CPL_TYPE_INT, &STACK_NITER),
               "Could not read parameter");

    /* Loop over one or two chips */
    for (chip = uves_chip_get_first(blue); 
     chip != UVES_CHIP_INVALID;
     chip = uves_chip_get_next(chip))
    {
       if(strcmp(PROCESS_CHIP,"REDU") == 0) {
	chip = uves_chip_get_next(chip);
      }
        raw_index = uves_chip_get_index(chip);

        uves_msg("Processing %s chip",
             uves_chip_tostring_upper(chip));

        uves_msg_debug("Binning = %dx%d", binx, biny);

        /* Process chip */
        uves_free_image(&master_bias);
        check( master_bias = uves_mbias_process_chip(raw_images[raw_index], 
                                                     raw_headers[raw_index],
                                                     product_header[raw_index],
                                                     binx, biny,
                                                     chip,CLEAN_TRAPS,
                                                     STACK_METHOD,
                                                     STACK_KLOW,
                                                     STACK_KHIGH,
                                                     STACK_NITER),
           "Error processing chip");
        
        
        cpl_free(product_filename);
        check( product_filename = uves_masterbias_filename(chip), 
           "Error getting filename");

        /* Finished. Calculate QC parameters and save */
        uves_msg("Calculating QC parameters");
            uves_qclog_delete(&qclog[0]);
            qclog[0] = uves_qclog_init(raw_headers[raw_index][0], chip);
            check(uves_mbias_qclog(raw_images[raw_index],
                                   raw_headers[raw_index],
                                   chip,
                                   master_bias,
                   /* binx,biny, */
                   qclog[0]),"error computing qclog");

        /* Insert into frame set */
        uves_msg("Saving product...");
        
        check( uves_frameset_insert(frames,
                    master_bias,
                    CPL_FRAME_GROUP_PRODUCT,
                    CPL_FRAME_TYPE_IMAGE,
                    CPL_FRAME_LEVEL_INTERMEDIATE,
                    product_filename,
                    UVES_MASTER_BIAS(chip),
                    raw_headers[raw_index][0], /* First frame */
                    product_header[raw_index],
                    NULL,
                    parameters,
                    recipe_id,
                    PACKAGE "/" PACKAGE_VERSION,qclog,
                    starttime, true, 
                    UVES_ALL_STATS),
           "Could not add master bias %s to frameset", product_filename);
            uves_qclog_delete(&qclog[0]);
        uves_msg("Master bias '%s' added to frameset", product_filename);

        if(strcmp(PROCESS_CHIP,"REDL") == 0) {
	chip = uves_chip_get_next(chip);
      }


        } /* For each chip */

  cleanup:
    /* Input */
    if (raw_images[0] != NULL)
    {

        for (i = 0; i < cpl_imagelist_get_size(raw_images[0]); i++) 
        {
            if (raw_headers[0] != NULL) uves_free_propertylist(&raw_headers[0][i]);
            if (raw_headers[1] != NULL) uves_free_propertylist(&raw_headers[1][i]);
        }
        cpl_free(raw_headers[0]); raw_headers[0] = NULL;
        cpl_free(raw_headers[1]); raw_headers[1] = NULL;
    }
    uves_free_imagelist(&raw_images[0]);
    uves_free_imagelist(&raw_images[1]);
    /* Output */

    uves_qclog_delete(&qclog[0]);
    uves_free_image(&master_bias);
    uves_free_propertylist(&product_header[0]);
    uves_free_propertylist(&product_header[1]);
    cpl_free(product_filename);
    uves_free_stats(&mbias_stats);
    
    return;
}



static int
count_good(const cpl_image *image)
{
    return 
        cpl_image_get_size_x(image) * cpl_image_get_size_y(image) - 
        cpl_image_count_rejected(image);
}
/*----------------------------------------------------------------------------*/
/**
  @brief    Reject outlier pixels
  @param    image       image with pixels
  @param    min         reject pixels below this value
  @param    max         reject pixels abouve this value
 */
/*----------------------------------------------------------------------------*/
static void
reject_lo_hi(cpl_image *image, double min, double max)
{
  cpl_mask *mask_lo = NULL;
  cpl_mask *mask_hi = NULL;

  mask_lo = cpl_mask_threshold_image_create(image, -DBL_MAX, min);
  mask_hi = cpl_mask_threshold_image_create(image, max, DBL_MAX);
  assure_mem( mask_lo );
  assure_mem( mask_hi );

  cpl_mask_or(mask_lo, mask_hi);

  cpl_image_reject_from_mask(image, mask_lo);
  
  cleanup:
  uves_free_mask(&mask_lo);
  uves_free_mask(&mask_hi);
  return;
}

/**
@brief computes QC log
@param raw_imgs list of input images
@param raw_headers list of input headers
@param chip chip ID
@param mbia master bias
@param qclog table with QC log parameters
@return QC parameters on a given product 
*/

static void uves_mbias_qclog(const cpl_imagelist* raw_imgs,
                             uves_propertylist **raw_headers,
                 enum uves_chip chip,
                 const cpl_image* mbia,
                 /* int sx_pix, Size of X bin in pix 
                 int sy_pix, Size of Y bin in pix */
                            cpl_table* qclog
                            ) 
{
  int nx_pix= 0;  /* No of X pix */
  int ny_pix= 0;  /* No of Y pix */

  int sample_x= 100; /* X size of sampling window in pix */
  int sample_y= 100; /* Y size of sampling window in pix */
  int x_cent_s= 0;   /* X sampling window starting point */
  int x_cent_e= 0;   /* X sampling window ending point */
  int y_cent_s= 0;   /* Y sampling window starting point */
  int y_cent_e= 0;   /* Y sampling window ending point */





  double upp_threshold= 0.0;
  double low_threshold= 0.0;
  double extra=0.1;
  double qc_ron_master= 0.0;

  double master_median=0.0;
  int pn= 0;


  double min=0.0;
  double max=0.0;
  double struct_col=0.0;
  double struct_row=0.0;

  double time_s=+9999999.0;
  double time_e=-9999999.0;
  int nraw=0;
  double qc_duty_cycle=0.;
  double exposure_time=0;
  int i=0;
  char key_name[80];

  const cpl_image* rbia=NULL;
  cpl_image* tima=NULL;
  cpl_image* avg_col=NULL;
  cpl_image* avg_row=NULL;
  
  uves_qclog_add_string(qclog,
                        "QC TEST1 ID",
                        "Test-on-Master-Bias",
                        "Name of QC test",
                        "%s");

  uves_msg("Computing duty cycle...");

  /* Set mbias exposure time to average of inputs */
  exposure_time = 0;
  nraw = cpl_imagelist_get_size(raw_imgs);
  check_nomsg(uves_qclog_add_int(qclog,
                        "PRO DATANCOM",
                        nraw,
                        "Number of frames combined",
                        "%d"));


  for (i = 0; i < nraw; i++)
    {
      check( exposure_time = uves_pfits_get_mjdobs(raw_headers[i]),
         "Error reading exposure time");
      if(exposure_time >= time_e) time_e = exposure_time;
      if(exposure_time <= time_s) time_s = exposure_time;
    }
  if(nraw > 1) {
    qc_duty_cycle = (time_e-time_s)/ (nraw-1);
  }
  else
      {
      qc_duty_cycle = 0;
      }
 
  check_nomsg(uves_qclog_add_double(qclog,
                        "QC DUTYCYCL",
                        qc_duty_cycle,
                        "Time to store a frame",
                        "%.5e"));

  /* The following is not really used in MIDAS so we comment
  strcpy(date,uves_pfits_get_tpl_start(plist));
  */


  /* CONVERT FROM MIDAS
  nx_pix = m$value({mbia},NPIX(1));
  ny_pix = m$value({mbia},NPIX(2));
  */

  nx_pix = cpl_image_get_size_x(mbia);
  ny_pix = cpl_image_get_size_y(mbia);


  x_cent_s = (nx_pix - sample_x)/2;
  x_cent_e = (nx_pix + sample_x)/2;
  y_cent_s = (ny_pix - sample_y)/2;
  y_cent_e = (ny_pix + sample_y)/2;
 
  
  check_nomsg(upp_threshold = 
          cpl_image_get_median_window(mbia,
                      x_cent_s,
                      y_cent_s,
                      x_cent_e,
                      y_cent_e)*(1 + extra));
  check_nomsg(low_threshold = 
          cpl_image_get_median_window(mbia,
                      x_cent_s,
                      y_cent_s,
                      x_cent_e,
                      y_cent_e)*(1 - extra));
  
  /* convert from MIDAS
  pn = {uves_portid({PATHID})};
  */
  check_nomsg(pn = PORT_ID(chip));
  uves_msg_debug("Port number = %d", pn);

  rbia = cpl_imagelist_get_const(raw_imgs,0);
  check_nomsg(uves_mbias_qc_ron_raw(rbia, chip,
                    x_cent_s,x_cent_e,y_cent_s,y_cent_e,qclog));


  /* convert from  MIDAS
  stat/ima {mbia} + bins=1 exc={low_threshold},{upp_threshold};
  */
  check_nomsg(tima=cpl_image_duplicate(mbia));

  check_nomsg( reject_lo_hi(tima, low_threshold, upp_threshold) );
  if (count_good(tima) >= 2)
      {
          check_nomsg(master_median = cpl_image_get_median(tima));
          check_nomsg(qc_ron_master = cpl_image_get_stdev(tima));
      }
  else
      {
          master_median = -1;
          qc_ron_master = -1;
          uves_msg_warning("Only %d good pixels in image. Setting QC parameters to -1",
                           count_good(tima));
      }
  uves_free_image(&tima);

  check_nomsg(uves_qclog_add_double(qclog,
                        "PRO DATAMED",
                        master_median,
                        "Median of pixel values",
                        "%7.3f"));

  sprintf(key_name, "QC OUT%d RON MASTER", pn);
  check_nomsg(uves_qclog_add_double(qclog,
                        key_name,
                        qc_ron_master,
                        "Read noise frame in ADU",
                        "%8.4f"));

  /* ==========================
   * Calculates Bias struct
   * ==========================
   */

 
  /*
   * in case of RED frame cuts out values greater than 300.
   * as the frame can be affected by this local operation we 
   * do the calculation on a copy of the original frame
   */

  
  check_nomsg(tima=cpl_image_duplicate(mbia));
  if (chip != UVES_CHIP_BLUE) {
    /*
    replace/ima {mbia} {tmpfrm} 300,>=300.;
    */
    check_nomsg(cpl_image_threshold(tima,
                                    -DBL_MAX,300,
                                    -DBL_MAX,300));
  }


  check_nomsg(avg_col = cpl_image_collapse_create(tima,1));
  check_nomsg(cpl_image_divide_scalar(avg_col,cpl_image_get_size_x(tima)));

  /* restricts statistics to +/- 2 ADU around mean */
  min = cpl_image_get_mean(avg_col) - 2;
  max = cpl_image_get_mean(avg_col) + 2; 

  /* replace with MIDAS
  stat/ima avg_col + exc={min},{max};
  */
  check_nomsg( reject_lo_hi(avg_col, min, max) );
  if (count_good(avg_col) >= 2)
      {
          check_nomsg(struct_col = cpl_image_get_stdev(avg_col));
      }
  else
      {
          struct_col = -1;
          uves_msg_warning("Only %d good pixels in image. Setting QC parameter to -1",
                           count_good(avg_col));
      }

  sprintf(key_name,"%s%d%s","QC OUT",pn," STRUCTY");
  check_nomsg(uves_qclog_add_double(qclog,
                        key_name,
                        struct_col,
                        "structure in Y (bias slope)",
                        "%8.4f"));



  check_nomsg(avg_row = cpl_image_collapse_create(tima,0));
  check_nomsg(cpl_image_divide_scalar(avg_row,cpl_image_get_size_y(tima)));

  /* restricts statistics to +/- 2 ADU around mean */
  min = cpl_image_get_mean(avg_row) - 2;
  max = cpl_image_get_mean(avg_row) + 2;
 
  /* replace with MIDAS
  stat/ima avg_row + exc={min},{max};
  */
  check_nomsg( reject_lo_hi(avg_row, min, max) );
  if (count_good(avg_row) >= 2)
      {
          check_nomsg(struct_row = cpl_image_get_stdev(avg_row));
      }
  else
      {
          struct_row = -1;
          uves_msg_warning("Only %d good pixels in image. Setting QC parameter to -1",
                           count_good(avg_row));
      }

  
  sprintf(key_name,"%s%d%s","QC OUT",pn," STRUCTX");
  check_nomsg(uves_qclog_add_double(qclog,
                        key_name,
                        struct_row,
                        "structure in X (bias slope)",
                        "%8.4f"));

 

 
 cleanup:
  uves_free_image(&avg_col);
  uves_free_image(&avg_row);
  uves_free_image(&tima);
  
  return;

}
/**
@brief computes Read Out Noise on raw frame
@param rbia raw bias frame
@param mbia master bias frame
@param raw_index index of frame in set
@param x_cent_s  start X of region to evaluate
@param x_cent_e  end   X of region to evaluate
@param y_cent_s  start Y of region to evaluate
@param y_cent_e  end   Y of region to evaluate
@param qclog table with QC parameters
@return Read Out Noise on raw frame

*/
static void
uves_mbias_qc_ron_raw(const cpl_image* rbia,
                      enum uves_chip chip,
                      const int x_cent_s,
                      const int x_cent_e,
                      const int y_cent_s,
                      const int y_cent_e,
                      cpl_table* qclog)
{

  double qc_ron_raw=0.0;
  double upp_threshold=0.0;
  double low_threshold=0.0;
  double extra=0.1;
  char key_name[80];
  int pn=0;
  cpl_image* tima=NULL;

  /* replace with MIDAS
  date   = "{{mbia},ESO.TPL.START}";
  store/frame infrm {incat} 1;
  */
  check_nomsg(upp_threshold = 
          cpl_image_get_median_window(rbia,
                      x_cent_s,
                      y_cent_s,
                      x_cent_e,
                      y_cent_e)*(1 + extra));
  
  check_nomsg(low_threshold = 
          cpl_image_get_median_window(rbia,
                      x_cent_s,
                      y_cent_s,
                      x_cent_e,
                      y_cent_e)*(1 - extra));
  


  /* replace from MIDAS
  stat/ima {rbia} + bins=1 exc={low_treshold},{upp_treshold};
  */
  check_nomsg(tima=cpl_image_duplicate(rbia));

  check_nomsg( reject_lo_hi(tima, low_threshold, upp_threshold) );
  if (count_good(tima) >= 2)
      {
          check_nomsg(qc_ron_raw = cpl_image_get_stdev(tima));
      }
  else
      {
          qc_ron_raw = -1;
          uves_msg_warning("Only %d good pixels in image. Setting QC parameter to -1",
                           count_good(tima));
      }


  /* replace from MIDAS
  pn = {uves_portid({PATHID})};
  */
  check_nomsg(pn = PORT_ID(chip));


  sprintf(key_name,"%s%d%s","QC OUT",pn," RON RAW");
  check_nomsg(uves_qclog_add_double(qclog,
                        key_name,
                        qc_ron_raw,
                        "Read noise frame in ADU",
                        "%8.4f"));

 cleanup:
  uves_free_image(&tima);
  return;
}

/**@}*/