vircam_sky.c

/* $Id: vircam_sky.c,v 1.10 2013-10-15 16:51:47 jim Exp $
 *
 * This file is part of the VIRCAM Pipeline
 * 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: 2013-10-15 16:51:47 $
 * $Revision: 1.10 $
 * $Name: not supported by cvs2svn $
 */
 
/* Includes */
 
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
 
#include <cpl.h>
#include <math.h>
#include <string.h>
 
#include <casu_utils.h>
#include <casu_mods.h>
#include <casu_mask.h>
#include <casu_stats.h>
#include <casu_wcsutils.h>
 
#include "vircam_sky.h"
#include "vircam_mods.h"
#include "vircam_pfits.h"
 
typedef struct {
    casu_fits **in;
    int      n;
    int      nalloc;
    int      status;
    casu_fits *xsky;
} xsky_struct;
 

 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
extern int vircam_tilesky(casu_fits **inlist, int nfiles, casu_mask *mask, 
                          casu_fits **skyout, int *status) { 
    int i,njit,nx,ny,ji,ns;
    unsigned char *inbpm,*rejmask,*rejplus;
    cpl_mask *cplmask;
    xsky_struct *xskys;
    cpl_propertylist *drs;
    cpl_image *outim;
    casu_fits **xsky_fits;
    const char *fctid = "vircam_tilesky";
 
    /* Inherited status */
 
    *skyout = NULL;
    if (*status != CASU_OK) 
        return(*status);
 
    /* If there aren't any images, then get out of here */
 
    if (nfiles == 0) {
        cpl_msg_error(fctid,"Sky correction impossible. No science frames");
        return(CASU_FATAL);
    }
 
    /* Wrap the input mask into a cpl_mask structure and use it to set the
       internal mask for each of the input images */
 
    inbpm = casu_mask_get_data(mask);
    nx = casu_mask_get_size_x(mask);
    ny = casu_mask_get_size_y(mask);
    cplmask = cpl_mask_wrap((cpl_size)nx,(cpl_size)ny,(cpl_binary *)inbpm);
    for (i = 0; i < nfiles; i++)
        cpl_image_reject_from_mask(casu_fits_get_image(inlist[i]),cplmask);
    cpl_mask_unwrap(cplmask);
 
    /* How many jitter steps are there? */
 
    (void)vircam_pfits_get_njsteps(casu_fits_get_phu(inlist[0]),&njit);
 
    /* Allocate some workspce */
 
    xskys = cpl_malloc(njit*sizeof(xsky_struct));
    for (i = 0; i < njit; i++) {
        xskys[i].n = 0;
        xskys[i].in = cpl_malloc(8*sizeof(casu_fits *));
        xskys[i].nalloc = 8;
        xskys[i].xsky = NULL;
        xskys[i].status = CASU_OK;
    }
    xsky_fits = cpl_malloc(njit*sizeof(casu_fits *));
 
    /* Now loop through all the files and put them into the relevant 
       structure depending upon where they are in the jitter pattern.
       NB: We don't have to check on the status of the individual
       frames because only good frames should have been passed in by
       the calling routine. */
 
    for (i = 0; i < nfiles; i++) {
        (void)vircam_pfits_get_jitteri(casu_fits_get_phu(inlist[i]),&ji);
        ji--;
        if (xskys[ji].n == xskys[ji].nalloc) {
            xskys[ji].in = cpl_realloc(xskys[ji].in,
                                       (xskys[ji].nalloc+8)*sizeof(casu_fits *));
            xskys[ji].nalloc += 8;
        }
        xskys[ji].in[xskys[ji].n] = inlist[i];
        xskys[ji].n += 1;
    }
    
    /* Now loop through each of these groups and form a sky file with
       rejection */
 
    ns = 0;
    for (i = 0; i < njit; i++) {
        *status = CASU_OK;
        if (xskys[i].n == 0) {
            xskys[i].xsky = NULL;
            continue;
        } else if (xskys[i].n == 1) {
            xskys[i].xsky = NULL;
            xsky_fits[ns++] = xskys[i].in[0];
            continue;
        }
        outim = NULL;
        (void)casu_imcombine(xskys[i].in,NULL,xskys[i].n,1,1,0,1.0,
                             "EXPTIME",&outim,NULL,&rejmask,
                             &rejplus,&drs,status);
        freespace(rejmask);
        freespace(rejplus);
        freepropertylist(drs);
        if (*status == CASU_OK) {
            xskys[i].xsky = casu_fits_wrap(outim,xskys[i].in[0],NULL,NULL);
            xsky_fits[ns++] = xskys[i].xsky;
        } else {
            freeimage(outim);
        }
    }
 
    /* Now combine all the intermediate skies to form the final sky */
 
    if (ns != 0) {
        outim = NULL;
        (void)casu_imcombine(xsky_fits,NULL,ns,1,1,0,1.0,"EXPTIME",&outim,
                             NULL,&rejmask,&rejplus,&drs,status);
        freespace(rejmask);
        freespace(rejplus);
        freepropertylist(drs);
        if (*status == CASU_OK) {
            *skyout = casu_fits_wrap(outim,inlist[0],NULL,NULL);
            (void)casu_inpaint(*skyout,64,status);
        } else {
            freeimage(outim);
            *skyout = NULL;
        }
    } else {
        *skyout = NULL;
    }
 
    /* Tidy and exit */
 
    freespace(xsky_fits);
    for (i = 0; i < njit; i++) {
        freespace(xskys[i].in);
        freefits(xskys[i].xsky);
    }
    freespace(xskys);
    return(*status);
}
 
    
 
/*
 
$Log: not supported by cvs2svn $
Revision 1.9  2012/01/27 12:25:10  jim
Fixed some casts
 
Revision 1.8  2012/01/15 17:40:09  jim
Minor modifications to take into accout the changes in cpl API for v6
 
Revision 1.7  2010/06/30 12:42:00  jim
A few fixes to stop compiler compaints
 
Revision 1.6  2009/02/23 10:46:26  jim
Modified tilesky to try and get rid of some of the low level glow around
bright and extended objects
 
Revision 1.5  2009/02/20 11:01:13  jim
Added vircam_tilesky. Commented out vircam_skycombine as it's not being
used at present and it's generating a compiler error. May delete it
completely later
 
Revision 1.4  2008/11/25 06:23:09  jim
Added some routine prologues
 
Revision 1.3  2008/11/21 10:12:36  jim
Fixed typo
 
Revision 1.2  2008/10/24 10:57:19  jim
Fixed bug in pawsky_mask so that the bad pixels in the confidence map are
folded into the object mask
 
Revision 1.1  2008/10/13 08:13:21  jim
New entry
 
*/