imcore_conf.c

/* $Id: imcore_conf.c,v 1.4 2015/08/12 11:16:55 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/08/12 11:16:55 $
 * $Revision: 1.4 $
 * $Name:  $
 */
 
 
 
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
 
#include <cpl.h>
#include "casu_fits.h"
#include "casu_wcsutils.h"
 
#include "ap.h"
#include "util.h"
#include "imcore.h"
#include "floatmath.h"
#include "imcore_version.h"
 
#define FATAL_ERR(_a) {freetable(tab); cpl_msg_error(fctid,_a); tidy(); return(CASU_FATAL);}
 
#define NW 5
 
static float *smoothed = NULL;
static float *smoothedc = NULL;
static unsigned char *mflag = NULL;
static float *indata = NULL;
static int *confdata = NULL;
static float *confsqrt = NULL;
static ap_t ap;
static int freeconf = 0;
 
static float weights[NW*NW];
static float weightc[NW*NW];
static long nx;
static long ny;
 
static void crweights(float);
static void convolve(int);
static void tidy(void);

 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
extern int imcore_conf(casu_fits *infile, casu_fits *conf, int ipix, 
                       float threshold, int icrowd, float rcore, int nbsize, 
                       int cattype, float filtfwhm, float gain,
                       casu_tfits **outcat) {
 
    int i,retval,mulpix,j,nw2,status,nobjects,imcore_xcol,imcore_ycol;
    float fconst,nullval,skymed,skysig,thresh,xintmin,offset;
    float *current,isatbc,isat,junk,*currentc;
    long npix,nxc,nyc,npts;
    cpl_image *map,*cmap;
    cpl_propertylist *plist,*extra;
    cpl_table *tab;
    char card[64];
    const char *fctid = "imcore_conf";
 
    /* Initialise output */
 
    *outcat = NULL;
 
    /* Useful constants */
 
    fconst = CPL_MATH_LOG2E;
    nullval = 0.0;
    nobjects = 0;
 
    /* Open input image */
 
    tab = NULL;
    map = casu_fits_get_image(infile);
    if ((indata = cpl_image_get_data_float(map)) == NULL) 
        FATAL_ERR("Error getting image data");
    nx = (long)cpl_image_get_size_x(map);
    ny = (long)cpl_image_get_size_y(map);
    npts = nx*ny;
 
    /* Open the associated confidence map, if it exists */
 
    if (conf != NULL) {
        cmap = casu_fits_get_image(conf);
        if ((confdata = cpl_image_get_data(cmap)) == NULL)
            FATAL_ERR("Error getting confidence map data");
        nxc = (long)cpl_image_get_size_x(cmap);
        nyc = (long)cpl_image_get_size_y(cmap);
        if ((nx != nxc) || (ny != nyc))
            FATAL_ERR("Input image and confidence dimensions don't match");
        freeconf = 0;
    } else {
        confdata = cpl_malloc(npts*sizeof(*confdata));
        for (i = 0; i < npts; i++) 
            confdata[i] = 100;
        freeconf = 1;
        cmap = NULL;
    }
    confsqrt = cpl_malloc(npts*sizeof(*confsqrt));
    for (i = 0; i < npts; i++)
        confsqrt[i] = sqrt(0.01*(float)confdata[i]);
    
    /* Get mflag array for flagging saturated pixels */
 
    npix = nx*ny;
    mflag = cpl_calloc(npix,sizeof(*mflag));
 
    /* Open the ap structure and define some stuff in it */
 
    ap.lsiz = nx;
    ap.csiz = ny;
    ap.inframe = map;
    ap.conframe = cmap;
    ap.xtnum = casu_fits_get_nexten(infile);
    imcore_apinit(&ap);
    ap.indata = indata;
    ap.confdata = confdata;
    ap.multiply = 1;
    ap.ipnop = ipix;
    ap.mflag = mflag;
    ap.rcore = rcore;
    ap.filtfwhm = filtfwhm;
    ap.icrowd = icrowd;
    ap.fconst = fconst;
 
    /* Open the output catalogue FITS table */
 
    imcore_tabinit(&ap,&imcore_xcol,&imcore_ycol,cattype,&tab);
 
    /* Set up the data flags */
 
    for (i = 0; i < npix; i++) 
        if (confdata[i] == 0)
            mflag[i] = MF_ZEROCONF;
        else if (indata[i] < STUPID_VALUE)
            mflag[i] = MF_STUPID_VALUE;
        else 
            mflag[i] = MF_CLEANPIX;
 
    /* Compute a saturation level before background correction. */
 
    retval = imcore_backstats(&ap,nullval,1,&skymed,&skysig,&isatbc);
    if (retval != CASU_OK) 
        FATAL_ERR("Error calculating saturation level");
 
    /* Flag up regions where the value is above the saturation level
       determined above. */
 
    for (i = 0; i < npix ; i++)
        if (mflag[i] == MF_CLEANPIX && indata[i] > isatbc)
            mflag[i] = MF_SATURATED;
 
    /* Compute the background variation and remove it from the data*/
 
    retval = imcore_background(&ap,nbsize,nullval);
    if (retval != CASU_OK) 
        FATAL_ERR("Error calculating background");
 
    /* Compute a saturation level after background correction. */
 
    retval = imcore_backstats(&ap,nullval,1,&skymed,&skysig,&isat);
    if (retval != CASU_OK) 
        FATAL_ERR("Error calculating saturation");
 
 
    /* Compute background statistics */
 
    retval = imcore_backstats(&ap,nullval,0,&skymed,&skysig,&junk);
    if (retval != CASU_OK) 
        FATAL_ERR("Error calculating background stats");
 
    /* Get the propertly list for the input file and add some info*/
 
    plist = casu_fits_get_ehu(infile);
    cpl_propertylist_update_float(plist,"ESO DRS SKYLEVEL",skymed);
    cpl_propertylist_set_comment(plist,"ESO DRS SKYLEVEL",
                                 "[adu] Median sky brightness");
    cpl_propertylist_update_float(plist,"ESO DRS SKYNOISE",skysig);
    cpl_propertylist_set_comment(plist,"ESO DRS SKYNOISE",
                                 "[adu] Pixel noise at sky level");
 
    /* Take mean sky level out of data. Uncomment the bit below if
       we start working out the covariance in dribbled data. See cirdr
       version. */
 
    for (i = 0; i < nx*ny; i++) {
        indata[i] -= skymed;
        /* if (indata[i]*confsqrt[i] > 3.0*skysig &&  */
        /*     mflag[i] != MF_SATURATED && mflag[i] != MF_STUPID_VALUE)  */
        /*     mflag[i] = MF_3SIG; */
    }
    
    /* Work out isophotal detection threshold levels */
 
    thresh = threshold*skysig;
    
    /* Minimum intensity for consideration */
 
    xintmin = 1.5*thresh*((float)ipix);
 
    /* Minimum size for considering multiple images */
 
    mulpix = MAX(8,2*ipix);
 
    /* Actual areal profile levels: T, 2T, 4T, 8T,...but written wrt T
       i.e. threshold as a power of 2 */
 
    offset = logf(thresh)*fconst;
 
    /* Get a bit of workspace for buffers */
 
    smoothed = cpl_malloc(nx*sizeof(*smoothed));
    smoothedc = cpl_malloc(nx*sizeof(*smoothedc));
 
    /* Define a few things more things in ap structure */
 
    ap.mulpix = mulpix;
    ap.areal_offset = offset;
    ap.thresh = thresh;
    ap.xintmin = xintmin;
    ap.sigma = skysig;
    ap.background = skymed;
    ap.saturation = (float)isat;
 
    /* Set the weights */
 
    crweights(filtfwhm);
    nw2 = NW/2;
 
    /* Right, now for the extraction loop.  Begin by defining a group of
       three rows of data and confidence */
 
    for (j = nw2; j < ny-nw2; j++) {
        current = indata + j*nx;
        currentc = confsqrt + j*nx;
        convolve(j);
   
        /* Do the detection now */
 
        imcore_apline(&ap,current,currentc,smoothed,smoothedc,j,NULL);
 
        /* Make sure we're not overruning the stacks */
 
        if (ap.ibstack > (ap.maxbl - ap.lsiz))
            imcore_apfu(&ap);
        if (ap.ipstack > (ap.maxpa*3/4))
            for (i = 0; i < ap.maxpa*3/8; i++)
                imcore_apfu(&ap);
 
        /* See if there are any images to terminate */
 
        if (ap.ipstack > 1)
            imcore_terminate(&ap,cattype,gain,&nobjects,tab);
    }
 
    /* Post process. First truncate the cpl_table to the correct size and then
       work out an estimate of the seeing */
 
    cpl_table_set_size(tab,nobjects);
    retval = imcore_do_seeing(&ap,cattype,nobjects,tab);
    if (retval != CASU_OK)
        FATAL_ERR("Error working out seeing");
    imcore_tabclose(&ap,cattype);
 
    /* Create a property list with extra parameters. First QC parameters */
 
    extra = cpl_propertylist_duplicate(casu_fits_get_ehu(infile));
    cpl_propertylist_update_float(extra,"ESO QC SATURATION",ap.saturation);
    cpl_propertylist_set_comment(extra,"ESO QC SATURATION",
                                 "[adu] Saturation level");
    cpl_propertylist_update_float(extra,"ESO QC MEAN_SKY",ap.background);
    cpl_propertylist_set_comment(extra,"ESO QC MEAN_SKY",
                                 "[adu] Median sky brightness");
    cpl_propertylist_update_float(extra,"ESO QC SKY_NOISE",ap.sigma);
    cpl_propertylist_set_comment(extra,"ESO QC SKY_NOISE",
                                 "[adu] Pixel noise at sky level");
 
    /* Now DRS parameters */
 
    cpl_propertylist_update_float(extra,"ESO DRS THRESHOL",ap.thresh);
    cpl_propertylist_set_comment(extra,"ESO DRS THRESHOL",
                                 "[adu] Isophotal analysis threshold");
    cpl_propertylist_update_int(extra,"ESO DRS MINPIX",ap.ipnop);
    cpl_propertylist_set_comment(extra,"ESO DRS MINPIX",
                                 "[pixels] Minimum size for images");
    cpl_propertylist_update_int(extra,"ESO DRS CROWDED",ap.icrowd);
    cpl_propertylist_set_comment(extra,"ESO DRS CROWDED",
                                 "Crowded field analysis flag");
    cpl_propertylist_update_float(extra,"ESO DRS RCORE",ap.rcore);
    cpl_propertylist_set_comment(extra,"ESO DRS RCORE",
                                 "[pixels] Core radius for default profile fit");
    cpl_propertylist_update_float(extra,"ESO DRS SEEING",ap.fwhm);
    cpl_propertylist_set_comment(extra,"ESO DRS SEEING",
                                 "[pixels] Average FWHM");
    cpl_propertylist_update_float(extra,"ESO DRS FILTFWHM",ap.filtfwhm);
    cpl_propertylist_set_comment(extra,"ESO DRS FILTFWHM",
                                 "[pixels] FWHM of smoothing kernel");
    cpl_propertylist_update_int(extra,"ESO DRS XCOL",imcore_xcol);
    cpl_propertylist_set_comment(extra,"ESO DRS XCOL","Column for X position");
    cpl_propertylist_update_int(extra,"ESO DRS YCOL",imcore_ycol);
    cpl_propertylist_set_comment(extra,"ESO DRS YCOL","Column for Y position");
    cpl_propertylist_update_int(extra,"ESO DRS NXOUT",nx);
    cpl_propertylist_set_comment(extra,"ESO DRS NXOUT",
                                 "X Dimension of input image");
    cpl_propertylist_update_int(extra,"ESO DRS NYOUT",ny);
    cpl_propertylist_set_comment(extra,"ESO DRS NYOUT",
                                 "Y Dimension of input image");
    snprintf(card,64,"IMCORE version: %s",imcore_version);
    cpl_propertylist_append_string(extra,"HISTORY",card);
    
    /* Now wrap all this stuff up and send it back */
 
    plist = cpl_propertylist_duplicate(casu_fits_get_phu(infile));
    status = CASU_OK;
    (void)casu_tabwcs(extra,imcore_xcol,imcore_ycol,&status);
    *outcat = casu_tfits_wrap(tab,NULL,plist,extra);
 
    /* Tidy and exit */  
 
    tidy();
    return(CASU_OK);
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
static void crweights(float filtfwhm) {
    int i,j,nw2,n;
    double gsigsq,di,dj;
    float renorm;
 
    /* Get the kernel size */
 
    nw2 = NW/2;
    
    /* Set the normalisation constants */
 
    gsigsq = 1.0/(2.0*pow(MAX(1.0,(double)filtfwhm)/2.35,2.0));
    renorm = 0.0;
 
    /* Now work out the weights */
 
    n = -1;
    for (i = -nw2; i <= nw2; i++) {
        di = (double)i;
        di *= gsigsq*di;
        for (j = -nw2; j <= nw2; j++) {
            dj = (double)j;
            dj *= gsigsq*dj;
            n++;
            weights[n] = (float)exp(-(di + dj));
            renorm += weights[n];
        }
    }
 
    /* Now normalise the weights */
 
    n = -1;
    for (i = -nw2; i <= nw2; i++) {
        for (j = -nw2; j <= nw2; j++) {
            n++;
            weights[n] /= renorm;
            /* weightc[n] = 0.01*weights[n]; */
            weightc[n] = weights[n];
        }
    }
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
static void convolve(int ir) {
    int i,nw2,ix,jx,jy,n;
    float *idata,*cdata;
 
    /* Zero the summations */
 
    for (i = 0; i < nx; i++) {
        smoothed[i] = 0.0;
        smoothedc[i] = 0.0;
    }
 
    /* Now big is the smoothing kernel? */
 
    nw2 = NW/2;
 
    /* Now loop for each column */
 
    for (ix = nw2; ix < nx-nw2; ix++) {
        n = -1;
        for (jy = ir-nw2; jy <= ir+nw2; jy++) {
            idata = indata + jy*nx;
            cdata = confsqrt + jy*nx;
            for (jx = ix-nw2; jx <= ix+nw2; jx++) {
                n++;
                smoothed[ix] += weights[n]*idata[jx];
                smoothedc[ix] += weightc[n]*idata[jx]*cdata[jx];
            }
        }
    }
}
 
static void tidy(void) {
 
    if (freeconf) 
        freespace(confdata);
    freespace(confsqrt);
    freespace(smoothed);
    freespace(smoothedc);
    freespace(mflag);
    imcore_apclose(&ap);
}

 
/* 
 
$Log: imcore_conf.c,v $
Revision 1.4  2015/08/12 11:16:55  jim
Modified procedure names to protect namespace
 
Revision 1.3  2015/08/07 13:06:54  jim
Fixed copyright to ESO
 
Revision 1.2  2015/08/06 05:34:02  jim
Fixes to get rid of compiler moans
 
Revision 1.1.1.1  2015/06/12 10:44:32  jim
Initial import
 
Revision 1.5  2015/01/09 11:42:36  jim
Fixed routines to remove globals
 
Revision 1.4  2014/12/11 12:23:34  jim
new version
 
Revision 1.3  2014/04/09 09:09:51  jim
Detabbed
 
Revision 1.2  2014/03/26 15:25:19  jim
Modified for floating point confidence maps
 
Revision 1.1.1.1  2013/08/27 12:07:48  jim
Imported
 
 
*/