vircam_imdither.c

/* $Id: vircam_imdither.c,v 1.14 2007/10/25 17:34:00 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/10/25 17:34:00 $
 * $Revision: 1.14 $
 * $Name:  $
 */

/* Includes */

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

#include <math.h>
#include <cpl.h>
#include "vircam_mods.h"
#include "vircam_utils.h"
#include "vircam_fits.h"
#include "vircam_stats.h"
#include "vircam_pfits.h"

typedef struct {
    vir_fits  *fname;
    vir_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

static dstrct *fileptrs = NULL;
static float *odata = NULL;
static float *owdata = NULL;
static keeptabs *clipmon = NULL;

static float lsig;
static float hsig;
static float sumweight;
static int   nxo;
static int   nyo;

static void average(keeptabs *, float *, float *, float, float, float);
static keeptabs *clip_open( int);
static void clip_close(keeptabs **);
static void fileptrs_close(void);
static void skyest(float *, long, float, float *, float *);
static void tidy(void);

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
    
extern int vircam_imdither(vir_fits **inf, vir_fits **inconf, int nimages, 
                           int nconfs, float lthr, float hthr,
                           cpl_propertylist **p, cpl_image **out, 
                           cpl_image **outc, int *status) {

    int i,itx,iy,ccur,clast,ix,n,iline,icol,jy,jx,*ocdata;
    long npts,ielm,iloc,index_y,index;
    dstrct *dd;
    keeptabs *c;
    float minxoff,minyoff,expref,sky,skynoise,clip1,clip2,outdata;
    float outconf,avlev,avvar,renorm,exposure;
    double crpix1,crpix2;
    cpl_propertylist *ehu,*p2;
    const char *fctid = "vircam_imdither";
    char timestamp[25];

    /* Inherited status */

    *out = NULL;
    *outc = NULL;
    *p = NULL;
    if (*status != VIR_OK)
        return(*status);

    /* Is there any point in being here? */

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

    /* Initialise some global variables */

    lsig = lthr;
    hsig = hthr;

    /* Allocate file struct array and fill in some values */

    fileptrs = cpl_malloc(nimages*sizeof(dstrct));
    (void)vircam_pfits_get_exptime(vircam_fits_get_phu(inf[0]),&exposure);
    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(vircam_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(vircam_fits_get_image(inconf[0]));
        } else {
            dd->conf = inconf[i];
            dd->cdata = cpl_image_get_data_int(vircam_fits_get_image(inconf[i]));
        }
        ehu = vircam_fits_get_ehu(dd->fname);
        (void)vircam_pfits_get_jxoff(ehu,&(dd->xoff));
        (void)vircam_pfits_get_jyoff(ehu,&(dd->yoff));
        minxoff = min(dd->xoff,minxoff);
        minyoff = min(dd->yoff,minyoff);
        (void)vircam_pfits_get_exptime(ehu,&exposure);
        dd->expscale = exposure/expref;

        /* Now work out a background and background noise estimate */

        dd->nx = cpl_image_get_size_x(vircam_fits_get_image(dd->fname));
        dd->ny = cpl_image_get_size_y(vircam_fits_get_image(dd->fname));
        npts = dd->nx*dd->ny;
        skyest(dd->data,npts,3.0,&sky,&skynoise);
        dd->sky = sky;
        dd->noise = skynoise;
        
        /* Double check to make sure the confidence maps and images have the
           same dimensions */

        if (cpl_image_get_size_x(vircam_fits_get_image(dd->conf)) != dd->nx || 
            cpl_image_get_size_y(vircam_fits_get_image(dd->conf)) != dd->ny) {
            cpl_msg_error(fctid,"VIRCAM_IMDITHER: Image %s and Confidence map %s don't match",
                          vircam_fits_get_fullname(dd->fname),
                          vircam_fits_get_fullname(dd->conf));
            tidy();
            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;
    }

    /* 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(nxo,nyo,CPL_TYPE_FLOAT);

    /* If an output confidence map has been specified, then create it now. */

    if (nconfs != 0) 
        *outc = cpl_image_new(nxo,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);

    /* 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] = dd->cdata[ielm];
                c->weights[n] = dd->weight;
                c->iff[n] = (short int)i;
                n++;
            }
            c->outindex = nxo*iy + ix;
            c->n = n;
            average(c,&outdata,&outconf,clip1,clip2,0.0);
            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(c,&outdata,&outconf,clip1,clip2,3.0*avvar);
            odata[c->outindex] = outdata;
            if (owdata != NULL) 
                owdata[c->outindex] = outconf;
        }
    }
    
    /* Normalise the output confidence map */
    
    if (owdata != NULL) {
        skyest(owdata,npts,3.0,&sky,&skynoise);
        renorm = 100.0/sky;
        ocdata = cpl_image_get_data_int(*outc);
        for (i = 0; i < npts; i++) 
            ocdata[i] = max(0,min(110,((int)(owdata[i]*renorm + 0.5))));
    }

    /* Create the output propertylist with some provenance info */

    *p = cpl_propertylist_duplicate(vircam_fits_get_ehu(inf[0]));    
    vircam_prov(*p,inf,nimages);

    /* Add a timestamp to the propertylist */

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

    /* Update the WCS in the header to reflect the new offset */

    (void)vircam_pfits_get_crpix1(*p,&crpix1);
    (void)vircam_pfits_get_crpix2(*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();
    GOOD_STATUS
}

static void average(keeptabs *c, float *outdata, float *outconf, float cliplow,
                    float cliphigh, float extra) {
    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) {
    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 i;

    for (i = 0; i < 2*nxo; i++) {
        freespace((*c+i)->values);
        freespace((*c+i)->iff);
    }
    freespace(*c);
}

static void fileptrs_close(void) {


    freespace(fileptrs);
}    


static void skyest(float *data, long npts, 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,2,-1000.0,65535.0,skymed,skynoise);

    /* Clean up */

    freespace(bpm);
}

static void tidy(void) {

    freespace(owdata);
    clip_close(&clipmon);
    fileptrs_close();
}    

/*

$Log: vircam_imdither.c,v $
Revision 1.14  2007/10/25 17:34:00  jim
Modified to remove lint warnings

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

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

Revision 1.11  2006/11/27 12:05:49  jim
Changed call from cpl_propertylist_append to cpl_propertylist_update

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

Revision 1.9  2006/06/13 21:25:41  jim
Fixed bug in normalising the offsets

Revision 1.8  2006/06/09 22:25:06  jim
tidied up a few bugs

Revision 1.7  2006/06/09 11:26:26  jim
Small changes to keep lint happy

Revision 1.6  2006/06/08 14:53:18  jim
Modified a few things to keep lint happy

Revision 1.5  2006/06/06 13:05:52  jim
Adds VIR_TIME to primary header

Revision 1.4  2006/05/26 15:03:32  jim
Fixed bug where status variable address was being changed

Revision 1.3  2006/05/18 09:48:16  jim
fixed docs

Revision 1.2  2006/05/17 12:07:12  jim
Fixed error condition returns

Revision 1.1  2006/05/15 13:14:45  jim
new routine


*/

---