vircam_dark_combine.c

/* $Id: vircam_dark_combine.c,v 1.57 2007/10/19 09:25:09 jim Exp $
 *
 * This file is part of the VIRCAM Pipeline
 * Copyright (C) 2005 Cambridge Astronomy Survey Unit
 *
 * 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: 2007/10/19 09:25:09 $
 * $Revision: 1.57 $
 * $Name:  $
 */

/* Includes */

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

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

#include "vircam_utils.h"
#include "vircam_pfits.h"
#include "vircam_dfs.h"
#include "vircam_mods.h"
#include "vircam_stats.h"
#include "vircam_fits.h"
#include "vircam_mask.h"
#include "vircam_channel.h"
#include "vircam_dfs.h"
#include "vircam_paf.h"
#include "vircam_wcsutils.h"

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

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

/* Function prototypes */

static int vircam_dark_combine_create(cpl_plugin *) ;
static int vircam_dark_combine_exec(cpl_plugin *) ;
static int vircam_dark_combine_destroy(cpl_plugin *) ;
static int vircam_dark_combine(cpl_parameterlist *, cpl_frameset *) ;
static int vircam_dark_combine_save(cpl_frameset *framelist, 
                                    cpl_parameterlist *parlist);
static void vircam_dark_combine_dummy_products(void);
static void vircam_dark_combine_hotpix(void);
static void vircam_dark_combine_normal(int jext, float exptime);
static int vircam_dark_combine_lastbit(int jext, cpl_frameset *framelist,
                                       cpl_parameterlist *parlist);
static void vircam_dark_combine_init(void);
static void vircam_dark_combine_tidy(int level);

/* Static global variables */

static struct {

    /* Input */

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

    /* Output */

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

} vircam_dark_combine_config;


static struct {
    int               *labels;
    cpl_frameset      *darklist;
    vir_fits          **darks;
    int               ndarks;
    vir_fits          **good;
    int               ngood;
    cpl_frame         *master_dark;
    vir_mask          *master_mask;
    cpl_frame         *chantab;
    cpl_image         *outimage;
    cpl_propertylist  *drs;
    unsigned char     *rejmask;
    unsigned char     *rejplus;
    vir_fits          *mdimage;
    cpl_image         *diffimg;
    cpl_table         *diffimstats;
    cpl_propertylist  *phupaf;
} ps;

static cpl_frame *product_frame_mean_dark = NULL;
static cpl_frame *product_frame_diffimg = NULL;
static cpl_frame *product_frame_diffimg_stats = NULL;
static int isfirst;
static int we_expect;
static int we_get;

static char vircam_dark_combine_description[] =
"vircam_dark_combine -- VIRCAM dark combine recipe.\n\n"
"Combine a list of dark frames into a mean dark frame. Optionally compare \n"
"the output frame to a master 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 Optional reference dark frame\n"
"    %-21s Optional master bad pixel map or\n"
"    %-21s Optional master confidence map\n"
"    %-21s Optional channel table or\n"
"    %-21s Optional initial channel table\n"
"If no master dark frame is made available, then no comparison will be done\n"
"This means there will be no output difference image. If a master dark is\n"
"available, but no channel table is, then a difference image will be formed\n"
"but no stats will be written."
"\n";

/* Function code */


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

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,vircam_dark_combine_description,
                   VIRCAM_DARK_RAW,VIRCAM_REF_DARK,VIRCAM_CAL_BPM,
                   VIRCAM_CAL_CONF,VIRCAM_CAL_CHANTAB,VIRCAM_CAL_CHANTAB_INIT);

    cpl_plugin_init(plugin,
                    CPL_PLUGIN_API,
                    VIRCAM_BINARY_VERSION,
                    CPL_PLUGIN_TYPE_RECIPE,
                    "vircam_dark_combine",
                    "VIRCAM dark combination recipe",
                    alldesc,
                    "Jim Lewis",
                    "jrl@ast.cam.ac.uk",
                    vircam_get_license(),
                    vircam_dark_combine_create,
                    vircam_dark_combine_exec,
                    vircam_dark_combine_destroy);

    cpl_pluginlist_append(list,plugin);

    return(0);
}

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

static int vircam_dark_combine_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_range("vircam.vircam_dark_combine.combtype",
                                CPL_TYPE_INT,
                                "1 == Median,\n 2 == Mean",
                                "vircam.vircam_dark_combine",
                                1,1,2);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"combtype");
    cpl_parameterlist_append(recipe->parameters,p);

    /* The requested scaling */

    p = cpl_parameter_new_range("vircam.vircam_dark_combine.scaletype",
                                CPL_TYPE_INT,
                                "0 == none,\n 1 == additive offset,\n 2 == multiplicative offset,\n 3 == exposure time scaling + additive offset",
                                "vircam.vircam_dark_combine",
                                1,0,3);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"scaletype");
    cpl_parameterlist_append(recipe->parameters,p);
    
    /* Extra rejection cycle */

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

    /* Rejection threshold */

    p = cpl_parameter_new_value("vircam.vircam_dark_combine.thresh",
                                CPL_TYPE_DOUBLE,
                                "Rejection threshold in sigma above background",
                                "vircam.vircam_dark_combine",5.0);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"thresh");
    cpl_parameterlist_append(recipe->parameters,p);

    /* How many cells to divide each data channel */

    p = cpl_parameter_new_enum("vircam.vircam_dark_combine.ncells",
                               CPL_TYPE_INT,
                               "Number of cells for data channel stats",
                               "vircam.vircam_dark_combine",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);     

    /* Extension number of input frames to use */

    p = cpl_parameter_new_range("vircam.vircam_dark_combine.extenum",
                                CPL_TYPE_INT,
                                "Extension number to be done, 0 == all",
                                "vircam.vircam_dark_combine",
                                1,0,16);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"ext");
    cpl_parameterlist_append(recipe->parameters,p);
        
    /* Get out of here */

    return(0);
}
    
    
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

static int vircam_dark_combine_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(vircam_dark_combine(recipe->parameters,recipe->frames));
}
                                
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

static int vircam_dark_combine_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);
}

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

static int vircam_dark_combine(cpl_parameterlist *parlist, 
                               cpl_frameset *framelist) {
    const char *fctid="vircam_dark_combine";
    const char *fname;
    int nlab,j,jst,jfn,retval,status,live,nx,ny;
    long i;
    float exptime;
    vir_fits *ff;
    cpl_parameter *p;
    cpl_propertylist *plist;

    /* 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 */

    vircam_dark_combine_init();
    we_expect |= MEANDARK;

    /* Get the parameters */

    p = cpl_parameterlist_find(parlist,"vircam.vircam_dark_combine.combtype");
    vircam_dark_combine_config.combtype = cpl_parameter_get_int(p);
    p = cpl_parameterlist_find(parlist,"vircam.vircam_dark_combine.scaletype");
    vircam_dark_combine_config.scaletype = cpl_parameter_get_int(p);
    p = cpl_parameterlist_find(parlist,"vircam.vircam_dark_combine.xrej");
    vircam_dark_combine_config.xrej = cpl_parameter_get_bool(p);
    p = cpl_parameterlist_find(parlist,"vircam.vircam_dark_combine.thresh");
    vircam_dark_combine_config.thresh = (float)cpl_parameter_get_double(p);
    p = cpl_parameterlist_find(parlist,"vircam.vircam_dark_combine.ncells");
    vircam_dark_combine_config.ncells = cpl_parameter_get_int(p);
    p = cpl_parameterlist_find(parlist,"vircam.vircam_dark_combine.extenum");
    vircam_dark_combine_config.extenum = cpl_parameter_get_int(p);

    /* Sort out raw from calib frames */

    if (vircam_dfs_set_groups(framelist) != VIR_OK) {
        cpl_msg_error(fctid,"Cannot identify RAW and CALIB frames");
        vircam_dark_combine_tidy(2);
        return(-1);
    }

    /* Get a list of the frame labels */ 

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

    /* Get the dark frames */

    if ((ps.darklist = vircam_frameset_subgroup(framelist,ps.labels,nlab,
                                                VIRCAM_DARK_RAW)) == NULL) {
        cpl_msg_error(fctid,"Cannot find dark frames in input frameset");
        vircam_dark_combine_tidy(2);
        return(-1);
    }
    ps.ndarks = cpl_frameset_get_size(ps.darklist);

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

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

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

    if ((ps.master_dark = vircam_frameset_subgroup_1(framelist,ps.labels,nlab,
                                                     VIRCAM_REF_DARK)) == NULL)
        cpl_msg_info(fctid,"No master dark found -- no difference image will be formed");
    else
        we_expect |= DIFFIMG;
        
    /* 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 = vircam_mask_define(framelist,ps.labels,nlab);

    /* Check to see if there is a channel table */

    if ((ps.chantab = vircam_frameset_subgroup_1(framelist,ps.labels,nlab,
                                                 VIRCAM_CAL_CHANTAB)) == NULL) {
        if ((ps.chantab = vircam_frameset_subgroup_1(framelist,ps.labels,nlab,
                                                     VIRCAM_CAL_CHANTAB_INIT)) == NULL)
            cpl_msg_info(fctid,"No channel table found -- no difference image stats will be done");
    } else if (we_expect & DIFFIMG)
        we_expect |= STATS_TAB;

    /* Now, how many image extensions do we want to do? If the extension
       number is zero, then we loop for all possible extensions. If it
       isn't then we just do the extension specified */

    vircam_exten_range(vircam_dark_combine_config.extenum,
                       (const cpl_frame *)cpl_frameset_get_frame(ps.darklist,0),
                       &jst,&jfn);
    if (jst == -1 || jfn == -1) {
        cpl_msg_error(fctid,"Unable to continue");
        vircam_dark_combine_tidy(2);
        return(-1);
    }

    /* Get some space for the good frames */

    ps.good = cpl_malloc(ps.ndarks*sizeof(vir_fits *));

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

    for (j = jst; j <= jfn; j++) {
        status = VIR_OK;
        we_get = 0;
        isfirst = (j == jst);

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

        ps.darks = vircam_fits_load_list(ps.darklist,CPL_TYPE_FLOAT,j);
        if (ps.darks == NULL) {
            cpl_msg_info(fctid,"Extension %d darks wouldn't load",j);
            retval = vircam_dark_combine_lastbit(j,framelist,parlist);
            if (retval != 0) 
                return(-1);
            continue;
        }

        /* Are any of these dark frames good? */

        ps.ngood = 0;
        for (i = 0; i < ps.ndarks; i++) {
            ff = ps.darks[i];
            vircam_pfits_get_detlive(vircam_fits_get_ehu(ff),&live);
            if (! live) {
                cpl_msg_info(fctid,"Detector flagged dead %s",
                             vircam_fits_get_fullname(ff));
                vircam_fits_set_error(ff,VIR_FATAL);
            } else {
                ps.good[ps.ngood] = ff;
                ps.ngood += 1;
            }
        }       

        /* If there are no good images, then signal that we need to 
           create some dummy products and move on */

        if (ps.ngood == 0) {
            cpl_msg_info(fctid,"All images flagged bad for this extension");
            retval = vircam_dark_combine_lastbit(j,framelist,parlist);
            if (retval != 0) 
                return(-1);
            continue;
        }
        
        /* Load the master mask extension */

        nx = cpl_image_get_size_x(vircam_fits_get_image(ps.good[0]));
        ny = cpl_image_get_size_y(vircam_fits_get_image(ps.good[0]));
        retval = vircam_mask_load(ps.master_mask,j,nx,ny);
        if (retval == VIR_FATAL) {
            cpl_msg_info(fctid,"Unable to load mask image %s[%d]",
                         vircam_mask_get_filename(ps.master_mask),j);
            cpl_msg_info(fctid,"Forcing all pixels to be good from now on");
            vircam_mask_force(ps.master_mask,nx,ny);
        }
            
        /* Call the combine module. If it fails, then signal that
           all products will be dummies */

        cpl_msg_info(fctid,"Doing combination for extension %d",j);
        (void)vircam_imcombine(ps.good,ps.ngood,
                               vircam_dark_combine_config.combtype,
                               vircam_dark_combine_config.scaletype,
                               vircam_dark_combine_config.xrej,
                               vircam_dark_combine_config.thresh,
                               &(ps.outimage),&(ps.rejmask),
                               &(ps.rejplus),&(ps.drs),&status);
        if (status == VIR_OK) {
            we_get |= MEANDARK;
            vircam_dark_combine_hotpix();
            vircam_dark_combine_normal(j,exptime);
        } 

        /* Create any dummies and save the products */

        retval = vircam_dark_combine_lastbit(j,framelist,parlist);
        if (retval != 0) 
            return(-1);
    }
    vircam_dark_combine_tidy(2);
    return(0);
}


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

static int vircam_dark_combine_save(cpl_frameset *framelist, 
                                    cpl_parameterlist *parlist) {
    cpl_propertylist *plist,*elist,*p,*pafprop;
    int status;
    const char *fctid = "vircam_dark_combine_save";
    const char *outfile = "darkcomb.fits";
    const char *outdiff = "darkdiff.fits";
    const char *outdimst = "darkdifftab.fits";
    const char *outfilepaf = "darkcomb";
    const char *outdiffpaf = "darkdiff";
    const char *recipeid = "vircam_dark_combine";

    /* 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);
        cpl_frame_set_tag(product_frame_mean_dark,VIRCAM_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 = vircam_fits_get_phu(ps.darks[0]);
        ps.phupaf = vircam_paf_phu_items(plist);
        vircam_dfs_set_product_primary_header(plist,product_frame_mean_dark,
                                              framelist,parlist,
                                              (char *)recipeid,
                                              "PRO-1.15");

        /* 'Save' the PHU image */                       

        if (cpl_image_save(NULL,outfile,CPL_BPP_8_UNSIGNED,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);

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

        if (we_expect & DIFFIMG) {
            product_frame_diffimg = cpl_frame_new();
            cpl_frame_set_filename(product_frame_diffimg,outdiff);
            cpl_frame_set_tag(product_frame_diffimg,VIRCAM_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 = vircam_fits_get_phu(ps.darks[0]);
            vircam_dfs_set_product_primary_header(plist,product_frame_diffimg,
                                                  framelist,parlist,
                                                  (char *)recipeid,
                                                  "PRO-1.15");

            /* 'Save' the PHU image */                   

            if (cpl_image_save(NULL,outdiff,CPL_BPP_8_UNSIGNED,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);
        }

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

        if (we_expect & STATS_TAB) {
            product_frame_diffimg_stats = cpl_frame_new();
            cpl_frame_set_filename(product_frame_diffimg_stats,outdimst);
            cpl_frame_set_tag(product_frame_diffimg_stats,
                              VIRCAM_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 = vircam_fits_get_phu(ps.darks[0]);
            vircam_dfs_set_product_primary_header(plist,
                                                  product_frame_diffimg_stats,
                                                  framelist,parlist,
                                                  (char *)recipeid,
                                                  "PRO-1.15");

            /* Fiddle with the extension header now */

            elist = vircam_fits_get_ehu(ps.darks[0]);
            p = cpl_propertylist_duplicate(elist);
            vircam_merge_propertylists(p,ps.drs);
            if (! (we_get & STATS_TAB))
                vircam_dummy_property(p);
            vircam_dfs_set_product_exten_header(p,product_frame_diffimg_stats,
                                                framelist,parlist,
                                                (char *)recipeid,
                                                "PRO-1.15");
            status = VIR_OK;
            vircam_removewcs(p,&status);

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

            if (cpl_table_save(ps.diffimstats,plist,p,outdimst,
                               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_frameset_insert(framelist,product_frame_diffimg_stats);
        }
    }

    /* Get the extension property list */

    plist = vircam_fits_get_ehu(ps.darks[0]);

    /* Fiddle with the header now */

    vircam_merge_propertylists(plist,ps.drs);
    p = cpl_propertylist_duplicate(plist);
    if (! (we_get & MEANDARK))
        vircam_dummy_property(p);
    vircam_dfs_set_product_exten_header(p,product_frame_mean_dark,
                                        framelist,parlist,
                                        (char *)recipeid,"PRO-1.15");
                
    /* Now save the mean dark image extension */

    cpl_propertylist_update_float(p,"ESO QC DARKMED",
                                  vircam_dark_combine_config.darkmed);
    cpl_propertylist_set_comment(p,"ESO QC DARKMED",
                                 "Median of mean dark frame");
    cpl_propertylist_update_float(p,"ESO QC DARKRMS",
                                  vircam_dark_combine_config.darkrms);
    cpl_propertylist_set_comment(p,"ESO QC DARKRMS",
                                 "RMS of mean dark frame");
    cpl_propertylist_update_float(p,"ESO QC PARTICLE_RATE",
                                  vircam_dark_combine_config.particle_rate);
    cpl_propertylist_set_comment(p,"ESO QC PARTICLE_RATE",
                                 "[N/(detector*sec)] Particle rate");
    cpl_propertylist_update_float(p,"ESO QC STRIPERMS",
                                  vircam_dark_combine_config.striperms);
    cpl_propertylist_set_comment(p,"ESO QC STRIPERMS","RMS of stripe pattern");
    cpl_propertylist_update_int(p,"ESO QC NHOTPIX",
                                vircam_dark_combine_config.nhot);
    cpl_propertylist_set_comment(p,"ESO QC NHOTPIX","Number of hot pixels");
    cpl_propertylist_update_float(p,"ESO QC HOTFRAC",
                                  vircam_dark_combine_config.hotfrac);
    cpl_propertylist_set_comment(p,"ESO QC HOTFRAC","Hot pixel fraction");
    if (cpl_image_save(ps.outimage,outfile,CPL_BPP_IEEE_FLOAT,p,
                       CPL_IO_EXTEND) != CPL_ERROR_NONE) {
        cpl_msg_error(fctid,"Cannot save product image extension");
        cpl_propertylist_delete(p);
        return(-1);
    }

    /* Write out PAF for mean image */

    pafprop = vircam_paf_req_items(p);
    vircam_merge_propertylists(pafprop,ps.phupaf);
    vircam_paf_append(pafprop,p,"ESO DET NDIT");
    if (vircam_paf_print((char *)outfilepaf,"VIRCAM/vircam_dark_combine",
                         "QC file",pafprop) != VIR_OK)
        cpl_msg_warning(fctid,"Unable to save PAF for mean dark");
    cpl_propertylist_delete(pafprop);
    cpl_propertylist_delete(p);

    /* Now save the dark difference image extension */

    if (we_expect & DIFFIMG) {
        p = cpl_propertylist_duplicate(plist);
        if (! (we_get & DIFFIMG))
            vircam_dummy_property(p);;
        cpl_propertylist_update_float(p,"ESO QC DARKDIFF_MED",
                                      vircam_dark_combine_config.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",
                                      vircam_dark_combine_config.darkdiff_rms);
        cpl_propertylist_set_comment(p,"ESO QC DARKDIFF_RMS",
                                     "RMS of dark difference image");
        vircam_dfs_set_product_exten_header(p,product_frame_diffimg,
                                            framelist,parlist,
                                            (char *)recipeid,
                                            "PRO-1.15");
        if (cpl_image_save(ps.diffimg,outdiff,CPL_BPP_IEEE_FLOAT,p,
                           CPL_IO_EXTEND) != CPL_ERROR_NONE) {
            cpl_msg_error(fctid,"Cannot save product image extension");
            cpl_propertylist_delete(p);
            return(-1);
        }
        
        /* Now write PAF for difference image */

        pafprop = vircam_paf_req_items(p);
        vircam_merge_propertylists(pafprop,ps.phupaf);
        if (vircam_paf_print((char *)outdiffpaf,"VIRCAM/vircam_dark_combine",
                             "QC file",pafprop) != VIR_OK)
            cpl_msg_warning(fctid,"Unable to save PAF for difference image");
        cpl_propertylist_delete(pafprop);
        cpl_propertylist_delete(p);
    }

    /* Now any further difference image stats tables */

    if (! isfirst && (we_expect & STATS_TAB)) {
        p = cpl_propertylist_duplicate(plist);
        if (! (we_get & STATS_TAB))
            vircam_dummy_property(p);
        vircam_dfs_set_product_exten_header(p,product_frame_diffimg_stats,
                                            framelist,parlist,
                                            (char *)recipeid,
                                            "PRO-1.15");
        status = VIR_OK;
        vircam_removewcs(p,&status);
        if (cpl_table_save(ps.diffimstats,NULL,p,outdimst,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);
}

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

static void vircam_dark_combine_dummy_products(void) {

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

    if (we_get == we_expect)
        return;

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

    if (! (we_get & MEANDARK)) {
        ps.outimage = vircam_dummy_image(ps.darks[0]);

        /* Set up the QC parameters */
    
        vircam_dark_combine_config.particle_rate = 0.0;
        vircam_dark_combine_config.darkmed = 0.0;
        vircam_dark_combine_config.darkrms = 0.0;
        vircam_dark_combine_config.nhot = 0;
        vircam_dark_combine_config.hotfrac = 0.0;
        vircam_dark_combine_config.striperms = 0.0;
    }

    /* Do the difference image */

    if ((we_expect & DIFFIMG) && ! (we_get & DIFFIMG)) {
        vircam_dark_combine_config.darkdiff_med = 0.0;
        vircam_dark_combine_config.darkdiff_rms = 0.0;
        ps.diffimg = vircam_dummy_image(ps.darks[0]);
    }

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

    if ((we_expect & STATS_TAB) && ! (we_get & STATS_TAB))
        ps.diffimstats = vircam_create_diffimg_stats(0);

    return;
}

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

static void vircam_dark_combine_hotpix(void) {
    int i,nx,ny,status,nh,nhot,j;
    long npts;
    cpl_image *im;
    unsigned char *bpm;
    float med,mad,lowcut,highcut,*data;
    vir_fits *f;

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

    im = vircam_fits_get_image(ps.good[0]);
    nx = cpl_image_get_size_x(im);
    ny = 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 and
       destripe it. */

    for (i = 0; i < ps.ngood; i++) {
        f = vircam_fits_duplicate(ps.good[i]);
        im = vircam_fits_get_image(f);
        cpl_image_subtract(im,ps.outimage);
        status = VIR_OK;
        vircam_destripe(f,NULL,&status);
        
        /* 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);
        vircam_medmad(data,NULL,npts,&med,&mad);
        lowcut = med - 1.48*mad*vircam_dark_combine_config.thresh;
        highcut = med + 1.48*mad*vircam_dark_combine_config.thresh;
        for (j = 0; j < npts; j++) 
            if (data[j] > highcut || data[j] < lowcut)
                bpm[j] += 1;

        /* Get rid of temporary image */

        vircam_fits_delete(f);
    }

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

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

    /* Clean up... */

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

    vircam_dark_combine_config.nhot = nhot;
    vircam_dark_combine_config.hotfrac = (float)nhot/(float)npts;
}


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

static void vircam_dark_combine_normal(int jext, float exptime) {
    int nx,ny,ndiff,ncells,status;
    long npi,i;
    unsigned char *bpm;
    float med,sig,*idata,grms,gdiff;
    const char *fctid="vircam_dark_combine_normal";
    cpl_table *ctable;
    cpl_image *im;
    cpl_propertylist *p;
    vir_fits *f;

    /* Load up the bad pixel mask */

    nx = cpl_image_get_size_x(ps.outimage);
    ny = cpl_image_get_size_y(ps.outimage);
    npi = nx*ny;
    vircam_dark_combine_config.particle_rate = 0;
    bpm = vircam_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];
    vircam_dark_combine_config.particle_rate = 
        (float)ndiff/(exptime*(float)(ps.ndarks));

    /* Work out the RMS of the mean dark frame */

    idata = cpl_image_get_data(ps.outimage);
    vircam_medmad(idata,bpm,npi,&med,&sig);
    sig *= 1.48;
    vircam_dark_combine_config.darkmed = med;
    vircam_dark_combine_config.darkrms = sig;

    /* Work out the RMS of the stripe pattern */

    im = cpl_image_duplicate(ps.outimage);
    f = vircam_fits_wrap(im,ps.good[0],NULL,NULL);
    status = VIR_OK;
    (void)vircam_destripe(f,ps.master_mask,&status);
    vircam_dark_combine_config.striperms = 
        cpl_propertylist_get_float(vircam_fits_get_ehu(f),"ESO DRS STRIPERMS");
    vircam_fits_delete(f);

    /* Load up the master dark */

    if (ps.master_dark != NULL) {
        ps.mdimage = vircam_fits_load(ps.master_dark,CPL_TYPE_FLOAT,jext);
        if (ps.mdimage == NULL) 
            cpl_msg_info(fctid,"Master dark extension %d won't load",jext);
        else if (vircam_is_dummy(vircam_fits_get_ehu(ps.mdimage))) {
            cpl_msg_info(fctid,"Master dark extension %d is a dummy!",jext);
            freefits(ps.mdimage);
        }
    } else 
        ps.mdimage = NULL;

    /* Load up the channel table */

    if (ps.chantab != NULL) {
        ctable = cpl_table_load(cpl_frame_get_filename(ps.chantab),jext,0);
        if (ctable == NULL) {
            cpl_error_reset();
            cpl_msg_info(fctid,"Channel table extension %d won't load",jext);
        } else if (vircam_chantab_verify(ctable) != VIR_OK) {
            cpl_msg_info(fctid,"Channel table extension %d has errors",jext);
            freetable(ctable);
        } else { 
            p = cpl_propertylist_load(cpl_frame_get_filename(ps.chantab),jext);
            if (vircam_is_dummy(p)) {
                cpl_msg_info(fctid,"Channel table extensions %d is a dummy",
                             jext);
                freetable(ctable);
            }
            freepropertylist(p);
        }
    } else 
        ctable = NULL;

    /* Form the difference image. NB: the difference image routine
       copes if the input mean image and or the channel tables are
       null. Thus if either or both are 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 */

    vircam_dark_combine_config.darkdiff_med = 0.0;
    vircam_dark_combine_config.darkdiff_rms = 0.0;
    ncells = vircam_dark_combine_config.ncells;
    vircam_difference_image(vircam_fits_get_image(ps.mdimage),
                            ps.outimage,bpm,ctable,ncells,1,
                            &gdiff,&grms,&(ps.diffimg),
                            &(ps.diffimstats));
    vircam_mask_clear(ps.master_mask);
    vircam_dark_combine_config.darkdiff_med = gdiff;
    vircam_dark_combine_config.darkdiff_rms = grms;
    freetable(ctable);
    if (ps.diffimg != NULL)
        we_get |= DIFFIMG;
    if (ps.diffimstats != NULL)
        we_get |= STATS_TAB;
    return;
}

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

static int vircam_dark_combine_lastbit(int jext, cpl_frameset *framelist,
                                       cpl_parameterlist *parlist) {
    int retval;
    const char *fctid="vircam_dark_combine_lastbit";

    /* Make whatever dummy products you need */

    vircam_dark_combine_dummy_products();

    /* Save everything */

    cpl_msg_info(fctid,"Saving products for extension %d",jext);
    retval = vircam_dark_combine_save(framelist,parlist);
    if (retval != 0) {
        vircam_dark_combine_tidy(2);
        return(-1);
    }

    /* Free some stuff up */

    vircam_dark_combine_tidy(1);
    return(0);
}

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

static void vircam_dark_combine_init(void) {
    ps.labels = NULL;
    ps.darklist = NULL;
    ps.darks = NULL;
    ps.ndarks = 0;
    ps.good = NULL;
    ps.ngood = 0;
    ps.master_dark = NULL;
    ps.master_mask = NULL;
    ps.chantab = NULL;
    ps.outimage = NULL;
    ps.drs = NULL;
    ps.rejmask = NULL;
    ps.rejplus = NULL;
    ps.mdimage = NULL;
    ps.diffimg = NULL;
    ps.diffimstats = NULL;
    ps.phupaf = NULL;
    we_expect = 0;
    we_get = 0;
}

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

static void vircam_dark_combine_tidy(int level) {

    freeimage(ps.outimage);
    freefitslist(ps.darks,ps.ndarks);
    freepropertylist(ps.drs);
    freefits(ps.mdimage);
    freeimage(ps.diffimg);
    freetable(ps.diffimstats);
    freespace(ps.rejmask);
    freespace(ps.rejplus);
    if (level == 1)
        return;
    freespace(ps.labels);
    freeframeset(ps.darklist);
    freespace(ps.good);
    freeframe(ps.master_dark);
    freemask(ps.master_mask);
    freeframe(ps.chantab);
    freepropertylist(ps.phupaf);
}

/*

$Log: vircam_dark_combine.c,v $
Revision 1.57  2007/10/19 09:25:09  jim
Fixed problems with missing includes

Revision 1.56  2007/10/15 12:53:26  jim
Modified for compatibiliity with cpl_4.0

Revision 1.55  2007/07/18 15:35:41  jim
Added better error handling for missing or corrupt mask extensions

Revision 1.54  2007/07/09 13:21:37  jim
Modified to use new vircam_exten_range and to fix comment on particle_rate
QC parameter

Revision 1.53  2007/04/30 09:40:17  jim
Added more stuff to paf files

Revision 1.52  2007/04/04 10:36:18  jim
Modified to use new dfs tags

Revision 1.51  2007/03/29 12:19:38  jim
Little changes to improve documentation

Revision 1.50  2007/03/02 12:37:16  jim
Removed WCS stuff from table headers

Revision 1.49  2007/03/01 12:41:48  jim
Modified slightly after code checking

Revision 1.48  2007/02/25 06:26:35  jim
Plugged a few memory leaks

Revision 1.47  2007/02/19 10:03:02  jim
Fixed small memory leak

Revision 1.46  2007/02/15 12:17:33  jim
Fixed typo

Revision 1.45  2007/02/15 11:54:09  jim
Modified to make a distinction between initial channel table and one that
has the proper linearity information

Revision 1.44  2007/02/15 06:59:37  jim
Added ability to write QC paf files

Revision 1.43  2007/02/09 14:49:05  jim
Added QC parameter NHOTPIX and HOTFRAC and routine vircam_dark_combine_hotpix

Revision 1.42  2007/02/07 10:12:39  jim
Removed calls to vircam_ndit_correct as this is now no longer necessary

Revision 1.41  2007/02/06 13:11:11  jim
Fixed entry for PRO dictionary in cpl_dfs_set_product_header

Revision 1.40  2007/02/05 14:14:05  jim
Input master frame is now tagged as REFERENCE. QC removed from stats table
headers

Revision 1.39  2006/12/13 13:14:38  jim
Fixed badly scaled sigma

Revision 1.38  2006/12/12 11:30:13  jim
Added QC STRIPERMS calculation

Revision 1.37  2006/11/27 12:13:21  jim
Swapped calls to cpl_propertylist_append to cpl_propertylist_update

Revision 1.36  2006/09/29 11:19:30  jim
changed aliases on parameter names

Revision 1.35  2006/09/09 16:49:39  jim
Header comment update

Revision 1.34  2006/08/27 20:30:02  jim
Major mods to structure of the main processing routine to deal with missing
and dummy frames. Deals better with lower level failures too

Revision 1.33  2006/06/20 19:07:00  jim
Corrects for ndit != 1

Revision 1.32  2006/06/15 09:58:57  jim
Minor changes to docs

Revision 1.31  2006/06/06 13:01:39  jim
Fixed so that the QC parameters go into the correct headers

Revision 1.30  2006/05/17 14:43:58  jim
Fixed problem in save routine which messed up the PRO CATG keywords

Revision 1.29  2006/05/16 13:58:47  jim
Fixed memory leaks that occur from not closing images at the end of
the image extension loop

Revision 1.28  2006/05/04 11:53:14  jim
Fixed the way the _save routine works to be more consistent with the
standard CPL way of doing things

Revision 1.27  2006/04/27 09:46:01  jim
Modified DFS frame types to conform to new dictionary

Revision 1.26  2006/04/25 13:45:56  jim
Fixed to adhere to new calling sequence for vircam_dfs routines

Revision 1.25  2006/04/24 13:46:35  jim
A bit more error trapping in case fits structures can't be loaded

Revision 1.24  2006/03/22 12:13:51  jim
Modified to use new vircam_mask capability

Revision 1.23  2006/03/15 10:43:40  jim
Fixed a few things

Revision 1.22  2006/03/08 14:32:35  jim
Lots of little mods

Revision 1.21  2006/03/03 14:29:06  jim
Now calls routines with vir_fits.

Revision 1.19  2006/02/27 14:05:07  jim
Fixed screwup

Revision 1.18  2006/02/27 13:51:17  jim
new routine

Revision 1.17  2006/02/22 10:01:38  jim
Modified to use new version of vircam_imcombine

Revision 1.16  2006/02/18 11:50:43  jim
Modified the way the dfs product keywords are written using the vircam
routines, rather than the cpl routine that doesn't understand image
extensions

Revision 1.15  2006/01/23 10:35:21  jim
Now allows both BPM or CPM to be used for masking

Revision 1.14  2005/12/14 22:19:11  jim
fixed docs

Revision 1.13  2005/12/12 14:16:20  jim
Fixed typo that caused compilation error

Revision 1.12  2005/12/09 09:47:57  jim
Many changes to add more documentation

Revision 1.11  2005/12/02 10:45:37  jim
The tags used in the sof are now written to the description string in the
constructor. This is so that if they change in the vircam_dfs.h file, they
aren't then hardcopied into each of the recipes...

Revision 1.10  2005/12/01 16:25:48  jim
Made the routine a bit more forgiving if certain master calibration files
were missing. Now does as much as it can with the info it has

Revision 1.9  2005/11/25 09:56:14  jim
Tidied up some more documentation

Revision 1.8  2005/11/23 14:57:40  jim
A bit of tidying in response to splint messages

Revision 1.7  2005/11/08 12:47:44  jim
Made garbage collection a little better

Revision 1.6  2005/11/07 13:14:41  jim
Added better garbage collection and fixed a few bugs

Revision 1.5  2005/11/03 15:16:28  jim
Lots of changes mainly to strengthen error reporting

Revision 1.4  2005/10/14 13:22:12  jim
Added lots of QC checking and diagnostics, so it pretty much does what the
docs say it should

Revision 1.3  2005/09/20 15:07:45  jim
Fixed a few bugs and added a few things

Revision 1.2  2005/08/09 10:24:37  jim
Replaced dodgy calls to cpl_msg_err with correct cpl_msg_error

Revision 1.1.1.1  2005/08/05 08:29:09  jim
Initial import


*/

---