casu_pawsky_minus.c

/* $Id: casu_pawsky_minus.c,v 1.4 2015/11/25 10:26:31 jim Exp $
 *
 * This file is part of the CASU Pipeline utilities
 * Copyright (C) 2015 European Southern Observatory
 *
 * 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: 2015/11/25 10:26:31 $
 * $Revision: 1.4 $
 * $Name:  $
 */
 
/* Includes */
 
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
 
#include <cpl.h>
#include <math.h>
#include <string.h>
 
#include "casu_sky.h"
#include "casu_mask.h"
#include "catalogue/casu_fits.h"
#include "casu_mods.h"
#include "catalogue/casu_utils.h"
#include "casu_stats.h"
#include "catalogue/casu_wcsutils.h"
 
#define DATAMIN -1000.0
#define DATAMAX 65535.0
 
static void masksky_zeros(float **datas, unsigned char **masks, int nfiles,
                          int npts, float **zeros);
static void combine(int nfiles, int stat, int npts, float **datas, float **vars,
                    unsigned char **masks, float **skyout, float **skyv,
                    unsigned char **skyout_bpm);
static void combine_mult(int nfiles, int stat, int npts, float **datas, 
                         float **vars, unsigned char **masks, float *skyout,
                         float *skyvar);
static void domed(float *buf, int n, float *val);

 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
extern int casu_pawsky_minus(casu_fits **infiles, casu_fits **invar,
                             casu_fits *conf, casu_fits *objmaskfits, 
                             int nfiles, casu_fits **skyout, casu_fits **skyvar,
                             int *status) {
    cpl_wcs *wcsmask,*wcsimg;
    unsigned char **masks,*sky_bpm;
    cpl_image *im,*skyim,*skyimv;
    double *xin = NULL, *yin = NULL, *xout = NULL, *yout = NULL, *dx = NULL, *dy = NULL, ddx, ddy;
    int nx,ny,npts,ind,i,j,kind,ix,iy,jx,jy,*opm = NULL;
    int nx_mask = -1; /* return value for cpl_image_get_size_x on error */
    int nfiles2, *confdata = NULL;
    float **datas,*sky,*skyv,*zeros,val,sig,**vars;
    cpl_propertylist *plist;
    casu_fits **infiles2,**invar2;
    const char *fctid = "casu_pawsky_minus";
 
    /* Inherited status */
 
    *skyout = NULL;
    *skyvar = NULL;
    if (*status != CASU_OK) 
        return(*status);
 
    /* If there aren't any images, then get out of here */
 
    if (nfiles == 0) {
        cpl_msg_error(fctid,"Sky correction impossible. No science frames");
        return(CASU_FATAL);
    }
 
    /* Check to see which files are good */
 
    infiles2 = cpl_malloc(nfiles*sizeof(casu_fits *));
    if (invar != NULL) 
        invar2 = cpl_malloc(nfiles*sizeof(casu_fits *));
    else
        invar2 = NULL;
    nfiles2 = 0;
    for (i = 0; i < nfiles; i++) {
        if (casu_fits_get_status(infiles[i]) == CASU_OK) 
            infiles2[nfiles2] = infiles[i];
        if (invar != NULL) 
            invar2[nfiles2] = invar[i];
        nfiles2++;
    }
 
    /* If no good images are present then wrap a dummy image and return it */
 
    if (nfiles2 == 0) {
        skyim = casu_dummy_image(infiles[0]);
        *skyout = casu_fits_wrap(skyim,infiles[0],NULL,NULL);
        casu_dummy_property(casu_fits_get_ehu(*skyout));
        casu_fits_set_status(*skyout,CASU_FATAL);
        if (invar != NULL) {
            skyimv = casu_dummy_image(infiles[0]);
            *skyvar = casu_fits_wrap(skyimv,invar[0],NULL,NULL);
            casu_dummy_property(casu_fits_get_ehu(*skyvar));
        } else {
            *skyvar = NULL;
        }
        cpl_msg_warning(fctid,"No good images in input list");
        freespace(infiles2);
        *status = CASU_WARN;
        return(CASU_WARN);
    }
 
    /* Get space for arrays */
 
    datas = cpl_malloc(nfiles2*sizeof(float *));
    if (invar != NULL)
        vars = cpl_malloc(nfiles2*sizeof(float *));
    else 
        vars = NULL;
    masks = cpl_malloc(nfiles2*sizeof(unsigned char *));
    im = casu_fits_get_image(infiles2[0]);
    nx = cpl_image_get_size_x(im);
    ny = cpl_image_get_size_y(im);
    npts = nx*ny;
 
    /* Info about object mask if it exists */
 
    if (objmaskfits != NULL) {
        nx_mask = cpl_image_get_size_x(casu_fits_get_image(objmaskfits));
        wcsmask = cpl_wcs_new_from_propertylist(casu_fits_get_ehu(objmaskfits));
        opm = cpl_image_get_data_int(casu_fits_get_image(objmaskfits));
 
        /* We need this space so that we can work out which part of the 
           input mask is relevant for each input image */
 
        xin = cpl_malloc(npts*sizeof(double));
        yin = cpl_malloc(npts*sizeof(double));
        xout = cpl_malloc(npts*sizeof(double));
        yout = cpl_malloc(npts*sizeof(double));
        dx = cpl_malloc(nfiles2*sizeof(double));
        dy = cpl_malloc(nfiles2*sizeof(double));
 
        /* Initialise the input xy arrays */
 
        ind = 0;
        for (j = 0; j < ny; j++) {
            for (i = 0; i < nx; i++) {
                xin[ind] = (double)(i+1);
                yin[ind++] = (double)(j+1);
            }
        }
 
        /* Define an output grid for the first image and then offsets for
           the subsequent images */
 
        ddx = 1.0;
        ddy = 1.0;
        for (i = 0; i < nfiles2; i++) {
            wcsimg = cpl_wcs_new_from_propertylist(casu_fits_get_ehu(infiles2[i]));
            if (i == 0) 
                casu_xytoxy_list(wcsimg,wcsmask,npts,xin,yin,xout,yout);
            casu_xytoxy_list(wcsimg,wcsmask,1,&ddx,&ddy,dx+i,dy+i);
            dx[i] -= xout[0];
            dy[i] -= yout[0];
            cpl_wcs_delete(wcsimg);
        }
        cpl_wcs_delete(wcsmask);
    }
 
    /* Loop for each file and get a reference for the data array. Get
       some space for a mask that will cover both the zeroed pixels in the
       confidence maps and any object mask that might be used */
 
    for (i = 0; i < nfiles2; i++) {
        im = casu_fits_get_image(infiles2[i]);
        datas[i] = cpl_image_get_data_float(im);
        masks[i] = cpl_calloc(npts,sizeof(unsigned char));
        if (i == 0) 
            confdata = cpl_image_get_data_int(casu_fits_get_image(conf));
        if (invar != NULL) 
            vars[i] = cpl_image_get_data_float(casu_fits_get_image(invar2[i]));
 
        /* Now loop for each pixel. If the confidence map is zero here then
           mark this as a bad pixel. If not, then look at the master mask
           and see which is its closest pixel to the current one. NB: the 
           -0.5 in the index is a combination of +0.5 to do a nint and a -1
           to take account of the fact that arrays start from index zero. */
 
        for (jy = 0; jy < ny; jy++) {
            for (jx = 0; jx < nx; jx++) {
                ind = jy*nx + jx;
                if (confdata != NULL && confdata[ind] == 0) {
                    masks[i][ind] = 1;
                } else if (objmaskfits != NULL) {
                    kind = (int)(yout[ind] + dy[i] - 0.5)*nx_mask + 
                        (int)(xout[ind] + dx[i] - 0.5);
                    if(opm != NULL) {
                        masks[i][ind] = opm[kind];
                    }
                }
            }
        }
    }
 
    /* Do an intermediate tidy */
 
    if (objmaskfits != NULL) {
        freespace(xin);
        freespace(yin);
        freespace(xout);
        freespace(yout);
        freespace(dx);
        freespace(dy);
    }
 
    /* Offset the DC level */
    
    masksky_zeros(datas,masks,nfiles2,npts,&zeros);
 
    /* Right, combine these now */
 
    combine(nfiles2,0,npts,datas,vars,masks,&sky,&skyv,&sky_bpm);
    
    /* Wrap the sky and reject the pixels that had no sky contribution */
 
    skyim = cpl_image_wrap_float(nx,ny,sky);
    for (i = 0; i < npts; i++) {
        iy = i/nx + 1;
        ix = npts - (iy-1)*nx;
        if (sky_bpm[i])
            cpl_image_reject(skyim,(cpl_size)ix,(cpl_size)iy);
    }
    freespace(sky_bpm);
    *skyout = casu_fits_wrap(skyim,infiles2[0],NULL,NULL);
    plist = casu_fits_get_ehu(*skyout);
    if (objmaskfits != NULL) {
        cpl_propertylist_update_string(plist,"ESO DRS MASKUSED",
                                       casu_fits_get_filename(objmaskfits));
        cpl_propertylist_set_comment(plist,"ESO DRS MASKUSED",
                                     "Object masked used to make sky");
    }
    cpl_propertylist_update_string(plist,"ESO DRS SKYALGO","pawsky_minus");
    cpl_propertylist_set_comment(plist,"ESO DRS SKYALGO",
                                 "Sky estimation algorithm");
    casu_prov(plist,infiles2,nfiles2,1);
    if (vars != NULL) {
        skyimv = cpl_image_wrap_float(nx,ny,skyv);
        *skyvar = casu_fits_wrap(skyimv,invar2[0],NULL,NULL);
    } else {
        *skyvar = NULL;
    }
 
    /* Fill in the bits with no sky contribution */
 
    casu_inpaint(*skyout,64,status);
 
    /* Do the multiple combine and subtraction */
 
    combine_mult(nfiles2,0,npts,datas,vars,masks,sky,skyv);
 
    /* Now add the DC level back onto the data */
 
    for (i = 0; i < nfiles2; i++) {
        casu_qmedsig(datas[i],masks[i],(long)npts,3.0,1,DATAMIN,DATAMAX,&val,
                       &sig);
        val = zeros[i] - val;
        for (j = 0; j < npts; j++) 
            datas[i][j] += val;
    }
    freespace(zeros);
 
    /* Add a little something to the header and start cleaning up */
 
    for (i = 0; i < nfiles2; i++) {
        plist = casu_fits_get_ehu(infiles2[i]);
        cpl_propertylist_update_string(plist,"ESO DRS SKYSUB",
                                       "Done with pawsky_minus");
        freespace(masks[i]);
    }
    freespace(masks);
    freespace(datas);
    freespace(vars);
    freespace(infiles2);
    freespace(invar2);
    return(CASU_OK);
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
static void masksky_zeros(float **datas, unsigned char **masks, int nfiles,
                          int npts, float **zeros) {
    int i,j;
    float sig,off,medval;
 
    /* Get some space for the zeros */
 
    *zeros = cpl_malloc(nfiles*sizeof(float));
 
    /* Loop for each input image and get the background median */
 
    for (i = 0; i < nfiles; i++) 
        casu_qmedsig(datas[i],masks[i],(long)npts,3.0,1,DATAMIN,
                     DATAMAX,*zeros+i,&sig);
 
    /* Get the median value of the array */
 
    medval = casu_med(*zeros,NULL,nfiles);
 
    /* Now work out the offset for each image and subtract it */
 
    for (i = 0; i < nfiles; i++) {
        off = medval - (*zeros)[i];
        for (j = 0; j < npts; j++)
            datas[i][j] += off;
    }
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
static void combine(int nfiles, int stat, int npts, float **datas, float **vars,
                    unsigned char **masks, float **skyout, float **skyv,
                    unsigned char **skyout_bpm) {
    int i,j,n;
    float *buf,sum;
 
    /* Get workspace for output arrays */
 
    *skyout = cpl_malloc(npts*sizeof(float));
    if (vars != NULL) 
        *skyv = cpl_malloc(npts*sizeof(float));
    else
        *skyv = NULL;
    *skyout_bpm = cpl_malloc(npts*sizeof(unsigned char));
    buf = cpl_malloc(nfiles*sizeof(float));
 
    /* Loop now for each input pixel and form the median */
 
    for (j = 0; j < npts; j++) {
        n = 0;
        sum = 0.0;
        for (i = 0; i < nfiles; i++) {
            if (masks[i][j] == 0) {
                buf[n] = datas[i][j];
                if (vars != NULL) 
                    sum += vars[i][j];
                n++;
            }
        }
        if (n == 0) {
            (*skyout)[j] = 0.0;
            if (vars != NULL)
                (*skyv)[j] = 0.0;
            (*skyout_bpm)[j] = 1;
        } else {
            if (vars != NULL)
                sum /= powf((float)n,2.0);
            if (stat == 0) {
                domed(buf,n,(*skyout)+j);
                if (vars != NULL)
                    (*skyv)[j] = sum*CPL_MATH_PI_2;
            } else {
                (*skyout)[j] = casu_mean(buf,NULL,n);
                if (vars != NULL) 
                    (*skyv)[j] = sum;
            }
            (*skyout_bpm)[j] = 0;
        }
    }
 
    /* Free some workspace and get out of here */
 
    freespace(buf);
    return;
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
static void combine_mult(int nfiles, int stat, int npts, float **datas, 
                         float **vars, unsigned char **masks, float *skyout,
                         float *skyvar) {
    int i,j,k,n;
    float *buf,*bufv = NULL,*val,*valv = NULL;
 
    /* Get some memory */
 
    buf = cpl_malloc(nfiles*sizeof(float));
    val = cpl_malloc(nfiles*sizeof(float));
    if (vars != NULL) {
        bufv = cpl_malloc(nfiles*sizeof(float));
        valv = cpl_malloc(nfiles*sizeof(float));
    }
 
    /* Loop now for each input pixel */
 
    for (j = 0; j < npts; j++) {
 
        /* Loop for each file. If the current pixel in the current file
           was flagged, then the sky for that pixel will be the same as
           the general sky we did before. If not, then remove that file's
           pixel and recalculate the sky */
 
        for (i = 0; i < nfiles; i++) {
            if (masks[i][j] != 0) {
                val[i] = skyout[j];
                if (vars != NULL) 
                    valv[i] = skyvar[j];
            } else {
                n = 0;
                for (k = 0; k < nfiles; k++) {
                    if (k == i || masks[k][j] != 0)
                        continue;
                    buf[n] = datas[k][j];
                    if (vars != NULL && vars[i] != NULL)
                        bufv[n++] = vars[k][j];
                }
                if (n == 0) {
                    val[i] = skyout[j];
                    if (vars != NULL)
                        valv[i] = skyvar[j];
                } else {
                    if (stat == 0) {
                        domed(buf,n,val+i);
                        if (vars != NULL) 
                            valv[i] = CPL_MATH_PI_2*casu_mean(bufv,NULL,n);
                    } else {
                        val[i] = casu_mean(buf,NULL,n);
                        if (vars != NULL) 
                            valv[i] = casu_mean(bufv,NULL,n);
                    }
                }
            }
        }
 
        /* Now correct these pixels in the input images */
 
        for (i = 0; i < nfiles; i++) {
            datas[i][j] -= val[i];
            if (vars != NULL) 
                vars[i][j] += valv[i];
        }
    }
 
    /* Free some workspace and get out of here */
 
    freespace(buf);
    freespace(val);
    if (vars != NULL) {
        freespace(bufv);
        freespace(valv);
    }
    return;
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
static void domed(float *buf, int n, float *val) {
    int is_even,n2;
 
    is_even = ! (n & 1);
    casu_sort(&buf,n,1);
    n2 = n/2;
    if (is_even) {
        *val = 0.5*(buf[n2-1] + buf[n2]);
    } else {
        if (n <= 5) 
            *val = buf[n2];
        else
            *val = 0.5*buf[n2] + 0.25*(buf[n2-1] + buf[n2+1]);
    }
}

 
/*
 
$Log: casu_pawsky_minus.c,v $
Revision 1.4  2015/11/25 10:26:31  jim
replaced some hardcoded numbers with defines
 
Revision 1.3  2015/09/11 09:27:53  jim
Fixed some problems with the docs
 
Revision 1.2  2015/08/07 13:06:54  jim
Fixed copyright to ESO
 
Revision 1.1.1.1  2015/06/12 10:44:32  jim
Initial import
 
Revision 1.9  2015/03/12 09:16:51  jim
Modified to remove some compiler moans
 
Revision 1.8  2015/01/29 11:52:36  jim
Modified comments + fixed error in variance propagation
 
Revision 1.7  2015/01/09 12:13:15  jim
*** empty log message ***
 
Revision 1.6  2014/12/11 12:23:33  jim
new version
 
Revision 1.5  2014/04/09 11:08:21  jim
Get rid of a couple of compiler moans
 
Revision 1.4  2014/03/26 15:55:40  jim
Uses floating point confidence
 
Revision 1.3  2013/11/21 09:38:14  jim
detabbed
 
Revision 1.2  2013-10-24 09:25:21  jim
traps for dummy input files and tries to remove them from the analysis
 
Revision 1.1.1.1  2013-08-27 12:07:48  jim
Imported
 
 
*/