/* $Id: vimos_ima_bias.c,v 1.9 2015/09/11 09:48:16 jim Exp $
 *
 * This file is part of the VIMOS Pipeline
 * Copyright (C) 2015 European Southern Observatory
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 */

/*
 * $Author: jim $
 * $Date: 2015/09/11 09:48:16 $
 * $Revision: 1.9 $
 * $Name:  $
 */

/* Includes */

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

#include <stdio.h>
#include <cpl.h>
#include <math.h>
#include <string.h>

#include "vmutils.h"
#include "vimos_imaging_utils.h"
#include "vimos_pfits.h"
#include "vimos_dfs.h"
#include "casu_utils.h"
#include "casu_stats.h"
#include "casu_fits.h"
#include "casu_wcsutils.h"
#include "casu_mods.h"

/* Define values for bit mask that flags dummy results */

#define MEANBIAS    1
#define DIFFIMG     2
#define STATS_TAB   4

/* Structure definitions */

typedef struct {

    /* Input */

    int         combtype;
    int         xrej;
    float       thresh;
    int         ncells;
    int         prettynames;

    /* Output */

    float       biasmed;
    float       biasrms;
    float       biasdiff_med;
    float       biasdiff_rms;
    float       underscan_med;
    float       underscan_rms;
    float       overscan_med;
    float       overscan_rms;

    /* Bitmasks used to make sure we get all the products we're expecting */

    int         we_expect;
    int         we_get;
} configstruct;

typedef struct {
    cpl_size          *labels;
    cpl_size          *labels2;
    cpl_frameset      *biaslist;
    casu_fits          **biases;
    int               nbiases;
    cpl_frame         *ref_bias;
    cpl_image         *outimage;
    cpl_propertylist  *drs;
    casu_fits          *rbimage;
    cpl_image         *diffimg;
    cpl_table         *diffimstats;
} memstruct;

/* Function prototypes */

static int vimos_ima_bias_create(cpl_plugin *);
static int vimos_ima_bias_exec(cpl_plugin *);
static int vimos_ima_bias_destroy(cpl_plugin *);
static int vimos_ima_bias(cpl_parameterlist *, cpl_frameset *) ;
static int vimos_ima_bias_save(cpl_frameset *framelist, 
                                   cpl_parameterlist *parlist, char *cname,
                                   configstruct *cs, memstruct *ps, int isfirst,
                                   cpl_frame **product_frame_mean_bias,
                                   cpl_frame **product_frame_diffimg,
                                   cpl_frame **product_frame_diffimg_stats);
static void vimos_ima_bias_dummy_products(configstruct *cs, memstruct *ps);
static void vimos_ima_bias_normal(int jext, configstruct *cs,
                                      memstruct *ps, char *vimos_names[]);
static int vimos_ima_bias_lastbit(cpl_frameset *framelist,
                                      cpl_parameterlist *parlist, 
                                      configstruct *cs, memstruct *ps,
                                      char *detname, int isfirst,
                                      cpl_frame **product_frame_mean_bias,
                                      cpl_frame **product_frame_diffimg,
                                      cpl_frame **product_frame_diffimg_stats);
static void vimos_ima_bias_init(memstruct *ps);
static void vimos_ima_bias_tidy(memstruct *ps, int level);


static char vimos_ima_bias_description[] =
"vimos_ima_bias -- VIMOS bias combine recipe.\n\n"
"Combine a list of bias frames into a mean bias frame. Optionally compare \n"
"the output frame to a reference bias frame\n\n"
"The program accepts the following files in the SOF:\n\n"
"    Tag                   Description\n"
"    -----------------------------------------------------------------------\n"
"    %-21s A list of raw bias images\n"
"    %-21s Optional reference bias frame\n"
"If no reference bias frame is made available, then no comparison will\n"
"be done. This means there will be no output difference image or stats table. "
"\n";

/**@{*/

/**

    \ingroup recipelist
    \defgroup vimos_ima_bias vimos_ima_bias
    \brief Combine a list of bias frames to form a mean bias frame

    \par Name: 
        vimos_ima_bias
    \par Purpose: 
        Combine a list of bias frames to form a mean bias frame. 
    \par Description: 
        A list of bias frames is combined with rejection to form a mean bias. 
        If a reference bias is supplied, then a difference image is formed 
        between it and the combined result from the current frame list. 
        This difference image can be useful for looking at the evolution 
        of the structure of the bias. The difference image is broken into 
        lots of little cells and the median value of the difference image 
        as well as the RMS in each cell is written to a difference image 
        statistics table.
    \par Language:
        C
    \par Parameters:
        - \b combtype (string): Determines the type of combination that is 
          done to form the output map. Can take the following values:
            - median: The output pixels are medians of the input pixels
            - mean: The output pixels are means of the input pixels
        - \b xrej (int): If set, then an extra rejection cycle will be run. 
        - \b thresh (float): The rejection threshold in numbers of sigmas.
        - \b ncells (int): If a difference image statistics table is being 
          done, then this is the number of cells in which to divide each
          image. The value should be a power of 2, up to 64.
        - \b prettynames (bool): If set then a descriptive file name will
          be used for the output products. Otherwise the standard ESO 
          'predictable' name will be used.
    \par Input File Types:
        The following list gives the file types that can appear in the SOF
        file. The word in bold is the DO CATG value.
        - \b BIAS (required): A list of raw bias images for combining
        - \b REFERENCE_BIAS (optional): A library reference bias frame.
             If this is given, then a difference image will be formed and some 
             basic statistics run on it.
    \par Output Products:
        The following list gives the output data products that are generated
        by this recipe. The word in bold gives the PRO CATG keyword value for
        each product:
        - The output mean/median bias frame, formed by either a mean or median
          combination of the input raw frames with rejection (\b MASTER_BIAS).
        - The output difference image, formed by subtracting the input reference
          bias frame from the current mean/median bias frame. This is only 
          done if a reference bias frame is specified from the outset
          (\b DIFFIMG_BIAS).
        - The output difference image statistics table. 
          (\b DIFFIMG_STATS_BIAS).
    \par Output QC Parameters:
        - \b BIASMED
            The median value of the combined bias frame.
        - \b BIASRMS
            The RMS value of the mean bias frame.
        - \b BIAS_DIFFMED
            The median of the bias difference image
        - \b BIAS_DIFFRMS
            The RMS of the BIAS difference image
        - \b USCAN_MED
            The median of the underscan region of the combined frame
        - \b USCAN_RMS
            The RMS of the underscan region of the combined frame
        - \b OSCAN_MED
            The median of the overscan region of the combined frame
        - \b OSCAN_RMS
            The RMS of the overscan region of the combined frame
    \par Notes
        None.
    \par Fatal Error Conditions:
        - NULL input frameset
        - Input frameset headers incorrect meaning that RAW and CALIB frame
          can't be distinguished
        - No bias frames in the input frameset
        - Unable to save data products
        - Unable to load any of the frames in the SOF
    \par Non-Fatal Error Conditions:
        - No reference bias. No difference image formed.
    \par Conditions Leading To Dummy Products:
        - The detector for the current image extension is flagged dead
        - Combination routine failed
        - Reference bias frame image extension won't load or is flagged as 
          a dummy
    \par Author:
        Jim Lewis, CASU
    \par Code Reference: 
        vimos_ima_bias.c

*/


/* Function code */


/*---------------------------------------------------------------------------*/
/**
  @brief    Build the list of available plugins, for this module.
  @param    list    the plugin list
  @return   0 if everything is ok

  This function is exported.
 */
/*---------------------------------------------------------------------------*/

int cpl_plugin_get_info(cpl_pluginlist *list) {
    cpl_recipe  *recipe = cpl_calloc(1,sizeof(*recipe));
    cpl_plugin  *plugin = &recipe->interface;
    char alldesc[SZ_ALLDESC];
    (void)snprintf(alldesc,SZ_ALLDESC,vimos_ima_bias_description,
                   VIMOS_BIAS_RAW,VIMOS_REF_BIAS);

    cpl_plugin_init(plugin,
                    CPL_PLUGIN_API,
                    VIMOS_BINARY_VERSION,
                    CPL_PLUGIN_TYPE_RECIPE,
                    "vimos_ima_bias",
                    "Bias combination for imaging",
                    alldesc,
                    "Jim Lewis",
                    "jrl@ast.cam.ac.uk",
                    vimos_get_license(),
                    vimos_ima_bias_create,
                    vimos_ima_bias_exec,
                    vimos_ima_bias_destroy);

    cpl_pluginlist_append(list,plugin);

    return(0);
}

/*---------------------------------------------------------------------------*/
/**
  @brief    Setup the recipe options
  @param    plugin  the plugin
  @return   0 if everything is ok

  Create the recipe instance and make it available to the application using the
  interface.
 */
/*---------------------------------------------------------------------------*/

static int vimos_ima_bias_create(cpl_plugin *plugin) {
    cpl_recipe      *recipe;
    cpl_parameter   *p;

    /* Get the recipe out of the plugin */

    if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE)
        recipe = (cpl_recipe *)plugin;
    else 
        return(-1);

    /* Create the parameters list in the cpl_recipe object */

    recipe->parameters = cpl_parameterlist_new();

    /* Fill in the parameters. First the combination type */

    p = cpl_parameter_new_enum("vimos.vimos_ima_bias.combtype",
                               CPL_TYPE_STRING,
                               "Combination algorithm",
                               "vimos.vimos_ima_bias",
                               "median",2,"median","mean");
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"combtype");
    cpl_parameterlist_append(recipe->parameters,p);

    /* Extra rejection cycle */

    p = cpl_parameter_new_value("vimos.vimos_ima_bias.xrej",
                                CPL_TYPE_BOOL,
                                "True if using extra rejection cycle",
                                "vimos.vimos_ima_bias",
                                TRUE);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"xrej");
    cpl_parameterlist_append(recipe->parameters,p);

    /* Rejection threshold */

    p = cpl_parameter_new_range("vimos.vimos_ima_bias.thresh",
                                CPL_TYPE_DOUBLE,
                                "Rejection threshold in sigma above background",
                                "vimos.vimos_ima_bias",5.0,1.0e-6,1.0e10);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"thresh");
    cpl_parameterlist_append(recipe->parameters,p);

    /* How many cells to divide each image */

    p = cpl_parameter_new_enum("vimos.vimos_ima_bias.ncells",
                               CPL_TYPE_INT,
                               "Number of cells for diff image stats",
                               "vimos.vimos_ima_bias",64,7,1,2,4,8,
                               16,32,64);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"ncells");
    cpl_parameterlist_append(recipe->parameters,p);     

    /* Flag to use pretty output product file names */

    p = cpl_parameter_new_value("vimos.vimos_ima_bias.prettynames",
                                CPL_TYPE_BOOL,
                                "Use pretty output file names?",
                                "vimos.vimos_ima_bias",FALSE);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"prettynames");
    cpl_parameterlist_append(recipe->parameters,p);     

    /* Get out of here */

    return(0);
}
    
    
/*---------------------------------------------------------------------------*/
/**
  @brief    Execute the plugin instance given by the interface
  @param    plugin  the plugin
  @return   0 if everything is ok
 */
/*---------------------------------------------------------------------------*/

static int vimos_ima_bias_exec(cpl_plugin *plugin) {
    cpl_recipe  *recipe;

    /* Get the recipe out of the plugin */

    if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE)
        recipe = (cpl_recipe *)plugin;
    else 
        return(-1);

    return(vimos_ima_bias(recipe->parameters,recipe->frames));
}
                                
/*---------------------------------------------------------------------------*/
/**
  @brief    Destroy what has been created by the 'create' function
  @param    plugin  the plugin
  @return   0 if everything is ok
 */
/*---------------------------------------------------------------------------*/

static int vimos_ima_bias_destroy(cpl_plugin *plugin) {
    cpl_recipe *recipe ;

    /* Get the recipe out of the plugin */

    if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE)
        recipe = (cpl_recipe *)plugin;
    else 
        return(-1);

    cpl_parameterlist_delete(recipe->parameters);
    return(0);
}

/*---------------------------------------------------------------------------*/
/**
  @brief    The recipe data reduction part is implemented here 
  @param    parlist     the parameters list
  @param    framelist   the frames list
  @return   0 if everything is ok
 */
/*---------------------------------------------------------------------------*/

static int vimos_ima_bias(cpl_parameterlist *parlist, 
                              cpl_frameset *framelist) {
    const char *fctid="vimos_ima_bias";
    int j,retval,status,isfirst;
    cpl_size nlab,nlab2;
    unsigned char *rejmask,*rejplus;
    cpl_parameter *p;
    cpl_frameset *chiplist;
    cpl_frame *product_frame_mean_bias,*product_frame_diffimg;
    cpl_frame *product_frame_diffimg_stats,*testfrm;
    cpl_image *junk;
    configstruct cs;
    memstruct ps;
    char *vimos_names[VIMOS_NEXTN];

    /* Check validity of input frameset */

    if (framelist == NULL || cpl_frameset_get_size(framelist) <= 0) {
        cpl_msg_error(fctid,"Input framelist NULL or has no input data");
        return(-1);
    }

    /* Initialise some things */

    cs.we_expect = 0;
    cs.we_get = 0;
    vimos_ima_bias_init(&ps);
    cs.we_expect |= MEANBIAS;

    /* Get the parameters */

    p = cpl_parameterlist_find(parlist,"vimos.vimos_ima_bias.combtype");
    cs.combtype = (strcmp(cpl_parameter_get_string(p),"median") ? 2 : 1);
    p = cpl_parameterlist_find(parlist,"vimos.vimos_ima_bias.xrej");
    cs.xrej = cpl_parameter_get_bool(p);
    p = cpl_parameterlist_find(parlist,"vimos.vimos_ima_bias.thresh");
    cs.thresh = (float)cpl_parameter_get_double(p);
    p = cpl_parameterlist_find(parlist,"vimos.vimos_ima_bias.ncells");
    cs.ncells = cpl_parameter_get_int(p);
    p = cpl_parameterlist_find(parlist,"vimos.vimos_ima_bias.prettynames");
    cs.prettynames = (cpl_parameter_get_bool(p) ? 1 : 0);

    /* Sort out raw from calib frames */

    if (vimos_dfs_set_groups(framelist) != CASU_OK) {
        cpl_msg_error(fctid,"Cannot identify RAW and CALIB frames");
        vimos_ima_bias_tidy(&ps,0);
        return(-1);
    }

    /* Get a list of the frame labels */ 

    if ((ps.labels = cpl_frameset_labelise(framelist,casu_compare_tags,
                                           &nlab)) == NULL) {
        cpl_msg_error(fctid,"Cannot labelise the input frames");
        vimos_ima_bias_tidy(&ps,0);
        return(-1);
    }

    /* Get all the bias frames frames */

    if ((ps.biaslist = casu_frameset_subgroup(framelist,ps.labels,nlab,
                                              VIMOS_BIAS_RAW)) == NULL) {
        cpl_msg_error(fctid,"Cannot find bias frames in input frameset");
        vimos_ima_bias_tidy(&ps,0);
        return(-1);
    }

    /* Now labelise the input frameset so that we separate them by the
       detector name */

    if ((ps.labels2 = cpl_frameset_labelise(ps.biaslist,vimos_compare_names,
                                            &nlab2)) == NULL) {
        cpl_msg_error(fctid,"Cannot labelise the bias frames");
        vimos_ima_bias_tidy(&ps,0);
        return(-1);
    }

    /* See what chip set this is */

    if (vimos_load_names(cpl_frameset_get_position(ps.biaslist,0),
            vimos_names) != CASU_OK) {
        cpl_msg_error(fctid,"Cannot get the name set");
        vimos_ima_bias_tidy(&ps,0);
        return(-1);
    }

    /* Check to see if there is a reference bias frame */

    if ((ps.ref_bias = casu_frameset_subgroup_1(framelist,ps.labels,nlab,
                                                VIMOS_REF_BIAS)) == NULL) {
        cpl_msg_info(fctid,"No reference bias found -- no difference image will be formed");
    } else {
        cs.we_expect |= DIFFIMG;
        cs.we_expect |= STATS_TAB;
    }
        
    /* Now loop for all the detectors... */

    for (j = 1; j <= VIMOS_NEXTN; j++) {
        status = CASU_OK;
        cs.we_get = 0;
        isfirst = (j == 1);

        /* Get the subset of images in the input frameset that belong to
           each detector */

        if ((chiplist = vimos_get_det_frameset(ps.biaslist,ps.labels2,nlab2,
                                               vimos_names[j-1])) == NULL) {
            cpl_msg_info(fctid,"No biases in input list for detector %s\n",
                         vimos_names[j-1]);
            retval = vimos_ima_bias_lastbit(framelist,parlist,&cs,&ps,
                                                vimos_names[j-1],isfirst,
                                                &product_frame_mean_bias,
                                                &product_frame_diffimg,
                                                &product_frame_diffimg_stats);
            if (retval != 0) 
                return(-1);
            continue;
        }

        /* Load up the images. If they won't load then signal a major error,
           create some dummy products and save them. */

        testfrm = cpl_frameset_get_position(chiplist,0);
        ps.biases = casu_fits_load_list(chiplist,CPL_TYPE_FLOAT,
                                        cpl_frame_get_nextensions(testfrm));
        ps.nbiases = cpl_frameset_get_size(chiplist);
        if (ps.biases == NULL) {
            cpl_msg_info(fctid,"Detector %s biases wouldn't load",
                         vimos_names[j-1]);
            retval = vimos_ima_bias_lastbit(framelist,parlist,&cs,&ps,
                                                vimos_names[j-1],isfirst,
                                                &product_frame_mean_bias,
                                                &product_frame_diffimg,
                                                &product_frame_diffimg_stats);
            if (retval != 0) 
                return(-1);
            continue;
        }

        /* Call the combine module. If it fails, then signal that
           all products will be dummies */

        cpl_msg_info(fctid,"Doing combination for detector %s",
                     vimos_names[j-1]);
        (void)casu_imcombine(ps.biases,NULL,ps.nbiases,cs.combtype,0,cs.xrej,
                             cs.thresh,"",&(ps.outimage),&junk,&rejmask,
                             &rejplus,&(ps.drs),&status);
        freespace(rejmask);
        freespace(rejplus);
        if (status == CASU_OK) {
            cs.we_get |= MEANBIAS;
            vimos_ima_bias_normal(j,&cs,&ps,vimos_names);
        } 

        /* Create any dummies and save the products */

        retval = vimos_ima_bias_lastbit(framelist,parlist,&cs,&ps,
                                            vimos_names[j-1],isfirst,
                                            &product_frame_mean_bias,
                                            &product_frame_diffimg,
                                            &product_frame_diffimg_stats);

        freeframeset(chiplist);
        if (retval != 0) 
            return(-1);
    }
    vimos_ima_bias_tidy(&ps,0);
    return(0);
}

/*---------------------------------------------------------------------------*/
/**
  @brief    The recipe data products are saved here
  @param    framelist    The input frame list
  @param    parlist      The input recipe parameter list
  @param    cname        The name of the current detector
  @param    cs           The configuration structure defined in the main routine
  @param    ps           The memory structure defined in the main routine
  @param    isfirst      Set if this is the first extension to be saved
  @param    product_frame_mean_bias   The frame for the mean bias product
  @param    product_frame_diffimg     The frame for the difference image product
  @param    product_frame_diffimg_stats The frame for the diffimg stats table
  @return   0 if everything is ok
 */
/*---------------------------------------------------------------------------*/

static int vimos_ima_bias_save(cpl_frameset *framelist, 
                                   cpl_parameterlist *parlist, char *cname,
                                   configstruct *cs, memstruct *ps, int isfirst,
                                   cpl_frame **product_frame_mean_bias,
                                   cpl_frame **product_frame_diffimg,
                                   cpl_frame **product_frame_diffimg_stats) {
    cpl_propertylist *plist,*elist,*p;
    int status,i,night;
    const char *fctid = "vimos_ima_bias_save";
    const char *esoout[] = {"biascomb.fits","biasdiff.fits","biasdifftab.fits"};
    const char *prettypfx[] = {"bias","biasdiff","biasdifftab"};
    char outfile[3][BUFSIZ],dateobs[81];
    const char *recipeid = "vimos_ima_bias";

    /* Get info to make pretty names if requested. */

    for (i = 0; i < 3; i++) 
        strcpy(outfile[i],esoout[i]);
    if (cs->prettynames) {
        plist = casu_fits_get_phu(ps->biases[0]);
        elist = casu_fits_get_ehu(ps->biases[0]);
        if (vimos_pfits_get_dateobs(plist,dateobs) != CASU_OK &&
            vimos_pfits_get_dateobs(elist,dateobs) != CASU_OK) {
            cpl_msg_warning(fctid,"Missing header information. Reverting to predictable names");
        } else {
            night = casu_night_from_dateobs(dateobs);
            for (i = 0; i < 3; i++) 
                (void)sprintf(outfile[i],"%s_%8d.fits",prettypfx[i],night);
        }
    }

    /* If we need to make a PHU then do that now. */

    if (isfirst) {

        /* Create a new product frame object and define some tags */

        *product_frame_mean_bias = cpl_frame_new();
        cpl_frame_set_filename(*product_frame_mean_bias,outfile[0]);
        cpl_frame_set_tag(*product_frame_mean_bias,VIMOS_PRO_BIAS);
        cpl_frame_set_type(*product_frame_mean_bias,CPL_FRAME_TYPE_IMAGE);
        cpl_frame_set_group(*product_frame_mean_bias,CPL_FRAME_GROUP_PRODUCT);
        cpl_frame_set_level(*product_frame_mean_bias,CPL_FRAME_LEVEL_FINAL);

        /* Base the header on the first image in the input framelist */

        plist = cpl_propertylist_duplicate(casu_fits_get_ehu(ps->biases[0]));
        vimos_dfs_set_product_primary_header(plist,*product_frame_mean_bias,
                                             framelist,parlist,
                                             (char *)recipeid,"PRO-1.15",
                                             NULL,0);

        /* 'Save' the PHU image */                       

        if (cpl_image_save(NULL,outfile[0],CPL_TYPE_UCHAR,plist,
                           CPL_IO_DEFAULT) != CPL_ERROR_NONE) {
            cpl_msg_error(fctid,"Cannot save product PHU");
            cpl_frame_delete(*product_frame_mean_bias);
            return(-1);
        }
        cpl_frameset_insert(framelist,*product_frame_mean_bias);
        cpl_propertylist_delete(plist);

        /* Create a new product frame object for the difference image */

        if (cs->we_expect & DIFFIMG) {
            *product_frame_diffimg = cpl_frame_new();
            cpl_frame_set_filename(*product_frame_diffimg,outfile[1]);
            cpl_frame_set_tag(*product_frame_diffimg,VIMOS_PRO_DIFFIMG_BIAS);
            cpl_frame_set_type(*product_frame_diffimg,CPL_FRAME_TYPE_IMAGE);
            cpl_frame_set_group(*product_frame_diffimg,CPL_FRAME_GROUP_PRODUCT);
            cpl_frame_set_level(*product_frame_diffimg,CPL_FRAME_LEVEL_FINAL);

            /* Base the header on the first image in the input framelist */

            plist = cpl_propertylist_duplicate(casu_fits_get_ehu(ps->biases[0]));
            vimos_dfs_set_product_primary_header(plist,*product_frame_diffimg,
                                                 framelist,parlist,
                                                 (char *)recipeid,
                                                 "PRO-1.15",NULL,0);

            /* 'Save' the PHU image */                   

            if (cpl_image_save(NULL,outfile[1],CPL_TYPE_UCHAR,plist,
                               CPL_IO_DEFAULT) != CPL_ERROR_NONE) {
                cpl_msg_error(fctid,"Cannot save product PHU");
                cpl_frame_delete(*product_frame_diffimg);
                return(-1);
            }
            cpl_frameset_insert(framelist,*product_frame_diffimg);
            cpl_propertylist_delete(plist);
        }

        /* Create a new product frame object for the difference image stats 
           table */

        if (cs->we_expect & STATS_TAB) {
            *product_frame_diffimg_stats = cpl_frame_new();
            cpl_frame_set_filename(*product_frame_diffimg_stats,outfile[2]);
            cpl_frame_set_tag(*product_frame_diffimg_stats,
                              VIMOS_PRO_DIFFIMG_BIAS_STATS);
            cpl_frame_set_type(*product_frame_diffimg_stats,
                               CPL_FRAME_TYPE_TABLE);
            cpl_frame_set_group(*product_frame_diffimg_stats,
                                CPL_FRAME_GROUP_PRODUCT);
            cpl_frame_set_level(*product_frame_diffimg_stats,
                                CPL_FRAME_LEVEL_FINAL);

            /* Base the header on the first image in the input framelist */

            plist = cpl_propertylist_duplicate(casu_fits_get_ehu(ps->biases[0]));
            vimos_dfs_set_product_primary_header(plist,
                                                 *product_frame_diffimg_stats,
                                                 framelist,parlist,
                                                 (char *)recipeid,
                                                 "PRO-1.15",NULL,0);

            /* Fiddle with the extension header now */

            elist = casu_fits_get_ehu(ps->biases[0]);
            p = cpl_propertylist_duplicate(elist);
            casu_merge_propertylists(p,ps->drs);
            if (! (cs->we_get & STATS_TAB))
                casu_dummy_property(p);
            cpl_propertylist_update_string(p,"EXTNAME",cname);
            vimos_dfs_set_product_exten_header(p,*product_frame_diffimg_stats,
                                               framelist,parlist,
                                               (char *)recipeid,
                                               "PRO-1.15",NULL);
            status = CASU_OK;
            casu_removewcs(p,&status);

            /* And finally save the difference image stats table */

            if (cpl_table_save(ps->diffimstats,plist,p,outfile[2],
                               CPL_IO_DEFAULT) != CPL_ERROR_NONE) {
                cpl_msg_error(fctid,"Cannot save product table extension");
                cpl_frame_delete(*product_frame_diffimg_stats);
                cpl_propertylist_delete(p);
                return(-1);
            }
            cpl_propertylist_delete(plist);
            cpl_propertylist_delete(p);
            cpl_frameset_insert(framelist,*product_frame_diffimg_stats);
        }
    }

    /* Get the extension property list */

    plist = cpl_propertylist_duplicate(casu_fits_get_ehu(ps->biases[0]));
    cpl_propertylist_update_int(plist,"ESO PRO DATANCOM",ps->nbiases);

    /* Fiddle with the header now */

    casu_merge_propertylists(plist,ps->drs);
    p = cpl_propertylist_duplicate(plist);
    if (! (cs->we_get & MEANBIAS))
        casu_dummy_property(p);
    cpl_propertylist_update_string(p,"EXTNAME",cname);
    vimos_dfs_set_product_exten_header(p,*product_frame_mean_bias,
                                       framelist,parlist,
                                       (char *)recipeid,"PRO-1.15",NULL);
                
    /* Now save the mean bias image extension */

    cpl_propertylist_update_float(p,"ESO QC BIASMED",cs->biasmed);
    cpl_propertylist_set_comment(p,"ESO QC BIASMED",
                                 "Median of mean bias frame");
    cpl_propertylist_update_float(p,"ESO QC BIASRMS",cs->biasrms);
    cpl_propertylist_set_comment(p,"ESO QC BIASRMS","RMS of mean bias frame");
    cpl_propertylist_update_float(p,"ESO QC USCAN_MED",cs->underscan_med);
    cpl_propertylist_set_comment(p,"ESO QC USCAN_MED",
                                 "Median of mean bias frame underscan");
    cpl_propertylist_update_float(p,"ESO QC USCAN_RMS",cs->underscan_rms);
    cpl_propertylist_set_comment(p,"ESO QC USCAN_RMS",
                                 "RMS of mean bias frame underscan");
    cpl_propertylist_update_float(p,"ESO QC OSCAN_MED",cs->overscan_med);
    cpl_propertylist_set_comment(p,"ESO QC OSCAN_MED",
                                 "Median of mean bias frame overscan");
    cpl_propertylist_update_float(p,"ESO QC OSCAN_RMS",cs->overscan_rms);
    cpl_propertylist_set_comment(p,"ESO QC OSCAN_RMS",
                                 "RMS of mean bias frame overscan");
    if (cpl_image_save(ps->outimage,outfile[0],CPL_TYPE_FLOAT,p,
                       CPL_IO_EXTEND) != CPL_ERROR_NONE) {
        cpl_msg_error(fctid,"Cannot save product image extension");
        cpl_propertylist_delete(p);
        return(-1);
    }
    cpl_propertylist_delete(p);

    /* Now save the bias difference image extension */

    if (cs->we_expect & DIFFIMG) {
        p = cpl_propertylist_duplicate(plist);
        if (! (cs->we_get & DIFFIMG))
            casu_dummy_property(p);;
        cpl_propertylist_update_float(p,"ESO QC BIAS_DIFFMED",cs->biasdiff_med);
        cpl_propertylist_set_comment(p,"ESO QC BIAS_DIFFMED",
                                     "Median of bias difference image");
        cpl_propertylist_update_float(p,"ESO QC BIAS_DIFFRMS",cs->biasdiff_rms);
        cpl_propertylist_set_comment(p,"ESO QC BIAS_DIFFRMS",
                                     "RMS of bias difference image");
        cpl_propertylist_update_string(p,"EXTNAME",cname);
        vimos_dfs_set_product_exten_header(p,*product_frame_diffimg,framelist,
                                           parlist,(char *)recipeid,"PRO-1.15",
                                           NULL);
        if (cpl_image_save(ps->diffimg,outfile[1],CPL_TYPE_FLOAT,p,
                           CPL_IO_EXTEND) != CPL_ERROR_NONE) {
            cpl_msg_error(fctid,"Cannot save product image extension");
            cpl_propertylist_delete(p);
            return(-1);
        }
        cpl_propertylist_delete(p);
    }

    /* Now any further difference image stats tables */

    if (! isfirst && (cs->we_expect & STATS_TAB)) {
        p = cpl_propertylist_duplicate(plist);
        if (! (cs->we_get & STATS_TAB))
            casu_dummy_property(p);
        cpl_propertylist_update_string(p,"EXTNAME",cname);
        vimos_dfs_set_product_exten_header(p,*product_frame_diffimg_stats,
                                           framelist,parlist,
                                           (char *)recipeid,
                                           "PRO-1.15",NULL);
        status = CASU_OK;
        casu_removewcs(p,&status);
        if (cpl_table_save(ps->diffimstats,NULL,p,outfile[2],CPL_IO_EXTEND)
                           != CPL_ERROR_NONE) {
            cpl_msg_error(fctid,"Cannot save product table extension");
            cpl_propertylist_delete(p);
            return(-1);
        }       
        cpl_propertylist_delete(p);
    }
    cpl_propertylist_delete(plist);

    return(0);
}

/*---------------------------------------------------------------------------*/
/**
  @brief    Fill undefined products with dummy products
  @param    cs           The configuration structure defined in the main routine
  @param    ps           The memory structure defined in the main routine
  @return   Nothing
 */
/*---------------------------------------------------------------------------*/

static void vimos_ima_bias_dummy_products(configstruct *cs, memstruct *ps) {

    /* See if you even need to be here */

    if (cs->we_get == cs->we_expect)
        return;

    /* We always expect a mean frame. If we don't have one, then create
       a dummy */

    if (! (cs->we_get & MEANBIAS)) {
        ps->outimage = casu_dummy_image(ps->biases[0]);

        /* Set up the QC parameters */
    
        cs->biasmed = 0.0;
        cs->biasrms = 0.0;
        cs->underscan_med = 0.0;
        cs->underscan_rms = 0.0;
        cs->overscan_med = 0.0;
        cs->overscan_rms = 0.0;
    }

    /* Do the difference image */

    if ((cs->we_expect & DIFFIMG) && ! (cs->we_get & DIFFIMG)) {
        cs->biasdiff_med = 0.0;
        cs->biasdiff_rms = 0.0;
        ps->diffimg = casu_dummy_image(ps->biases[0]);
    }

    /* If a difference image stats table is required, then do that now */

    if ((cs->we_expect & STATS_TAB) && ! (cs->we_get & STATS_TAB))
        ps->diffimstats = vimos_create_diffimg_stats(0);

    return;
}

/*---------------------------------------------------------------------------*/
/**
  @brief    Do stats on bias frame and create diff image and stats table
  @param    jext         the extension number
  @param    cs           The configuration structure defined in the main routine
  @param    ps           The memory structure defined in the main routine
  @params   vimos_names  The names of the detectors in the chip set
 */
/*---------------------------------------------------------------------------*/

static void vimos_ima_bias_normal(int jext, configstruct *cs, 
                                      memstruct *ps, char *vimos_names[]) {
    int nx,ny,ncells,i,j,n,k;
    long npi;
    float med,sig,*idata,grms,gdiff,*work;
    const char *fctid="vimos_ima_bias_normal";

    /* Work out the median and RMS of the mean bias frame */

    nx = cpl_image_get_size_x(ps->outimage);
    ny = cpl_image_get_size_y(ps->outimage);
    npi = nx*ny;
    idata = cpl_image_get_data(ps->outimage);
    casu_medmad(idata,NULL,npi,&med,&sig);
    sig *= 1.48;
    cs->biasmed = med;
    cs->biasrms = sig;

    /* Work out the median and RMS of the underscan and overscan regions */

    npi = (underlim[1] - underlim[0] + 1)*(underlim[3] - underlim[2] + 1);
    work = cpl_malloc(npi*sizeof(float));
    n = 0;
    for (j = underlim[2]; j <= underlim[3]; j++) {
        for (i = underlim[0]; i <= underlim[1]; i++) {
            k = (j-1)*nx + i - 1;
            work[n++] = idata[k];
        }
    }
    casu_medmad(work,NULL,n,&med,&sig);
    sig *= 1.48;
    cs->underscan_med = med;
    cs->underscan_rms = sig;
    freespace(work);
    npi = (overlim[1] - overlim[0] + 1)*(overlim[3] - overlim[2] + 1);
    work = cpl_malloc(npi*sizeof(float));
    n = 0;
    for (j = overlim[2]; j <= overlim[3]; j++) {
        for (i = overlim[0]; i <= overlim[1]; i++) {
            k = (j-1)*nx + i - 1;
            work[n++] = idata[k];
        }
    }
    casu_medmad(work,NULL,n,&med,&sig);
    sig *= 1.48;
    cs->overscan_med = med;
    cs->overscan_rms = sig;
    freespace(work);

    /* Load up the reference bias */

    if (ps->ref_bias != NULL) {
        ps->rbimage = casu_fits_load(ps->ref_bias,CPL_TYPE_FLOAT,jext);
        if (ps->rbimage == NULL) 
            cpl_msg_info(fctid,
                         "Reference bias extension %" CPL_SIZE_FORMAT " (%s) won't load",
                         (cpl_size)jext,vimos_names[jext-1]);
        else if (casu_is_dummy(casu_fits_get_ehu(ps->rbimage))) {
            cpl_msg_info(fctid,
                         "Reference bias extension %" CPL_SIZE_FORMAT " (%s) is a dummy!",
                         (cpl_size)jext,vimos_names[jext-1]);
            freefits(ps->rbimage);
        }
    } else 
        ps->rbimage = NULL;

    /* Form the difference image. NB: the difference image routine copes
       if the input mean image is NULL. Thus if it is null because of a 
       failure to load then the routine will do as much as it can and return
       allowing you to fill in the rest with dummy products */

    cs->biasdiff_med = 0.0;
    cs->biasdiff_rms = 0.0;
    ncells = cs->ncells;
    vimos_difference_image(casu_fits_get_image(ps->rbimage),ps->outimage,
                           NULL,ncells,1,&gdiff,&grms,&(ps->diffimg),
                           &(ps->diffimstats));
    cs->biasdiff_med = gdiff;
    cs->biasdiff_rms = grms;
    if (ps->diffimg != NULL)
        cs->we_get |= DIFFIMG;
    if (ps->diffimstats != NULL)
        cs->we_get |= STATS_TAB;
    return;
}

/*---------------------------------------------------------------------------*/
/**
  @brief    Make any dummies that are required and then save all products
  @param    framelist    the input frame list
  @param    parlist      the input recipe parameter list
  @param    cs           The configuration structure defined in the main routine
  @param    ps           The memory structure defined in the main routine
  @param    name         The name of the current detector
  @param    isfirst      Set if this is the first extension to be saved
  @param    product_frame_mean_bias   The frame for the mean bias product
  @param    product_frame_diffimg     The frame for the difference image product
  @param    product_frame_diffimg_stats The frame for the diffimg stats table
  @return   0 if everything is ok
 */
/*---------------------------------------------------------------------------*/

static int vimos_ima_bias_lastbit(cpl_frameset *framelist,
                                      cpl_parameterlist *parlist, 
                                      configstruct *cs, memstruct *ps,
                                      char *name, int isfirst,
                                      cpl_frame **product_frame_mean_bias,
                                      cpl_frame **product_frame_diffimg,
                                      cpl_frame **product_frame_diffimg_stats) {
    int retval;
    const char *fctid="vimos_ima_bias_lastbit";

    /* Make whatever dummy products you need */

    vimos_ima_bias_dummy_products(cs,ps);

    /* Save everything */

    cpl_msg_info(fctid,"Saving products for detector %s",name);
    retval = vimos_ima_bias_save(framelist,parlist,name,cs,ps,isfirst,
                                     product_frame_mean_bias,
                                     product_frame_diffimg,
                                     product_frame_diffimg_stats);
    if (retval != 0) {
        vimos_ima_bias_tidy(ps,0);
        return(-1);
    }

    /* Free some stuff up */

    vimos_ima_bias_tidy(ps,1);
    return(0);
}

/*---------------------------------------------------------------------------*/
/**
  @brief    Initialise the pointers in the memory structure
  @param    ps           The memory structure defined in the main routine
  @return   Nothing
 */
/*---------------------------------------------------------------------------*/

static void vimos_ima_bias_init(memstruct *ps) {
    ps->labels = NULL;
    ps->labels2 = NULL;
    ps->biaslist = NULL;
    ps->biases = NULL;
    ps->nbiases = 0;
    ps->ref_bias = NULL;
    ps->outimage = NULL;
    ps->drs = NULL;
    ps->rbimage = NULL;
    ps->diffimg = NULL;
    ps->diffimstats = NULL;
}

/*---------------------------------------------------------------------------*/
/**
  @brief    Free any allocated workspace in the memory structure
  @param    ps           The memory structure defined in the main routine
  @param    level        The release level. 0 == everything, 1 == just
                         what is needed for the current detector
 */
/*---------------------------------------------------------------------------*/

static void vimos_ima_bias_tidy(memstruct *ps, int level) {

    freeimage(ps->outimage);
    freefitslist(ps->biases,ps->nbiases);
    freepropertylist(ps->drs);
    freefits(ps->rbimage);
    freeimage(ps->diffimg);
    freetable(ps->diffimstats);
    if (level == 1)
        return;
    freespace(ps->labels);
    freespace(ps->labels2);
    freeframeset(ps->biaslist);
    freeframe(ps->ref_bias);
}

/**@}*/

/*

$Log: vimos_ima_bias.c,v $
Revision 1.9  2015/09/11 09:48:16  jim
changed declaration of some parameters to ranges where appropriate

Revision 1.8  2015/08/07 13:07:15  jim
Fixed copyright to ESO

Revision 1.7  2015/05/01 09:02:22  jim
Modified to pass detector name information differently so that Mac compiler
doesn't gripe

Revision 1.6  2015/01/29 12:10:22  jim
Modified comments, some command line arguments are now strings where it
makes sense to do so, some support routines moved to vimos_utils.c

Revision 1.5  2014/12/11 12:27:28  jim
new version

Revision 1.4  2013/11/21 09:39:08  jim
detabbed

Revision 1.3  2013/10/25 09:31:22  jim
Fixed docs

Revision 1.2  2013-10-24 09:16:28  jim
Changes made to remove compiler warnings

Revision 1.1.1.1  2013-10-24 08:36:12  jim
Initial import


*/