vircam_interleave.c

/* $Id: vircam_interleave.c,v 1.17 2008/07/10 13:01:35 jim Exp $
 *
 * This file is part of the VIRCAM Pipeline
 * Copyright (C) 2006 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: 2008/07/10 13:01:35 $
 * $Revision: 1.17 $
 * $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_wcsutils.h"

/* Function prototypes */

static int vircam_interleave_create(cpl_plugin *) ;
static int vircam_interleave_exec(cpl_plugin *) ;
static int vircam_interleave_destroy(cpl_plugin *) ;
static int vircam_interleave_test(cpl_parameterlist *, cpl_frameset *) ;
static int vircam_interleave_save(cpl_frameset *framelist, 
                                  cpl_parameterlist *parlist);
static void vircam_interleave_init(void);
static void vircam_interleave_tidy(void);

/* Static global variables */

static struct {

    /* Input */

    int         extenum;

} vircam_interleave_config;


static struct {
    int              *labels;
    cpl_frameset     *imagelist;
    vir_fits         **images;
    cpl_frameset     *conflist;
    vir_fits         **confs;
    int              nimages;
    int              nconfs;
    cpl_image        *outimage;
    cpl_image        *outconf;
    cpl_propertylist *plist;
} ps;

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

static char vircam_interleave_description[] =
"vircam_interleave -- VIRCAM test interleave recipe.\n\n"
"Interleave a list of frames into an output frame.\n\n"
"The program accepts the following files in the SOF:\n\n"
"    Tag                   Description\n"
"    -----------------------------------------------------------------------\n"
"    %-21s A list of images\n"
"    %-21s A list of confidence maps\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_interleave_description,
                   VIRCAM_TEST_SCIENCE_RAW,VIRCAM_CAL_CONF);

    cpl_plugin_init(plugin,
                    CPL_PLUGIN_API,
                    VIRCAM_BINARY_VERSION,
                    CPL_PLUGIN_TYPE_RECIPE,
                    "vircam_interleave",
                    "VIRCAM interleaf test recipe [test]",
                    alldesc,
                    "Jim Lewis",
                    "jrl@ast.cam.ac.uk",
                    vircam_get_license(),
                    vircam_interleave_create,
                    vircam_interleave_exec,
                    vircam_interleave_destroy);

    cpl_pluginlist_append(list,plugin);

    return(0);
}

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

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

    /* Extension number of input frames to use */

    p = cpl_parameter_new_range("vircam.vircam_interleave.extenum",
                                CPL_TYPE_INT,
                                "Extension number to be done, 0 == all",
                                "vircam.vircam_interleave",
                                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_interleave_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_interleave_test(recipe->parameters,recipe->frames));
}
                                
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

static int vircam_interleave_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_interleave_test(cpl_parameterlist *parlist, 
                                  cpl_frameset *framelist) {
    const char *fctid="vircam_interleave";
    int nlab,j,jst,jfn,retval,status,i,nstep,*dims;
    long npts;
    float val;
    double refx,refy,refra,refdec,x,y;
    cpl_parameter *p;
    cpl_image *image;
    cpl_array *a;
    cpl_wcs *wcs;
    

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

    /* Get the parameters */

    p = cpl_parameterlist_find(parlist,"vircam.vircam_interleave.extenum");
    vircam_interleave_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_interleave_tidy();
        return(-1);
    }

    /* Get the frames frames */

    if ((ps.labels = cpl_frameset_labelise(framelist,vircam_compare_tags,
                                           &nlab)) == NULL) {
        cpl_msg_error(fctid,"Cannot labelise the input frames");
        vircam_interleave_tidy();
        return(-1);
    }
    if ((ps.imagelist = vircam_frameset_subgroup(framelist,ps.labels,nlab,
                                                VIRCAM_TEST_SCIENCE_RAW)) == NULL) {
        cpl_msg_error(fctid,"Cannot get images in input frameset");
        vircam_interleave_tidy();
        return(-1);
    }
    ps.nimages = cpl_frameset_get_size(ps.imagelist);
    nstep = (int)sqrt((float)ps.nimages);
    if ((ps.conflist = vircam_frameset_subgroup(framelist,ps.labels,nlab,
                                                VIRCAM_CAL_CONF)) == NULL) {
        cpl_msg_error(fctid,"Cannot get confidence maps in input frameset");
        vircam_interleave_tidy();
        return(-1);
    }
    ps.nconfs = cpl_frameset_get_size(ps.conflist);

    /* 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_interleave_config.extenum,
                       (const cpl_frame *)cpl_frameset_get_frame(ps.imagelist,0),
                       &jst,&jfn);
    if (jst == -1 || jfn == -1) {
        cpl_msg_error(fctid,"Unable to continue");
        vircam_interleave_tidy();
        return(-1);
    }

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

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

        /* Load the images */

        ps.images = vircam_fits_load_list(ps.imagelist,CPL_TYPE_FLOAT,j);
        ps.confs = vircam_fits_load_list(ps.conflist,CPL_TYPE_INT,j);

        /* Get some information that you need for the interleaving. Start
           by getting the background level add it to the extension property
           list */

        for (i = 0; i < ps.nimages; i++) {
            image = vircam_fits_get_image(ps.images[i]);
            npts = vircam_getnpts(image);
            val = vircam_med(cpl_image_get_data_float(image),NULL,npts);
            cpl_propertylist_update_float(vircam_fits_get_ehu(ps.images[i]),
                                          "ESO DRS BACKMED",val);
            wcs = cpl_wcs_new_from_propertylist(vircam_fits_get_ehu(ps.images[i]));
            if (i == 0) {
                a = cpl_wcs_get_image_dims(wcs);
                dims = cpl_array_get_data_int(a);
                refx = 0.5*(double)dims[0];
                refy = 0.5*(double)dims[1];
                vircam_xytoradec(wcs,refx,refy,&refra,&refdec);
            }
            vircam_radectoxy(wcs,refra,refdec,&x,&y);
            x = refx - x;
            y = refy - y;
            cpl_propertylist_update_double(vircam_fits_get_ehu(ps.images[i]),
                                           "ESO DRS XOFFMICRO",x);
            cpl_propertylist_update_double(vircam_fits_get_ehu(ps.images[i]),
                                           "ESO DRS YOFFMICRO",y);
            cpl_wcs_delete(wcs);
        }

        /* Call the interleaf module */

        cpl_msg_info(fctid,"Doing interleaving for extension %d\n",j);
        (void)vircam_interleave(ps.images,ps.nimages,ps.confs,ps.nconfs,
                                nstep,&(ps.plist),&(ps.outimage),&(ps.outconf),
                                &status);
        if (status != VIR_OK) {
            vircam_interleave_tidy();
            return(-1);
        }

        /* Save everything */

        cpl_msg_info(fctid,"Saving combined image extension %d\n",j);
        retval = vircam_interleave_save(framelist,parlist);
        if (retval != 0) {
            vircam_interleave_tidy();
            return(-1);
        }
        freefitslist(ps.images,ps.nimages);
        freefitslist(ps.confs,ps.nconfs);
        freeimage(ps.outimage);
        freeimage(ps.outconf);
        freepropertylist(ps.plist);
    }
    vircam_interleave_tidy();
    return(0);
}


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

static int vircam_interleave_save(cpl_frameset *framelist, 
                                  cpl_parameterlist *parlist) {
    cpl_propertylist *plist;
    const char *recipeid = "vircam_interleave";
    const char *fctid = "vircam_interleave_save";
    const char *outfile = "comb.fits";
    const char *outconf = "combconf.fits";

    /* 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 = cpl_frame_new();
        cpl_frame_set_filename(product_frame,outfile);
        cpl_frame_set_tag(product_frame,VIRCAM_PRO_INTER_TEST);
        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 up header for phu */

        plist = vircam_fits_get_phu(ps.images[0]);
        vircam_dfs_set_product_primary_header(plist,product_frame,framelist,
                                              parlist,(char *)recipeid,
                                              "?Dictionary?");

        /* 'Save' the PHU interleaved 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);
            return(-1);
        }
        cpl_frameset_insert(framelist,product_frame);

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

        product_conf = cpl_frame_new();
        cpl_frame_set_filename(product_conf,outconf);
        cpl_frame_set_tag(product_conf,VIRCAM_PRO_CONF_TEST);
        cpl_frame_set_type(product_conf,CPL_FRAME_TYPE_IMAGE);
        cpl_frame_set_group(product_conf,CPL_FRAME_GROUP_PRODUCT);
        cpl_frame_set_level(product_conf,CPL_FRAME_LEVEL_FINAL);

        /* 'Save' the PHU confidence map image */                        

        if (cpl_image_save(NULL,outconf,CPL_BPP_8_UNSIGNED,plist,
                           CPL_IO_DEFAULT) != CPL_ERROR_NONE) {
            cpl_msg_error(fctid,"Cannot save product PHU");
            cpl_frame_delete(product_conf);
            return(-1);
        }
        cpl_frameset_insert(framelist,product_conf);
    }

    /* Get the extension property list */

    plist = ps.plist;

    /* Fiddle with the header now */

    vircam_dfs_set_product_exten_header(plist,product_frame,framelist,
                                        parlist,(char *)recipeid,
                                        "?Dictionary?");
                
    /* Now save the interleaved image extension */

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

    /* And the confidence map */

    if (cpl_image_save(ps.outconf,outconf,CPL_BPP_16_SIGNED,plist,
                       CPL_IO_EXTEND) != CPL_ERROR_NONE) {
        cpl_msg_error(fctid,"Cannot save confidence map image extension");
        return(-1);
    }

    /* Get out of here */

    return(0);
}

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

static void vircam_interleave_init(void) {
    ps.labels = NULL;
    ps.imagelist = NULL;
    ps.images = NULL;
    ps.conflist = NULL;
    ps.confs = NULL;
    ps.outimage = NULL;
    ps.outconf = NULL;
    ps.plist = NULL;
}

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

static void vircam_interleave_tidy(void) {
    freespace(ps.labels);
    freeframeset(ps.imagelist);
    freefitslist(ps.images,ps.nimages);
    freeframeset(ps.conflist);
    freefitslist(ps.confs,ps.nconfs);
    freeimage(ps.outimage);
    freeimage(ps.outconf);
    freepropertylist(ps.plist);
}

/*

$Log: vircam_interleave.c,v $
Revision 1.17  2008/07/10 13:01:35  jim
Modified to use v4.2 version of cpl_wcs

Revision 1.16  2008/06/20 11:13:02  jim
Fixed dodgy call to cpl_wcs_get_image_dims

Revision 1.15  2008/05/06 08:40:43  jim
Modified to use cpl_wcs interface

Revision 1.14  2007/10/25 19:38:22  jim
modified to keep lint happy

Revision 1.13  2007/10/15 12:53:55  jim
Modified for compatibility with cpl_4.0

Revision 1.12  2007/07/09 13:22:09  jim
Modified to use new version of vircam_exten_range

Revision 1.11  2007/05/02 12:53:11  jim
typo fixes in docs

Revision 1.10  2007/04/13 12:27:39  jim
Added some extra docs

Revision 1.9  2007/04/04 10:36:29  jim
Modified to use new dfs tags

Revision 1.8  2007/03/02 12:38:33  jim
Fixed small memory leak

Revision 1.7  2007/03/01 12:42:59  jim
Modified slightly after code checking

Revision 1.6  2006/07/11 14:59:09  jim
Fixed offsets

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

Revision 1.4  2006/05/15 12:55:42  jim
Fixed a few typos

Revision 1.3  2006/05/04 11:53:42  jim
Fixed _save routine so that it's more consistent with the standard CPL
way of doing things

Revision 1.2  2006/04/27 14:22:05  jim
Fixed docs

Revision 1.1  2006/04/24 10:42:45  jim
New routine


*/

---