cr2res_calib.c

/*
 * This file is part of the CR2RES 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
 */
 
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
 
/*-----------------------------------------------------------------------------
                                   Includes
 -----------------------------------------------------------------------------*/
 
#include <math.h>
#include <float.h>
#include <cpl.h>
#include "cr2res_calib.h"
#include "cr2res_bpm.h"
#include "cr2res_pfits.h"
#include "cr2res_io.h"
#include "cr2res_detlin.h"
#include "cr2res_utils.h"
 
/*-----------------------------------------------------------------------------
                                Functions prototypes
 -----------------------------------------------------------------------------*/
 
int cr2res_add_shotnoise(hdrl_image * in, double dit, int ndit, int ndsamples,
        int chip);
int cr2res_calib_subtract_interorder_column(hdrl_image * in,
    const hdrl_image * flat, cpl_size degree);
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/

 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
hdrl_imagelist * cr2res_calib_imagelist(
        const hdrl_imagelist    *   in,
        int                         chip,
        int                         clean_bad,
        int                         subtract_nolight_rows,
        int                         subtract_interorder_column,
        int                         cosmics_corr,
        const cpl_frame         *   flat,
        const cpl_frame         *   dark,
        const cpl_frame         *   bpm,
        const cpl_frame         *   detlin,
        const cpl_vector        *   dits,
        const cpl_vector        *   ndits,
        const cpl_vector        *   ndsamples)
{
    hdrl_imagelist      *   out ;
    double                  dit ;
    int                     ndit ;
    int                     ndsamp ;
    cpl_size                i ;
 
    /* Check Inputs */
    if (in == NULL) return NULL ;
 
    /* Initialise */
    dit = 0.0 ;
    ndit = 1 ;
    ndsamp = -1 ;
 
    /* Create calibrated image list */
    out = hdrl_imagelist_new() ;
 
    /* Loop on the images */
    for (i = 0; i < hdrl_imagelist_get_size(in); i++) {
        const hdrl_image *cur_ima;
        hdrl_image *cur_ima_calib;
        cur_ima = hdrl_imagelist_get(in, i) ;
        //if (dark != NULL)  dit = cpl_vector_get(dits, i) ;
        if (dits != NULL)  dit = cpl_vector_get(dits, i) ;
        if (ndits != NULL) ndit = (int)cpl_vector_get(ndits, i) ;
        if (ndsamples != NULL) ndsamp = (int)cpl_vector_get(ndsamples, i) ;
 
        /* Calibrate */
        cpl_msg_info(__func__, 
                "Apply calibrations for image #%"CPL_SIZE_FORMAT, i+1) ;
        cpl_msg_indent_more() ;
        if ((cur_ima_calib = cr2res_calib_image(cur_ima, chip, clean_bad, 
                        subtract_nolight_rows, subtract_interorder_column,
                        cosmics_corr, flat, dark, bpm, 
                        detlin, dit, ndit, ndsamp)) == NULL) {
            cpl_msg_error(__func__, "Failed to Calibrate the Data") ;
            hdrl_imagelist_delete(out) ;
            cpl_msg_indent_less() ;
            return NULL ;
        } else {
            /* All the calibrated image in the list */
            hdrl_imagelist_set(out, cur_ima_calib, i);
        }
        cpl_msg_indent_less() ;
    }
    return out ;
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
hdrl_image * cr2res_calib_image(
        const hdrl_image    *   in,
        int                     chip,
        int                     clean_bad,
        int                     subtract_nolight_rows,
        int                     subtract_interorder_column,
        int                     cosmics_corr,
        const cpl_frame     *   flat,
        const cpl_frame     *   dark,
        const cpl_frame     *   bpm,
        const cpl_frame     *   detlin,
        double                  dit,
        int                     ndit,
        int                     ndsamples)
{
    hdrl_image          *   out ;
    hdrl_image          *   calib ;
    hdrl_imagelist      *   calib_list ;
 
    /* Test entries */
    if (in == NULL) return NULL ;
    if (chip < 1 || chip > CR2RES_NB_DETECTORS) return NULL ;
 
    /* Create out image */
    out = hdrl_image_duplicate(in) ;
 
    /* Clean the bad pixels */
    if (bpm != NULL) {
        cpl_msg_info(__func__, "Correct the bad pixels") ;
        if (cr2res_bpm_set_and_correct_image(hdrl_image_get_image(out),
                    cpl_frame_get_filename(bpm), chip, clean_bad) != 0) {
            cpl_msg_error(__func__, "Cannot clean the bad pixels");
            hdrl_image_delete(out);
            return NULL ;
        }
    }
 
    /* Apply the non linearity correction */
    if (detlin != NULL) {
        cpl_msg_info(__func__, "Correct for the Non-Linearity") ;
        if ((calib_list = cr2res_io_load_DETLIN_COEFFS(
                        cpl_frame_get_filename(detlin), chip)) == NULL) {
            cpl_msg_error(__func__, "Cannot load the detlin") ;
            hdrl_image_delete(out);
            return NULL ;
        }
        /* Detlin correction */
        if (cr2res_detlin_correct(out, calib_list)) {
            hdrl_imagelist_delete(calib_list) ;
            hdrl_image_delete(out);
            cpl_msg_error(__func__, "Cannot correct for the Non-Linearity") ;
            return NULL ;
        }
    }
 
    /* Add shot-noise */
    cpl_msg_info(__func__, "Add shot-noise") ;
    if (cr2res_add_shotnoise(out, dit, ndit, ndsamples, chip)){
        cpl_msg_error(__func__, "Cannot add shot-noise") ;
        if(detlin != NULL)
            hdrl_imagelist_delete(calib_list) ;
        hdrl_image_delete(out);
        return NULL ;
    }
 
    /* Apply the dark */
    if (dark != NULL) {
        cpl_propertylist *plist;
        double dark_dit;
        cpl_msg_info(__func__, "Correct for the dark") ;
 
        /* Load the dark */
        if ((calib = cr2res_io_load_MASTER_DARK(cpl_frame_get_filename(dark), 
                        chip)) == NULL) {
            cpl_msg_error(__func__, "Cannot load the dark") ;
            hdrl_image_delete(out);
            return NULL ;
        }
        if (detlin != NULL) {
            cpl_msg_indent_more() ;
            cpl_msg_info(__func__, "Correct DARK for Non-Linearity") ;
            if (cr2res_detlin_correct(calib, calib_list)) {
                hdrl_imagelist_delete(calib_list) ;
                hdrl_image_delete(calib) ;
                hdrl_image_delete(out);
                cpl_msg_error(__func__,"Cannot correct DARK for Non-Linearity");
                cpl_msg_indent_less() ;
                return NULL ;
            }
            cpl_msg_indent_less() ;
        }
 
        /* Get the dark DIT */
        plist = cpl_propertylist_load(cpl_frame_get_filename(dark), 0);
        dark_dit = cr2res_pfits_get_dit(plist) ;
        cpl_propertylist_delete(plist) ;
        /* Multiply the dark by dit/dark_dit */
        hdrl_value hdrl_dit_corr = {dit/dark_dit, 0.0};
        hdrl_image_mul_scalar(calib, hdrl_dit_corr) ;
 
        /* Subtract the dark */
        if (hdrl_image_sub_image(out, calib) != CPL_ERROR_NONE) {
            cpl_msg_error(__func__, "Cannot apply the dark") ;
            hdrl_image_delete(calib) ;
            hdrl_image_delete(out);
            return NULL ;
        }
        hdrl_image_delete(calib) ;
    }
    if (detlin != NULL) {
       hdrl_imagelist_delete(calib_list) ;
    }
 
    /* Subtract residual bias/dark from vignetted rows at bottom */
    if (subtract_nolight_rows) {
        cpl_size i;
        cpl_image *img_tmp;
        cpl_msg_info(__func__, "Subtract median of bottom rows");
 
        img_tmp = cpl_image_collapse_median_create(
                hdrl_image_get_image(out), 0, CR2RES_NB_BPM_EDGEPIX, 
                CR2RES_DETECTOR_SIZE-CR2RES_NB_BPM_VIGN_BOTTOM);
        calib = hdrl_image_new(CR2RES_DETECTOR_SIZE,CR2RES_DETECTOR_SIZE);
        for (i=1; i<=CR2RES_DETECTOR_SIZE; i++) {
            hdrl_image_insert(calib, img_tmp, NULL, 1, i);
        }
        hdrl_image_sub_image(out, calib);
        hdrl_image_delete(calib) ;
        cpl_image_delete(img_tmp);
    }
 
    /* Apply the interorder correction */
    if (subtract_interorder_column) {
        if (flat == NULL) {
            cpl_msg_info(__func__, 
                "Skip subtracting fit to inter-order pixels (no flat-field)");
        } else {
            cpl_msg_info(__func__, "Subtract fit to inter-order pixels");
            if ((calib = cr2res_io_load_MASTER_FLAT(
                            cpl_frame_get_filename(flat), chip)) == NULL) {
                cpl_msg_error(__func__, "Cannot load the flat field") ;
                hdrl_image_delete(out);
                return NULL ;
            }
            if (cr2res_calib_subtract_interorder_column(out, calib, 0)){
                cpl_msg_error(__func__, "Could not subtract inter-order fit");
            }
            hdrl_image_delete(calib);
        }
    }
 
    /* Apply the flatfield */
    if (flat != NULL) {
        cpl_msg_info(__func__, "Correct for the flat field") ;
        /* Load the flat */
        if ((calib = cr2res_io_load_MASTER_FLAT(
                        cpl_frame_get_filename(flat), chip)) == NULL) {
            cpl_msg_error(__func__, "Cannot load the flat field") ;
            hdrl_image_delete(out);
            return NULL ;
        }
 
        /* Divide */
        if (hdrl_image_div_image(out, calib) != CPL_ERROR_NONE) {
            cpl_msg_error(__func__, "Cannot apply the flat field") ;
            hdrl_image_delete(calib) ;
            hdrl_image_delete(out);
            return NULL ;
        }
        hdrl_image_delete(calib) ;
    }
 
 
    /* Comics correction */
    if (cosmics_corr) {
        cpl_msg_info(__func__, "Apply the cosmics corrections") ;
        hdrl_parameter * cospar = hdrl_lacosmic_parameter_create(12.0,2.5,1);
        if (hdrl_lacosmic_parameter_verify(cospar)!=CPL_ERROR_NONE) 
            cpl_msg_warning(__func__,"invalid cospar");
        cpl_mask* cosmask = hdrl_lacosmic_edgedetect(out,cospar);
        cpl_mask_or(hdrl_image_get_mask(out), cosmask);
        cpl_mask_delete(cosmask);
        hdrl_parameter_delete(cospar);
    }
 
    return out ;
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
static double cr2res_utr_shot_factor(int ndsamples)
{
    double n ;
    if (ndsamples < 2) return 1.0 ;
    n = (double)ndsamples ;
    return 6.0 * (n * n + 1.0) / (5.0 * n * (n + 1.0)) ;
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
int cr2res_add_shotnoise(hdrl_image * in, double dit, int ndit, int ndsamples,
        int chip)
{
    cpl_image   *   error ;
    cpl_image   *   var ;
    double          gain, alpha, read_noise, min_rn ;
    const double    min_dit = 1.427 ;
    const double    lim_dit = 50.0 ;
    const double    lim_rn = 6.0 ;
 
    if (in == NULL) return -1 ;
    error = hdrl_image_get_error(in) ;
    if (hdrl_image_get_image_const(in) == NULL || error == NULL) {
        cpl_msg_error(__func__, "Broken input image") ;
        return -1 ;
    }
    if (ndit < 1) {
        cpl_msg_error(__func__, "Invalid NDIT: %d", ndit) ;
        return -1 ;
    }
 
    if      (chip == 1) gain = CR2RES_GAIN_CHIP1 ;
    else if (chip == 2) gain = CR2RES_GAIN_CHIP2 ;
    else if (chip == 3) gain = CR2RES_GAIN_CHIP3 ;
    else {
        cpl_msg_error(__func__, "Unknown detector") ;
        return -1 ;
    }
 
    /* Effective read noise per DIT in e-, from the CRIRES user manual */
    min_rn = (chip == 1) ? 11.0 : 12.0 ;
    if (dit <= min_dit)         read_noise = min_rn ;
    else if (dit >= lim_dit)    read_noise = lim_rn ;
    else read_noise = min_rn + (lim_rn-min_rn)*(dit-min_dit)/(lim_dit-min_dit);
 
    alpha = cr2res_utr_shot_factor(ndsamples) ;
 
    cpl_msg_debug(__func__, "chip:%d, gain:%g, ndit:%d, ndsamples:%d, "
            "alpha:%g, read_noise:%g", chip, gain, ndit, ndsamples, alpha,
            read_noise) ;
 
    if ((var = cpl_image_duplicate(hdrl_image_get_image_const(in))) == NULL) {
        cpl_msg_error(__func__, "Cannot duplicate the input image") ;
        return -1 ;
    }
    cpl_image_threshold(var, 0.0, DBL_MAX, 0.0, DBL_MAX) ;
    cpl_image_multiply_scalar(var, alpha / (gain * ndit)) ;
    cpl_image_add_scalar(var, read_noise * read_noise / (gain * gain * ndit)) ;
 
    /* The BPM should not be stored in the error image */
    cpl_image_accept_all(var) ;
 
    if (cpl_image_power(error, 2.0) != CPL_ERROR_NONE ||
            cpl_image_add(error, var) != CPL_ERROR_NONE ||
            cpl_image_power(error, 0.5) != CPL_ERROR_NONE) {
        cpl_msg_error(__func__, "Cannot combine the errors") ;
        cpl_image_delete(var) ;
        return -1 ;
    }
    cpl_image_delete(var) ;
    return 0 ;
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
int cr2res_calib_subtract_interorder_column(hdrl_image * in, 
                                    const hdrl_image * flat, cpl_size degree)
{
    // Define Parameters
    cpl_image * img;
    const cpl_image * flat_img;
    cpl_polynomial * poly;
    cpl_vector * tmp;
    cpl_size ncolumns, nrow;
    cpl_size i, j;
    int badpix;
    int badcols = 0;
    double value, median, mad;
 
    // Check Input 
    if (in == NULL || flat == NULL) return -1;
    if (hdrl_image_get_size_x(in) != hdrl_image_get_size_x(flat)) return -1;
    if (hdrl_image_get_size_y(in) != hdrl_image_get_size_y(flat)) return -1;
 
    // Loop over the columns
    img = hdrl_image_get_image(in);
    flat_img = hdrl_image_get_image_const(flat);
    ncolumns = cpl_image_get_size_x(img);
    nrow = cpl_image_get_size_y(img);
 
 
    for (i = CR2RES_NB_BPM_EDGEPIX + 1; i < ncolumns - CR2RES_NB_BPM_EDGEPIX;
         i++) {
        cpl_size npixel;
        cpl_matrix *px;
        cpl_vector *py;
        // Fill vectors and matrix
        px = cpl_matrix_new(1, nrow);
        py = cpl_vector_new(nrow);
        npixel = 0;
        for ( j = CR2RES_NB_BPM_VIGN_BOTTOM + 5; j < nrow -20 ; j++)
        {
            // Filter pixels that are rejected in the flat (i.e. between orders)
            // but not rejected by the image (i.e. not bad pixels)
            if (cpl_image_is_rejected(flat_img, i, j) 
                    && !cpl_image_is_rejected(img, i, j) &&
                    isnan(cpl_image_get(img, i, j, &badpix)) == 0){
                cpl_matrix_set(px, 0, npixel, j);
                cpl_vector_set(py, npixel, cpl_image_get(img, i, j, &badpix));
                npixel++;
            }
        }
        cpl_msg_debug(__func__, "Column %lld has %lld pixels for background",
                i,npixel);
        if (npixel < 10){
            badcols += 1;
            cpl_matrix_delete(px);
            cpl_vector_delete(py);
            continue;
        }
        cpl_matrix_set_size(px, 1, npixel);
        cpl_vector_set_size(py, npixel);
 
        // Reject outliers and filter
        median = cr2res_vector_get_mad(py, &mad);
        for (j=0;j<npixel;j++){
            // Filtering for the range: median - 1*median to median + 1*median
            if ((cpl_vector_get(py,j) > median*2) || (cpl_vector_get(py,j) < 0.0) )
                cpl_vector_set(py,j,median);
        }
        //tmp = cpl_vector_filter_median_create(py,50);
        tmp = cpl_vector_filter_lowpass_create(py,CPL_LOWPASS_LINEAR,100);
        cpl_vector_delete(py);
        py=tmp;
        if (cpl_msg_get_level() == CPL_MSG_DEBUG){
            tmp = cpl_vector_wrap(npixel, cpl_matrix_get_data(px));
            cpl_vector_save(tmp, "debug_background.fits",
                         CPL_TYPE_DOUBLE, NULL, CPL_IO_CREATE);
            cpl_vector_save(py, "debug_background.fits",
                         CPL_TYPE_DOUBLE, NULL, CPL_IO_EXTEND);
            cpl_vector_unwrap(tmp);
        }
 
        // Fit polynomial
        poly = cpl_polynomial_new(1);
        if (cpl_polynomial_fit(poly, px, NULL, py, NULL, CPL_FALSE, NULL,
                                                                    &degree)
                != CPL_ERROR_NONE) {
            badcols += 1;
            cpl_msg_warning(__func__,
            "Could not fit the background of column %"CPL_SIZE_FORMAT, i);
            cpl_error_reset();
            cpl_matrix_delete(px);
            cpl_vector_delete(py);
            cpl_polynomial_delete(poly);
            continue;
        }
 
        // Subtract column
        for ( j = 1; j < nrow + 1; j++)
        {
            value = cpl_image_get(img, i, j, &badpix);
            // Just skip this pixel if it is bad
            if (badpix) continue; 
            value -= cpl_polynomial_eval_1d(poly, j, NULL);
            cpl_image_set(img, i, j, value);
        }
 
        // Clean up
        cpl_matrix_delete(px);
        cpl_vector_delete(py);
        cpl_polynomial_delete(poly);
    }
 
    if (badcols){
        cpl_msg_warning(__func__, 
        "%d columns could not be background-subtracted from inter-order gaps",
            badcols);
    }
    return 0;
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
void cr2res_calib_add_var(
        hdrl_image      *   ima,
        const cpl_image *   var)
{
    cpl_image * err = hdrl_image_get_error(ima) ;
    cpl_image_power(err, 2.0) ;
    cpl_image_add(err, var) ;
    cpl_image_threshold(err, 0.0, DBL_MAX, 0.0, DBL_MAX) ;
    cpl_image_power(err, 0.5) ;
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
cpl_image * cr2res_calib_flat_relerr2(const hdrl_image * mflat)
{
    cpl_image * relerr2 ;
 
    if (mflat == NULL) return NULL ;
    relerr2 = cpl_image_divide_create(hdrl_image_get_error_const(mflat),
            hdrl_image_get_image_const(mflat)) ;
    if (relerr2 == NULL) return NULL ;
    cpl_image_power(relerr2, 2.0) ;
    return relerr2 ;
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
void cr2res_calib_flat_var(
        hdrl_image      *   ima,
        const cpl_image *   relerr2,
        double              sign)
{
    cpl_image * term = cpl_image_power_create(
            hdrl_image_get_image_const(ima), 2.0) ;
    cpl_image_multiply(term, relerr2) ;
    cpl_image_multiply_scalar(term, sign) ;
    cr2res_calib_add_var(ima, term) ;
    cpl_image_delete(term) ;
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
int cr2res_calib_dark_var(
        hdrl_imagelist      *   list,
        const cpl_frame     *   master_dark_frame,
        int                     chip,
        const hdrl_image    *   mflat,
        const cpl_vector    *   dits,
        double                  sign)
{
    hdrl_image          *   mdark ;
    cpl_propertylist    *   plist ;
    cpl_image           *   var ;
    double                  dark_dit ;
    cpl_size                i ;
 
    if (list == NULL || master_dark_frame == NULL || dits == NULL) return -1 ;
    if (cpl_vector_get_size(dits) != hdrl_imagelist_get_size(list)) return -1;
 
    if ((mdark = cr2res_io_load_MASTER_DARK(
                    cpl_frame_get_filename(master_dark_frame), chip)) == NULL)
        return -1 ;
    plist = cpl_propertylist_load(cpl_frame_get_filename(master_dark_frame),0);
    dark_dit = cr2res_pfits_get_dit(plist) ;
    cpl_propertylist_delete(plist) ;
    if (cpl_error_get_code() || dark_dit == 0.0) {
        hdrl_image_delete(mdark) ;
        return -1 ;
    }
 
    var = cpl_image_power_create(hdrl_image_get_error_const(mdark), 2.0) ;
    hdrl_image_delete(mdark) ;
    if (mflat != NULL) {
        cpl_image * flat2 = cpl_image_power_create(
                hdrl_image_get_image_const(mflat), 2.0) ;
        cpl_image_divide(var, flat2) ;
        cpl_image_delete(flat2) ;
    }
    /* Keep the error image BPM untouched */
    cpl_image_fill_rejected(var, 0.0) ;
    cpl_image_accept_all(var) ;
 
    for (i=0 ; i<hdrl_imagelist_get_size(list) ; i++) {
        double      scale = cpl_vector_get(dits, i) / dark_dit ;
        cpl_image * term = cpl_image_multiply_scalar_create(var,
                sign * scale * scale) ;
        cr2res_calib_add_var(hdrl_imagelist_get(list, i), term) ;
        cpl_image_delete(term) ;
    }
    cpl_image_delete(var) ;
    return 0 ;
}