vircam_imcombine.c

/* $Id: vircam_imcombine.c,v 1.24 2007/05/14 15:26:03 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/14 15:26:03 $
 * $Revision: 1.24 $
 * $Name:  $
 */

/* Includes */

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

#include <cpl.h>

#include "vircam_utils.h"
#include "vircam_pfits.h"
#include "vircam_stats.h"
#include "vircam_fits.h"
#include "vircam_mods.h"

/* Macro definitions */

#define FATAL_ERR(_f,_a) {cpl_msg_error(_f,_a); tidy(); *status = VIR_FATAL; return(*status);}
#define WARN_ERR(_f,_a) {cpl_msg_error(_f,_a); tidy(); *status = VIR_WARN; return(VIR_WARN);}
#define INFO_ERR(_f,_a) {cpl_msg_info(_f,_a);}
#define MEDIANCALC 1
#define MEANCALC 2
#define SZBUF 1024

/* Definition of litestruct which is used to hold useful information about
   each of the input frames and their associated images */

typedef struct {
    vir_fits         *frame;
    float            exptime;
    float            expfudge;
    float            skylevel;
    float            skynoise;
    float            skyfudge;
    float            skyrenorm;
} litestruct;

/* Static Global variables */

static litestruct *fileptrs = NULL;
static float **datas = NULL;
static float *odata;
static long npts;
static int nf;
static float oskylevel;
static long nx;
static long ny;
static unsigned char *rmask = NULL;
static unsigned char *rplus = NULL;

/* Static subroutine prototypes */

static void skyest(float *data, float thresh, float *skymed, float *skynoise);
static void medcalc(float, float, int);
static void meancalc(float, float, int);
static void xclip_med(float thresh, int scaletype);
static void xclip_mean(float thresh, int scaletype);

static void tidy(void);

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

extern int vircam_imcombine(vir_fits **fset, int nfits, int combtype, 
                            int scaletype, int xrej, float thresh,
                            cpl_image **outimage, unsigned char **rejmask,
                            unsigned char **rejplus, cpl_propertylist **drs,
                            int *status) {
    int i,k,j;
    const char *ic_fctid = "vircam_imcombine";
    char msg[SZBUF];
    float sumsky,sumsig,exp1,exp2,expfudge,skylevel,skynoise,oskynoise;
    float *dat;
    litestruct *ff;
    cpl_propertylist *plist_p;

    /* Inherited status */

    *rejmask = NULL;
    *rejplus = NULL;
    *drs = NULL;
    *outimage = NULL;
    if (*status != VIR_OK) 
        return(*status);

    /* Check that there are any files in the first place...*/ 

    nf = nfits;
    if (nf == 0) 
        WARN_ERR(ic_fctid,"No files to combine")

    /* Get some file structures */

    fileptrs = cpl_calloc(nf,sizeof(litestruct));

    /* Get workspace for convenience array */

    datas = cpl_malloc(nf*sizeof(float*));
    npts = vircam_getnpts(vircam_fits_get_image(fset[0]));
    for (k = 0; k < nf; k++)
        datas[k] = cpl_malloc(npts*sizeof(float));

    /* Get pointers to the data arrays */

    for (k = 0; k < nf; k++) {

        dat = cpl_image_get_data_float(vircam_fits_get_image(fset[k]));
        if (dat == NULL) {
            snprintf(msg,SZBUF,"Failed to load data from extension %d in %s",
                     vircam_fits_get_nexten(fset[k]),
                     vircam_fits_get_filename(fset[k]));
            FATAL_ERR(ic_fctid,msg)
        }
        for (i = 0; i < npts; i++)
            datas[k][i] = dat[i];
    }

    /* Open each file in turn and fill in the necessary information. Start with
       the file name...*/

    sumsky = 0.0;
    sumsig = 0.0;
    for (i = 0; i < nf; i++) {
        ff = fileptrs + i;
        ff->frame = fset[i];

        /* If this is the first frame, then keep the size of the data
           array for future reference */

        if (i == 0) {
            nx = cpl_image_get_size_x(vircam_fits_get_image(fset[0]));
            ny = cpl_image_get_size_y(vircam_fits_get_image(fset[0]));
            npts = nx*ny;
        }

        /* Get the header and the exposure time */

        plist_p = vircam_fits_get_phu(ff->frame);
        if (plist_p == NULL) {
            snprintf(msg,SZBUF,"Failed to load primary property list %s",
                     vircam_fits_get_filename(ff->frame));
            INFO_ERR(ic_fctid,msg)
            exp2 = 1.0;
        } else {
            if (vircam_pfits_get_exptime(plist_p,&exp2) != VIR_OK) {
                exp2 = 1.0;
                snprintf(msg,SZBUF,"Failed to get exposure time from %s",
                         vircam_fits_get_filename(ff->frame));
                INFO_ERR(ic_fctid,msg)
            }
        }

        /* Set a few properties */

        ff->exptime = exp2;
        exp1 = fileptrs->exptime;
        expfudge = exp1/exp2;
        ff->expfudge = expfudge;

        /* If scaling by relative exposure time, then do it now. NB: This
           isn't necessary for the first file as all the others are scaled
           relative to it */

        if (scaletype == 3 && i > 0) {
            for (j = 0; j < npts; j++) 
                datas[i][j] *= ff->expfudge;
        }

        /* Get the background estimate and noise */

        skyest(datas[i],thresh,&skylevel,&skynoise);
        ff->skylevel = skylevel;
        ff->skynoise = skynoise;
        sumsky += skylevel;
        sumsig += skynoise;
    }

    /* Work out average background and noise. Then create background zeropoint
       or scale factor, depending upon which was requested in the call */

    sumsky /= (float)nf;
    sumsig /= (float)nf;
    switch (scaletype) {
    case 1:
        for (i = 0; i < nf; i++)
            (fileptrs+i)->skyfudge = sumsky - (fileptrs+i)->skylevel;
        break;
    case 2:
        for (i = 0; i < nf; i++)
            (fileptrs+i)->skyfudge = sumsky/(fileptrs+i)->skylevel;
        break;      
    case 3:
        for (i = 0; i < nf; i++)
            (fileptrs+i)->skyfudge = sumsky - (fileptrs+i)->skylevel;
        break;
    default:
        for (i = 0; i < nf; i++) 
            (fileptrs+i)->skyfudge = 0.0;
        break;
    }

    /* Open an output image based on the first frame in the frameset */

    *outimage = cpl_image_new(nx,ny,CPL_TYPE_FLOAT);
    odata = cpl_image_get_data_float(*outimage);
    if (*outimage == NULL || odata == NULL) 
        FATAL_ERR(ic_fctid,"Couldn't create output image")
    *rejmask = cpl_calloc(npts,sizeof(*rejmask));
    rmask = *rejmask;
    *rejplus = cpl_calloc(npts,sizeof(*rejplus));
    rplus = *rejplus;

    /* Now do the averaging/medianing */

    switch (combtype) {
    case MEDIANCALC:
        medcalc(thresh,sumsig,scaletype);
        break;
    case MEANCALC:
        meancalc(thresh,sumsig,scaletype);
        break;
    }

    /* Do the extra clipping here if you want it */

    if (xrej) {
        
        /* First get sky background and sigma from output data */

        skyest(odata,thresh,&oskylevel,&oskynoise);

        /* Now loop for all the files subtract off the mean frame (suitably
           scaled and zero pointed depending on what was done in the first 
           place) */

        for (i = 0; i < nf; i++) {
            ff = fileptrs + i;
            ff->skyrenorm = ff->skylevel/oskylevel;
            switch (scaletype) {
            case 1:
                for (k = 0; k < npts; k++) 
                    datas[i][k] -= (odata[k] - oskylevel);
                break;
            case 2:
                for (k = 0; k < npts; k++)
                    datas[i][k] -= (odata[k] - oskylevel)*ff->skyrenorm;
                break;
            case 3:
                for (k = 0; k < npts; k++) 
                    datas[i][k] -= (odata[k] - oskylevel);
                break;
            case 0:
                for (k = 0; k < npts; k++) 
                    datas[i][k] -= (odata[k] - oskylevel);
                break;
            }

            /* Re-estimate the noise for this image */

            skyest(datas[i],thresh,&skylevel,&skynoise);
            ff->skynoise = skynoise;
        }

        /* Now do the extra clip... */

        switch (combtype) {
        case MEDIANCALC:
            xclip_med(thresh,scaletype);
            break;
        case MEANCALC:
            xclip_mean(thresh,scaletype);
            break;
        }
    }

    /* Write provenance keywords */

    *drs = cpl_propertylist_new();
    vircam_prov(*drs,fset,nfits);

    /* Right, tidy and get out of here */

    tidy();
    return(VIR_OK);
}

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

static void xclip_med(float thresh, int scaletype) {
    int nf1,nf2,nfm,nrejmax,is_even,k,is_even2,nrej,nremain,nm,nmm,nplus;
    int nminus;
    float **work,*dork,value,cliplev;
    long i;
    int j;
    litestruct *ff;

    /* Set up a few useful variables */

    nf1 = nf/2 - 1;
    nf2 = nf1 + 1;
    nfm = (nf + 1)/2 - 1;
    nrejmax = nf/2;
    is_even = !(nf & 1);

    /* Get some workspace */

    work = cpl_malloc(3*sizeof(float *));
    for (i = 0; i < 3; i++)
        work[i] = cpl_malloc(nf*sizeof(float));
    dork = cpl_malloc(nf*sizeof(float));

    /* Loop for each input pixel now... */

    for (i = 0; i < npts; i++) {

        /* Scale or shift data */

        switch (scaletype) {
        case 0:
            for (k = 0; k < nf; k++) {
                ff = fileptrs + k;
                work[0][k] = datas[k][i];
                work[1][k] = ff->skynoise;
                work[2][k] = datas[k][i] + odata[i] - oskylevel;
            }
            break;
        case 1:
            for (k = 0; k < nf; k++) {
                ff = fileptrs + k;
                work[0][k] = datas[k][i] + ff->skyfudge;
                work[1][k] = ff->skynoise;
                work[2][k] = datas[k][i] + odata[i] - oskylevel + ff->skyfudge;
            }
            break;
        case 2:
            for (k = 0; k < nf; k++) {
                ff = fileptrs + k;
                work[0][k] = datas[k][i]*ff->skyfudge;
                work[1][k] = ff->skynoise*ff->skyfudge;
                work[2][k] = (datas[k][i] + odata[i]*ff->skyrenorm - 
                              ff->skylevel)*ff->skyfudge;
            }
            break;
        case 3:
            for (k = 0; k < nf; k++) {
                ff = fileptrs + k;
                work[0][k] = datas[k][i] + ff->skyfudge;
                work[1][k] = ff->skynoise;
                work[2][k] = datas[k][i] + odata[i] - oskylevel + ff->skyfudge;
            }
            break;
        }       

        /* Sort and get a first pass median */

        vircam_sort(work,nf,3);
        if (is_even)
            value = 0.5*(work[0][nf1] + work[0][nf2]);
        else 
            if (nf < 5) 
                value = work[0][nfm];
            else 
                value = 0.25*(work[0][nfm-1] + work[0][nfm+1]) + 0.5*work[0][nfm];
        
        /* Do clipping */

        nplus = 0;
        cliplev = value + thresh*work[1][nf-1];
        while (nplus < nrejmax && work[0][nf-nplus-1] > cliplev)
            nplus++;
        nminus = 0;
        cliplev = value - thresh*work[1][nf-1];
        while ((nplus+nminus) < nrejmax && work[0][nminus] < cliplev)
            nminus++;
        nrej = nplus + nminus;

        /* If there were any clipped out, the re-estimate the value */

        if (nrej > 0) {
            nremain = nf - nrej;
            if (nremain != 0) {
                nm = nremain/2 - 1;
                for (j = 0; j < nremain; j++) 
                    dork[j] = work[2][j+nminus];
                nmm = (nremain + 1)/2 - 1;
                is_even2 = !(nremain & 1);
                vircam_sort(&dork,nm,1);
                if (is_even2) 
                    value = 0.5*(dork[nm] + dork[nm+1]);
                else 
                    if (nremain < 3) 
                        value = dork[nmm];
                    else 
                        value = 0.5*dork[nmm] + 0.25*(dork[nmm-1] + dork[nmm+1]);
            }
        
            /* Store the result away */

            odata[i] = value;
            rmask[i] = min(255,nrej);
            rplus[i] = min(255,nplus);
        } else {
            rmask[i] = 0;
            rplus[i] = 0;
        }
    }

    /* Ditch workspace and get out of here */

    for (i = 0; i < 3; i++)
        cpl_free(work[i]);
    cpl_free(work);
    cpl_free(dork);
}

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

static void xclip_mean(float thresh, int scaletype) {
    int k,nf2,nrej,nplus;
    float *work[3],value,cliplev1,cliplev2,value2;
    long i;
    litestruct *ff;

    /* Get some workspace */

    for (i = 0; i < 3; i++)
        work[i] = cpl_malloc(nf*sizeof(float));

    /* Loop for each input pixel now... */

    for (i = 0; i < npts; i++) {

        /* Scale or shift data */

        switch (scaletype) {
        case 0:
            for (k = 0; k < nf; k++) {
                ff = fileptrs + k;
                work[0][k] = datas[k][i];
                work[1][k] = ff->skynoise;
                work[2][k] = datas[k][i] + odata[i] - oskylevel;
            }
            break;
        case 1:
            for (k = 0; k < nf; k++) {
                ff = fileptrs + k;
                work[0][k] = datas[k][i] + ff->skyfudge;
                work[1][k] = ff->skynoise;
                work[2][k] = datas[k][i] + odata[i] - oskylevel + ff->skyfudge;
            }
            break;
        case 2:
            for (k = 0; k < nf; k++) {
                ff = fileptrs + k;
                work[0][k] = datas[k][i]*ff->skyfudge;
                work[1][k] = ff->skynoise*ff->skyfudge;
                work[2][k] = (datas[k][i] + odata[i]*ff->skyrenorm - 
                           ff->skylevel)*ff->skyfudge;
            }
            break;
        case 3:
            for (k = 0; k < nf; k++) {
                ff = fileptrs + k;
                work[0][k] = datas[k][i] + ff->skyfudge;
                work[1][k] = ff->skynoise;
                work[2][k] = datas[k][i] + odata[i] - oskylevel + ff->skyfudge;
            }
            break;
        }       

        /* Get a first pass mean */

        value = 0.0;
        for (k = 0; k < nf; k++) 
            value += work[0][k];
        value /= (float)nf;
        
        /* Do upper clipping */

        value2 = 0.0;
        nplus = 0;
        nf2 = 0;
        for (k = 0; k < nf; k++) {
            cliplev1 = value - thresh*work[1][k];
            cliplev2 = value + thresh*work[1][k];
            if (work[0][k] >= cliplev1 && work[0][k] <= cliplev2) {
                value2 += work[2][k];
                nf2++;
            } else if (work[0][k] > cliplev2) {
                nplus++;
            }
        }

        /* If there were any clipped out, the re-estimate the value */

        if (nf2 < nf) {
            if (nf2 != 0)
                value = value2/(float)nf2;
        
            /* Store the result away */

            odata[i] = value;
        }
        nrej = nf - nf2;
        rmask[i] = min(255,nrej);       
        rplus[i] = min(255,nplus);
    }

    /* Ditch workspace and get out of here */

    for (k = 0; k < 3; k++)
        cpl_free(work[k]);
}


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

static void medcalc(float thresh, float avskynoise, int scaletype) {
    int nf1,nf2,nfm,nrejmax,is_even,nrej,nremain,nm,nmm,is_even2,k,nminus;
    int nplus;
    long i;
    float value,cliplev,*work;

    /* Set up a few useful variables */

    nf1 = nf/2 - 1;
    nf2 = nf1 + 1;
    nfm = (nf + 1)/2 - 1;
    nrejmax = nf/2;
    is_even = !(nf & 1);

    /* Get a workspace */

    work = cpl_malloc(nf*sizeof(*work));

    /* Ok, loop for each pixel... */

    for (i = 0; i < npts; i++) {

        /* Scale or shift data */

        switch (scaletype) {
        case 0:
            for (k = 0; k < nf; k++) 
                work[k] = datas[k][i];
            break;
        case 1:
            for (k = 0; k < nf; k++) 
                work[k] = datas[k][i] + (fileptrs+k)->skyfudge;
            break;
        case 2:
            for (k = 0; k < nf; k++) 
                work[k] = datas[k][i]*(fileptrs+k)->skyfudge;
            break;
        case 3:
            for (k = 0; k < nf; k++) 
                work[k] = datas[k][i] + (fileptrs+k)->skyfudge;
            break;
        }
        
        /* Sort data and get the median */

        vircam_sort(&work,nf,1);
        if (is_even)
            value = 0.5*(work[nf1] + work[nf2]);
        else 
            if (nf < 5) 
                value = work[nfm];
            else 
                value = 0.25*(work[nfm-1] + work[nfm+1]) + 0.5*work[nfm];

        /* Enter a rejection loop */

        nplus = 0;
        cliplev = value + thresh*avskynoise;
        while (nplus < nrejmax && work[nf-nplus-1] > cliplev)
            nplus++;
        nminus = 0;
        cliplev = value - thresh*avskynoise;
        while ((nplus+nminus) < nrejmax && work[nminus] < cliplev)
            nminus++;
        nrej = nplus + nminus;

        /* If you've clipped any, then recalculate the median...*/

        if (nrej > 0) {
            nremain = nf - nrej;
            nm = nremain/2 - 1 + nminus;
            nmm = (nremain + 1)/2 - 1 + nminus;
            is_even2 = !(nremain & 1);
            if (is_even2) 
                value = 0.5*(work[nm] + work[nm+1]);
            else 
                if (nremain < 3) 
                    value = work[nmm];
                else 
                    value = 0.5*work[nmm] + 0.25*(work[nmm-1] + work[nmm+1]);
        }

        /* Store the result away */

        odata[i] = value;
        rmask[i] = min(255,nrej);
        rplus[i] = min(255,nplus);
    }

    /* Get rid of workspace */

    cpl_free(work);
}

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

static void meancalc(float thresh, float avskynoise, int scaletype) {
    int nf2,k,nrej,nplus;
    long i;
    float *work,value,cliplev1,cliplev2,value2;

    /* Get a vector for workspace */

    work = cpl_malloc(nf*sizeof(*work));

    /* Ok, loop for each pixel... */

    for (i = 0; i < npts; i++) {

        /* Scale or shift data */

        switch (scaletype) {
        case 0:
            for (k = 0; k < nf; k++) 
                work[k] = datas[k][i];
            break;
        case 1:
            for (k = 0; k < nf; k++) 
                work[k] = datas[k][i] + (fileptrs+k)->skyfudge;
            break;
        case 2:
            for (k = 0; k < nf; k++) 
                work[k] = datas[k][i]*(fileptrs+k)->skyfudge;
            break;
        case 3:
            for (k = 0; k < nf; k++) 
                work[k] = datas[k][i] + (fileptrs+k)->skyfudge;
            break;
        }
        
        /* Get the mean */
        
        value = 0.0;
        for (k = 0; k < nf; k++) 
            value += work[k];
        value /= (float)nf;

        /* Enter a rejection loop  */

        value2 = 0.0;
        nf2 = 0;
        nplus = 0;
        cliplev1 = value - thresh*avskynoise;
        cliplev2 = value + thresh*avskynoise;
        for (k = 0; k < nf; k++) {
            if (work[k] >= cliplev1 && work[k] <= cliplev2) {
                value2 += work[k];
                nf2 += 1;
            } else if (work[k] > cliplev2) {
                nplus++;
            }
        }

        /* If you've clipped any, then recalculate the mean...*/

        if (nf2 < nf && nf2 != 0)
            value = value2/(float)nf2;
        nrej = nf - nf2;

        /* Store the result away */

        odata[i] = value;
        rmask[i] = min(255,nrej);
        rplus[i] = min(255,nplus);
    }

    /* Get rid of workspace */

    cpl_free(work);
}
            
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

static void skyest(float *data, float thresh, float *skymed, float *skynoise) {
    unsigned char *bpm;

    /* Get a dummy bad pixel mask */

    bpm = cpl_calloc(npts,sizeof(*bpm));

    /* Get the stats */

    vircam_qmedsig(data,bpm,npts,thresh,3,-65535.0,65535.0,skymed,
                   skynoise);

    /* Clean up */

    cpl_free(bpm);
}

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

static void tidy(void) {
    int i;

    /* Free up work space associated with file structures */

    freespace(fileptrs);
    for (i = 0; i < nf; i++)
        freespace(datas[i]);
    freespace(datas);
}

/*

$Log: vircam_imcombine.c,v $
Revision 1.24  2007/05/14 15:26:03  jim
Fixed bug in multiplicative scaling

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

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

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

Revision 1.20  2006/09/08 09:21:37  jim
Fixed bug in median routine

Revision 1.19  2006/08/21 09:08:10  jim
Fixed a few subtle bugs in the median routines

Revision 1.18  2006/07/04 09:19:05  jim
replaced all sprintf statements with snprintf

Revision 1.17  2006/05/08 10:56:14  jim
Fixed bug where exposure time scaling was not happening to the data copy

Revision 1.16  2006/04/20 11:26:45  jim
Fixed so that original data isn't modified in the input vir_fits list.

Revision 1.15  2006/03/23 21:18:48  jim
Minor changes mainly to comment headers

Revision 1.14  2006/03/22 13:58:31  jim
Cosmetic fixes to keep lint happy

Revision 1.13  2006/03/08 14:32:21  jim
Lots of little modifications

Revision 1.12  2006/03/03 14:29:46  jim
Modified definition of vir_fits and channel table

Revision 1.11  2006/03/01 10:31:28  jim
Now uses new vir_fits objects

Revision 1.10  2006/02/22 10:09:09  jim
Added status variable to call

Revision 1.9  2006/01/23 10:30:49  jim
Mainly documentation mods

Revision 1.8  2005/12/14 22:17:33  jim
Updated docs

Revision 1.7  2005/11/29 11:54:19  jim
Modified call in skyest to allow for larger negative numbers as it is
likely that the reset frames will be negative.

Revision 1.6  2005/11/25 09:56:15  jim
Tidied up some more documentation

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

Revision 1.4  2005/11/07 13:15:16  jim
Fixed lots of bugs and added some error checking

Revision 1.3  2005/11/03 13:28:49  jim
All sorts of changes to tighten up error handling

Revision 1.2  2005/10/14 13:17:54  jim
Now spews out a mask of rejected pixels and rejected positive pixels
(indicative of cosmic ray hits)

Revision 1.1.1.1  2005/08/05 08:29:09  jim
Initial import


*/

---