/* $Id: vimos_ima_dark.c,v 1.12 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.12 $
 * $Name:  $
 */

/* Includes */

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

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

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

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

#define MEANDARK    1
#define DIFFIMG     2
#define STATS_TAB   4

/* Static global variables */

typedef struct {

    /* Input */

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

    /* Output */

    float       particle_rate;
    float       darkmed;
    float       darkrms;
    float       darkdiff_med;
    float       darkdiff_rms;
    int         nhot;
    float       hotfrac;

    /* 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      *darklist;
    casu_fits         **darks;
    int               ndarks;
    cpl_frame         *master_bias;
    cpl_frame         *ref_dark;
    casu_mask         *master_mask;
    cpl_image         *outimage;
    cpl_propertylist  *drs;
    casu_fits         *rdimage;
    cpl_image         *diffimg;
    cpl_table         *diffimstats;
    unsigned char     *rejmask;
    unsigned char     *rejplus;
    casu_fits         *mbias;
} memstruct;

/* Function prototypes */

static int vimos_ima_dark_create(cpl_plugin *);
static int vimos_ima_dark_exec(cpl_plugin *);
static int vimos_ima_dark_destroy(cpl_plugin *);
static int vimos_ima_dark(cpl_parameterlist *, cpl_frameset *) ;
static int vimos_ima_dark_save(cpl_frameset *framelist, 
                                   cpl_parameterlist *parlist, char *cname,
                                   configstruct *cs, memstruct *ps, int isfirst,
                                   cpl_frame **product_frame_mean_dark,
                                   cpl_frame **product_frame_diffimg,
                                   cpl_frame **product_frame_diffimg_stats);
static void vimos_ima_dark_dummy_products(configstruct *cs, memstruct *ps);
static void vimos_ima_dark_hotpix(configstruct *cs, memstruct *ps);
static void vimos_ima_dark_normal(int jext, float exptime, 
                                      configstruct *cs, memstruct *ps);
static int vimos_ima_dark_lastbit(cpl_frameset *framelist,
                                      cpl_parameterlist *parlist,
                                      char *chipname, configstruct *cs,
                                      memstruct *ps, int isfirst,
                                      cpl_frame **product_frame_mean_dark,
                                      cpl_frame **product_frame_diffimg,
                                      cpl_frame **product_frame_diffimg_stats);
static void vimos_ima_dark_init(memstruct *ps);
static void vimos_ima_dark_tidy(memstruct *ps, int level);

static char vimos_ima_dark_description[] =
"vimos_ima_dark -- VIMOS dark combine recipe.\n\n"
"Combine a list of dark frames into a mean dark frame. Optionally compare \n"
"the output frame to a reference dark frame\n\n"
"The program accepts the following files in the SOF:\n\n"
"    Tag                   Description\n"
"    -----------------------------------------------------------------------\n"
"    %-21s A list of raw dark images\n"
"    %-21s A master bias frame\n"
"    %-21s Optional reference dark frame\n"
"    %-21s Optional master bad pixel map or\n"
"    %-21s Optional master confidence map\n"
"If no reference dark frame is made available, then no comparison will be \n"
"done. This means there will be no output difference image and no stats table"
"If neither a master bad pixel map nor confidence map is provided then the"
"stats will be done assuming all pixels are good"
"\n";

/**@{*/

/**

    \ingroup recipelist
    \defgroup vimos_ima_dark vimos_ima_dark
    \brief Combine a list of darks to form a mean dark frame

    \par Name: 
        vimos_ima_dark
    \par Purpose: 
        Combine a list of darks to form a mean dark frame. 
    \par Description: 
        A list of darks exposures is combined with rejection to form a mean 
        dark. If a reference dark 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 dark current with time. The difference 
        image is broken into lots of little cells. 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 scaletype (string): Determines how the input data are scaled or 
          offset before they are combined. Can take the following values:
            - none: No scaling of offsetting
            - additive: All input frames are biassed additively to bring their
              backgrounds to a common mean.
            - multiplicative: All input frames are scaled multiplicatively to 
              bring their backgrounds to a common mean.
            - exptime: All input frames are scaled to a uniform exposure time 
              and then additively corrected to bring their backgrounds to a 
              common mean.
        - \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 DARK (required): A list of raw dark images for combining
        - \b MASTER_BIAS (required): A master bias frame
        - \b REFERENCE_DARK (optional): A library reference dark frame.
             If this is given, then a difference image will be formed and some 
             basic statistics run on it.
        - \b MASTER_BPM or \b MASTER_CONF (optional): If either is given, then
             it will be used to mask out bad pixels during any statistical 
             analysis that is done.
    \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 dark frame, formed by either a mean or median
          combination of the input raw frames with rejection (\b MASTER_DARK).
        - The output difference image, formed by subtracting the input reference
          dark frame from the current mean/median dark frame. This is only 
          done if a reference dark frame is specified from the outset
          (\b DIFFIMG_DARK).
        - The output difference image statistics table. 
          (\b DIFFIMG_STATS_DARK).
    \par Output QC Parameters:
        - \b DARKMED
            The median value of the combined dark frame.
        - \b DARKRMS
            The RMS value of the mean dark frame.
        - \b DARKDIFF_MED
            The median of the dark difference image
        - \b DARKDIFF_RMS
            The RMS of the dark difference image
        - \b PARTICLE_RATE
            The average number of cosmic ray hits per pixel per second 
            per frame.
        - \b NHOTPIX
            The number of hot pixels on the mean frame
        - \b HOTFRAC
            The fraction of pixels on the mean frame that are hot
    \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 dark frames in the input frameset
        - Unable to save data products
        - Image extensions can't be loaded
    \par Non-Fatal Error Conditions:
        - Exposure time missing from header. Assumed to be 1 second.
        - No reference dark. No difference image formed.
        - No reference bad pixel or confidence map. All pixels considered to be 
          good.
    \par Conditions Leading To Dummy Products:
        - The detector for the current image extension is flagged dead
        - Combination routine failed
        - Reference dark frame image extension is flagged as a dummy
    \par Author:
        Jim Lewis, CASU
    \par Code Reference: 
        vimos_ima_dark.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_dark_description,
                   VIMOS_DARK_RAW,VIMOS_CAL_BIAS,VIMOS_REF_DARK,VIMOS_CAL_BPM,
                   VIMOS_CAL_CONF);

    cpl_plugin_init(plugin,
                    CPL_PLUGIN_API,
                    VIMOS_BINARY_VERSION,
                    CPL_PLUGIN_TYPE_RECIPE,
                    "vimos_ima_dark",
                    "Dark combination for imaging",
                    alldesc,
                    "Jim Lewis",
                    "jrl@ast.cam.ac.uk",
                    vimos_get_license(),
                    vimos_ima_dark_create,
                    vimos_ima_dark_exec,
                    vimos_ima_dark_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_dark_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_dark.combtype",
                               CPL_TYPE_STRING,
                               "Combination algorithm",
                               "vimos.vimos_ima_dark",
                               "median",2,"median","mean");
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"combtype");
    cpl_parameterlist_append(recipe->parameters,p);

    /* The requested scaling */

    p = cpl_parameter_new_enum("vimos.vimos_ima_dark.scaletype",
                               CPL_TYPE_STRING,
                               "Scaling algorithm",
                               "vimos.vimos_ima_dark",
                               "exptime",4,"none","additive",
                               "multiplicative","exptime");
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"scaletype");
    cpl_parameterlist_append(recipe->parameters,p);
    
    /* Extra rejection cycle */

    p = cpl_parameter_new_value("vimos.vimos_ima_dark.xrej",
                                CPL_TYPE_BOOL,
                                "True if using extra rejection cycle",
                                "vimos.vimos_ima_dark",
                                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_dark.thresh",
                                CPL_TYPE_DOUBLE,
                                "Rejection threshold in sigma above background",
                                "vimos.vimos_ima_dark",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_dark.ncells",
                               CPL_TYPE_INT,
                               "Number of cells for diff image stats",
                               "vimos.vimos_ima_dark",8,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_dark.prettynames",
                                CPL_TYPE_BOOL,
                                "Use pretty output file names?",
                                "vimos.vimos_ima_dark",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_dark_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_dark(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_dark_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_dark(cpl_parameterlist *parlist, 
                              cpl_frameset *framelist) {
    const char *fctid="vimos_ima_dark";
    const char *fname;
    int j,retval,status,nx,ny,isfirst;
    cpl_size nlab,nlab2;
    long i;
    float exptime;
    cpl_parameter *p;
    cpl_propertylist *plist;
    cpl_frameset *chiplist;
    cpl_frame *fr;
    cpl_image *junk;
    configstruct cs;
    memstruct ps;
    cpl_frame *product_frame_mean_dark,*product_frame_diffimg;
    cpl_frame *product_frame_diffimg_stats,*testfrm;
    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_dark_init(&ps);
    cs.we_expect |= MEANDARK;

    /* Get the parameters */

    p = cpl_parameterlist_find(parlist,"vimos.vimos_ima_dark.combtype");
    cs.combtype = (strcmp(cpl_parameter_get_string(p),"median") ? 2 : 1);
    p = cpl_parameterlist_find(parlist,"vimos.vimos_ima_dark.scaletype");
    if (! strcmp(cpl_parameter_get_string(p),"none")) 
        cs.scaletype = 0;
    else if (! strcmp(cpl_parameter_get_string(p),"additive")) 
        cs.scaletype = 1;
    else if (! strcmp(cpl_parameter_get_string(p),"multiplicative")) 
        cs.scaletype = 2;
    else if (! strcmp(cpl_parameter_get_string(p),"exptime")) 
        cs.scaletype = 3;
    p = cpl_parameterlist_find(parlist,"vimos.vimos_ima_dark.xrej");
    cs.xrej = cpl_parameter_get_bool(p);
    p = cpl_parameterlist_find(parlist,"vimos.vimos_ima_dark.thresh");
    cs.thresh = (float)cpl_parameter_get_double(p);
    p = cpl_parameterlist_find(parlist,"vimos.vimos_ima_dark.ncells");
    cs.ncells = cpl_parameter_get_int(p);
    p = cpl_parameterlist_find(parlist,"vimos.vimos_ima_dark.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_dark_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_dark_tidy(&ps,0);
        return(-1);
    }

    /* Get the dark frames */

    if ((ps.darklist = casu_frameset_subgroup(framelist,ps.labels,nlab,
                                              VIMOS_DARK_RAW)) == NULL) {
        cpl_msg_error(fctid,"Cannot find dark frames in input frameset");
        vimos_ima_dark_tidy(&ps,0);
        return(-1);
    }
    ps.ndarks = cpl_frameset_get_size(ps.darklist);

    /* Get the exposure time from the first dark frame in the list */

    fr = cpl_frameset_get_position(ps.darklist,0);
    fname = cpl_frame_get_filename(fr);
    plist = cpl_propertylist_load(fname,cpl_frame_get_nextensions(fr));
    if (vimos_pfits_get_exptime(plist,&exptime) != CASU_OK) {
        cpl_msg_warning(fctid,"Unable to get exposure time for %s",fname);
        exptime = 1.0;
    }
    cpl_propertylist_delete(plist);

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

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

    /* See what chip set this is */

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

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

    if ((ps.master_bias = casu_frameset_subgroup_1(framelist,ps.labels,nlab,
                                                   VIMOS_CAL_BIAS)) == NULL) {
        cpl_msg_error(fctid,"No master bias found");
        vimos_ima_dark_tidy(&ps,0);
        return(-1);
    }
        
    /* Check to see if there is a reference dark frame */

    if ((ps.ref_dark = casu_frameset_subgroup_1(framelist,ps.labels,nlab,
                                                VIMOS_REF_DARK)) == NULL) {
        cpl_msg_info(fctid,"No master dark found -- no difference image will be formed");
    } else {
        cs.we_expect |= DIFFIMG;
        cs.we_expect |= STATS_TAB;
    }
        
    /* Check to see if there is a master bad pixel map. If there isn't one 
       then look for a confidence map */

    ps.master_mask = casu_mask_define(framelist,ps.labels,nlab,VIMOS_CAL_CONF,
                                      VIMOS_CAL_BPM);

    /* Now loop for all the extensions... */

    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.darklist,ps.labels2,nlab2,
                                               vimos_names[j-1])) == NULL) {
            cpl_msg_info(fctid,"No darks in input list for detector %s\n",
                         vimos_names[j-1]);
            retval = vimos_ima_dark_lastbit(framelist,parlist,
                                                vimos_names[j-1],&cs,&ps,
                                                isfirst,
                                                &product_frame_mean_dark,
                                                &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.darks = casu_fits_load_list(chiplist,CPL_TYPE_FLOAT,
                                       cpl_frame_get_nextensions(testfrm));
        ps.ndarks = cpl_frameset_get_size(chiplist);
        freeframeset(chiplist);
        if (ps.darks == NULL) {
            cpl_msg_info(fctid,"Detector %s darks wouldn't load",
                         vimos_names[j-1]);
            retval = vimos_ima_dark_lastbit(framelist,parlist,
                                                vimos_names[j-1],&cs,&ps,
                                                isfirst,
                                                &product_frame_mean_dark,
                                                &product_frame_diffimg,
                                                &product_frame_diffimg_stats);
            if (retval != 0) 
                return(-1);
            continue;
        }

        /* Right, we want to bias correct, so we need to load the mean 
           bias and make sure it isn't a dummy */

        ps.mbias = casu_fits_load(ps.master_bias,CPL_TYPE_FLOAT,j);
        if (ps.mbias == NULL) {
            cpl_msg_info(fctid,
                         "Can't load master bias for detector %s",
                         vimos_names[j-1]);
            retval = vimos_ima_dark_lastbit(framelist,parlist,
                                                vimos_names[j-1],&cs,&ps,
                                                isfirst,
                                                &product_frame_mean_dark,
                                                &product_frame_diffimg,
                                                &product_frame_diffimg_stats);
            if (retval != 0)
                return(-1);
            continue;
        } else if (casu_is_dummy(casu_fits_get_ehu(ps.mbias))) {
            cpl_msg_info(fctid,
                         "Master bias for detector %s is a dummy",
                         vimos_names[j-1]);
            retval = vimos_ima_dark_lastbit(framelist,parlist,
                                                vimos_names[j-1],&cs,&ps,
                                                isfirst,
                                                &product_frame_mean_dark,
                                                &product_frame_diffimg,
                                                &product_frame_diffimg_stats);
            if (retval != 0)
                return(-1);
            continue;
        }

        /* Load the master mask extension */

        nx = (int)cpl_image_get_size_x(casu_fits_get_image(ps.darks[0]));
        ny = (int)cpl_image_get_size_y(casu_fits_get_image(ps.darks[0]));
        retval = casu_mask_load(ps.master_mask,j,nx,ny);
        if (retval == CASU_FATAL) {
            cpl_msg_info(fctid,
                         "Unable to load mask image %s[%" CPL_SIZE_FORMAT "]",
                         casu_mask_get_filename(ps.master_mask),(cpl_size)j);
            cpl_msg_info(fctid,"Forcing all pixels to be good from now on");
            casu_mask_force(ps.master_mask,nx,ny);
        }

        /* Loop for each image and dark correct */

        cpl_msg_info(fctid,"Bias correcting detector %s",vimos_names[j-1]);
        for (i = 0; i < ps.ndarks; i++) 
            vimos_biascor((ps.darks)[i],ps.mbias,1,1,&status);
            
        /* 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.darks,NULL,ps.ndarks,cs.combtype,cs.scaletype,
                             cs.xrej,cs.thresh,"EXPTIME",&(ps.outimage),&junk,
                             &(ps.rejmask),&(ps.rejplus),&(ps.drs),&status);
        if (status == CASU_OK) {
            cs.we_get |= MEANDARK;
            vimos_ima_dark_hotpix(&cs,&ps);
            vimos_ima_dark_normal(j,exptime,&cs,&ps);
        } 

        /* Create any dummies and save the products */

        retval = vimos_ima_dark_lastbit(framelist,parlist,
                                            vimos_names[j-1],&cs,&ps,isfirst,
                                            &product_frame_mean_dark,
                                            &product_frame_diffimg,
                                            &product_frame_diffimg_stats);
        if (retval != 0) 
            return(-1);
    }
    vimos_ima_dark_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_dark   The frame for the mean dark 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_dark_save(cpl_frameset *framelist, 
                                   cpl_parameterlist *parlist, char *cname,
                                   configstruct *cs, memstruct *ps, int isfirst,
                                   cpl_frame **product_frame_mean_dark,
                                   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_dark_save";
    const char *esoout[] = {"darkcomb.fits","darkdiff.fits","darkdifftab.fits"};
    const char *prettypfx[] = {"dark","darkdiff","darkdifftab"};
    char outfile[3][BUFSIZ],dateobs[81];
    const char *recipeid = "vimos_ima_dark";

    /* 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->darks[0]);
        elist = casu_fits_get_ehu(ps->darks[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_dark = cpl_frame_new();
        cpl_frame_set_filename(*product_frame_mean_dark,outfile[0]);
        cpl_frame_set_tag(*product_frame_mean_dark,VIMOS_PRO_DARK);
        cpl_frame_set_type(*product_frame_mean_dark,CPL_FRAME_TYPE_IMAGE);
        cpl_frame_set_group(*product_frame_mean_dark,CPL_FRAME_GROUP_PRODUCT);
        cpl_frame_set_level(*product_frame_mean_dark,CPL_FRAME_LEVEL_FINAL);

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

        plist = cpl_propertylist_duplicate(casu_fits_get_ehu(ps->darks[0]));
        vimos_dfs_set_product_primary_header(plist,*product_frame_mean_dark,
                                             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_dark);
            return(-1);
        }
        cpl_frameset_insert(framelist,*product_frame_mean_dark);
        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_DARK);
            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->darks[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_DARK_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->darks[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->darks[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(p);
            cpl_propertylist_delete(plist);
            cpl_frameset_insert(framelist,*product_frame_diffimg_stats);
        }
    }

    /* Get the extension property list */

    plist = casu_fits_get_ehu(ps->darks[0]);
    cpl_propertylist_update_int(plist,"ESO PRO DATANCOM",ps->ndarks);

    /* Fiddle with the header now */

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

    cpl_propertylist_update_float(p,"ESO QC DARKMED",cs->darkmed);
    cpl_propertylist_set_comment(p,"ESO QC DARKMED",
                                 "Median of mean dark frame");
    cpl_propertylist_update_float(p,"ESO QC DARKRMS",cs->darkrms);
    cpl_propertylist_set_comment(p,"ESO QC DARKRMS","RMS of mean dark frame");
    cpl_propertylist_update_float(p,"ESO QC PARTICLE_RATE",cs->particle_rate);
    cpl_propertylist_set_comment(p,"ESO QC PARTICLE_RATE",
                                 "[N/(detector*sec)] Particle rate");
    cpl_propertylist_update_int(p,"ESO QC NHOTPIX",cs->nhot);
    cpl_propertylist_set_comment(p,"ESO QC NHOTPIX","Number of hot pixels");
    cpl_propertylist_update_float(p,"ESO QC HOTFRAC",cs->hotfrac);
    cpl_propertylist_set_comment(p,"ESO QC HOTFRAC","Hot pixel fraction");
    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 dark 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 DARKDIFF_MED",cs->darkdiff_med);
        cpl_propertylist_set_comment(p,"ESO QC DARKDIFF_MED",
                                     "Median of dark difference image");
        cpl_propertylist_update_float(p,"ESO QC DARKDIFF_RMS",cs->darkdiff_rms);
        cpl_propertylist_set_comment(p,"ESO QC DARKDIFF_RMS",
                                     "RMS of dark 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);
    }

    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_dark_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 & MEANDARK)) {
        ps->outimage = casu_dummy_image(ps->darks[0]);

        /* Set up the QC parameters */
    
        cs->particle_rate = 0.0;
        cs->darkmed = 0.0;
        cs->darkrms = 0.0;
        cs->nhot = 0;
        cs->hotfrac = 0.0;
    }

    /* Do the difference image */

    if ((cs->we_expect & DIFFIMG) && ! (cs->we_get & DIFFIMG)) {
        cs->darkdiff_med = 0.0;
        cs->darkdiff_rms = 0.0;
        ps->diffimg = casu_dummy_image(ps->darks[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    Work out the number of hot pixels
  @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_dark_hotpix(configstruct *cs, memstruct *ps) {
    int i,nx,ny,nh,nhot,j;
    long npts;
    cpl_image *im;
    unsigned char *bpm;
    float med,mad,lowcut,highcut,*data;
    casu_fits *f;

    /* Get some workspace to hold the bad pixel mask */

    im = casu_fits_get_image(ps->darks[0]);
    nx = (int)cpl_image_get_size_x(im);
    ny = (int)cpl_image_get_size_y(im);
    npts = (long)(nx*ny);
    bpm = cpl_calloc(npts,sizeof(*bpm));

    /* Create a difference image for each of the good frames */

    for (i = 0; i < ps->ndarks; i++) {
        f = casu_fits_duplicate(ps->darks[i]);
        im = casu_fits_get_image(f);
        cpl_image_subtract(im,ps->outimage);
        
        /* Work out the stats of the difference image. Define a lower and
           upper cut. NB: a lower cut is needed since we are doing stats
           on a difference image and hot pixels will appear as either 
           bright or dark. Dead pixels will probably correct properly and
           hence shouldn't be flagged using this procedure. */

        data = cpl_image_get_data_float(im);
        casu_medmad(data,NULL,npts,&med,&mad);
        lowcut = med - 1.48*mad*cs->thresh;
        highcut = med + 1.48*mad*cs->thresh;
        for (j = 0; j < npts; j++) 
            if (data[j] > highcut || data[j] < lowcut)
                bpm[j] += 1;

        /* Get rid of temporary image */

        casu_fits_delete(f);
    }

    /* Define a pixel as hot so long as it is discordant on at least half of 
       the frames */

    nh = (ps->ndarks + 1)/2;
    nhot = 0;
    for (j = 0; j < npts; j++)
        if (bpm[j] >= nh)
            nhot++;

    /* Clean up... */

    cpl_free(bpm);
        
    /* Set QC parameters */

    cs->nhot = nhot;
    cs->hotfrac = (float)nhot/(float)npts;
}


/*---------------------------------------------------------------------------*/
/**
  @brief    Normalise the dark frame and create diff image and stats table
  @param    jext         the extension number
  @param    exptime      the exposure time
  @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_dark_normal(int jext, float exptime,
                                      configstruct *cs, memstruct *ps) {
    int nx,ny,ndiff,ncells;
    long npi,i;
    unsigned char *bpm;
    float med,sig,*idata,grms,gdiff;
    const char *fctid="vimos_ima_dark_normal";

    /* Load up the bad pixel mask */

    nx = (int)cpl_image_get_size_x(ps->outimage);
    ny = (int)cpl_image_get_size_y(ps->outimage);
    npi = nx*ny;
    cs->particle_rate = 0;
    bpm = casu_mask_get_data(ps->master_mask);

    /* Now find out how many 'good' pixels were rejected for
       being too high during the combination phase */

    ndiff = 0;
    for (i = 0; i < npi; i++)
        if ((ps->rejplus)[i] > 0 && bpm[i] == 0)
            ndiff += (ps->rejplus)[i];
    cs->particle_rate = (float)ndiff/(exptime*(float)(ps->ndarks));

    /* Divide by the exposure time and work out the RMS of the mean 
       dark frame */

    cpl_image_divide_scalar(ps->outimage,(double)exptime);
    idata = cpl_image_get_data(ps->outimage);
    casu_medmad(idata,bpm,npi,&med,&sig);
    sig *= 1.48;
    cs->darkmed = med;
    cs->darkrms = sig;

    /* Load up the master dark */

    if (ps->ref_dark != NULL) {
        ps->rdimage = casu_fits_load(ps->ref_dark,CPL_TYPE_FLOAT,jext);
        if (ps->rdimage == NULL) 
            cpl_msg_info(fctid,
                         "Master dark extension %" CPL_SIZE_FORMAT " won't load",
                         (cpl_size)jext);
        else if (casu_is_dummy(casu_fits_get_ehu(ps->rdimage))) {
            cpl_msg_info(fctid,
                         "Master dark extension %" CPL_SIZE_FORMAT " is a dummy!",
                         (cpl_size)jext);
            freefits(ps->rdimage);
        }
    } else 
        ps->rdimage = NULL;

    /* Form the difference image. NB: the difference image routine
       copes if the input mean image is null and will return nulls */

    cs->darkdiff_med = 0.0;
    cs->darkdiff_rms = 0.0;
    ncells = cs->ncells;
    vimos_difference_image(casu_fits_get_image(ps->rdimage),ps->outimage,bpm,
                           ncells,1,&gdiff,&grms,&(ps->diffimg),
                           &(ps->diffimstats));
    casu_mask_clear(ps->master_mask);
    cs->darkdiff_med = gdiff;
    cs->darkdiff_rms = grms;
    if (ps->diffimg != NULL)
        cs->we_get |= DIFFIMG;
    if (ps->diffimstats != NULL)
        cs->we_get |= STATS_TAB;
    return;
}

/*---------------------------------------------------------------------------*/
/**
  @brief    Make dummys and save
  @param    framelist    the input frame list
  @param    parlist      the input recipe parameter list
  @param    chipname     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_dark   The frame for the mean dark 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_dark_lastbit(cpl_frameset *framelist,
                                      cpl_parameterlist *parlist,
                                      char *chipname, configstruct *cs,
                                      memstruct *ps, int isfirst,
                                      cpl_frame **product_frame_mean_dark,
                                      cpl_frame **product_frame_diffimg,
                                      cpl_frame **product_frame_diffimg_stats) {
    int retval;
    const char *fctid="vimos_ima_dark_lastbit";

    /* Make whatever dummy products you need */

    vimos_ima_dark_dummy_products(cs,ps);

    /* Save everything */

    cpl_msg_info(fctid,"Saving products for detector %s",chipname);
    retval = vimos_ima_dark_save(framelist,parlist,chipname,cs,ps,isfirst,
                                     product_frame_mean_dark, 
                                     product_frame_diffimg, 
                                     product_frame_diffimg_stats);
    if (retval != 0) {
        vimos_ima_dark_tidy(ps,2);
        return(-1);
    }

    /* Free some stuff up */

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

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

static void vimos_ima_dark_init(memstruct *ps) {
    ps->labels = NULL;
    ps->labels2 = NULL;
    ps->darklist = NULL;
    ps->darks = NULL;
    ps->ndarks = 0;
    ps->ref_dark = NULL;
    ps->master_mask = NULL;
    ps->master_bias = NULL;
    ps->outimage = NULL;
    ps->drs = NULL;
    ps->rdimage = NULL;
    ps->diffimg = NULL;
    ps->diffimstats = NULL;
    ps->rejmask = NULL;
    ps->rejplus = NULL;
    ps->mbias = 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_dark_tidy(memstruct *ps, int level) {

    freeimage(ps->outimage);
    freefitslist(ps->darks,ps->ndarks);
    freepropertylist(ps->drs);
    freefits(ps->rdimage);
    freeimage(ps->diffimg);
    freetable(ps->diffimstats);
    freespace(ps->rejplus);
    freespace(ps->rejmask);
    freefits(ps->mbias);
    if (level == 1)
        return;
    freespace(ps->labels);
    freespace(ps->labels2)
    freeframeset(ps->darklist);
    freeframe(ps->ref_dark);
    freemask(ps->master_mask);
    freeframe(ps->master_bias);
}

/**@}*/

/*

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

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

Revision 1.10  2015/06/25 13:08:25  jim
Fixed bug in propertylist loading

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

Revision 1.8  2015/04/30 12:23:35  jim
Fixed manpage info

Revision 1.7  2015/01/29 12:10:23  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.6  2014/12/11 12:27:28  jim
new version

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

Revision 1.4  2013/11/04 05:55:51  jim
Forgot to divide by exposure time!

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


*/