vircam_dome_flat_combine.c

/* $Id: vircam_dome_flat_combine.c,v 1.45 2007/11/26 09:59:06 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/11/26 09:59:06 $
 * $Revision: 1.45 $
 * $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_pfits.h"
#include "vircam_dfs.h"
#include "vircam_mods.h"
#include "vircam_stats.h"
#include "vircam_fits.h"
#include "vircam_tfits.h"
#include "vircam_channel.h"
#include "vircam_paf.h"
#include "vircam_wcsutils.h"

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

#define MEANDOME    1
#define RATIMG      2
#define STATS_TAB   4

/* Function prototypes */

static int vircam_dome_flat_combine_create(cpl_plugin *) ;
static int vircam_dome_flat_combine_exec(cpl_plugin *) ;
static int vircam_dome_flat_combine_destroy(cpl_plugin *) ;
static int vircam_dome_flat_combine(cpl_parameterlist *, cpl_frameset *) ;
static int vircam_dome_flat_combine_save(cpl_frameset *framelist, 
                                         cpl_parameterlist *parlist);
static void vircam_dome_flat_combine_dummy_products(void);
static void vircam_dome_flat_combine_normal(int jext);
static int vircam_dome_flat_combine_lastbit(int jext, cpl_frameset *framelist,
                                            cpl_parameterlist *parlist);
static void vircam_dome_flat_combine_init(void);
static void vircam_dome_flat_combine_tidy(int level);

/* Static global variables */

static struct {

    /* Input */

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

    /* Output */

    float       flatrms;
    float       flatratio_med;
    float       flatratio_rms;

} vircam_dome_flat_combine_config;


static struct {
    int              *labels;
    cpl_frameset     *domelist;
    cpl_frame        *master_dark;
    cpl_frame        *master_dome_flat;
    cpl_frame        *chantab;
    vir_fits         **good;
    int              ngood;
    vir_mask         *master_mask;
    cpl_image        *outimage;
    vir_fits         **domes;
    int              ndomes;
    cpl_propertylist *drs;
    unsigned char    *rejmask;
    unsigned char    *rejplus;
    vir_tfits        *ctable;
    vir_fits         *mfimage;
    vir_fits         *mdark;
    cpl_image        *ratioimg;
    cpl_table        *ratioimstats;
    cpl_propertylist *phupaf;
} ps;

static int isfirst;
static cpl_frame *product_frame_mean_dome = NULL;
static cpl_frame *product_frame_ratioimg = NULL;
static cpl_frame *product_frame_ratioimg_stats = NULL;
static int we_expect;
static int we_get;

static char vircam_dome_flat_combine_description[] =
"vircam_dome_flat_combine -- VIRCAM dome flat combine recipe.\n\n"
"Combine a list of dome flat frames into a mean frame. Optionally compare \n"
"the output frame to a master dome flat frame\n\n"
"The program accepts the following files in the SOF:\n\n"
"    Tag                   Description\n"
"    -----------------------------------------------------------------------\n"
"    %-21s A list of raw dome flat images\n"
"    %-21s A master dark frame\n"
"    %-21s Optional reference dome flat frame\n"
"    %-21s Optional channel table or\n"
"    %-21s Optional initial channel table\n"
"    %-21s Optional master bad pixel map or\n"
"    %-21s Optional master confidence map\n"
"If no master dome flat is made available, then no comparison will be done\n"
"This means there will be no output ratio image. If a master dome is\n"
"available, but no channel table is, then a ratio 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_dome_flat_combine_description,
                   VIRCAM_DOME_RAW,VIRCAM_CAL_DARK,VIRCAM_REF_DOME_FLAT,
                   VIRCAM_CAL_CHANTAB,VIRCAM_CAL_CHANTAB_INIT,VIRCAM_CAL_BPM,
                   VIRCAM_CAL_CONF);

    cpl_plugin_init(plugin,
                    CPL_PLUGIN_API,
                    VIRCAM_BINARY_VERSION,
                    CPL_PLUGIN_TYPE_RECIPE,
                    "vircam_dome_flat_combine",
                    "VIRCAM dome flat combination recipe",
                    alldesc,
                    "Jim Lewis",
                    "jrl@ast.cam.ac.uk",
                    vircam_get_license(),
                    vircam_dome_flat_combine_create,
                    vircam_dome_flat_combine_exec,
                    vircam_dome_flat_combine_destroy);

    cpl_pluginlist_append(list,plugin);

    return(0);
}

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

static int vircam_dome_flat_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();

    /* Lower threshold for rejecting underexposed images */

    p = cpl_parameter_new_value("vircam.vircam_dome_flat_combine.lthr",
                                CPL_TYPE_DOUBLE,
                                "Low rejection threshold for underexpsed images",
                                "vircam.vircam_dome_flat_combine",
                                0.0);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"lthr");
    cpl_parameterlist_append(recipe->parameters,p);

    /* Upper threshold for rejecting overexposed images */

    p = cpl_parameter_new_value("vircam.vircam_dome_flat_combine.hthr",
                                CPL_TYPE_DOUBLE,
                                "High rejection threshold for overexposed images",
                                "vircam.vircam_dome_flat_combine",
                                65535.0);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"hthr");
    cpl_parameterlist_append(recipe->parameters,p);

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

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

static int vircam_dome_flat_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_dome_flat_combine(cpl_parameterlist *parlist, 
                                    cpl_frameset *framelist) {
    const char *fctid="vircam_dome_flat_combine";
    int nlab,j,jst,jfn,retval,status,live,nx,ny,ndit;
    long i;
    cpl_parameter *p;
    cpl_propertylist *pp;
    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 some things */

    vircam_dome_flat_combine_init();
    we_expect = MEANDOME;

    /* Get the parameters */

    p = cpl_parameterlist_find(parlist,"vircam.vircam_dome_flat_combine.lthr");
    vircam_dome_flat_combine_config.lthr = (float)cpl_parameter_get_double(p);
    p = cpl_parameterlist_find(parlist,"vircam.vircam_dome_flat_combine.hthr");
    vircam_dome_flat_combine_config.hthr = (float)cpl_parameter_get_double(p);
    p = cpl_parameterlist_find(parlist,"vircam.vircam_dome_flat_combine.combtype");
    vircam_dome_flat_combine_config.combtype = cpl_parameter_get_int(p);
    p = cpl_parameterlist_find(parlist,"vircam.vircam_dome_flat_combine.scaletype");
    vircam_dome_flat_combine_config.scaletype = cpl_parameter_get_int(p);
    p = cpl_parameterlist_find(parlist,"vircam.vircam_dome_flat_combine.xrej");
    vircam_dome_flat_combine_config.xrej = cpl_parameter_get_bool(p);
    p = cpl_parameterlist_find(parlist,"vircam.vircam_dome_flat_combine.thresh");
    vircam_dome_flat_combine_config.thresh = (float)cpl_parameter_get_double(p);
    p = cpl_parameterlist_find(parlist,"vircam.vircam_dome_flat_combine.ncells");
    vircam_dome_flat_combine_config.ncells = cpl_parameter_get_int(p);
    p = cpl_parameterlist_find(parlist,"vircam.vircam_dome_flat_combine.extenum");
    vircam_dome_flat_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_dome_flat_combine_tidy(2);
        return(-1);
    }

    /* Get the dome frames */

    if ((ps.labels = cpl_frameset_labelise(framelist,vircam_compare_tags,
                                           &nlab)) == NULL) {
        cpl_msg_error(fctid,"Cannot labelise the input frames");
        vircam_dome_flat_combine_tidy(2);
        return(-1);
    }
    if ((ps.domelist = vircam_frameset_subgroup(framelist,ps.labels,nlab,
                                                VIRCAM_DOME_RAW)) == NULL) {
        cpl_msg_error(fctid,"Cannot find dome frames in input frameset");
        vircam_dome_flat_combine_tidy(2);
        return(-1);
    }
    ps.ndomes = cpl_frameset_get_size(ps.domelist);

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

    if ((ps.master_dark = vircam_frameset_subgroup_1(framelist,ps.labels,nlab,
                                                     VIRCAM_CAL_DARK)) == NULL) {
        cpl_msg_error(fctid,"No master dark found");
        vircam_dome_flat_combine_tidy(2);
        return(-1);
    }
        
    /* Check to see if there is a master dome flat frame */

    if ((ps.master_dome_flat = vircam_frameset_subgroup_1(framelist,ps.labels,
                                                          nlab,VIRCAM_REF_DOME_FLAT)) == NULL)
        cpl_msg_info(fctid,"No master dome flat found -- no ratio image will be formed");
    else 
        we_expect |= RATIMG;
        
    /* 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 ratio image stats and no linearisation will be done");
        } else {
            cpl_msg_info(fctid,"Channel table is labelled INIT -- no linearisation will be done");
            if (we_expect & RATIMG)
                we_expect |= STATS_TAB;
        }
    } else if (we_expect & RATIMG) {
        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_dome_flat_combine_config.extenum,
                       (const cpl_frame *)cpl_frameset_get_frame(ps.domelist,0),
                       &jst,&jfn);
    if (jst == -1 || jfn == -1) {
        cpl_msg_error(fctid,"Unable to continue");
        vircam_dome_flat_combine_tidy(2);
        return(-1);
    }

    /* Get the number of DITs */

    pp = cpl_propertylist_load(cpl_frame_get_filename(cpl_frameset_get_frame(ps.domelist,0)),0);
    if (vircam_pfits_get_ndit(pp,&ndit) != VIR_OK) {
        cpl_msg_error(fctid,"No value for NDIT available");
        freepropertylist(pp);
        vircam_dome_flat_combine_tidy(2);
        return(-1);
    }
    cpl_propertylist_delete(pp);

    /* Get some space for the good frames */

    ps.good = cpl_malloc(ps.ndomes*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 the images and the master dark. */

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

        /* Are any of these dome flats any good? */
        
        ps.ngood = 0;
        for (i = 0; i < ps.ndomes; i++) {
            ff = ps.domes[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 
           dummy products and move on */

        if (ps.ngood == 0) {
            cpl_msg_info(fctid,"All images flagged bad for this extension");
            retval = vircam_dome_flat_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);
        }

        /* Sort out the images that are either over or under exposed */

        vircam_overexp(ps.good,&(ps.ngood),
                       vircam_dome_flat_combine_config.lthr,
                       vircam_dome_flat_combine_config.hthr,0);

        /* Check to see how many are left. If there aren't any, then
           signal a major error */

        if (ps.ngood == 0) {
            cpl_msg_info(fctid,"All images either under or overexposed");
            retval = vircam_dome_flat_combine_lastbit(j,framelist,parlist);
            if (retval != 0)
                return(-1);
            continue;
        }

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

        ps.mdark = vircam_fits_load(ps.master_dark,CPL_TYPE_FLOAT,j);
        if (ps.mdark == NULL) {
            cpl_msg_info(fctid,"Can't load master dark for extension %d",j);
            retval = vircam_dome_flat_combine_lastbit(j,framelist,parlist);
            if (retval != 0)
                return(-1);
            continue;
        } else if (vircam_is_dummy(vircam_fits_get_ehu(ps.mdark))) {
            cpl_msg_info(fctid,"Can't master dark extension %d is a dummy",j);
            retval = vircam_dome_flat_combine_lastbit(j,framelist,parlist);
            if (retval != 0)
                return(-1);
            continue;
        }

        /* Loop for each image and dark correct */

        cpl_msg_info(fctid,"Dark correcting extension %d\n",j);
        for (i = 0; i < ps.ngood; i++)
            vircam_darkcor((ps.good)[i],ps.mdark,1.0,&status);

        /* We need to load the channel table (if it exists) for linearisation
           and for the ratio image stats table */

        if (ps.chantab != NULL) {
            ps.ctable = vircam_tfits_load(ps.chantab,j);
            if (ps.ctable == NULL) {
                cpl_msg_info(fctid,"Channel table extension %d won't load",j);
            } else if (vircam_chantab_verify(vircam_tfits_get_table(ps.ctable)) != VIR_OK) {
                cpl_msg_info(fctid,"Channel table extension %d has errors",j);
                freetfits(ps.ctable);
            } else { 
                pp = cpl_propertylist_load(cpl_frame_get_filename(ps.chantab),
                                           j);
                if (vircam_is_dummy(pp)) {
                    cpl_msg_info(fctid,
                                 "Channel table extensions %d is a dummy",j);
                    freetfits(ps.ctable);
                }
                freepropertylist(pp);
            }
        } else 
            ps.ctable = NULL;
                
        /* Loop for each of the input images and linearise it if there
           is a channel table */

        if (ps.ctable != NULL) {
            cpl_msg_info(fctid,"Linearising extension %d\n",j);
            for (i = 0; i < ps.ngood; i++)
                (void)vircam_lincor((ps.good)[i],ps.ctable,1,ndit,&status);
        }

        /* Call the combine module */

        cpl_msg_info(fctid,"Doing combination for extension %d\n",j);
        (void)vircam_imcombine(ps.good,ps.ngood,
                               vircam_dome_flat_combine_config.combtype,
                               vircam_dome_flat_combine_config.scaletype,
                               vircam_dome_flat_combine_config.xrej,
                               vircam_dome_flat_combine_config.thresh,
                               &(ps.outimage),&(ps.rejmask),&(ps.rejplus),
                               &(ps.drs),&status);

        /* If these correction and combination routines failed at any stage
           then get out of here */

        if (status == VIR_OK) {
            we_get |= MEANDOME;
            vircam_dome_flat_combine_normal(j);
        } else {
            cpl_msg_info(fctid,"A processing step failed");
        }

        /* Create any dummies and save the products */
        
        retval = vircam_dome_flat_combine_lastbit(j,framelist,parlist);
        if (retval != 0)
            return(-1);

    }
    vircam_dome_flat_combine_tidy(2);
    return(0);
}


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

static int vircam_dome_flat_combine_save(cpl_frameset *framelist, 
                                         cpl_parameterlist *parlist) {
    cpl_propertylist *plist,*elist,*p,*paf;
    int status;
    const char *fctid = "vircam_dome_flat_combine_save";
    const char *outfile = "domecomb.fits";
    const char *outdiff = "domeratio.fits";
    const char *outdimst = "domeratiotab.fits";
    const char *outpaf = "domecomb";
    const char *outdpaf = "domeratio";
    const char *recipeid = "vircam_dome_flat_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_dome = cpl_frame_new();
        cpl_frame_set_filename(product_frame_mean_dome,outfile);
        cpl_frame_set_tag(product_frame_mean_dome,VIRCAM_PRO_DOME_FLAT);
        cpl_frame_set_type(product_frame_mean_dome,CPL_FRAME_TYPE_IMAGE);
        cpl_frame_set_group(product_frame_mean_dome,CPL_FRAME_GROUP_PRODUCT);
        cpl_frame_set_level(product_frame_mean_dome,CPL_FRAME_LEVEL_FINAL);

        /* Set up the PHU header */

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

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

        if (we_expect & RATIMG) {
            product_frame_ratioimg = cpl_frame_new();
            cpl_frame_set_filename(product_frame_ratioimg,outdiff);
            cpl_frame_set_tag(product_frame_ratioimg,
                              VIRCAM_PRO_RATIOIMG_DOME_FLAT);
            cpl_frame_set_type(product_frame_ratioimg,CPL_FRAME_TYPE_IMAGE);
            cpl_frame_set_group(product_frame_ratioimg,
                                CPL_FRAME_GROUP_PRODUCT);
            cpl_frame_set_level(product_frame_ratioimg,CPL_FRAME_LEVEL_FINAL);

            /* Set up the PHU header */

            plist = vircam_fits_get_phu(ps.domes[0]);
            vircam_dfs_set_product_primary_header(plist,product_frame_ratioimg,
                                                  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_ratioimg);
                return(-1);
            }
            cpl_frameset_insert(framelist,product_frame_ratioimg);
        }

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

        if (we_expect & STATS_TAB) {
            product_frame_ratioimg_stats = cpl_frame_new();
            cpl_frame_set_filename(product_frame_ratioimg_stats,outdimst);
            cpl_frame_set_tag(product_frame_ratioimg_stats,
                              VIRCAM_PRO_RATIOIMG_DOME_FLAT_STATS);
            cpl_frame_set_type(product_frame_ratioimg_stats,
                               CPL_FRAME_TYPE_TABLE);
            cpl_frame_set_group(product_frame_ratioimg_stats,
                                CPL_FRAME_GROUP_PRODUCT);
            cpl_frame_set_level(product_frame_ratioimg_stats,
                                CPL_FRAME_LEVEL_FINAL);

            /* Set up the PHU header */

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

            /* Fiddle with the extension header now */

            elist = vircam_fits_get_ehu(ps.domes[0]);
            p = cpl_propertylist_duplicate(elist);
            vircam_merge_propertylists(p,ps.drs);
            vircam_paf_append(p,vircam_fits_get_phu(ps.domes[0]),
                              "ESO INS FILT1 NAME");
            if (! (we_get & STATS_TAB)) 
                vircam_dummy_property(p);
            vircam_dfs_set_product_exten_header(p,product_frame_ratioimg_stats,
                                                framelist,parlist,
                                                (char *)recipeid,
                                                "PRO-1.15");
            status = VIR_OK;
            vircam_removewcs(p,&status);

            /* And finally the difference image stats table */

            if (cpl_table_save(ps.ratioimstats,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_ratioimg_stats);
        }
    }

    /* Get the extension property list */

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

    /* Fiddle with the header now */

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

    cpl_propertylist_update_float(p,"ESO QC FLATRMS",
                                  vircam_dome_flat_combine_config.flatrms);
    cpl_propertylist_set_comment(p,"ESO QC FLATRMS","RMS of output flat");
    if (cpl_image_save(ps.outimage,outfile,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 the PAF */

    paf = vircam_paf_req_items(p);
    vircam_merge_propertylists(paf,ps.phupaf);
    vircam_paf_append(paf,vircam_fits_get_phu(ps.domes[0]),
                      "ESO INS FILT1 NAME");
    if (vircam_paf_print((char *)outpaf,"VIRCAM/vircam_dome_flat_combine",
                         "QC file",paf) != VIR_OK)
        cpl_msg_warning(fctid,"Unable to save PAF for mean dome");
    cpl_propertylist_delete(paf);    
    cpl_propertylist_delete(p);

    /* Now save the dome ratio image extension */

    if (we_expect & RATIMG) {
        p = cpl_propertylist_duplicate(plist);
        if (! (we_get & RATIMG))
            vircam_dummy_property(p);
        cpl_propertylist_update_float(p,"ESO QC FLATRATIO_MED",
                                      vircam_dome_flat_combine_config.flatratio_med);
        cpl_propertylist_set_comment(p,"ESO QC FLATRATIO_MED",
                                     "Median of ratio map");
        cpl_propertylist_update_float(p,"ESO QC FLATRATIO_RMS",
                                      vircam_dome_flat_combine_config.flatratio_rms);
        cpl_propertylist_set_comment(p,"ESO QC FLATRATIO_RMS",
                                     "RMS of ratio map");
        vircam_dfs_set_product_exten_header(p,product_frame_ratioimg,
                                            framelist,parlist,
                                            (char *)recipeid,
                                            "PRO-1.15");
        if (cpl_image_save(ps.ratioimg,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 the PAF */

        paf = vircam_paf_req_items(p);
        vircam_merge_propertylists(paf,ps.phupaf);
        vircam_paf_append(paf,vircam_fits_get_phu(ps.domes[0]),
                          "ESO INS FILT1 NAME");
        if (vircam_paf_print((char *)outdpaf,"VIRCAM/vircam_dome_flat_combine",
                             "QC file",paf) != VIR_OK)
            cpl_msg_warning(fctid,"Unable to save PAF for ratio image");
        cpl_propertylist_delete(paf);    
        cpl_propertylist_delete(p);
    }

    /* Now any further ratio 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(plist,product_frame_ratioimg_stats,
                                            framelist,parlist,
                                            (char *)recipeid,
                                            "PRO-1.15");
        status = VIR_OK;
        vircam_removewcs(p,&status);
        if (cpl_table_save(ps.ratioimstats,NULL,p,outdimst,CPL_IO_EXTEND)
                           != CPL_ERROR_NONE) {
            cpl_propertylist_delete(p);
            cpl_msg_error(fctid,"Cannot save product table extension");
            return(-1);
        }       
        cpl_propertylist_delete(p);
    }

    /* Get out of here */

    return(0);
}

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

static void vircam_dome_flat_combine_dummy_products(void) {

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

    if (we_get == we_expect)
        return;

    /* First an output combined dome  frame */

    if (! (we_get & MEANDOME)) {
        ps.outimage = vircam_dummy_image(ps.domes[0]);
        vircam_dome_flat_combine_config.flatrms = 0.0;
    }

    /* Do a ratio image */

    if ((we_expect & RATIMG) && ! (we_get & RATIMG)) {
        vircam_dome_flat_combine_config.flatratio_med = 0.0;
        vircam_dome_flat_combine_config.flatratio_rms = 0.0;
        ps.ratioimg = vircam_dummy_image(ps.domes[0]);
    }

    /* If a ratio image stats table is required, then do that now */
   
    if ((we_expect & STATS_TAB) && ! (we_get & STATS_TAB)) 
        ps.ratioimstats = vircam_create_diffimg_stats(0);

    return;
}

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

static void vircam_dome_flat_combine_normal(int jext) {
    int nx,ny,ncells;
    long npi;
    unsigned char *bpm;
    float *idata,med,sig,gdiff,grms;
    const char *fctid="vircam_dome_flat_combine_normal";

    /* Load up a 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 dome frame */

    idata = cpl_image_get_data(ps.outimage);
    vircam_medsig(idata,bpm,npi,&med,&sig);

    /* Divide through by the median */

    cpl_image_divide_scalar(ps.outimage,med);
    vircam_medmad(idata,bpm,npi,&med,&sig);
    sig *= 1.48;
    vircam_dome_flat_combine_config.flatrms = sig;

    /* Load up the master dome flat */

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


    /* Create a ratio image. NB: the difference image routine copes if the 
       input mean image 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_dome_flat_combine_config.flatratio_med = 0.0;
    vircam_dome_flat_combine_config.flatratio_rms = 0.0;
    ncells = vircam_dome_flat_combine_config.ncells;
    vircam_difference_image(vircam_fits_get_image(ps.mfimage),ps.outimage,bpm,
                            vircam_tfits_get_table(ps.ctable),ncells,2,
                            &gdiff,&grms,&(ps.ratioimg),
                            &(ps.ratioimstats));
    vircam_mask_clear(ps.master_mask);
    vircam_dome_flat_combine_config.flatratio_med = gdiff;
    vircam_dome_flat_combine_config.flatratio_rms = grms;
    if (ps.ratioimg != NULL)
        we_get |= RATIMG;
    if (ps.ratioimstats != NULL)
        we_get |= STATS_TAB;
    return;
}

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

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

    /* Make whatever dummy products you need */

    vircam_dome_flat_combine_dummy_products();

    /* Save everything */

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

    /* Free some stuff up */

    vircam_dome_flat_combine_tidy(1);
    return(0);
}

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

static void vircam_dome_flat_combine_init(void) {
    ps.labels = NULL;
    ps.domelist = NULL;
    ps.domes = NULL;
    ps.good = NULL;
    ps.master_dark = NULL;
    ps.master_dome_flat = NULL;
    ps.master_mask = NULL;
    ps.chantab = NULL;
    ps.ctable = NULL;
    ps.outimage = NULL;
    ps.drs = NULL;
    ps.rejmask = NULL;
    ps.rejplus = NULL;
    ps.mfimage = NULL;
    ps.ratioimg = NULL;
    ps.ratioimstats = NULL;
    ps.phupaf = NULL;
}

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

static void vircam_dome_flat_combine_tidy(int level) {
    freeimage(ps.outimage);
    freefitslist(ps.domes,ps.ndomes);
    freepropertylist(ps.drs);
    freespace(ps.rejmask);
    freespace(ps.rejplus);
    freetfits(ps.ctable);
    freefits(ps.mfimage);
    freefits(ps.mdark);
    freeimage(ps.ratioimg);
    freetable(ps.ratioimstats);
    if (level == 1)
        return;

    freespace(ps.labels);
    freeframeset(ps.domelist);
    freeframe(ps.master_dark);
    freeframe(ps.master_dome_flat);
    freeframe(ps.chantab);
    freespace(ps.good);
    freemask(ps.master_mask);
    freepropertylist(ps.phupaf);

}

/*

$Log: vircam_dome_flat_combine.c,v $
Revision 1.45  2007/11/26 09:59:06  jim
Recipe now takes ndit into account when doing linearity correction

Revision 1.44  2007/10/19 09:25:09  jim
Fixed problems with missing includes

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Revision 1.28  2006/12/13 13:26:09  jim
Fxied bad sigma estimate

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

Revision 1.26  2006/11/10 09:21:48  jim
Fixed bug in save routine to get the correct propertylist written to product
header. Also fixed memory deallocation error in tidy routine

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

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

Revision 1.23  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.22  2006/07/11 14:55:12  jim
Now checks for zeros in the output flat and replaces them

Revision 1.21  2006/06/20 19:07:01  jim
Corrects for ndit != 1

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

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

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

Revision 1.17  2006/05/26 19:38:20  jim
Fixed recipe so that the output flat is normalised as it should be

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

Revision 1.15  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.14  2006/05/10 14:30:13  jim
deleted ps.domes twice

Revision 1.13  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.12  2006/04/27 09:34:48  jim
Fixed data type bug in parameter lthr and hthr. Modified headers to change the
dfs keywords to the latest values. Removed confidence map part.

Revision 1.11  2006/04/25 14:01:37  jim
fixed bug in docs

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

Revision 1.9  2006/04/20 11:32:10  jim
Added new parameters to cull out over or under exposed images

Revision 1.8  2006/03/23 21:18:45  jim
Minor changes mainly to comment headers

Revision 1.7  2006/03/22 14:02:51  jim
cosmetic changes to keep lint happy

Revision 1.6  2006/03/22 13:19:59  jim
Added confidence map generation

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

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

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

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

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


*/

---