vircam_defringe.c

/* $Id: vircam_defringe.c,v 1.6 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.6 $
 * $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"

#define NTERMAX 15
#define ITERMAX 7
#define OGAP 0.01
#define GAPMIN 0.25
#define CLIPLEV 4.0

static int vircam_defringe_1(vir_fits **infiles, int nimages, vir_fits *fringe,
                             vir_mask *mask, int nbsize, float *scaleout,
                             int *status);

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

extern int vircam_defringe(vir_fits **infiles, int nimages, vir_fits **fringes,
                           int nfringes, vir_mask *mask, int nbsize,
                           int *status) {
    float *sig_before,*scale,*idata,med,mad,ratio;
    unsigned char *bpm;
    int i,j,npts;
    char pname1[64],comment1[64],pname2[64],comment2[64];
    cpl_propertylist *ehu;

    /* Inherited status */

    if (*status != VIR_OK)
        return(*status);

    /* Get some workspace to use for working out the background sigmas */

    sig_before = cpl_malloc(nimages*sizeof(float));
    scale = cpl_malloc(nimages*sizeof(float));

    /* Now work out the background sigma of each image */

    bpm = vircam_mask_get_data(mask);
    npts = cpl_image_get_size_x(vircam_fits_get_image(infiles[0]))*
        cpl_image_get_size_y(vircam_fits_get_image(infiles[0]));
    for (i = 0; i < nimages; i++) {
        idata = cpl_image_get_data_float(vircam_fits_get_image(infiles[i]));
        vircam_medmad(idata,bpm,(long)npts,&med,&mad);
        sig_before[i] = 1.48*mad;
    }

    /* Now do the defringing */

    for (i = 0; i < nfringes; i++) {
        (void)vircam_defringe_1(infiles,nimages,fringes[i],mask,nbsize,scale,
                                status);

        /* Create scale factor and fringe name information for headers */

        (void)sprintf(pname1,"ESO DRS FRINGE%d",i);
        (void)sprintf(comment1,"Fringe frame # %d",i);
        (void)sprintf(pname2,"ESO DRS FRNGSC%d",i);
        (void)sprintf(comment2,"Fringe scale # %d",i);
        
        /* Now loop for each image and write these to the headers */
        
        for (j = 0; j < nimages; j++) {
            ehu = vircam_fits_get_ehu(infiles[j]);
            cpl_propertylist_update_string(ehu,pname1,
                                           vircam_fits_get_fullname(fringes[i]));
            cpl_propertylist_set_comment(ehu,pname1,comment1);
            cpl_propertylist_update_float(ehu,pname2,scale[i]);
            cpl_propertylist_set_comment(ehu,pname2,comment2);
        }
    }
        
    /* Now work out the final background sigma and add the ratio of the
       before to after to the QC header of each image */

    for (i = 0; i < nimages; i++) {
        ehu = vircam_fits_get_ehu(infiles[i]);
        idata = cpl_image_get_data_float(vircam_fits_get_image(infiles[i]));
        vircam_medmad(idata,bpm,(long)npts,&med,&mad);
        ratio = sig_before[i]/(1.48*mad);
        cpl_propertylist_update_float(ehu,"ESO QC FRINGE_RATIO",ratio);
        cpl_propertylist_set_comment(ehu,"ESO QC FRINGE_RATIO",
                                     "Ratio RMS before to after defringing");
    }
    freespace(sig_before);
    freespace(scale);
    GOOD_STATUS
}

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

static int vircam_defringe_1(vir_fits **infiles, int nimages, vir_fits *fringe,
                             vir_mask *mask, int nbsize, float *scaleout,
                             int *status) {
    cpl_image *frback,*imback,*fringefit,*imfit;
    float frmed,immed,*frdata,*frorig,*wptr,*imdata,*data,scaleth,scalemin;
    float scaleprev,spreadmin,spreadfbest,gap,offset,clipmax,clip,spreadf;
    float scalefound,scale,scalelist[3],diff,spreadlist[3],a,b,c;
    long npts,ntot;
    int i,iter,nter,k,j;
    vir_fits *im;
    unsigned char *bpm;

    /* Inherited status */

    if (*status != VIR_OK)
        return(*status);

    /* Prepare the fringe frame by creating a background mask */

    vircam_backmap(fringe,mask,nbsize,&frback,&frmed);

    /* Create a fringe frame image that is background corrected and normalised
       to zero median */

    fringefit = cpl_image_subtract_create(vircam_fits_get_image(fringe),
                                          frback);
    cpl_image_subtract_scalar(fringefit,frmed);
    frdata = cpl_image_get_data_float(fringefit);
    npts = cpl_image_get_size_x(fringefit)*cpl_image_get_size_y(fringefit);
    frorig = cpl_image_get_data_float(vircam_fits_get_image(fringe));

    /* Get some workspace */

    wptr = cpl_malloc(npts*sizeof(float));

    /* Get BPM */

    bpm = vircam_mask_get_data(mask);

    /* Now loop for each of the input images */

    for (i = 0; i < nimages; i++) {
        im = infiles[i];

        /* Create a background map and correct for it */

        vircam_backmap(im,mask,nbsize,&imback,&immed);
        imdata = cpl_image_get_data_float(vircam_fits_get_image(im));
        imfit = cpl_image_subtract_create(vircam_fits_get_image(im),imback);
        cpl_image_subtract_scalar(imfit,immed);
        data = cpl_image_get_data_float(imfit);

        /* The overall medians are used as a first guess of the scaling 
           between the two images */

        scaleth = immed/frmed;

        /* Set up some values for tracking the goodness of fit */

        scalemin = 0.0;
        scaleprev = 1.0e6;
        spreadmin = 1.0e3;
        spreadfbest = 1.0e6;
        iter = 0;
        nter = 0;
        gap = scaleth;
        offset = 0.5;
        clipmax = 1.0e3;

        /* Begin the iteration loop */

        while (nter < NTERMAX && iter < ITERMAX && (fabs(offset)*gap > OGAP ||
                                                    gap > GAPMIN)) {
            iter++;
            nter++;

            /* Clip levels */

            clip = min(clipmax,CLIPLEV*1.48*spreadmin);

            /* Speed up convergence if sitting on top of a solution. 
               Slow it down if outside range */

            if (fabs(scalemin - scaleprev) < 0.5*gap) 
                iter += 1;
            scaleprev = scalemin;
            if (fabs(offset) > 0.9)
                iter -= 1;
            gap = scaleth*pow(2.0,(double)(1-iter));

            /* Initialise a few things */

            spreadf = 1.0e6;
            scalefound = 2.0;
            
            /* Do three calculations -- just below, at and just above
               the current best guess scale factor */

            for (k = 0; k < 3; k++) {
                scale = scalemin + gap*(float)(k-1);
                scalelist[k] = scale;
                ntot = 0;
            
                /* Do a scaled subtraction of the fringe from the data */

                for (j = 0; j < npts; j++) {
                    if (bpm[j] == 0) {
                        diff = fabs(data[j] - scale*frdata[j]);
                        if (diff < clip)
                            wptr[ntot++] = diff;
                    }
                }

                /* Find the MAD */

                spreadlist[k] = vircam_med(wptr,NULL,ntot);
                if (spreadlist[k] < spreadf) {
                    spreadf = spreadlist[k]; 
                    scalefound = scale;
                }
            }

            /* Right, how have we done on this iteration? If the spread
               has started to increase then this is the last iteration */

            if (spreadf > spreadfbest) {
                nter = NTERMAX + 1;

            /* Otherwise interpolate to find the best solution */

            } else {
                a = spreadlist[1];
                b = 0.5*(spreadlist[2] - spreadlist[0]);
                c = 0.5*(spreadlist[2] + spreadlist[0] - 2.0*spreadlist[1]);
                offset = max(min((-0.5*b/c),1.0),-1.0);
                spreadmin = a + b*offset + c*offset*offset;
                scalemin = scalelist[1] + offset*gap;
                
                /* Make sure we're not going for a maximum instead of
                   a minimum */

                if (spreadmin > spreadf) {
                    spreadmin = spreadf;
                    scalemin = scalefound;
                }
            }
          
            /* Define the best spread found so far */

            spreadfbest = min(spreadfbest,spreadf);

        } /* End of iteration loop */

        /* Trap for no refinement */

        if (iter == 0)
            scalemin = scaleth;
            
        /* Subtract the fringe frame now with the defined scale factor */

        for (j = 0; j < npts; j++) 
            imdata[j] -= scalemin*(frorig[j] - frmed);
        scaleout[i] = scalemin;

        /* Tidy */

        freeimage(imfit);
        freeimage(imback);
    }

    /* Do a bit more tidying */

    freeimage(fringefit);
    freeimage(frback);
    freespace(wptr);
    GOOD_STATUS
}

/*

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

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

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

Revision 1.3  2006/12/06 13:03:25  jim
Oops, left some diagnostic prints in...

Revision 1.2  2006/12/06 12:59:14  jim
Fixed a small bug in defringe_1

Revision 1.1  2006/12/01 14:09:40  jim
Initial entry (not yet debugged)


*/

---