imcore_opm.c

/* $Id: imcore_opm.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 <string.h>
 
#include "ap.h"
#include "util.h"
#include "imcore.h"
#include "floatmath.h"
 
#include "casu_fits.h"
 
#define FATAL_ERR(_a) {freetable(tab); cpl_msg_error(fctid,_a); tidy(); return(CASU_FATAL);}
 
#define NW 5
 
static float *indata = NULL;
static int *confdata = NULL;
static float *confsqrt = NULL;
static float *smoothed = NULL;
static float *smoothedc = NULL;
static unsigned char *mflag = NULL;
static ap_t ap;
static int freeconf = 0;
static float *incopy = NULL;
static int *ccopy = NULL;
 
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_opm(casu_fits *infile, casu_fits *conf, int ipix, 
                      float threshold, int nbsize, float filtfwhm, 
                      int niter) {
 
    int i,retval,j,nw2,iter,nclip,nxc,nyc,cattype,nobj,xcol,ycol;
    float fconst,nullval,skymed,skysig,thresh,xintmin,offset;
    float isat,isatbc,*current,junk,*currentc;
    long npix;
    cpl_image *map,*cmap;
    cpl_propertylist *plist;
    cpl_binary *opm;
    cpl_table *tab;
    const char *fctid = "imcore_opm";
 
    /* Useful constants */
 
    fconst = CPL_MATH_LOG2E;
    nullval = 0.0;
    cattype = 4;
 
    /* 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);
    npix = 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(npix*sizeof(*confdata));
        for (i = 0; i < npix; i++) 
            confdata[i] = 100;
        freeconf = 1;
        cmap = NULL;
    }
    confsqrt = cpl_malloc(npix*sizeof(*confsqrt));
    for (i = 0; i < npix; i++) 
        confsqrt[i] = sqrt(0.01*(float)confdata[i]);
 
    /* Make a copy of each */
 
    incopy = cpl_malloc(npix*sizeof(*incopy));
    ccopy = cpl_malloc(npix*sizeof(*ccopy));
    memmove(incopy,indata,npix*sizeof(*incopy));
    memmove(ccopy,confdata,npix*sizeof(*ccopy));
 
    /* Get mflag array for flagging saturated pixels */
 
    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.fconst = fconst;
    ap.filtfwhm = filtfwhm;
 
    /* Open the output catalogue FITS table */
 
    imcore_tabinit(&ap,&xcol,&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 saturated pixels */
 
    for (i = 0; i < npix ; i++) 
        if (mflag[i] == MF_CLEANPIX && indata[i] > isatbc)
            mflag[i] = MF_SATURATED;
 
    /* Get a bit of workspace for buffers */
 
    smoothed = cpl_malloc(nx*sizeof(*smoothed));
    smoothedc = cpl_malloc(nx*sizeof(*smoothedc));
 
    /* Set the weights */
 
    crweights(filtfwhm);
    nw2 = NW/2;
 
    /* Iteration loop */
 
    for (iter = 0; iter < niter; iter++) {
 
        /* 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 */
 
        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 */
 
        for (i = 0; i < nx*ny; i++)
            indata[i] -= skymed;
 
        /* Work out isophotal detection threshold levels */
 
        thresh = threshold*skysig;
 
        /* Minimum intensity for consideration */
 
        xintmin = 1.5*thresh*((float)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;
 
        /* Define a few things more things in ap structure */
 
        ap.areal_offset = offset;
        ap.thresh = thresh;
        ap.xintmin = xintmin;
        ap.sigma = skysig;
        ap.background = skymed;
        ap.saturation = isat;
 
        /* 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,1.0,&nobj,tab);
        }
    
        /* Restore input data to its former glory. Update confidence map */
    
        memmove(indata,incopy,npix*sizeof(*indata));
        nclip = 0;
        opm = cpl_mask_get_data(ap.opmask);
        for (i = 0; i < npix; i++) {
            if (opm[i]) {
                confdata[i] = 0;
                opm[i] = 0;
                nclip++;
            }
        }
        if (ap.backmap.bvals != NULL) {
            for (i = 0; i < ap.backmap.nby; i++) 
                freespace(ap.backmap.bvals[i]);
            freespace(ap.backmap.bvals);
        }
        if (nclip == 0)
            break;
    }
    opm = cpl_mask_get_data(ap.opmask);
    for (i = 0; i < npix; i++)
        opm[i] = (confdata[i] == 0);
    memmove(confdata,ccopy,npix*sizeof(*confdata));
    retval = imcore_tabclose(&ap,cattype);
 
    /* 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);
    freespace(incopy);
    freespace(ccopy);
    imcore_apclose(&ap);
}
 
/* 
 
$Log: imcore_opm.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
 
 
*/