sofi_utils.c

/* $Id: sofi_utils.c,v 1.17 2013-03-12 08:04:54 llundin Exp $
 *
 * This file is part of the SOFI Pipeline
 * Copyright (C) 2002,2003 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., 51 Franklin St, Fifth Floor, Boston, MA  02111-1307  USA
 */
 
/*
 * $Author: llundin $
 * $Date: 2013-03-12 08:04:54 $
 * $Revision: 1.17 $
 * $Name: not supported by cvs2svn $
 */
 
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
 
/*-----------------------------------------------------------------------------
                                   Includes
 -----------------------------------------------------------------------------*/
 
#include <math.h>
#include <string.h>
#include <cpl.h>
 
#include "sofi_utils.h"
#include "sofi_pfits.h"
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/

 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
const char * sofi_get_license(void)
{
    const char  *   sofi_license = 
        "This file is part of the SOFI Instrument Pipeline\n"
        "Copyright (C) 2002,2003 European Southern Observatory\n"
        "\n"
        "This program is free software; you can redistribute it and/or modify\n"
        "it under the terms of the GNU General Public License as published by\n"
        "the Free Software Foundation; either version 2 of the License, or\n"
        "(at your option) any later version.\n"
        "\n"
        "This program is distributed in the hope that it will be useful,\n"
        "but WITHOUT ANY WARRANTY; without even the implied warranty of\n"
        "MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the\n"
        "GNU General Public License for more details.\n"
        "\n"
        "You should have received a copy of the GNU General Public License\n"
        "along with this program; if not, write to the Free Software\n"
        "Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, \n"
        "MA  02111-1307  USA" ;
    return sofi_license ;
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
int sofi_correct_crosstalk_list(cpl_imagelist * ilist)
{
    cpl_image           *   ima ;
    int                     nima ;
    int                     i ;
 
    /* Test entries */
    if (ilist == NULL) return -1 ;
    ima  = cpl_imagelist_get(ilist, 0) ;
    if (cpl_image_get_type(ima) != CPL_TYPE_FLOAT) return -1 ;
 
    /* Initialise */
    nima = cpl_imagelist_get_size(ilist) ;
 
    /* Loop on the images */
    for (i=0 ; i<nima ; i++) {
        ima = cpl_imagelist_get(ilist, i) ;
        if (sofi_correct_crosstalk(ima) == -1) {
            cpl_msg_error(cpl_func, "Cannot correct crosstalk in image %d", i+1) ;
            return -1 ;
        }
    }
    return 0 ;
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
int sofi_correct_crosstalk(cpl_image * ima)
{
    int                     nx, ny ;
    cpl_image           *   collapsed ;
    cpl_image           *   corr ;
    float               *   pcorr ;
    float               *   pcollapsed ;
    float               *   pima ;
    float                   val ;
    int                     i, j, k ;
 
    /* Test entries */
    if (ima == NULL) return -1 ;
    if (cpl_image_get_type(ima) != CPL_TYPE_FLOAT) return -1 ;
 
    /* Initialise */
    nx = cpl_image_get_size_x(ima) ;
    ny = cpl_image_get_size_y(ima) ;
 
    /* Collapse the input image */
    collapsed = cpl_image_new(1, ny, CPL_TYPE_FLOAT) ;
    pcollapsed = cpl_image_get_data_float(collapsed) ;
    pima = cpl_image_get_data_float(ima) ;
    for (j=0 ; j<ny ; j++) {
        val = 0.0 ;
        for (i=0 ; i<nx ; i++) val += pima[i+j*nx] ;
        pcollapsed[j] = val * 1.4e-5 ;
    }
 
    /* Upper part plus lower part */
    corr = cpl_image_new(1, ny, CPL_TYPE_FLOAT) ;
    pcorr = cpl_image_get_data_float(corr) ;
    for (j=0 ; j<ny ; j++) {
        if (j < ny/2)   k = j + ny/2 ;
        else            k = j - ny/2 ;
        pcorr[j] = pcollapsed[j] + pcollapsed[k] ;
    }
    cpl_image_delete(collapsed) ;
    
    /* Subtract to the input image */
    pima = cpl_image_get_data_float(ima) ;
    for (j=0 ; j<ny ; j++) {
        for (i=0 ; i<nx ; i++) {
            pima[i+j*nx] -= pcorr[j] ;
        }
    }
    cpl_image_delete(corr) ;
    return 0 ;
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
cpl_bivector * sofi_get_offsets(cpl_frameset * fset)
{
    cpl_bivector        *   offsets ;
    double              *   offsets_x ;
    double              *   offsets_y ;
    cpl_frame           *   frame ;
    cpl_propertylist    *   plist ;
    int                     nfiles ;
    int                     i ;
 
    /* Test entries */
    if (fset == NULL) return NULL ;
 
    /* Create the offsets bi vector */
    nfiles = cpl_frameset_get_size(fset) ;
    offsets = cpl_bivector_new(nfiles) ;
    offsets_x = cpl_bivector_get_x_data(offsets) ;
    offsets_y = cpl_bivector_get_y_data(offsets) ; 
    for (i=0 ; i<nfiles ; i++) {
        if (cpl_error_get_code()) {
            cpl_bivector_delete(offsets) ;
            return NULL ;
        }
        /* X and Y offsets */
        frame = cpl_frameset_get_position(fset, i) ;
        plist=cpl_propertylist_load(cpl_frame_get_filename(frame),0);
        offsets_x[i] = sofi_pfits_get_cumoffsetx(plist) ;
        offsets_y[i] = sofi_pfits_get_cumoffsety(plist) ;
        cpl_propertylist_delete(plist) ;
        if (cpl_error_get_code()) {
            cpl_msg_error(cpl_func, "Cannot get offsets from header") ;
            cpl_bivector_delete(offsets) ;
            return NULL ;
        }
    }
    /* Subtract the first offset to all offsets */
    for (i=1 ; i<nfiles ; i++) {
        offsets_x[i] -= offsets_x[0] ;
        offsets_y[i] -= offsets_y[0] ;
    }
    offsets_x[0] = offsets_y[0] = 0.00 ;
 
    return offsets ;
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
sofi_band sofi_get_bbfilter(char * f)
{
    if (!strcmp(f, "NB_2.195"))         return SOFI_BAND_KS ;
    if (!strcmp(f, "NB_2.09"))          return SOFI_BAND_KS ;
    if (!strcmp(f, "NB_1.187"))         return SOFI_BAND_J ;
    if (!strcmp(f, "NB_1.061"))         return SOFI_BAND_J ;
    if (!strcmp(f, "Z"))                return SOFI_BAND_J ;
    if (!strcmp(f, "Js"))               return SOFI_BAND_J ;
    if (!strcmp(f, "J"))                return SOFI_BAND_J ;
    if (!strcmp(f, "H"))                return SOFI_BAND_H ;
    if (!strcmp(f, "Ks"))               return SOFI_BAND_KS ;
    if (!strcmp(f, "NB_2.28"))          return SOFI_BAND_KS ;
    if (!strcmp(f, "NB_2.248"))         return SOFI_BAND_KS ;
    if (!strcmp(f, "NB_HeI_J"))         return SOFI_BAND_J ;
    if (!strcmp(f, "NB_1.215"))         return SOFI_BAND_J ;
    if (!strcmp(f, "NB_FeII_J"))        return SOFI_BAND_J ;
    if (!strcmp(f, "NB_Pbeta"))         return SOFI_BAND_J ;
    if (!strcmp(f, "NB_FeII_H"))        return SOFI_BAND_H ;
    if (!strcmp(f, "NB_1.71"))          return SOFI_BAND_H ;
    if (!strcmp(f, "NB_HeI_K"))         return SOFI_BAND_KS ;
    if (!strcmp(f, "NB_H2_S1"))         return SOFI_BAND_KS ;
    if (!strcmp(f, "NB_Brgamma"))       return SOFI_BAND_KS ;
    if (!strcmp(f, "NB_CO"))            return SOFI_BAND_KS ;
    return SOFI_BAND_UNKNOWN ;
}
 
/*-------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/
const char * sofi_std_band_name(sofi_band band)
{
    switch (band) {
        case SOFI_BAND_J:        return "J" ;
        case SOFI_BAND_JS:       return "Js" ;
        case SOFI_BAND_JBLOCK:   return "J+Block" ;
        case SOFI_BAND_H:        return "H" ;
        case SOFI_BAND_K:        return "K" ;
        case SOFI_BAND_KS:       return "Ks" ;
        case SOFI_BAND_L:        return "L" ;
        case SOFI_BAND_M:        return "M" ;
        case SOFI_BAND_LP:       return "Lp" ;
        case SOFI_BAND_MP:       return "Mp" ;
        case SOFI_BAND_Z:        return "Z" ;
        case SOFI_BAND_SZ:       return "SZ" ;
        case SOFI_BAND_SH:       return "SH" ;
        case SOFI_BAND_SK:       return "SK" ;
        case SOFI_BAND_SL:       return "SL" ;
        default:            return "Unknown" ;
    }
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
cpl_frameset * sofi_extract_frameset(
        const cpl_frameset  *   in,
        const char          *   tag)
{
    cpl_frameset        *   out ;
    const cpl_frame     *   cur_frame ;
    cpl_frame           *   loc_frame ;
    int                     nbframes, nbext ;
    int                     i ;
 
    /* Test entries */
    if (in == NULL) return NULL ;
    if (tag == NULL) return NULL ;
 
    /* Initialise */
    nbframes = cpl_frameset_get_size(in) ;
 
    /* Count the frames with the tag */
    if ((nbext = cpl_frameset_count_tags(in, tag)) == 0) return NULL ;
 
    /* Create the output frameset */
    out = cpl_frameset_new() ;
 
    /* Loop on the requested frames and store them in out */
    nbext = 0 ;
    for (i=0 ; i<nbframes ; i++) {
        cur_frame = cpl_frameset_get_position_const(in, i) ;
        if (!strcmp(cpl_frame_get_tag(cur_frame), tag)) {
            loc_frame = cpl_frame_duplicate(cur_frame) ;
            cpl_frameset_insert(out, loc_frame) ;
            nbext ++ ;
        }
    }
    return out ;
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
const char * sofi_extract_filename(
        const cpl_frameset  *   in,
        const char          *   tag)
{
    const cpl_frame *   cur_frame ;
 
    /* Get the frame  */
    if ((cur_frame = cpl_frameset_find_const(in, tag)) == NULL) return NULL ;
    return cpl_frame_get_filename(cur_frame) ;
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
int sofi_detlin_correct(
        cpl_imagelist       *   ilist,
        const char          *   detlin_a,
        const char          *   detlin_b,
        const char          *   detlin_c)
{
    cpl_image       *   ima ;
    cpl_image       *   imb ;
    cpl_image       *   imc ;
    float           *   pima ;
    float           *   pimb ;
    float           *   pimc ;
    float           *   pdata ;
    int                 nx, ny, ni ;
    float               val, val2, val3 ;
    int                 i, j ;
 
    /* Test entries */
    if (!ilist || !detlin_a || !detlin_b || !detlin_c) return -1 ;
 
    /* Load the 3 coeffs images */
    ima = cpl_image_load(detlin_a, CPL_TYPE_FLOAT, 0, 0) ;
    imb = cpl_image_load(detlin_b, CPL_TYPE_FLOAT, 0, 0) ;
    imc = cpl_image_load(detlin_c, CPL_TYPE_FLOAT, 0, 0) ;
    if (!ima || !imb || !imc) {
        cpl_msg_error(cpl_func, "Cannot load the detlin images") ;
        if (ima) cpl_image_delete(ima) ;
        if (imb) cpl_image_delete(imb) ;
        if (imc) cpl_image_delete(imc) ;
        return -1 ;
    }
    pima = cpl_image_get_data_float(ima) ;
    pimb = cpl_image_get_data_float(imb) ;
    pimc = cpl_image_get_data_float(imc) ;
 
    /* Test sizes */
    nx = cpl_image_get_size_x(cpl_imagelist_get(ilist, 0)) ;
    ny = cpl_image_get_size_y(cpl_imagelist_get(ilist, 0)) ;
    ni = cpl_imagelist_get_size(ilist) ;
    if ((cpl_image_get_size_x(ima) != nx) ||
            (cpl_image_get_size_x(imb) != nx) ||
            (cpl_image_get_size_x(imc) != nx) ||
            (cpl_image_get_size_y(ima) != ny) ||
            (cpl_image_get_size_y(imb) != ny) ||
            (cpl_image_get_size_y(imc) != ny)) {
        cpl_msg_error(cpl_func, "Incompatible sizes") ;
        cpl_image_delete(ima) ;
        cpl_image_delete(imb) ;
        cpl_image_delete(imc) ;
        return -1 ;
    }
 
    /* Loop on pixels */
    for (i=0 ; i<nx*ny ; i++) {
        if (fabs(pimc[i]) < 1e-5) {
            /* Correct this pixel in each plane */
            for (j=0 ; j<ni ; j++) {
                pdata = cpl_image_get_data_float(cpl_imagelist_get(ilist, j)) ;
                val = pdata[i] ;
                pdata[i] = val-pima[i] ;
            }
        } else if (fabs(pimb[i]) < 1e-3) {
            for (j=0 ; j<ni ; j++) {
                pdata = cpl_image_get_data_float(cpl_imagelist_get(ilist, j)) ;
                pdata[i] = 0.0 ;
            }
        } else {
            val2 = 2 * pimc[i] / (pimb[i] * pimb[i]) ;
            /* Correct this pixel in each plane */
            for (j=0 ; j<ni ; j++) {
                pdata = cpl_image_get_data_float(cpl_imagelist_get(ilist, j)) ;
                val = pdata[i] ;
                val3 = 1-2*val2*(pima[i]-val) ;
                if (val3 < 0.0) {
                    pdata[i] = val-pima[i] ;
                } else {
                    pdata[i]=((float)sqrt(val3)-1) / val2 ;
                }
            }
        }
    }
    /* Free and return */
    cpl_image_delete(ima) ;
    cpl_image_delete(imb) ;
    cpl_image_delete(imc) ;
    return 0 ;
}