vircam_dark_current.c

/* $Id: vircam_dark_current.c,v 1.46 2007/11/20 09:40:27 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/20 09:40:27 $
 * $Revision: 1.46 $
 * $Name:  $
 */

/* Includes */

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

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

#include "vircam_utils.h"
#include "vircam_mask.h"
#include "vircam_pfits.h"
#include "vircam_dfs.h"
#include "vircam_stats.h"
#include "vircam_paf.h"

/* Function prototypes */

static int vircam_dark_current_create(cpl_plugin *) ;
static int vircam_dark_current_exec(cpl_plugin *) ;
static int vircam_dark_current_destroy(cpl_plugin *) ;
static int vircam_dark_current(cpl_parameterlist *, cpl_frameset *) ;
static int vircam_dark_current_save(cpl_frameset *filelist,
                                    cpl_parameterlist *parlist);
static int vircam_dark_current_lastbit(int jext, cpl_frameset *framelist,
                                       cpl_parameterlist *parlist);
static void vircam_dark_current_init(void);
static void vircam_dark_current_tidy();

/* Static global variables */

static struct {

    /* Input */

    float       thresh;
    int         extenum;

    /* Output */

    float       mean_dark_current;
} vircam_dark_current_config;

static struct {
    int              *labels;
    cpl_frameset     *darklist;
    vir_mask         *master_mask;
    int              nframes;
    float            **data;
    vir_fits         **allfits;
    double           *subset;
    double           *exps;
    cpl_image        *outimage;
    vir_fits         **good;
    int              ngood;
    cpl_propertylist *phupaf;
} ps;

static cpl_frame *product_frame = NULL;
static int isfirst;
static int dummy;

static char vircam_dark_current_description[] =
"vircam_dark_current -- VIRCAM recipe for measuring dark current.\n"
"A list of dark frames is given. A robust estimate of the dark current\n"
"is calculated by fitting a median slope to the dark value vs exposure time\n"
"for each pixel. The output is to a dark current map which shows the dark\n"
"current in counts per second for each input pixel.\n\n"
"The program requires the following files in the SOF:\n\n"
"    Tag                   Description\n"
"    -----------------------------------------------------------------------\n"
"    %-21s A list of raw dark images with various exposure times\n"
"    %-21s Optional master bad pixel map or\n"
"    %-21s Optional master confidence map\n"
"\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_current_description,
                   VIRCAM_DARKCUR_RAW,VIRCAM_CAL_BPM,VIRCAM_CAL_CONF);

    cpl_plugin_init(plugin,
                    CPL_PLUGIN_API,
                    VIRCAM_BINARY_VERSION,
                    CPL_PLUGIN_TYPE_RECIPE,
                    "vircam_dark_current",
                    "VIRCAM recipe to determine detector dark current",
                    alldesc,
                    "Jim Lewis",
                    "jrl@ast.cam.ac.uk",
                    vircam_get_license(),
                    vircam_dark_current_create,
                    vircam_dark_current_exec,
                    vircam_dark_current_destroy);

    cpl_pluginlist_append(list,plugin);

    return(0);
}

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

static int vircam_dark_current_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 rejection threshold parameter */

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

    /* Extension number of input frames to use */

    p = cpl_parameter_new_range("vircam.vircam_dark_current.extenum",
                                CPL_TYPE_INT,
                                "Extension number to be done, 0 == all",
                                "vircam.vircam_dark_current",
                                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_current_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_current_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_current(recipe->parameters,recipe->frames));
}

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

static int vircam_dark_current(cpl_parameterlist *parlist, 
                               cpl_frameset *framelist) {
    int jst,jfn,nlab,i,j,nx,ny,n,retval,live;
    long npts;
    double intercept,slope,sig;
    const char *fctid = "vircam_dark_current";
    float *outdata,val;
    unsigned char *bpm;
    vir_fits *ff;
    cpl_frame *cur_frame;
    cpl_propertylist *plist;
    cpl_parameter *p;
    
    /* 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 variables */

    vircam_dark_current_init();

    /* Get the parameters */

    p = cpl_parameterlist_find(parlist,"vircam.vircam_dark_current.thresh");
    vircam_dark_current_config.thresh = (float)cpl_parameter_get_double(p);
    p = cpl_parameterlist_find(parlist,"vircam.vircam_dark_current.extenum");
    vircam_dark_current_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 dark frames. Make sure there are at least 2 of them */

    if ((ps.labels = cpl_frameset_labelise(framelist,vircam_compare_tags,
                                           &nlab)) == NULL) {
        cpl_msg_error(fctid,"Cannot labelise the input frameset");
        vircam_dark_current_tidy();
        return(-1);
    }
    if ((ps.darklist = vircam_frameset_subgroup(framelist,ps.labels,nlab,
                                                VIRCAM_DARKCUR_RAW)) == NULL) {
        cpl_msg_error(fctid,"Cannot find dark frames in input frameset");
        vircam_dark_current_tidy();
        return(-1);
    }
    ps.nframes = cpl_frameset_get_size(ps.darklist);
    if (ps.nframes < 2) {
        cpl_msg_error(fctid,"Dark frameset doesn't have enough frames");
        vircam_dark_current_tidy();
        return(-1);
    }

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

    /* Get some workspace for the data arrays */

    ps.data = cpl_malloc(ps.nframes*sizeof(float *));
    ps.subset = cpl_malloc(ps.nframes*sizeof(double));
    ps.exps = cpl_malloc(ps.nframes*sizeof(double));

    /* Fill in the exposure times */

    for (i = 0; i < ps.nframes; i++) {
        cur_frame = cpl_frameset_get_frame(ps.darklist,i);
        plist = cpl_propertylist_load(cpl_frame_get_filename(cur_frame),0);
        if (vircam_pfits_get_exptime(plist,&val) != VIR_OK) {
            cpl_msg_error(fctid,"Unable to get exposure time for %s",
                          cpl_frame_get_filename(cur_frame));
            return(-1);
        }
        ps.exps[i] = (double)val;
        cpl_propertylist_delete(plist);
    }

    /* 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_current_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_current_tidy();
        return(-1);
    }

    /* Get some space for the good frames */

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

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

    for (j = jst; j <= jfn; j++) {
        isfirst = (j == jst);
        dummy = 0;
        vircam_dark_current_config.mean_dark_current = 0.0;

        /* Load the image data from each frame. If there was an
           error loading, then just create a dummy output */

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

        /* Are any of these images good? */

        ps.ngood = 0;
        for (i = 0; i < ps.nframes; i++) {
            ff = ps.allfits[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 too few good images, then signal that we need to 
           create some dummy products and move on */

        if (ps.ngood < 2) {
            cpl_msg_warning(fctid,"Need at least 2 good images -- %d found",
                          ps.ngood);
            dummy = 1;
            retval = vircam_dark_current_lastbit(j,framelist,parlist);
            freefitslist(ps.allfits,ps.nframes);
            freeimage(ps.outimage);
            if (retval != 0) 
                return(-1);
            continue;
        }

        /* Get the data arrays */

        for (i = 0; i < ps.ngood; i++)
            ps.data[i] = cpl_image_get_data(vircam_fits_get_image(ps.allfits[i]));

        /* Load the BPM */

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

        /* Get an output image */

        ps.outimage = cpl_image_new(nx,ny,CPL_TYPE_FLOAT);
        outdata = cpl_image_get_data(ps.outimage);

        /* Now loop over all pixels and work out the slope */

        cpl_msg_info(fctid,"Doing dark current fits for extension %d",j);
        npts = (long)(nx*ny);
        for (n = 0; n < npts; n++) {
            if (bpm[n] != 0) {
                slope = 0.0;
            } else {
                for (i = 0; i < ps.ngood; i++)
                    ps.subset[i] = (double)(ps.data[i][n]);
                vircam_linfit(ps.ngood,ps.exps,ps.subset,&intercept,&slope,
                    &sig);
            }

            /* Store the slope away */

            outdata[n] = (float)slope;
        }

        /* Get the median value of the dark current */

        vircam_dark_current_config.mean_dark_current = 
            vircam_med(outdata,bpm,npts);

        /* Save the last part of the processing and saving */

        (void)vircam_dark_current_lastbit(j,framelist,parlist);

        /* Tidy up */

        freeimage(ps.outimage);
        vircam_mask_clear(ps.master_mask);
        freefitslist(ps.allfits,ps.nframes);
    }

    /* Tidy up */

    vircam_dark_current_tidy();

    return(0);
}

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

static int vircam_dark_current_save(cpl_frameset *framelist, 
                                    cpl_parameterlist *parlist) {
    cpl_propertylist *plist,*p;
    const char *fctid = "vircam_dark_current_save";
    const char *outfile = "darkcurrent.fits";
    const char *outpaf = "darkcurrent";
    const char *recipeid = "vircam_dark_current";
    float darkcur_med;

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

    darkcur_med = vircam_dark_current_config.mean_dark_current;
    if (isfirst) {

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

        product_frame = cpl_frame_new();
        cpl_frame_set_filename(product_frame,outfile);
        cpl_frame_set_tag(product_frame,VIRCAM_PRO_DARKCUR);
        cpl_frame_set_type(product_frame,CPL_FRAME_TYPE_IMAGE);
        cpl_frame_set_group(product_frame,CPL_FRAME_GROUP_PRODUCT);
        cpl_frame_set_level(product_frame,CPL_FRAME_LEVEL_FINAL);

        /* Set the PHU header */

        plist = vircam_fits_get_phu(ps.allfits[0]);
        vircam_dfs_set_product_primary_header(plist,product_frame,framelist,
                                              parlist,(char *)recipeid,
                                              "PRO-1.15");
        ps.phupaf = vircam_paf_phu_items(plist);

        /* '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_propertylist_delete(plist);
            return(-1);
        }
        cpl_frameset_insert(framelist,product_frame);
    }

    /* Get the header for the extension */

    plist = vircam_fits_get_ehu(ps.allfits[0]);
    vircam_dfs_set_product_exten_header(plist,product_frame,framelist,
                                        parlist,(char *)recipeid,
                                        "PRO-1.15");

    /* Add the mean dark current to the header as a QC parameter */

    cpl_propertylist_update_float(plist,"ESO QC DARKCURRENT",darkcur_med);
    cpl_propertylist_set_comment(plist,"ESO QC DARKCURRENT",
                                 "[ADU/s] Median dark current");
    if (dummy)
        vircam_dummy_property(plist);

    /* Now save the image */

    if (cpl_image_save(ps.outimage,outfile,CPL_BPP_IEEE_FLOAT,plist,
                       CPL_IO_EXTEND) != CPL_ERROR_NONE) {
        cpl_msg_error(fctid,"Cannot save product image extension");
        cpl_propertylist_delete(plist);
        return(-1);
    }

    /* Write a PAF now */

    p = vircam_paf_req_items(plist);
    vircam_merge_propertylists(p,ps.phupaf);
    if (vircam_paf_print((char *)outpaf,"VIRCAM/vircam_dark_current","QC file",
                         p) != VIR_OK)
        cpl_msg_warning(fctid,"Unable to write PAF\n");
    cpl_propertylist_delete(p);

    /* Get out of here */

    return(0);
}

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

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

    /* If this is a dummy result then create it now */

    if (dummy)
        ps.outimage = vircam_dummy_image(ps.allfits[0]);
                       
    /* Save the result */

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


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

static void vircam_dark_current_init(void) {
    ps.labels = NULL;
    ps.darklist = NULL;
    ps.master_mask = NULL;
    ps.data = NULL;
    ps.allfits = NULL;
    ps.subset = NULL;
    ps.exps = NULL;
    ps.outimage = NULL;
    ps.nframes = 0;
    ps.good = NULL;
    ps.phupaf = NULL;
}

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

static void vircam_dark_current_tidy(void) {

    freespace(ps.labels);
    freeframeset(ps.darklist);
    freemask(ps.master_mask);
    freespace(ps.data);
    freespace(ps.subset);
    freespace(ps.exps);
    freeimage(ps.outimage);
    freefitslist(ps.allfits,ps.nframes);
    ps.nframes = 0;
    freespace(ps.good);
    freepropertylist(ps.phupaf);
}

/*

$Log: vircam_dark_current.c,v $
Revision 1.46  2007/11/20 09:40:27  jim
changed values for linear fit to doubles

Revision 1.45  2007/11/14 14:47:53  jim
vircam_linfit now works only with doubles

Revision 1.44  2007/10/25 18:39:22  jim
Altered to remove some lint messages

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

Revision 1.42  2007/09/06 21:37:53  jim
fixed call to vircam_dfs_setup_product_ routines to use the full input
frameset

Revision 1.41  2007/08/23 09:01:34  jim
Error when there aren't enough frames is now just a warning

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

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

Revision 1.38  2007/06/13 08:11:27  jim
Modified docs to reflect changes in DFS tags

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

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

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

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

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:11  jim
Fixed entry for PRO dictionary in cpl_dfs_set_product_header

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

Revision 1.29  2006/11/10 09:19:47  jim
fixed typo

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

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

Revision 1.26  2006/09/04 23:02:14  jim
Modified to deal with det live issues. Also does a better job of dealing
with duff input

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

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

Revision 1.23  2006/06/09 11:32:59  jim
A few more minor fixes for lint

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

Revision 1.21  2006/05/17 11:15:38  jim
plugged memory leak

Revision 1.20  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.19  2006/04/27 09:46:01  jim
Modified DFS frame types to conform to new dictionary

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

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

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

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

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

Revision 1.13  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.12  2006/01/23 10:35:55  jim
Now allows either an BPM or CPM to be used as a mask

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

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

Revision 1.9  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.8  2005/12/01 16:25:06  jim
Fixed default output file extension

Revision 1.7  2005/11/25 09:37:10  jim
Fitting now done by vircam_linfit

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

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

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

Revision 1.3  2005/08/09 11:09:39  jim
Replaced dodgy call to cpl_framelist_delete with correct cpl_frameset_delete

Revision 1.2  2005/08/09 10:24:38  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


*/

---