vircam_interleave.c

/* $Id: vircam_interleave.c,v 1.16 2007/05/15 08:54:32 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/05/15 08:54:32 $
 * $Revision: 1.16 $
 * $Name:  $
 */

/* Includes */

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

#include <cpl.h>
#include "vircam_utils.h"
#include "vircam_wcsutils.h"
#include "vircam_fits.h"
#include "vircam_mods.h"

static float *odata = NULL;
static int *ocdata = NULL;
static cpl_propertylist *pp = NULL;
static int outok = 0;
static float *xoff = NULL;
static float *yoff = NULL;
static float *backmeds = NULL;

static void tidy(void);

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

extern int vircam_interleave(vir_fits **infiles, int ninputs, 
                             vir_fits **inconfs, int nconfs, int nsteps, 
                             cpl_propertylist **p, cpl_image **outimage, 
                             cpl_image **outconf, int *status) {
    const char *fctid = "vircam_interleave";
    char timestamp[25];
    int i,j,k,nxo,nyo,itx,nx,ny,jindex,joindex;
    int iindex,outindex,*icdata,ival;
    float minoff,minzero,*idata,fval;
    double offs[2],fnm;
    cpl_image *inimg,*icimg;
    cpl_propertylist *ehu,*p2;

    /* Inherited status */

    *outimage = NULL;
    *outconf = NULL;
    *p = NULL;
    if (*status != VIR_OK)
        return(*status);

    /* Look to see how many files there are */

    if (ninputs == 0) {
        cpl_msg_error(fctid,"No input files to interleave");
        tidy();
        FATAL_ERROR
    }

    /* Do we have any input confidence maps? If we don't then just get
       on with life. If there aren't the same number of confidence maps
       as input files, then just use the first one. Otherwise use all
       that are there */

    if (inconfs == NULL) {
        nconfs = 0;
    } else {
        if (nconfs != ninputs) 
            nconfs = 1;
    }

    /* Get the propertylist from the first file in the frameset. We need
       this because we need to reset the WCS for the output frame */

    pp = cpl_propertylist_duplicate(vircam_fits_get_ehu(infiles[0]));

    /* Get the information we need from the extension propertylists. Force
       the offsets to be the exact fraction the should be. */

    xoff = cpl_malloc(ninputs*sizeof(float));
    yoff = cpl_malloc(ninputs*sizeof(float));
    backmeds = cpl_malloc(ninputs*sizeof(float));
    fnm = (double)nsteps;
    for (i = 0; i < ninputs; i++) {
        ehu = vircam_fits_get_ehu(infiles[i]);
        xoff[i] = (float)cpl_propertylist_get_double(ehu,"ESO DRS XOFFMICRO");
        fval = xoff[i] - (int)xoff[i];
        ival = vircam_nint(fval*fnm);
        xoff[i] = (int)xoff[i] + (float)ival/fnm;
        yoff[i] = (float)cpl_propertylist_get_double(ehu,"ESO DRS YOFFMICRO");
        fval = yoff[i] - (int)yoff[i];
        ival = vircam_nint(fval*fnm);
        yoff[i] = (int)yoff[i] + (float)ival/fnm;
        backmeds[i] = cpl_propertylist_get_float(ehu,"ESO DRS BACKMED");
    }

    /* Normalise the x and y offsets so that they have a minimum of zero */

    minoff = xoff[0];
    for (i = 1; i < ninputs; i++)
        minoff = min(minoff,xoff[i]);
    for (i = 0; i < ninputs; i++)
        xoff[i] -= minoff;
    offs[0] = (double)minoff;
    minoff = yoff[0];
    for (i = 1; i < ninputs; i++)
        minoff = min(minoff,yoff[i]);
    for (i = 0; i < ninputs; i++)
        yoff[i] -= minoff;
    offs[1] = (double)minoff;

    /* Normalise the background offsets to a minimum of zero */

    minzero = backmeds[0];
    for (i = 1; i < ninputs; i++)
        minzero = min(minzero,backmeds[i]);
    for (i = 0; i < ninputs; i++)
        backmeds[i] -= minzero;

    /* How big does the output image need to be? */

    nxo = 0;
    nyo = 0;
    for (i = 0; i < ninputs; i++) {
        itx = (int)((float)nsteps*((float)cpl_image_get_size_x(vircam_fits_get_image(infiles[i])) + xoff[i]) + 0.5) - 1;
        nxo = max(nxo,itx);
        itx = (int)((float)nsteps*((float)cpl_image_get_size_y(vircam_fits_get_image(infiles[i])) + yoff[i]) + 0.5) - 1;
        nyo = max(nyo,itx);
    }

    /* Right, get some memory for the output array(s) */

    odata = cpl_calloc(nxo*nyo,sizeof(*odata));
    if (nconfs > 0) 
        ocdata = cpl_calloc(nxo*nyo,sizeof(*ocdata));

    /* Right, now loop through all the images and deposit the pixels into
       the correct location in the output image. First get the data for the
       current image and (if used) the current confidence map */

    for (i = 0; i < ninputs; i++) {
        inimg = vircam_fits_get_image(infiles[i]);
        if (nconfs == 1) 
            icimg = vircam_fits_get_image(inconfs[0]);
        else
            icimg = vircam_fits_get_image(inconfs[i]);
        idata = cpl_image_get_data_float(inimg);
        if (idata == NULL) {
            cpl_msg_error(fctid,"Unable to map data for image %d",i);
            tidy();
            FATAL_ERROR
        }
        if (nconfs == 0) {
            icdata = NULL;
        } else {
            icdata = cpl_image_get_data_int(icimg);
            if (icdata == NULL) {
                cpl_msg_error(fctid,
                              "Unable to map data for confidence map %d",i);
                tidy();
                FATAL_ERROR
            }
        }

        /* Now distribute the data */

        nx = cpl_image_get_size_x(inimg);
        ny = cpl_image_get_size_y(inimg);
        for (j = 0; j < ny; j++) {
            jindex = j*nx;
            joindex = nxo*((int)(((float)j + yoff[i])*fnm + 0.5));
            for (k = 0; k < nx; k++) {
                iindex = jindex + k;
                outindex = joindex + (int)(((float)k + xoff[i])*fnm + 0.5);
                if (outindex < 0 || outindex >= nxo*nyo) {
                    cpl_msg_error(fctid,"Programming error -- output index out of bounds %d\n",outindex);
                    tidy();
                    FATAL_ERROR
                }
                odata[outindex] = idata[iindex] - backmeds[i];
                if (ocdata != NULL)
                    ocdata[outindex] = icdata[iindex];
            }
        }
    }


    /* Update the propertylist to rescale the CD matrix and the CRPIX vector.
       The latter must also be offset by the minimum shift introduced during
       the interleaving */

    fnm = 1.0/fnm;
    if (vircam_crpixshift(pp,fnm,offs) != VIR_OK) {
        tidy();
        FATAL_ERROR
    }
    if (vircam_rescalecd(pp,fnm) != VIR_OK) {
        tidy();
        FATAL_ERROR
    }

    /* Add the provenance information to the output header */

    vircam_prov(pp,infiles,ninputs);

    /* Wrap the data in a cpl_image */

    *outimage = cpl_image_wrap_float(nxo,nyo,odata);
    *outconf = cpl_image_wrap_int(nxo,nyo,ocdata);

    /* Add a timestamp to the primary propertylist */

    vircam_timestamp(timestamp,25);
    p2 = vircam_fits_get_phu(infiles[0]);
    cpl_propertylist_update_string(p2,"VIR_TIME",timestamp);
    cpl_propertylist_set_comment(p2,"VIR_TIME",
                                 "Timestamp for matching to conf map");

    /* Get out of here */

    *p = pp;
    outok = 1;
    tidy();
    GOOD_STATUS
}

static void tidy(void) {
    if (! outok) {
        freespace(odata);
        freespace(ocdata);
        freepropertylist(pp)
    }
    freespace(xoff);
    freespace(yoff);
    freespace(backmeds);
}

/*

$Log: vircam_interleave.c,v $
Revision 1.16  2007/05/15 08:54:32  jim
Decreased size of output map by 1 pixel

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

Revision 1.14  2007/03/22 10:58:46  jim
Force offsets to be the exact fraction that they should be

Revision 1.13  2007/03/01 12:42:41  jim
Modified slightly after code checking

Revision 1.12  2006/11/27 12:04:53  jim
Changed call from cpl_propertylist_append to cpl_propertylist_update

Revision 1.11  2006/10/02 13:47:33  jim
Added missing .h file to include list

Revision 1.10  2006/07/11 14:54:41  jim
The output propertylist is now a duplicate to avoid possible clash when
freeing up memory

Revision 1.9  2006/07/07 09:33:20  jim
Changed datatype of x,y offset properties

Revision 1.8  2006/06/06 13:06:23  jim
Adds VIR_TIME to primary header

Revision 1.7  2006/05/15 12:58:53  jim
fixed a small typo in docs

Revision 1.6  2006/03/23 21:18:50  jim
Minor changes mainly to comment headers

Revision 1.5  2006/03/22 13:58:32  jim
Cosmetic fixes to keep lint happy

Revision 1.4  2006/03/08 14:32:22  jim
Lots of little modifications

Revision 1.3  2006/03/01 10:31:29  jim
Now uses new vir_fits objects

Revision 1.2  2006/02/22 14:09:02  jim
Fixed missing entry in docs

Revision 1.1  2006/02/18 11:53:40  jim
new file


*/

---