casu_imdither.c

/* $Id: casu_imdither.c,v 1.5 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.5 $
 * $Name:  $
 */
 
/* Includes */
 
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
 
#include <math.h>
#include <cpl.h>
#include "casu_mods.h"
#include "catalogue/casu_utils.h"
#include "catalogue/casu_fits.h"
#include "casu_stats.h"
 
typedef struct {
    casu_fits  *fname;
    casu_fits  *conf;
    float     xoff;
    float     yoff;
    int       ixoff;
    int       iyoff;
    int       nx;
    int       ny;
    float     sky;
    float     skydiff;
    float     noise;
    float     expscale;
    float     weight;
    float     *data;
    int       *cdata;
    int       ndata;
} dstrct;
 
typedef struct {
    float         *values;
    float         *confs;
    float         *weights;
    short int     *iff;
    int           n;
    long          outindex;
    unsigned char clipped;
} keeptabs;
 
#define NROWSBUF 512
#define MINDATA -1000.0
#define MAXDATA 65535.0
 
static void average(dstrct *fileptrs, keeptabs *c, float *outdata, 
                    float *outconf, float cliplow, float cliphigh, float lthr, 
                    float hthr, float extra, float sumweight);
static keeptabs *clip_open(int nimages, int nxo);
static void clip_close(keeptabs **c, int nxo);
static void skyest(float *data, int *cdata, long npts, float thresh, 
                   float dmin, float dmax, float *skymed, float *skynoise);
static void tidy(int nxo, dstrct *fileptrs, keeptabs *clipmon, float *owdata);

 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
    
extern int casu_imdither(casu_fits **inf, casu_fits **inconf, int nimages, 
                         int nconfs, float lthr, float hthr,
                         cpl_propertylist **p, const char *expkey,
                         cpl_image **out, cpl_image **outc, int *status) {
 
    int i,itx,iy,ccur,clast,ix,n,iline,icol,jy,jx,j,nxo,nyo,*ocdata;
    long npts,ielm,iloc,index_y,index;
    dstrct *fileptrs=NULL,*dd;
    keeptabs *clipmon=NULL,*c;
    float minxoff,minyoff,expref,sky,skynoise,clip1,clip2,outdata;
    float outconf,avlev,avvar,renorm,exposure,sumweight,*odata,*owdata=NULL;
    double crpix1,crpix2;
    cpl_propertylist *ehu,*p2,*phu;
    const char *fctid = "casu_imdither";
    char timestamp[25];
 
    /* Inherited status */
 
    *out = NULL;
    *outc = NULL;
    *p = NULL;
    if (*status != CASU_OK)
        return(*status);
 
    /* Is there any point in being here? */
 
    if (nimages == 0) {
        cpl_msg_error(fctid,"No input files to combine");
        tidy(0,fileptrs,clipmon,owdata);
        FATAL_ERROR
    }
 
    /* Allocate file struct array and fill in some values */
 
    fileptrs = cpl_malloc(nimages*sizeof(dstrct));
    if (cpl_propertylist_has(casu_fits_get_phu(inf[0]),expkey)) {
        exposure = (float)cpl_propertylist_get_double(casu_fits_get_phu(inf[0]),
                                                      expkey);
    } else {
        exposure = 0.5;
    }
    expref = max(0.5,exposure);
    minxoff = 1.0e10;
    minyoff = 1.0e10;
    for (i = 0; i < nimages; i++) {
        dd = fileptrs + i;
        dd->fname = inf[i];
        dd->data = cpl_image_get_data_float(casu_fits_get_image(inf[i]));
        if (nconfs == 0) {
            dd->conf = NULL;
        } else if (nconfs == 1) {
            dd->conf = inconf[0];
            dd->cdata = cpl_image_get_data_int(casu_fits_get_image(inconf[0]));
        } else {
            dd->conf = inconf[i];
            dd->cdata = cpl_image_get_data_int(casu_fits_get_image(inconf[i]));
        }
        phu = casu_fits_get_phu(dd->fname);
        ehu = casu_fits_get_ehu(dd->fname);
        dd->xoff = cpl_propertylist_get_float(ehu,"ESO DRS XOFFDITHER");
        dd->yoff = cpl_propertylist_get_float(ehu,"ESO DRS YOFFDITHER");
        minxoff = min(dd->xoff,minxoff);
        minyoff = min(dd->yoff,minyoff);
        if (cpl_propertylist_has(phu,expkey)) {
            exposure = (float)cpl_propertylist_get_double(phu,expkey);
            exposure = max(0.5,exposure);
        } else {
            exposure = 0.5;
        }
        dd->expscale = exposure/expref;
 
        /* Now work out a background and background noise estimate */
 
        dd->nx = (int)cpl_image_get_size_x(casu_fits_get_image(dd->fname));
        dd->ny = (int)cpl_image_get_size_y(casu_fits_get_image(dd->fname));
        npts = dd->nx*dd->ny;
        skyest(dd->data,dd->cdata,npts,3.0,MINDATA,MAXDATA,&sky,&skynoise);
        dd->sky = sky;
        dd->noise = skynoise;
        
        /* Double check to make sure the confidence maps and images have the
           same dimensions */
 
        if ((int)cpl_image_get_size_x(casu_fits_get_image(dd->conf)) != dd->nx || 
            (int)cpl_image_get_size_y(casu_fits_get_image(dd->conf)) != dd->ny) {
            cpl_msg_error(fctid,"Image %s and Confidence map %s don't match",
                          casu_fits_get_fullname(dd->fname),
                          casu_fits_get_fullname(dd->conf));
            tidy(0,fileptrs,clipmon,owdata);
            FATAL_ERROR
        }
    }
 
    /* Redo the offsets so that they are all positive. */
 
    for (i = 0; i < nimages; i++) {
        dd = fileptrs + i;
        dd->xoff -= minxoff;
        dd->yoff -= minyoff;
        dd->ixoff = (int)(dd->xoff + 0.5);
        dd->iyoff = (int)(dd->yoff + 0.5);
    }
 
    /* Redo the zero point offsets so that they are all relative to
       the first image in the list. Make sure to divide by the relative
       exposure time first! Set up weights*/
 
    fileptrs->sky /= fileptrs->expscale;
    fileptrs->skydiff = 0.0;
    fileptrs->weight = 1.0;
    sumweight = 1.0;
    for (i = 1; i < nimages; i++) {
        dd = fileptrs + i;
        dd->sky /= dd->expscale;
        dd->skydiff = fileptrs->sky - dd->sky;
        dd->noise /= (float)sqrt((double)dd->expscale);
        dd->weight = (float)(pow((double)fileptrs->noise,2.0)/pow((double)dd->noise,2.0)); 
        sumweight += dd->weight;
    }
 
    /* Scale data (don't do image 0 since that has 1 scale and 0 offset) */
 
    for (i = 1; i < nimages; i++) {
        dd = fileptrs + i;
        npts = dd->nx*dd->ny;
        for (j = 0; j < npts; j++) 
            dd->data[j] = dd->data[j]/dd->expscale + dd->skydiff;
    }
 
    /* Set up clipping levels */
 
    clip1 = fileptrs->sky - lthr*fileptrs->noise;
    clip2 = fileptrs->sky + hthr*fileptrs->noise;
 
    /* Open the output file. First of all work out how big the output map
       needs to be. Then create it based on the first image in the list */
 
    nxo = 0;
    nyo = 0;
    for (i = 0; i < nimages; i++) {
        dd = fileptrs + i;
        itx = dd->nx + dd->ixoff;
        nxo = max(nxo,itx);
        itx = dd->ny + dd->iyoff;
        nyo = max(nyo,itx);
    }
 
    /* Create the output image */
 
    *out = cpl_image_new((cpl_size)nxo,(cpl_size)nyo,CPL_TYPE_FLOAT);
 
    /* If an output confidence map has been specified, then create it now. */
 
    if (nconfs != 0) 
        *outc = cpl_image_new((cpl_size)nxo,(cpl_size)nyo,CPL_TYPE_INT);
    else 
        *outc = NULL;
 
    /* Get the data arrays for the output images */
 
    npts = nxo*nyo;
    odata = cpl_image_get_data_float(*out);
    if (*outc != NULL) 
        owdata = cpl_malloc(npts*sizeof(float));
    clipmon = clip_open(nimages,nxo);
 
    /* Right, now try and do the work. Start by deciding whether for a given
       output line an input line is able to contribute */
 
    for (iy = 0; iy < nyo; iy++) {
        ccur = (iy % 2)*nxo;
        clast = nxo - ccur;
        for (ix = 0; ix < nxo; ix++) {
            c = clipmon + ccur + ix;
            c->n = 0;
            c->clipped = 0;
            n = 0;
            for (i = 0; i < nimages; i++) {             
                dd = fileptrs + i;
                iline = iy - dd->iyoff;
                if (iline < 0 || iline >= dd->ny)
                    continue;
                icol = ix - dd->ixoff;
                if (icol < 0 || icol >= dd->nx)
                    continue;
 
                /* Load up any input data for this pixel from the current 
                   image */
                
                ielm = dd->nx*iline + icol;
                c->values[n] = dd->data[ielm];
                c->confs[n] = (float)dd->cdata[ielm];
                c->weights[n] = dd->weight;
                c->iff[n] = (short int)i;
                n++;
            }
            c->outindex = nxo*iy + ix;
            c->n = n;
            average(fileptrs,c,&outdata,&outconf,clip1,clip2,lthr,hthr,0.0,
                    sumweight);
            odata[c->outindex] = outdata;
            if (owdata != NULL) 
                owdata[c->outindex] = outconf;
        }
 
        /* If we're away from the edges, have a look and see which pixels in
           the previous row had clipping. Evaluate whether that clipping was
           really justified or not */
 
        if (iy < 2)
            continue;
        for (ix = 1; ix < nxo-1; ix++) {
            c = clipmon + clast + ix;
            if (! c->clipped)
                continue;
 
            /* If it was clipped, then evaluate the amount of 'noise' there 
               is spatially */
 
            iloc = c->outindex;
            avlev = 0.0;
            for (jy = -1; jy <= 1; jy++) {
                index_y = iloc + jy*nxo;
                for (jx = -1; jx <= 1; jx++) {
                    index = index_y + jx;
                    avlev += odata[index];
                }
            }
            avlev /= 9.0;
            avvar = 0.0;
            for (jy = -1; jy <= 1; jy++) {
                index_y = iloc + jy*nxo;
                for (jx = -1; jx <= 1; jx++) {
                    index = index_y + jx;
                    avvar += fabs(odata[index] - avlev);
                }
            }
            avvar /= 9.0;
 
            /* If the average level in this cell is below the upper clip level
               or the mean absolute deviation is smaller than the poisson 
               noise in the cell, then the clip was probably justified. */
 
            if (avlev < clip2 || avvar <= 2.0*fileptrs->noise)
                continue;
 
            /* Otherwise, create new clip levels and redo the average */
 
            average(fileptrs,c,&outdata,&outconf,clip1,clip2,lthr,hthr,
                    3.0*avvar,sumweight);
            odata[c->outindex] = outdata;
            if (owdata != NULL) 
                owdata[c->outindex] = outconf;
        }
    }
    
    /* Normalise the output confidence map */
    
    if (owdata != NULL) {
        skyest(owdata,NULL,npts,3.0,1.0,MAXDATA,&sky,&skynoise);
        renorm = 100.0/sky;
        ocdata = cpl_image_get_data_int(*outc);
        for (i = 0; i < npts; i++) 
            ocdata[i] = max(0,min(1000,(int)(owdata[i]*renorm + 0.5)));
    }
 
    /* Create the output propertylist with some provenance info */
 
    *p = cpl_propertylist_duplicate(casu_fits_get_ehu(inf[0]));    
    casu_prov(*p,inf,nimages,1);
 
    /* Add a timestamp to the propertylist */
 
    casu_timestamp(timestamp,25);
    p2 = casu_fits_get_phu(inf[0]);
    cpl_propertylist_update_string(p2,"ESO CASU_TIME",timestamp);
    cpl_propertylist_set_comment(p2,"ESO CASU_TIME",
                                 "Timestamp for matching to conf map");
 
    /* Update the WCS in the header to reflect the new offset */
 
    if (cpl_propertylist_has(*p,"CRPIX1") && 
        cpl_propertylist_has(*p,"CRPIX2")) {
        crpix1 = cpl_propertylist_get_double(*p,"CRPIX1");
        crpix2 = cpl_propertylist_get_double(*p,"CRPIX2");
        crpix1 += fileptrs->xoff;
        crpix2 += fileptrs->yoff;
        cpl_propertylist_update_double(*p,"CRPIX1",crpix1);
        cpl_propertylist_update_double(*p,"CRPIX2",crpix2);
    }
 
    /* Get out of here */
 
    tidy(nxo,fileptrs,clipmon,owdata);
    GOOD_STATUS
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
static void average(dstrct *fileptrs, keeptabs *c, float *outdata, 
                    float *outconf, float cliplow, float cliphigh, float lsig, 
                    float hsig, float extra, float sumweight) {
    int i,imin,imax;
    float valuemax,valuemin,cwmin,cwmax,sum,cnumb,numb,cw,cv,reflev,noise;
    float sky,clipval;
    
    /* If there aren't any pixels defined for this (kind of a funny state
       to be in, but never mind), give it some nominal value, which is the
       sky background of the first input image. Flag it with zero confidence */
    
    if (c->n <= 0) {
        *outdata = fileptrs->sky;
        *outconf = 0.0;
        return;
    }
    
    /* Initialise a few things (avoid boring compiler errors about 
       uninitialised variables */
    
    valuemin = 1.0e10;
    valuemax = -1.0e10;
    cwmin = 1.0e10;
    cwmax = -1.0e10;
    imin = 0;
    imax = 0;
    sum = 0.0;
    cnumb = 0.0;
    numb = 0.0;
    
    /* Now loop through all the data for this point, keeping track of the 
       min and max */
    
    for (i = 0; i < c->n; i++) {
        cw = c->weights[i]*c->confs[i];
        cv = c->values[i];
        sum += cv*cw;
        cnumb +=cw;
        numb += c->confs[i];
        if (cv < valuemin) {
            valuemin = cv;
            cwmin = cw;
            imin = i;
        } 
        if (cv > valuemax) {
            valuemax = cv;
            cwmax = cw;
            imax = i;
        }
    }
    if (cnumb > 0.0)
        *outdata = sum/cnumb;
    else
        *outdata = fileptrs->sky;
 
    /* See if we need to clip. Look at bright one first */
    
    if (valuemax > cliphigh && numb > 150.0 && cnumb > 150.0) {
        reflev = (sum - valuemax*cwmax)/(cnumb - cwmax);
        noise = (fileptrs+c->iff[imax])->noise;
        sky = (fileptrs+c->iff[imax])->sky;
        clipval = reflev + hsig*noise*max(1.0,reflev)/max(1.0,sky) + extra;
        if (valuemax > clipval) {
            sum -= valuemax*cwmax;
            cnumb -= cwmax;
            *outdata = reflev;
            c->clipped = 1;
        }
    }
    
    /* Now look at the lowest value */
    
    if (valuemin < cliplow && numb > 150.0 && cnumb > 150.0) {
        reflev = (sum - valuemin*cwmin)/(cnumb - cwmin);
        noise = (fileptrs+c->iff[imin])->noise;
        clipval = reflev - lsig*noise;
        if (valuemin < clipval) {
            cnumb -= cwmin;
            *outdata = reflev;
        }
    }
    
    /* Do the output confidence */
    
    *outconf = cnumb/sumweight;
}
        
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
    
static keeptabs *clip_open(int nimages, int nxo) {
    keeptabs *c;
    int i;
    short int *iff;
    float *dptr;
 
    c = cpl_malloc(2*nxo*sizeof(keeptabs));
    for (i = 0; i < 2*nxo; i++) {
        dptr = cpl_malloc(3*nimages*sizeof(*dptr));
        iff = cpl_malloc(nimages*sizeof(*iff));
        (c+i)->values = dptr;
        (c+i)->confs = dptr + nimages;
        (c+i)->weights = dptr + 2*nimages;
        (c+i)->iff = iff;
        (c+i)->n = 0;
        (c+i)->outindex = -1;
        (c+i)->clipped = 0;
    }
    return(c);
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
    
static void clip_close(keeptabs **c, int nxo) {
    int i;
 
    for (i = 0; i < 2*nxo; i++) {
        freespace((*c+i)->values);
        freespace((*c+i)->iff);
    }
    freespace(*c);
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
    
static void skyest(float *data, int *cdata, long npts, float thresh, 
                   float dmin, float dmax, float *skymed, float *skynoise) {
    unsigned char *bpm;
    int i;
 
    /* Get a dummy bad pixel mask */
 
    bpm = cpl_calloc(npts,sizeof(*bpm));
    if (cdata != NULL) {
        for (i = 0; i < npts; i++) 
            bpm[i] = (cdata[i] == 0);
    }
 
    /* Get the stats */
 
    casu_qmedsig(data,bpm,npts,thresh,2,dmin,dmax,skymed,skynoise);
 
    /* Clean up */
 
    freespace(bpm);
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
static void tidy(int nxo, dstrct *fileptrs, keeptabs *clipmon, float *owdata) {
 
    if(owdata != NULL) {cpl_free(owdata);}
    //freespace(owdata);
    clip_close(&clipmon,nxo);
    if(fileptrs != NULL) {cpl_free(fileptrs);}
    //freespace(fileptrs);
}    

    
/*
 
$Log: casu_imdither.c,v $
Revision 1.5  2015/11/25 10:26:31  jim
replaced some hardcoded numbers with defines
 
Revision 1.4  2015/11/18 20:06:12  jim
Fixed CASU_TIME keyword
 
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.7  2015/01/29 11:51:56  jim
modified comments
 
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/04/09 09:09:51  jim
Detabbed
 
Revision 1.3  2014/03/26 15:45:45  jim
Modified to remove globals and to allow for floating point confidence maps
 
Revision 1.2  2013/11/21 09:38:14  jim
detabbed
 
Revision 1.1.1.1  2013-08-27 12:07:48  jim
Imported
 
 
*/