vircam_genbpm.c

/* $Id: vircam_genbpm.c,v 1.1 2007/11/14 10:43:35 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/11/14 10:43:35 $
 * $Revision: 1.1 $
 * $Name:  $
 */

/* Includes */

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

#include <cpl.h>

#include <math.h>

#include "vircam_mods.h"
#include "vircam_utils.h"
#include "vircam_stats.h"

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

extern int vircam_genbpm(vir_fits **flatlist, int nflatlist, float lthr,
                         float hthr, cpl_array **bpm_array, int *nbad,
                         float *badfrac, int *status) {
    cpl_image *master_img,*im;
    unsigned char *rejmask,*rejplus;
    cpl_propertylist *drs;
    int npi,i,status2,*bpm,k,nbmax;
    float *idata,med,sig,low,high;
    char *fctid = "vircam_genbpm";

    /* Inherited status */

    *nbad = 0;
    *badfrac = 0.0;
    *bpm_array = NULL;
    if (*status != VIR_OK)
        return(*status);

    /* Combine all of the dome flats into a master */

    status2 = VIR_OK;
    (void)vircam_imcombine(flatlist,nflatlist,1,3,1,5.0,&master_img,&rejmask,
                           &rejplus,&drs,&status2);
    freespace(rejmask);
    freespace(rejplus);
    freepropertylist(drs);
    if (status2 != VIR_OK) {
        cpl_msg_error(fctid,"Flat combination failed");
        *status = VIR_FATAL;
        return(*status);
    }

    /* Divide the resulting image by its median */

    idata = cpl_image_get_data_float(master_img);
    npi = vircam_getnpts(master_img);
    vircam_medsig(idata,NULL,npi,&med,&sig);
    cpl_image_divide_scalar(master_img,(double)med);
    for (i = 0; i < npi; i++)
        if (idata[i] == 0.0)
            idata[i] = 1.0;

    /* Create and zero the rejection mask */

    *bpm_array = cpl_array_new(npi,CPL_TYPE_INT);
    bpm = cpl_array_get_data_int(*bpm_array);
    for (i = 0; i < npi; i++)
        bpm[i] = 0;

    /* Now loop for all input images and divide each by the master */

    for (i = 0; i < nflatlist; i++) {
        im = cpl_image_duplicate(vircam_fits_get_image(flatlist[i]));
        cpl_image_divide(im,master_img);
        idata = cpl_image_get_data_float(im);
        vircam_medmad(idata,NULL,npi,&med,&sig);
        sig *= 1.48;

        /* Divide the resulting image by its median */

        cpl_image_divide_scalar(im,med);

        /* Get the standard deviation of the image */

        low = 1.0 - lthr*sig/med;
        high = 1.0 + hthr*sig/med;

        /* Now define the bad pixels */

        for (k = 0; k < npi; k++)
            if (idata[k] < low || idata[k] > high)
                bpm[k] += 1;
        cpl_image_delete(im);
    }
    cpl_image_delete(master_img);

    /* Go through the rejection mask and if a pixel has been marked bad 
       more than a set number of times, then it is defined as bad */

    nbmax = max(2,nflatlist/4);
    for (i = 0; i < npi; i++) {
        if (bpm[i] >= nbmax) {
            bpm[i] = 1;
            (*nbad)++;
        } else
            bpm[i] = 0;
    }
     *badfrac = (float)(*nbad)/(float)npi;

    /* Get out of here */

    return(VIR_OK);
}

/*

$Log: vircam_genbpm.c,v $
Revision 1.1  2007/11/14 10:43:35  jim
Initial version


*/

---