vircam_reset_combine.c

/* $Id: vircam_reset_combine.c,v 1.50 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.50 $
 * $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_mask.h"
#include "vircam_dfs.h"
#include "vircam_mods.h"
#include "vircam_stats.h"
#include "vircam_fits.h"
#include "vircam_pfits.h"
#include "vircam_channel.h"
#include "vircam_paf.h"
#include "vircam_wcsutils.h"

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

#define MEANRESET    1
#define DIFFIMG      2
#define STATS_TAB    4

/* Function prototypes */

static int vircam_reset_combine_create(cpl_plugin *) ;
static int vircam_reset_combine_exec(cpl_plugin *) ;
static int vircam_reset_combine_destroy(cpl_plugin *) ;
static int vircam_reset_combine(cpl_parameterlist *, cpl_frameset *) ;
static int vircam_reset_combine_save(cpl_frameset *framelist, 
                                     cpl_parameterlist *parlist);
static void vircam_reset_combine_dummy_products(void);
static void vircam_reset_combine_normal(int jext);
static int vircam_reset_combine_lastbit(int jext, cpl_frameset *framelist,
                                        cpl_parameterlist *parlist);
static void vircam_reset_combine_init(void);
static void vircam_reset_combine_tidy(int level);

/* Static global variables */

static struct {

    /* Input */

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

    /* Output */

    float       resetmed;
    float       resetrms;
    float       resetdiff_med;
    float       resetdiff_rms;

} vircam_reset_combine_config ;

static struct {
    int              *labels;
    cpl_frameset     *resetlist;
    vir_fits         **resets;
    int              nresets;
    vir_fits         **good;
    int              ngood;
    cpl_frame        *master_reset;
    vir_mask         *master_mask;
    cpl_frame        *chantab;
    cpl_image        *outimage;
    cpl_propertylist *drs;
    unsigned char    *rejmask;
    unsigned char    *rejplus;
    vir_fits         *mrimage;
    cpl_image        *diffimg;
    cpl_table        *diffimstats;
    cpl_propertylist *phupaf;
} ps;

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

static char vircam_reset_combine_description[] =
"vircam_reset_combine -- VIRCAM reset combine recipe.\n\n"
"Combine a list of reset frames into a mean reset frame. Optionally compare \n"
"the output frame to a master reset frame\n\n"
"The program requires the following files in the SOF:\n\n"
"    Tag                   Description\n"
"    -----------------------------------------------------------------------\n"
"    %-21s A list of raw reset images\n"
"    %-21s Optional reference reset 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"
"\n"
"If no master reset frame is made available, then no comparison will be done\n"
"This means there will be no output difference image. If a master reset 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_reset_combine_description,
                   VIRCAM_RESET_RAW,VIRCAM_REF_RESET,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_reset_combine",
                    "VIRCAM reset combination recipe",
                    alldesc,
                    "Jim Lewis",
                    "jrl@ast.cam.ac.uk",
                    vircam_get_license(),
                    vircam_reset_combine_create,
                    vircam_reset_combine_exec,
                    vircam_reset_combine_destroy);

    cpl_pluginlist_append(list,plugin);

    return(0);
}

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

static int vircam_reset_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_reset_combine.combtype",
                                CPL_TYPE_INT,
                                "1 == Median,\n 2 == Mean",
                                "vircam.vircam_reset_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_reset_combine.scaletype",
                                CPL_TYPE_INT,
                                "0 == none,\n 1 == additive offset,\n 2 == multiplicative offset,\n 3 == exposure time scaling + additive offset",
                                "vircam.vircam_reset_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_reset_combine.xrej",
                                CPL_TYPE_BOOL,
                                "True if using extra rejection cycle",
                                "vircam.vircam_reset_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_reset_combine.thresh",
                                CPL_TYPE_DOUBLE,
                                "Rejection threshold in sigma above background",
                                "vircam.vircam_reset_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_reset_combine.ncells",
                               CPL_TYPE_INT,
                               "Number of cells for data channel stats",
                               "vircam.vircam_reset_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_reset_combine.extenum",
                                CPL_TYPE_INT,
                                "Extension number to be done, 0 == all",
                                "vircam.vircam_reset_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_reset_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_reset_combine(recipe->parameters,recipe->frames));
}
                                
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

static int vircam_reset_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_reset_combine(cpl_parameterlist *parlist, 
                                cpl_frameset *framelist) {
    const char *fctid="vircam_reset_combine";
    int nlab,j,jst,jfn,retval,status,i,live,nx,ny;
    cpl_parameter *p;
    vir_fits *ff;

    /* 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\n");
        return(-1);
    }

    /* Initialise a few things */

    vircam_reset_combine_init();
    we_expect |= MEANRESET;

    /* Get the parameters */

    p = cpl_parameterlist_find(parlist,"vircam.vircam_reset_combine.combtype");
    vircam_reset_combine_config.combtype = cpl_parameter_get_int(p);
    p = cpl_parameterlist_find(parlist,"vircam.vircam_reset_combine.scaletype");
    vircam_reset_combine_config.scaletype = cpl_parameter_get_int(p);
    p = cpl_parameterlist_find(parlist,"vircam.vircam_reset_combine.xrej");
    vircam_reset_combine_config.xrej = cpl_parameter_get_bool(p);
    p = cpl_parameterlist_find(parlist,"vircam.vircam_reset_combine.thresh");
    vircam_reset_combine_config.thresh = (float)cpl_parameter_get_double(p);
    p = cpl_parameterlist_find(parlist,"vircam.vircam_reset_combine.ncells");
    vircam_reset_combine_config.ncells = cpl_parameter_get_int(p);
    p = cpl_parameterlist_find(parlist,"vircam.vircam_reset_combine.extenum");
    vircam_reset_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");
        return(-1);
    }

    /* Get the reset frames */

    if ((ps.labels = cpl_frameset_labelise(framelist,vircam_compare_tags,
                                           &nlab)) == NULL) {
        cpl_msg_error(fctid,"Cannot labelise the input frames");
        return(-1);
    }
    if ((ps.resetlist = vircam_frameset_subgroup(framelist,ps.labels,nlab,
                                                 VIRCAM_RESET_RAW)) == NULL) {
        cpl_msg_error(fctid,"Cannot find reset frames in input frameset");
        return(-1);
    }
    ps.nresets = cpl_frameset_get_size(ps.resetlist);
        
    /* Check to see if there is a master reset frame */

    if ((ps.master_reset = vircam_frameset_subgroup_1(framelist,ps.labels,nlab,
                                                      VIRCAM_REF_RESET)) == NULL) 
        cpl_msg_info(fctid,"No master reset 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_reset_combine_config.extenum,
                       (const cpl_frame *)cpl_frameset_get_frame(ps.resetlist,0),
                       &jst,&jfn);
    if (jst == -1 || jfn == -1) {
        cpl_msg_error(fctid,"Unable to continue");
        vircam_reset_combine_tidy(2);
        return(-1);
    }

    /* Get some space for the good frames */

    ps.good = cpl_malloc(ps.nresets*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 the signal a major error */

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

        /* Are any of these reset images good? */

        ps.ngood = 0;
        for (i = 0; i < ps.nresets; i++) {
            ff = ps.resets[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 wee 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_reset_combine_lastbit(j,framelist,parlist);
            if (retval != 0)
                return(-1);
            continue;
        }

        /* Load the mask */

        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]));
        if (vircam_mask_load(ps.master_mask,j,nx,ny) == 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\n",j);
        (void)vircam_imcombine(ps.good,ps.ngood,
                               vircam_reset_combine_config.combtype,
                               vircam_reset_combine_config.scaletype,
                               vircam_reset_combine_config.xrej,
                               vircam_reset_combine_config.thresh,
                               &(ps.outimage),&(ps.rejmask),
                               &(ps.rejplus),&(ps.drs),&status);
        if (status == VIR_OK) {
            we_get |= MEANRESET;
            vircam_reset_combine_normal(j);
        }

        /* Create any dummies and save products */

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

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

static int vircam_reset_combine_save(cpl_frameset *framelist, 
                                     cpl_parameterlist *parlist) {
    cpl_propertylist *plist,*elist,*p,*pafprop;
    int status;
    const char *fctid = "vircam_reset_combine_save";
    const char *outfile = "resetcomb.fits";
    const char *outdiff = "resetdiff.fits";
    const char *outdimst = "resetdifftab.fits";
    const char *outfilepaf = "resetcomb";
    const char *outdiffpaf = "resetdiff";
    const char *recipeid = "vircam_reset_combine";

    /* If we need to make a PHU then do that now based on the first frame
       in the input frame list */

    if (isfirst) {

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

        product_frame_mean_reset = cpl_frame_new();
        cpl_frame_set_filename(product_frame_mean_reset,outfile);
        cpl_frame_set_tag(product_frame_mean_reset,VIRCAM_PRO_RESET);
        cpl_frame_set_type(product_frame_mean_reset,CPL_FRAME_TYPE_IMAGE);
        cpl_frame_set_group(product_frame_mean_reset,CPL_FRAME_GROUP_PRODUCT);
        cpl_frame_set_level(product_frame_mean_reset,CPL_FRAME_LEVEL_FINAL);

        /* Set up product phu */

        plist = vircam_fits_get_phu(ps.resets[0]);
        ps.phupaf = vircam_paf_phu_items(plist);
        vircam_dfs_set_product_primary_header(plist,product_frame_mean_reset,
                                              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_reset);
            return(-1);
        }
        cpl_frameset_insert(framelist,product_frame_mean_reset);

        /* 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_RESET);
            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);

            /* Set up product phu */

            plist = vircam_fits_get_phu(ps.resets[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_RESET_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);

            /* Set up product phu */

            plist = vircam_fits_get_phu(ps.resets[0]);
            vircam_dfs_set_product_primary_header(plist,
                                                  product_frame_diffimg_stats,
                                                  framelist,parlist,
                                                  (char *)recipeid,
                                                  "PRO-1.15");

            /* Fiddle with the extension header */

            elist = vircam_fits_get_ehu(ps.resets[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);
            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_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.resets[0]);

    /* Fiddle with the header now */

    vircam_merge_propertylists(plist,ps.drs);
    p = cpl_propertylist_duplicate(plist);
    if (! (we_get & MEANRESET))
        vircam_dummy_property(p);
    vircam_dfs_set_product_exten_header(p,product_frame_mean_reset,framelist,
                                        parlist,(char *)recipeid,
                                        "PRO-1.15");

    /* Now save the reset image extension */

    cpl_propertylist_update_float(p,"ESO QC RESETMED",
                                  vircam_reset_combine_config.resetmed);
    cpl_propertylist_set_comment(p,"ESO QC RESETMED",
                                 "Median of mean reset frame");
    cpl_propertylist_update_float(p,"ESO QC RESETRMS",
                                  vircam_reset_combine_config.resetrms);
    cpl_propertylist_set_comment(p,"ESO QC RESETRMS",
                                 "RMS of mean reset frame");
    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);
    if (vircam_paf_print((char *)outfilepaf,"VIRCAM/vircam_reset_combine",
                         "QC file",pafprop) != VIR_OK)
        cpl_msg_warning(fctid,"Unable to save PAF for mean reset");
    cpl_propertylist_delete(pafprop);
    cpl_propertylist_delete(p);

    /* Now save the reset 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 RESETDIFF_MED",
                                      vircam_reset_combine_config.resetdiff_med);
        cpl_propertylist_set_comment(p,"ESO QC RESETDIFF_MED",
                                     "Median value of difference image");
        cpl_propertylist_update_float(p,"ESO QC RESETDIFF_RMS",
                                      vircam_reset_combine_config.resetdiff_rms);
        cpl_propertylist_set_comment(p,"ESO QC RESETDIFF_RMS",
                                     "RMS value of 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_propertylist_delete(p);
            cpl_msg_error(fctid,"Cannot save product image extension");
            return(-1);
        }
        /* Write out PAF for difference image */

        pafprop = vircam_paf_req_items(p);
        vircam_merge_propertylists(pafprop,ps.phupaf);
        if (vircam_paf_print((char *)outdiffpaf,"VIRCAM/vircam_reset_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_reset_combine_dummy_products(void) {

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

    if (we_get == we_expect)
        return;

    /* First an output combined reset frame */

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

        /* Set up the QC parameters */

        vircam_reset_combine_config.resetmed = 0.0;
        vircam_reset_combine_config.resetrms = 0.0;
    }

    /* Do the difference image */

    if ((we_expect & DIFFIMG) && ! (we_get & DIFFIMG)) {
        vircam_reset_combine_config.resetdiff_med = 0.0;
        vircam_reset_combine_config.resetdiff_rms = 0.0;

        /* Is a difference image required? If so then let's have it... */

        ps.diffimg = vircam_dummy_image(ps.resets[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_reset_combine_normal(int jext) {
    int nx,ny,ncells;
    long npi;
    unsigned char *bpm;
    float med,sig,*idata,grms,gdiff;
    const char *fctid="vircam_reset_combine_normal";
    cpl_table *ctable;
    cpl_propertylist *p;

    /* 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;
    bpm = vircam_mask_get_data(ps.master_mask);

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

    idata = cpl_image_get_data(ps.outimage);
    vircam_medmad(idata,bpm,npi,&med,&sig);
    sig *= 1.48;
    vircam_reset_combine_config.resetmed = med;
    vircam_reset_combine_config.resetrms = sig;

    /* Load up the master reset */

    if (ps.master_reset != NULL) {
        ps.mrimage = vircam_fits_load(ps.master_reset,CPL_TYPE_FLOAT,jext);
        if (ps.mrimage == NULL) 
            cpl_msg_info(fctid,"Master reset extension %d won't load",jext);
        else if (vircam_is_dummy(vircam_fits_get_ehu(ps.mrimage))) {
            cpl_msg_info(fctid,"Master reset extension %d is a dummy!",jext);
            freefits(ps.mrimage);
        }
    } else 
        ps.mrimage = 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_reset_combine_config.resetdiff_med = 0.0;
    vircam_reset_combine_config.resetdiff_rms = 0.0;
    ncells = vircam_reset_combine_config.ncells;
    vircam_difference_image(vircam_fits_get_image(ps.mrimage),
                            ps.outimage,bpm,ctable,ncells,1,
                            &gdiff,&grms,&(ps.diffimg),
                            &(ps.diffimstats));
    vircam_mask_clear(ps.master_mask);
    vircam_reset_combine_config.resetdiff_med = gdiff;
    vircam_reset_combine_config.resetdiff_rms = grms;
    freetable(ctable);
    if (ps.diffimg != NULL)
        we_get |= DIFFIMG;
    if (ps.diffimstats != NULL)
        we_get |= STATS_TAB;
    return;
}

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

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

    /* Make whatever dummy products you need */

    vircam_reset_combine_dummy_products();

    /* Save everything */

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

    /* Free some stuff up */

    vircam_reset_combine_tidy(1);
    return(0);
}

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

static void vircam_reset_combine_init(void) {
    ps.labels = NULL;
    ps.resetlist = NULL;
    ps.resets = NULL;
    ps.nresets = 0;
    ps.good = NULL;
    ps.master_reset = NULL;
    ps.master_mask = NULL;
    ps.chantab = NULL;
    ps.outimage = NULL;
    ps.drs = NULL;
    ps.rejmask = NULL;
    ps.rejplus = NULL;
    ps.mrimage = NULL;
    ps.diffimg = NULL;
    ps.diffimstats = NULL;
    ps.phupaf = NULL;
}

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

static void vircam_reset_combine_tidy(int level) {
    freeimage(ps.outimage);
    freefitslist(ps.resets,ps.nresets);
    freespace(ps.rejmask);
    freespace(ps.rejplus);
    freepropertylist(ps.drs);
    freefits(ps.mrimage);
    freeimage(ps.diffimg);
    freetable(ps.diffimstats);
    if (level == 1)
        return;
    freespace(ps.labels);
    freeframeset(ps.resetlist);
    freeframe(ps.master_reset);
    freemask(ps.master_mask);
    freeframe(ps.chantab);
    freespace(ps.good);
    freepropertylist(ps.phupaf);
}

/*

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

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

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

Revision 1.47  2007/07/09 13:21:56  jim
Modified to use new version of vircam_exten_range

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

Revision 1.45  2007/03/29 12:19:39  jim
Little changes to improve documentation

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

Revision 1.43  2007/03/01 12:41:49  jim
Modified slightly after code checking

Revision 1.42  2007/02/25 06:27:41  jim
plugged a few memory leaks

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

Revision 1.40  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.39  2007/02/15 06:59:38  jim
Added ability to write QC paf files

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

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

Revision 1.36  2007/01/09 11:39:02  jim
Moved free for ps.good in tidy routine to the correct place

Revision 1.35  2007/01/08 19:09:11  jim
Fixed memory leak

Revision 1.34  2006/12/13 13:19:52  jim
Fixed problem with bad sigma estimate

Revision 1.33  2006/12/08 11:39:27  jim
Fixed bug where we_expect didn't check to see if the difference image was
being produced before deciding whether or not a table would be produced.

Revision 1.32  2006/11/27 12:15:08  jim
changed calls to cpl_propertylist_append to cpl_propertylist_update

Revision 1.31  2006/09/29 11:19:31  jim
changed aliases on parameter names

Revision 1.30  2006/09/09 16:49:40  jim
Header comment update

Revision 1.29  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.28  2006/06/15 09:58:58  jim
Minor changes to docs

Revision 1.27  2006/06/09 11:26:25  jim
Small changes to keep lint happy

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

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

Revision 1.24  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.23  2006/05/09 09:27:06  jim
removed unecessary call to cpl_propertylist_delete

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

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

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

Revision 1.19  2006/03/23 21:18:46  jim
Minor changes mainly to comment headers

Revision 1.18  2006/03/22 12:13:52  jim
Modified to use new vircam_mask capability

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

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

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

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

Revision 1.12  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.11  2006/01/23 10:37:21  jim
Now allows either a BPM or a CPM to be used as a mask

Revision 1.10  2005/12/14 22:19:12  jim
fixed docs

Revision 1.9  2005/12/09 09:47:58  jim
Many changes to add more documentation

Revision 1.8  2005/12/02 10:45:38  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.7  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.6  2005/11/25 09:56:14  jim
Tidied up some more documentation

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

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

Revision 1.3  2005/11/07 13:13:43  jim
Added some docs and calls to vircam_getnpts

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

Revision 1.1  2005/09/29 08:58:25  jim
new routine



*/

---