irplib_mkmaster.c

/* $Id: irplib_mkmaster.c,v 1.6 2013-02-27 16:00:51 jtaylor Exp $
 *
 * This file is part of the irplib package
 * 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: jtaylor $
 * $Date: 2013-02-27 16:00:51 $
 * $Revision: 1.6 $
 * $Name: not supported by cvs2svn $
 */
 
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
 
/*-----------------------------------------------------------------------------
                                   Includes
 -----------------------------------------------------------------------------*/
 
#include <math.h>
#include <string.h>
#include "irplib_mkmaster.h"
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/

/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
/*-------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/
static cpl_vector *
irplib_imagelist_get_clean_mean_levels(const cpl_imagelist *iml,
                                       const double kappa,
                                       const int nclip,
                                       const double tolerance)
{
    int size = 0;
    int i = 0;
    cpl_vector *levels = NULL;
    double *pval = NULL;
    double mean = 0;
    double stdev = 0;
 
 
    cpl_error_ensure(iml != NULL, CPL_ERROR_NULL_INPUT, return (levels),
                     "Null input image list");
    cpl_error_ensure(kappa >= 0, CPL_ERROR_ILLEGAL_INPUT, return (levels),
                     "Must be kappa>0");
 
    size = cpl_imagelist_get_size(iml);
    levels = cpl_vector_new(size);
    pval = cpl_vector_get_data(levels);
 
    for (i = 0; i < size; i++) {
        const cpl_image *img = cpl_imagelist_get_const(iml, i);
        irplib_ksigma_clip(img, 1, 1, cpl_image_get_size_x(img),
                           cpl_image_get_size_y(img), nclip, kappa, tolerance,
                           &mean, &stdev);
        cpl_msg_info(cpl_func, "Ima %d mean level: %g", i + 1, mean);
        pval[i] = mean;
    }
 
 
    return levels;
}
 
/*-------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/
static cpl_error_code
irplib_imagelist_subtract_values(cpl_imagelist **iml, cpl_vector *values)
{
    int size = 0;
    int i = 0;
    double *pval = NULL;
 
    size = cpl_imagelist_get_size(*iml);
    pval = cpl_vector_get_data(values);
 
    for (i = 0; i < size; i++) {
        cpl_image *img = cpl_imagelist_get(*iml, i);
        cpl_image_subtract_scalar(img, pval[i]);
        cpl_imagelist_set(*iml, img, i);
    }
 
    return cpl_error_get_code();
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
static double
irplib_vector_ksigma(cpl_vector *values,
                     const double klow,
                     const double khigh,
                     int kiter)
{
    double mean /*= 0.0*/; /* Comment out to suppress cppcheck warning. */
    double sigma = 0.0;
    double *data = cpl_vector_get_data(values);
    int n = cpl_vector_get_size(values);
    int ngood = n;
    int i;
 
    /*
    * At first iteration the mean is taken as the median, and the
    * standard deviation relative to this value is computed.
    */
 
    mean = cpl_vector_get_median(values);
 
    for (i = 0; i < n; i++) {
        sigma += (mean - data[i]) * (mean - data[i]);
    }
    sigma = sqrt(sigma / (n - 1));
 
    while (kiter) {
        cpl_vector *accepted;
        int count = 0;
        for (i = 0; i < ngood; i++) {
            if (data[i] - mean < khigh * sigma &&
                mean - data[i] < klow * sigma) {
                data[count] = data[i];
                ++count;
            }
        }
 
        if (count == 0)  // This cannot happen at first iteration.
            break;       // So we can break: we have already computed a mean.
 
        /*
       * The mean must be computed even if no element was rejected
       * (count == ngood), because at first iteration median instead
       * of mean was computed.
       */
 
        accepted = cpl_vector_wrap(count, data);
        mean = cpl_vector_get_mean(accepted);
        if (count > 1) {
            sigma = cpl_vector_get_stdev(accepted);
        }
        cpl_vector_unwrap(accepted);
 
        if (count == ngood) {
            break;
        }
        ngood = count;
        --kiter;
    }
 
    return mean;
}
 

 
static cpl_image *
irplib_imagelist_ksigma_stack(const cpl_imagelist *imlist,
                              double klow,
                              double khigh,
                              int kiter)
{
    cpl_size ni, nx, ny, npix;
    cpl_image *out_ima = NULL;
    cpl_imagelist *loc_iml = NULL;
    double *pout_ima = NULL;
    cpl_image *image = NULL;
    const double **data = NULL;
    double *med = NULL;
    cpl_vector *time_line = NULL;
 
    double *ptime_line = NULL;
    cpl_size i, j;
    double mean_of_medians = 0;
 
    cpl_error_ensure(imlist != NULL, CPL_ERROR_NULL_INPUT, return (out_ima),
                     "Null input image list");
 
    ni = cpl_imagelist_get_size(imlist);
    loc_iml = cpl_imagelist_duplicate(imlist);
    image = cpl_imagelist_get(loc_iml, 0);
    nx = cpl_image_get_size_x(image);
    ny = cpl_image_get_size_y(image);
    npix = nx * ny;
 
    out_ima = cpl_image_new(nx, ny, CPL_TYPE_DOUBLE);
    pout_ima = cpl_image_get_data_double(out_ima);
 
    time_line = cpl_vector_new(ni);
 
    ptime_line = cpl_vector_get_data(time_line);
 
    cpl_error_ensure(ni < (cpl_size)(SIZE_MAX / sizeof(double *)),
                     CPL_ERROR_INCOMPATIBLE_INPUT, return (NULL),
                     "Image to big for allocation");
    data = cpl_calloc((size_t)ni, sizeof(double *));
    med = cpl_calloc((size_t)ni, sizeof(double));
 
    for (i = 0; i < ni; i++) {
        image = cpl_imagelist_get(loc_iml, i);
        med[i] = cpl_image_get_median(image);
        cpl_image_subtract_scalar(image, med[i]);
        data[i] = cpl_image_get_data_double(image);
        mean_of_medians += med[i];
    }
    mean_of_medians /= ni;
 
    for (i = 0; i < npix; i++) {
        for (j = 0; j < ni; j++) {
            ptime_line[j] = data[j][i];
        }
        pout_ima[i] = irplib_vector_ksigma(time_line, klow, khigh, kiter);
    }
 
    cpl_image_add_scalar(out_ima, mean_of_medians);
 
 
    cpl_free(data);
    cpl_free(med);
    cpl_vector_delete(time_line);
    cpl_imagelist_delete(loc_iml);
 
    return out_ima;
}
 
 
/*-------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/
cpl_image *
irplib_mkmaster_mean(cpl_imagelist *images,
                     const double kappa,
                     const int nclip,
                     const double tolerance,
                     const double klow,
                     const double khigh,
                     const int niter)
{
    cpl_image *master = NULL;
    cpl_vector *levels = NULL;
    double mean = 0;
    cpl_imagelist *iml = NULL;
 
    cpl_msg_info(cpl_func, "method mean");
    iml = cpl_imagelist_duplicate(images);
    levels =
        irplib_imagelist_get_clean_mean_levels(iml, kappa, nclip, tolerance);
    mean = cpl_vector_get_mean(levels);
    cpl_msg_info(cpl_func, "Master mean level: %g", mean);
 
    irplib_imagelist_subtract_values(&iml, levels);
 
    master = irplib_imagelist_ksigma_stack(iml, klow, khigh, niter);
    cpl_image_add_scalar(master, mean);
 
    cpl_vector_delete(levels);
    cpl_imagelist_delete(iml);
    return master;
}
 
/*-------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/
cpl_image *
irplib_mkmaster_median(cpl_imagelist *images,
                       const double kappa,
                       const int nclip,
                       const double tolerance)
{
    cpl_image *master = NULL;
    cpl_vector *levels = NULL;
    double mean = 0;
    cpl_imagelist *iml = NULL;
 
    cpl_msg_info(cpl_func, "method median");
    iml = cpl_imagelist_duplicate(images);
    levels =
        irplib_imagelist_get_clean_mean_levels(iml, kappa, nclip, tolerance);
 
    mean = cpl_vector_get_mean(levels);
    cpl_msg_info(cpl_func, "Master mean level: %g", mean);
    irplib_imagelist_subtract_values(&iml, levels);
 
    master = cpl_imagelist_collapse_median_create(iml);
 
    cpl_image_add_scalar(master, mean);
 
    cpl_vector_delete(levels);
    cpl_imagelist_delete(iml);
 
    return master;
}
 
/* Work in progress */
static cpl_error_code
irplib_mkmaster_dark_qc(const cpl_imagelist *raw_images,
                        cpl_imagelist *preproc_images,
                        const cpl_parameterlist *parameters,
                        const int pr_num_x,
                        const int pr_num_y,
                        const int pr_box_sx,
                        const int pr_box_sy,
                        const char *recipe_id,
                        cpl_table *qclog)
{
    cpl_ensure_code(qclog != NULL, CPL_ERROR_NULL_INPUT);
    cpl_ensure_code(recipe_id != NULL, CPL_ERROR_NULL_INPUT);
    cpl_ensure_code(parameters != NULL, CPL_ERROR_NULL_INPUT);
 
    if (pr_num_x != 0 && pr_num_y != 0 && pr_box_sx != 0 && pr_box_sy != 0) {
        int i;
        for (i = 0; i < cpl_imagelist_get_size(raw_images); i++) {
            cpl_image *current_dark =
                cpl_image_duplicate(cpl_imagelist_get_const(preproc_images, i));
            cpl_msg_info(cpl_func,
                         "Calculating QC parameters on raw dark frame %d", i);
            /* Here To be defined more general way to qc-log */
            /* UVES specific stuff: may be this function should not be put in irplib
      irplib_mdark_region_qc(current_dark, parameters, raw_images, recipe_id,qclog);
       */
            /* FIXME: still safe if irplib_mdark_region_qc is commented in? */
            cpl_image_delete(current_dark);
        }
    }
    return cpl_error_get_code();
}
 
/*-------------------------------------------------------------------------*/
/*-------------------------------------------------------------------------*/
static double
irplib_head_get_exptime(const cpl_propertylist *plist)
{
    double result = 0; /* Conversion from electrons to ADUs */
 
    result = cpl_propertylist_get_double(plist, "EXPTIME");
    cpl_ensure_code(result >= 0, CPL_ERROR_ILLEGAL_OUTPUT);
 
    return result;
}
 
/*-------------------------------------------------------------------------*/
/*-------------------------------------------------------------------------*/
static cpl_error_code
irplib_head_set_exptime(cpl_propertylist *plist, double exptime)
{
    cpl_propertylist_update_double(plist, "EXPTIME", exptime);
    cpl_propertylist_set_comment(plist, "EXPTIME", "Total integration time");
 
    return cpl_error_get_code();
}
 
static cpl_imagelist *
irplib_mkmaster_dark_fill_imagelist(const cpl_imagelist *raw_images,
                                    cpl_propertylist **raw_headers,
                                    const cpl_image *master_bias,
                                    double *mean_exptime)
{
    /* First process each input image and store the results in a
   new image list */
 
    cpl_imagelist *preproc_images = NULL;
    int i = 0;
    double min_exptime = 0;
    double max_exptime = 0;
 
    preproc_images = cpl_imagelist_new();
    for (i = 0; i < cpl_imagelist_get_size(raw_images); i++) {
        double exposure_time = 0.0;
        cpl_image *current_dark = NULL;
        const cpl_propertylist *current_header;
 
        current_dark =
            cpl_image_duplicate(cpl_imagelist_get_const(raw_images, i));
        current_header = raw_headers[i];
 
        /* Subtract master bias */
        if (master_bias != NULL) {
            cpl_msg_info(cpl_func, "Subtracting master bias");
            cpl_image_subtract(current_dark, master_bias);
        }
        else {
            cpl_msg_info(cpl_func, "Skipping bias subtraction");
        }
 
        exposure_time = irplib_head_get_exptime(current_header);
 
        /* Initialize/update min/max exposure time*/
        if (i == 0 || exposure_time < min_exptime) {
            min_exptime = exposure_time;
        }
        if (i == 0 || exposure_time > max_exptime) {
            max_exptime = exposure_time;
        }
 
        /* Do not normalize to unit exposure time */
        /*        If this is uncommented, then remember to also calculate the
     correct master dark exposure time below.
     irplib_msg("Normalizing from %f s to unit exposure time", exposure_time);
     check( cpl_image_divide_scalar(current_dark, exposure_time),
     "Error normalizing dark frame");   */
 
        /* Append to imagelist */
        cpl_imagelist_set(preproc_images, current_dark, i);
 
        /* Don't deallocate the image. It will be deallocated when
     the image list is deallocated */
        current_dark = NULL;
    }
 
 
    /* Check exposure times */
    cpl_msg_info(cpl_func,
                 "Exposure times range from %e s to %e s (%e %% variation)",
                 min_exptime, max_exptime,
                 100 * (max_exptime - min_exptime) / min_exptime);
 
    if ((max_exptime - min_exptime) / min_exptime > .001) {
        cpl_msg_warning(cpl_func, "Exposure times differ by %e %%",
                        100 * (max_exptime - min_exptime) / min_exptime);
    }
 
    /* compute correct exposure time */
    *mean_exptime = 0.5 * (max_exptime + min_exptime);
    return preproc_images;
}
 
 
cpl_image *
irplib_mdark_process_chip(const cpl_imagelist *raw_images,
                          cpl_propertylist **raw_headers,
                          const cpl_image *master_bias,
                          cpl_propertylist *mdark_header,
                          const cpl_parameterlist *parameters,
                          const char *recipe_id,
                          cpl_table *qclog,
                          const int do_qc,
                          const char *STACK_METHOD,
                          const double STACK_KLOW,
                          const double STACK_KHIGH,
                          const int STACK_NITER,
                          const int pr_num_x,
                          const int pr_num_y,
                          const int pr_box_sx,
                          const int pr_box_sy)
{
    cpl_image *master_dark = NULL; /* Result */
    cpl_image *current_dark = NULL;
    cpl_imagelist *preproc_images = NULL;
    double mean_exptime = 0;
 
    /* First process each input image and store the results in a
   new image list */
    preproc_images =
        irplib_mkmaster_dark_fill_imagelist(raw_images, raw_headers,
                                            master_bias, &mean_exptime);
    if (do_qc) {
        /* Here we should compute QC but a a better way to log it is TBD  */
        irplib_mkmaster_dark_qc(raw_images, preproc_images, parameters,
                                pr_num_x, pr_num_y, pr_box_sx, pr_box_sy,
                                recipe_id, qclog);
    }
    /* Get median stack of input darks */
    if (strcmp(STACK_METHOD, "MEDIAN") == 0) {
        cpl_msg_info(cpl_func, "Calculating stack median");
        master_dark = cpl_imagelist_collapse_median_create(preproc_images);
    }
    else {
        cpl_msg_info(cpl_func, "Calculating stack mean");
        master_dark = irplib_imagelist_ksigma_stack(preproc_images, STACK_KLOW,
                                                    STACK_KHIGH, STACK_NITER);
    }
    irplib_head_set_exptime(mdark_header, mean_exptime);
 
    cpl_image_delete(current_dark);
    cpl_imagelist_delete(preproc_images);
    if (cpl_error_get_code() != CPL_ERROR_NONE) {
        cpl_image_delete(master_dark);
        master_dark = NULL;
    }
 
    return master_dark;
}