hdrl_resample.c

Go to the documentation of this file.

/*
 * This file is part of the HDRL
 * Copyright (C) 2020 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  02110-1301  USA
 */
 
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
 
/*-----------------------------------------------------------------------------
                                    Includes
 -----------------------------------------------------------------------------*/
#include <cpl.h>
#include <string.h>
#include <math.h>
 
#include "hdrl_resample.h"
#include "hdrl_utils.h"
 
#ifndef _OPENMP
#define omp_get_max_threads() 1
#define omp_get_thread_num() 0
#else
#include <omp.h>
#endif
 
#include <sys/time.h>
 
/*-----------------------------------------------------------------------------
                                Defines
 -----------------------------------------------------------------------------*/
 
/* Not very elegant but it removes compiler warnings */
#if !defined(_POSIX_C_SOURCE)
#define _POSIX_C_SOURCE 200112L
#elif (_POSIX_C_SOURCE - 0) < 200112L
#undef _POSIX_C_SOURCE
#define _POSIX_C_SOURCE 200112L
#endif
 
 
/* maximum keyword length for FITS headers, including '\0' */
#define KEYWORD_LENGTH 81
 
/* Default field margin (in percent), if the user does not specify any. 5
 * percent is also used in the software package swarp */
#define FIELDMARGIN 5.
 
/* use bits 0-52 for the value (the pixel table row), this allows to convert
 * pixel tables with up to 9e15 pixels  into a pixel grid    */
#define PT_IDX_MASK 0x1FFFFFFFFFFFFFll
 
/* use bits 53-62 to store the thread ID, this allows parallelization with up to
 *  1024 cores */
#define XMAP_BITMASK 0x3FFll /* 1023 */
#define XMAP_LSHIFT 53ll
 
typedef struct {
  double crpix1, crpix2;         /* the reference point                       */
  double crval1, crval2;         /* the coordinate values at the ref. point   */
  double cd11, cd12, cd21, cd22; /* the four elements of the CDi_j matrix     */
  double cddet;                  /* the determinant of the CDi_j matrix       */
} hdrl_resample_smallwcs;
 
typedef struct {
  unsigned int npix; /* number of pixels in this grid point                   */
  cpl_size *pix;     /* the row number(s) in the pixel table                  */
} hdrl_resample_pixels_ext;
 
typedef struct {
   /* The pixel grid array, elements can be:
  0:        empty
  positive: row_number in the pixel table
  negative: -(i_ext+1) in the extension array,
  bits 53-62 contain the map index   */
  cpl_size *pix;
  cpl_size nx;             /* horizontal spatial size                         */
  cpl_size ny;             /* vertical spatial size                           */
  cpl_size nz;             /* size in dispersion direction                    */
  unsigned short nmaps;    /* number of extension maps                        */
  cpl_size *nxmap;         /* number of filled pixels in the extension maps   */
  cpl_size *nxalloc;       /* number of allocated pixels in the extension maps */
  hdrl_resample_pixels_ext **xmaps; /* the extension maps                              */
} hdrl_resample_pixgrid;
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
/*-----------------------------------------------------------------------------
                        RESAMPLE Parameter Definition
 -----------------------------------------------------------------------------*/
 
typedef struct {
  HDRL_PARAMETER_HEAD;
  hdrl_resample_outgrid method;
  double              delta_ra ;     /* step size in right ascension [deg] */
  double              delta_dec ;    /* step size in declination [deg] */
  double              delta_lambda ; /* step size in wavelength direction [m] */
  const cpl_wcs*      wcs ;          /* World Coordinate System */
  cpl_boolean         recalc_limits;
  double              ra_min;        /* Minimal Right ascension [deg] */
  double              ra_max;        /* Maximal Right ascension [deg] */
  double              dec_min;       /* Minimal Declination [deg] */
  double              dec_max;       /* Maximal Declination [deg] */
  double              lambda_min;    /* Minimal wavelength [m] */
  double              lambda_max;    /* Maximal wavelength [m] */
  double              fieldmargin;   /* Field margin to add [percent] */
} hdrl_resample_outgrid_parameter;
 
/* Parameter type */
static hdrl_parameter_typeobj hdrl_resample_outgrid_parameter_type = {
    HDRL_PARAMETER_RESAMPLE_OUTGRID,          /* type */
    (hdrl_alloc *)&cpl_malloc,                /* fp_alloc */
    (hdrl_free *)&cpl_free,                   /* fp_free */
    NULL,                                     /* fp_destroy */
    sizeof(hdrl_resample_outgrid_parameter),  /* obj_size */
};
 
typedef struct {
  HDRL_PARAMETER_HEAD;
  hdrl_resample_method method;
  int loop_distance;
  cpl_boolean use_errorweights;
  double drizzle_pix_frac_x;
  double drizzle_pix_frac_y;
  double drizzle_pix_frac_lambda;
  double renka_critical_radius;
  int lanczos_kernel_size;
} hdrl_resample_method_parameter;
 
/* Parameter type */
static hdrl_parameter_typeobj hdrl_resample_method_parameter_type = {
    HDRL_PARAMETER_RESAMPLE_METHOD,           /* type */
    (hdrl_alloc *)&cpl_malloc,                /* fp_alloc */
    (hdrl_free *)&cpl_free,                   /* fp_free */
    NULL,                                     /* fp_destroy */
    sizeof(hdrl_resample_method_parameter),   /* obj_size */
};
 
/*-----------------------------------------------------------------------------
                            Functions
 -----------------------------------------------------------------------------*/
 
static hdrl_resample_result *
hdrl_resample_cube(const cpl_table *ResTable,
             hdrl_resample_method_parameter *aParams_method,
             hdrl_resample_outgrid_parameter *aParams_outputgrid,
             hdrl_resample_pixgrid **aGrid);
 
static cpl_error_code
hdrl_resample_inputtable_verify(const cpl_table *ResTable);
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
static cpl_error_code
hdrl_resample_wcs_print(const cpl_wcs *wcs)
{
  cpl_ensure_code(wcs, CPL_ERROR_NULL_INPUT);
 
  const cpl_array *crval = cpl_wcs_get_crval(wcs);
  const cpl_array *crpix = cpl_wcs_get_crpix(wcs);
  const cpl_array *ctype = cpl_wcs_get_ctype(wcs);
  const cpl_array *cunit = cpl_wcs_get_cunit(wcs);
 
  const cpl_matrix *cd = cpl_wcs_get_cd(wcs);
  const cpl_array *dims = cpl_wcs_get_image_dims(wcs);
  int naxis = cpl_wcs_get_image_naxis(wcs);
 
  cpl_msg_debug(cpl_func, "NAXIS:  %d", naxis);
 
  int testerr = 0;
 
  cpl_msg_indent_more();
  /* Check NAXIS */
  for (cpl_size i = 0; i < naxis; i++) {
      cpl_msg_debug(cpl_func, "NAXIS%lld: %d", i + 1,
           cpl_array_get_int(dims, i, &testerr));
  }
  cpl_msg_indent_less();
 
  double cd11 = cpl_matrix_get(cd, 0, 0);
  double cd12 = cpl_matrix_get(cd, 0, 1);
  double cd21 = cpl_matrix_get(cd, 1, 0);
  double cd22 = cpl_matrix_get(cd, 1, 1);
  double crpix1 = cpl_array_get_double(crpix, 0, &testerr);
  double crpix2 = cpl_array_get_double(crpix, 1, &testerr);
  double crval1 = cpl_array_get_double(crval, 0, &testerr);
  double crval2 = cpl_array_get_double(crval, 1, &testerr);
 
  cpl_msg_debug(cpl_func, "1st and 2nd dimension");
  cpl_msg_indent_more();
  cpl_msg_debug(cpl_func, "CD1_1:  %g", cd11);
  cpl_msg_debug(cpl_func, "CD1_2:  %g", cd12);
  cpl_msg_debug(cpl_func, "CD2_1:  %g", cd21);
  cpl_msg_debug(cpl_func, "CD2_2:  %g", cd22);
 
  cpl_msg_debug(cpl_func, "CRPIX1: %g", crpix1);
  cpl_msg_debug(cpl_func, "CRPIX2: %g", crpix2);
  cpl_msg_debug(cpl_func, "CRVAL1: %g", crval1);
  cpl_msg_debug(cpl_func, "CRVAL2: %g", crval2);
  if (ctype) {
      cpl_msg_debug(cpl_func, "CTYPE1: %s", cpl_array_get_string(ctype, 0));
      cpl_msg_debug(cpl_func, "CTYPE2: %s", cpl_array_get_string(ctype, 1));
  }
 
  if (cunit) {
      cpl_msg_debug(cpl_func, "CUNIT1: %s", cpl_array_get_string(cunit, 0));
      cpl_msg_debug(cpl_func, "CUNIT2: %s", cpl_array_get_string(cunit, 1));
  }
  cpl_msg_indent_less();
 
  /* Is it a 3D cube or a 2D image */
  cpl_size cdncol = cpl_matrix_get_ncol(cd);
  if (cdncol == 3) {
 
      double cd13 = cpl_matrix_get(cd, 0, 2);
      double cd23 = cpl_matrix_get(cd, 1, 2);
      double cd31 = cpl_matrix_get(cd, 2, 0);
      double cd32 = cpl_matrix_get(cd, 2, 1);
      double cd33 = cpl_matrix_get(cd, 2, 2);
      double crval3 = cpl_array_get_double(crval, 2, &testerr);
      double crpix3 = cpl_array_get_double(crpix, 2, &testerr);
 
      cpl_msg_debug(cpl_func, "3rd dimension");
      cpl_msg_indent_more();
      cpl_msg_debug(cpl_func, "CD1_3:  %g", cd13);
      cpl_msg_debug(cpl_func, "CD2_3:  %g", cd23);
      cpl_msg_debug(cpl_func, "CD3_1:  %g", cd31);
      cpl_msg_debug(cpl_func, "CD3_2:  %g", cd32);
      cpl_msg_debug(cpl_func, "CD3_3:  %g", cd33);
 
      cpl_msg_debug(cpl_func, "CRPIX3: %g", crpix3);
      cpl_msg_debug(cpl_func, "CRVAL3: %g", crval3);
 
      if (ctype) {
          cpl_msg_debug(cpl_func, "CTYPE3: %s", cpl_array_get_string(ctype, 2));
      }
 
      if (cunit) {
          cpl_msg_debug(cpl_func, "CUNIT3: %s", cpl_array_get_string(cunit, 2));
      }
 
      cpl_msg_indent_less();
  }
 
  return cpl_error_get_code();
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
static cpl_error_code
hdrl_resample_outgrid_parameter_print(hdrl_resample_outgrid_parameter* p)
{
 
  cpl_ensure_code(p, CPL_ERROR_NULL_INPUT);
  /* Content of the outgrid parameter structure */
 
  cpl_msg_indent_more();
  cpl_msg_debug(cpl_func, "delta_ra:       %g", p->delta_ra);
  cpl_msg_debug(cpl_func, "delta_dec:      %g", p->delta_dec);
  cpl_msg_debug(cpl_func, "delta_lambda:   %g", p->delta_lambda);
  cpl_msg_debug(cpl_func, "ra_min:         %g", p->ra_min);
  cpl_msg_debug(cpl_func, "ra_max:         %g", p->ra_max);
  cpl_msg_debug(cpl_func, "dec_min:        %g", p->dec_min);
  cpl_msg_debug(cpl_func, "dec_max:        %g", p->dec_max);
  cpl_msg_debug(cpl_func, "lambda_min:     %g", p->lambda_min);
  cpl_msg_debug(cpl_func, "lambda_max:     %g", p->lambda_max);
  cpl_msg_debug(cpl_func, "fieldmargin:    %g", p->fieldmargin);
  cpl_msg_debug(cpl_func, "recalc_limits:  %d", p->recalc_limits);
 
  /* World Coordinate System */
  hdrl_resample_wcs_print(p->wcs);
  cpl_msg_indent_less();
  return cpl_error_get_code();
 
}
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
static cpl_error_code
hdrl_resample_method_parameter_print(hdrl_resample_method_parameter* p)
{
    if (hdrl_resample_parameter_method_verify((hdrl_parameter*)p) != CPL_ERROR_NONE) {
        return cpl_error_get_code();
    }
 
    /* Content of the outgrid parameter structure */
 
    cpl_msg_indent_more();
    if(p->method == HDRL_RESAMPLE_METHOD_NEAREST) {
        cpl_msg_debug(cpl_func, "method:                   %s", "NEAREST" );
    } else if (p->method == HDRL_RESAMPLE_METHOD_RENKA) {
        cpl_msg_debug(cpl_func, "method:                   %s", "RENKA" );
        cpl_msg_debug(cpl_func, "loop_distance:            %d", p->loop_distance);
        cpl_msg_debug(cpl_func, "use_errorweights:         %s",
                     p->use_errorweights == CPL_TRUE ? "TRUE" : "FALSE");
        cpl_msg_debug(cpl_func, "renka_critical_radius:    %g", p->renka_critical_radius);
    } else if (p->method == HDRL_RESAMPLE_METHOD_LINEAR) {
        cpl_msg_debug(cpl_func, "method:                   %s", "LINEAR" );
        cpl_msg_debug(cpl_func, "loop_distance:            %d", p->loop_distance);
        cpl_msg_debug(cpl_func, "use_errorweights:         %s",
                     p->use_errorweights == CPL_TRUE ? "TRUE" : "FALSE");
    } else if (p->method == HDRL_RESAMPLE_METHOD_QUADRATIC) {
        cpl_msg_debug(cpl_func, "method:                   %s", "QUADRATIC)" );
        cpl_msg_debug(cpl_func, "loop_distance:            %d", p->loop_distance);
        cpl_msg_debug(cpl_func, "use_errorweights:         %s",
                     p->use_errorweights == CPL_TRUE ? "TRUE" : "FALSE");
    } else if (p->method == HDRL_RESAMPLE_METHOD_DRIZZLE) {
        cpl_msg_debug(cpl_func, "method:                   %s", "DRIZZLE" );
        cpl_msg_debug(cpl_func, "loop_distance:            %d", p->loop_distance);
        cpl_msg_debug(cpl_func, "use_errorweights:         %s",
                     p->use_errorweights == CPL_TRUE ? "TRUE" : "FALSE");
        cpl_msg_debug(cpl_func, "drizzle_pix_frac_x:       %g", p->drizzle_pix_frac_x);
        cpl_msg_debug(cpl_func, "drizzle_pix_frac_y:       %g", p->drizzle_pix_frac_y);
        cpl_msg_debug(cpl_func, "drizzle_pix_frac_lambda:  %g", p->drizzle_pix_frac_lambda);
    } else if (p->method == HDRL_RESAMPLE_METHOD_LANCZOS) {
        cpl_msg_debug(cpl_func, "method:                   %s", "LANCZOS" );
        cpl_msg_debug(cpl_func, "loop_distance:            %d", p->loop_distance);
        cpl_msg_debug(cpl_func, "use_errorweights:         %s",
                     p->use_errorweights == CPL_TRUE ? "TRUE" : "FALSE");
        cpl_msg_debug(cpl_func, "lanczos_kernel_size:      %d", p->lanczos_kernel_size);
    }
 
    cpl_msg_indent_less();
    return cpl_error_get_code();
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
cpl_error_code
hdrl_wcs_xy_to_radec(const cpl_wcs *wcs, double x, double y, double *ra,
         double *dec) {
  cpl_ensure_code(wcs && ra && dec, CPL_ERROR_NULL_INPUT);
  double * xy = NULL;
  double * radec = NULL;
  cpl_matrix * from = NULL;
  cpl_matrix * to = NULL;
  cpl_array * status = NULL;
 
  /* Load up the information */
  int naxis = cpl_wcs_get_image_naxis(wcs);
  from = cpl_matrix_new(1, naxis);
  xy = cpl_matrix_get_data(from);
  xy[0] = x;
  xy[1] = y;
 
  /* Call the conversion routine */
 
  cpl_wcs_convert(wcs, from, &to, &status, CPL_WCS_PHYS2WORLD);
 
  /* Pass it back now */
 
  radec = cpl_matrix_get_data(to);
  *ra = radec[0];
  *dec = radec[1];
 
  /* Tidy and exit */
 
  cpl_matrix_delete(from);
  cpl_matrix_delete(to);
  cpl_array_delete(status);
  return cpl_error_get_code();
}
 
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
static double
hdrl_resample_pfits_get_crpix(const cpl_propertylist *aHeaders, unsigned int aAxis)
{
  cpl_errorstate prestate = cpl_errorstate_get();
  cpl_ensure(aHeaders, CPL_ERROR_NULL_INPUT, 0);
 
  char keyword[KEYWORD_LENGTH];
  snprintf(keyword, KEYWORD_LENGTH, "CRPIX%u", aAxis);
  const double value = cpl_propertylist_get_double(aHeaders, keyword);
  /* default to 0.0 as per FITS Standard v3.0 */
  cpl_ensure(cpl_errorstate_is_equal(prestate), cpl_error_get_code(), 0.0);
  return value;
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
static double
hdrl_resample_pfits_get_crval(const cpl_propertylist *aHeaders, unsigned int aAxis)
{
  cpl_errorstate prestate = cpl_errorstate_get();
  cpl_ensure(aHeaders, CPL_ERROR_NULL_INPUT, 0);
 
  char keyword[KEYWORD_LENGTH];
  snprintf(keyword, KEYWORD_LENGTH, "CRVAL%u", aAxis);
  const double value = cpl_propertylist_get_double(aHeaders, keyword);
  /* default to 0.0 as per FITS Standard v3.0 */
  cpl_ensure(cpl_errorstate_is_equal(prestate), cpl_error_get_code(), 0.0);
  return value;
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
static double
hdrl_resample_pfits_get_cd(const cpl_propertylist *aHeaders, unsigned int aAxisI,
          unsigned int aAxisJ)
{
  cpl_errorstate prestate = cpl_errorstate_get();
  cpl_ensure(aHeaders, CPL_ERROR_NULL_INPUT, 0);
 
  char keyword[KEYWORD_LENGTH];
  snprintf(keyword, KEYWORD_LENGTH, "CD%u_%u", aAxisI, aAxisJ);
  const double value = cpl_propertylist_get_double(aHeaders, keyword);
  /* default to 0.0 as per FITS Standard v3.0 */
  cpl_ensure(cpl_errorstate_is_equal(prestate), cpl_error_get_code(), 0.0);
  return value;
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
static hdrl_resample_smallwcs *
hdrl_resample_smallwcs_new(cpl_propertylist *aHeader)
{
 
  cpl_ensure(aHeader, CPL_ERROR_NULL_INPUT, NULL);
  hdrl_resample_smallwcs *wcs = cpl_calloc(1, sizeof(hdrl_resample_smallwcs));
 
  cpl_errorstate prestate = cpl_errorstate_get();
  wcs->crpix1 = hdrl_resample_pfits_get_crpix(aHeader, 1);
  wcs->crpix2 = hdrl_resample_pfits_get_crpix(aHeader, 2);
  wcs->crval1 = hdrl_resample_pfits_get_crval(aHeader, 1);
  wcs->crval2 = hdrl_resample_pfits_get_crval(aHeader, 2);
  if (!cpl_errorstate_is_equal(prestate)) {
      /* all these headers default to 0.0 following the FITS *
       * Standard, so we can ignore any errors set here      */
      cpl_errorstate_set(prestate);
  }
 
  prestate = cpl_errorstate_get();
  wcs->cd11 = hdrl_resample_pfits_get_cd(aHeader, 1, 1);
  wcs->cd22 = hdrl_resample_pfits_get_cd(aHeader, 2, 2);
  wcs->cd12 = hdrl_resample_pfits_get_cd(aHeader, 1, 2);
  wcs->cd21 = hdrl_resample_pfits_get_cd(aHeader, 2, 1);
  if (!cpl_errorstate_is_equal(prestate) &&
      wcs->cd11 == 0. && wcs->cd12 == 0. && wcs->cd21 == 0. && wcs->cd22 == 0.) {
      /* FITS Standard says to handle the CD matrix like the PC *
       * matrix in this case, with 1 for the diagonal elements  */
      wcs->cd11 = wcs->cd22 = wcs->cddet = 1.;
      cpl_errorstate_set(prestate); /* not a real error */
  }
  wcs->cddet = wcs->cd11 * wcs->cd22 - wcs->cd12 * wcs->cd21;
  if (wcs->cddet == 0.) {
      cpl_error_set(cpl_func, CPL_ERROR_SINGULAR_MATRIX);
  }
 
  return wcs;
} /* hdrl_resample_smallwcs_new() */
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
static void
hdrl_resample_pixgrid_delete(hdrl_resample_pixgrid *aGrid)
{
  if (!aGrid) {
      return;
  }
  cpl_free(aGrid->pix);
  aGrid->pix = NULL;
  /* clean up extension maps */
  unsigned short ix;
  for (ix = 0; ix < aGrid->nmaps; ix++) {
      cpl_size iext; /* index in this extension map */
      for (iext = 0; iext < aGrid->nxalloc[ix]; iext++) {
      cpl_free(aGrid->xmaps[ix][iext].pix);
      } /* for iext (all allocated pixels in this extension map) */
 
      cpl_free(aGrid->xmaps[ix]);
  } /* for ix (all extension maps) */
  cpl_free(aGrid->xmaps);
  aGrid->xmaps = NULL;
  cpl_free(aGrid->nxmap);
  aGrid->nxmap = NULL;
  cpl_free(aGrid->nxalloc);
  aGrid->nxalloc = NULL;
  cpl_free(aGrid);
} /* hdrl_resample_pixgrid_delete() */
 
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
hdrl_resample_result *
hdrl_resample_compute(const cpl_table *ResTable,
              hdrl_parameter *method,
              hdrl_parameter *outputgrid,
              const cpl_wcs* wcs)
{
  cpl_ensure(ResTable && method, CPL_ERROR_NULL_INPUT, NULL);
  cpl_ensure(wcs != NULL, CPL_ERROR_NULL_INPUT, NULL);
 
  if (hdrl_resample_inputtable_verify(ResTable)) {
        return NULL;
  }
  if (hdrl_resample_parameter_method_verify((hdrl_parameter *)method)) {
      return NULL;
  }
  if (hdrl_resample_parameter_outgrid_verify((hdrl_parameter *)outputgrid)) {
      return NULL;
  }
 
  if (hdrl_resample_inputtable_verify(ResTable)) {
      return NULL;
  }
 
  hdrl_resample_outgrid_parameter *aParams_outputgrid =
      (hdrl_resample_outgrid_parameter *) outputgrid;
  hdrl_resample_method_parameter *aParams_method =
      (hdrl_resample_method_parameter *) method;
 
  /*Assign the wcs*/
  aParams_outputgrid->wcs = wcs;
 
  /* Recalculate the limits if the user did not specify any */
  cpl_boolean recalc_limits = aParams_outputgrid->recalc_limits;
 
  if (recalc_limits == CPL_TRUE) {
      double ramin =  cpl_table_get_column_min(ResTable, HDRL_RESAMPLE_TABLE_RA);
      double ramax =  cpl_table_get_column_max(ResTable, HDRL_RESAMPLE_TABLE_RA);
      double decmin = cpl_table_get_column_min(ResTable, HDRL_RESAMPLE_TABLE_DEC);
      double decmax = cpl_table_get_column_max(ResTable, HDRL_RESAMPLE_TABLE_DEC);
      double lmin =   cpl_table_get_column_min(ResTable, HDRL_RESAMPLE_TABLE_LAMBDA);
      double lmax =   cpl_table_get_column_max(ResTable, HDRL_RESAMPLE_TABLE_LAMBDA);
 
      /* We have the rare case that the image spans over ra = 0.*/
      if(ramax - ramin > 180){
      const double *ra = cpl_table_get_data_double_const(ResTable,
                                 HDRL_RESAMPLE_TABLE_RA);
 
      /* set to extreme values for a start */
      ramin = 0.;
      ramax = 360.;
      cpl_size nrow = cpl_table_get_nrow(ResTable);
 
      for (cpl_size i = 0; i < nrow; i++) {
          if (ra[i] > ramin && ra[i] <= 180.) ramin = ra[i]; /* get the maximum */
          if (ra[i] < ramax && ra[i] >  180.) ramax = ra[i]; /* get the minimum */
      }
      }
 
      aParams_outputgrid->ra_min = ramin;
      aParams_outputgrid->ra_max = ramax;
      aParams_outputgrid->dec_min = decmin;
      aParams_outputgrid->dec_max = decmax;
      aParams_outputgrid->lambda_min = lmin;
      aParams_outputgrid->lambda_max = lmax;
  }
 
  cpl_msg_debug(cpl_func, "Content of the outgrid parameter structure "
      "hdrl_resample_outgrid_parameter when resampling starts:");
  hdrl_resample_outgrid_parameter_print(aParams_outputgrid);
 
  cpl_msg_debug(cpl_func, "Content of the method parameter structure "
      "hdrl_resample_method_parameter when resampling starts:");
  hdrl_resample_method_parameter_print(aParams_method);
 
  /* create cube and cast to generic pointer to save code duplication */
  hdrl_resample_result *cube = NULL;
  hdrl_resample_pixgrid *grid = NULL;
 
  cpl_msg_debug(cpl_func, "Resampling starts ...");
  cpl_msg_indent_more();
  cube = hdrl_resample_cube(ResTable,
                  (hdrl_resample_method_parameter *)aParams_method,
                  (hdrl_resample_outgrid_parameter *)aParams_outputgrid,
                  &grid);
  cpl_msg_indent_less();
  cpl_ensure(cube, CPL_ERROR_NULL_INPUT, NULL);
  cpl_ensure(cube->header, CPL_ERROR_NULL_INPUT, NULL);
  cpl_ensure(cube->himlist, CPL_ERROR_NULL_INPUT, NULL);
 
  /* Cleanup wcs for 2D / 3D case */
 
  if(hdrl_imagelist_get_size(cube->himlist) == 1) { /* 2D case */
      cpl_propertylist *header = cpl_propertylist_new();
      cpl_wcs* wcs_local = cpl_wcs_new_from_propertylist(cube->header);
      hdrl_wcs_to_propertylist(wcs_local, header, CPL_TRUE);
      cpl_wcs_delete(wcs_local);
      cpl_propertylist_delete(cube->header);
      cpl_propertylist_set_comment(header, "CTYPE1", "Gnomonic projection");
      cpl_propertylist_set_comment(header, "CTYPE2", "Gnomonic projection");
      cube->header = header;
  }
 
 
  hdrl_resample_pixgrid_delete(grid);
  return cube;
} /* hdrl_resample_compute() */
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
void
hdrl_resample_result_delete(hdrl_resample_result *aCube)
{
  /* if the cube does not exists at all, we don't need to do anything */
  if (!aCube) {
      return;
  }
 
  /* checks for the existence of the sub-images *
   * are done in the CPL functions              */
 
  hdrl_imagelist_delete(aCube->himlist);
  aCube->himlist = NULL;
  /* delete the FITS header, too */
  cpl_propertylist_delete(aCube->header);
  aCube->header = NULL;
 
  cpl_free(aCube);
} /* hdrl_resample_result_delete() */
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
static cpl_error_code
hdrl_resample_wcs_projplane_from_celestial(hdrl_resample_outgrid_parameter *
                  aParams_outputgrid, double aRA,
                  double aDEC, double *aX, double *aY)
{
  cpl_ensure_code(aParams_outputgrid && aX && aY, CPL_ERROR_NULL_INPUT);
  /* make sure that the header represents TAN */
  //const char * type1 = aParams_outputgrid->limits.ctype1;
  //const char * type2 = aParams_outputgrid->limits.ctype2;
  /*
  cpl_ensure_code(type1 && type2 && !strncmp(type1, "RA---TAN", 9) &&
                  !strncmp(type2, "DEC--TAN", 9),
                  CPL_ERROR_UNSUPPORTED_MODE);
   */
 
  int err = 0;
  const cpl_array *crval = cpl_wcs_get_crval(aParams_outputgrid->wcs);
  double crval1 = cpl_array_get_double(crval, 0, &err);
  double crval2 = cpl_array_get_double(crval, 1, &err);
 
  /* spherical coordinate shift/translation */
  double a = aRA / CPL_MATH_DEG_RAD,  /* RA in radians */
      d = aDEC / CPL_MATH_DEG_RAD, /* DEC in radians */
      /* alpha_p and delta_p in Paper II (in radians) */
      //ap = aParams_outputgrid->limits.crval1 / CPL_MATH_DEG_RAD,
      //dp = aParams_outputgrid->limits.crval2 / CPL_MATH_DEG_RAD,
      ap = crval1 / CPL_MATH_DEG_RAD,
      dp = crval2 / CPL_MATH_DEG_RAD,
      phi = atan2(-cos(d) * sin(a - ap),
          sin(d) * cos(dp) - cos(d) * sin(dp) * cos(a-ap))
          + 180 / CPL_MATH_DEG_RAD,
          theta = asin(sin(d) * sin(dp) + cos(d) * cos(dp) * cos(a-ap)),
          R_theta = CPL_MATH_DEG_RAD / tan(theta);
  /* spherical deprojection */
  *aX = R_theta * sin(phi),
      *aY = -R_theta * cos(phi);
 
  return CPL_ERROR_NONE;
} /* hdrl_resample_wcs_projplane_from_celestial() */
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
static inline void
hdrl_resample_wcs_pixel_from_celestial_fast(hdrl_resample_smallwcs *aWCS,
        double aRA, double aDEC, double *aX, double *aY)
{
  if(!aWCS) return;
  /* spherical coordinate shift/translation */
  /*Calabretta & Greisen 2002 A&A 395, 1077 (Paper II) */
  /* aRA is alpha in Paper II, aDEC is delta, aWCS->crval1 is alpha_p, *
   * aWCS->crval2 is delta_p, all of them in units of radians, eq (5), arg=atan2        */
 
  double phi = atan2(-cos(aDEC) * sin(aRA - aWCS->crval1),
             sin(aDEC) * cos(aWCS->crval2)
             - cos(aDEC) * sin(aWCS->crval2) * cos(aRA - aWCS->crval1))
                 + 180 / CPL_MATH_DEG_RAD,
             theta = asin(sin(aDEC) * sin(aWCS->crval2)
                  + cos(aDEC) * cos(aWCS->crval2) * cos(aRA - aWCS->crval1)),
                  R_theta = CPL_MATH_DEG_RAD / tan(theta);
  /* spherical deprojection */
  double x = R_theta * sin(phi),
      y = -R_theta * cos(phi);
  /* inverse linear transformation */
  *aX = (aWCS->cd22 * x - aWCS->cd12 * y) / aWCS->cddet + aWCS->crpix1;
  *aY = (aWCS->cd11 * y - aWCS->cd21 * x) / aWCS->cddet + aWCS->crpix2;
} /* hdrl_resample_wcs_pixel_from_celestial_fast() */
 
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
/* Forward declaration */
static cpl_error_code
hdrl_resample_wcs_get_scales(hdrl_resample_outgrid_parameter *aParams_outputgrid,
                 double *aXScale, double *aYScale);
 
static cpl_error_code
hdrl_resample_compute_size(hdrl_resample_outgrid_parameter *aParams_outputgrid,
                 int *aX, int *aY, int *aZ)
{
  cpl_ensure_code(aParams_outputgrid && aX && aY && aZ, CPL_ERROR_NULL_INPUT);
  const char func[] = "hdrl_resample_compute_size"; /* pretend to be in that fct... */
  double x1, y1, x2, y2;
 
  /* Fix for PIPE-12792: Correct delta_ra/delta_dec if they appear incorrect.
   * When delta values are extracted using only diagonal CD matrix elements
   * (ignoring rotation), they can be orders of magnitude too small, leading
   * to unreasonably large grid dimensions. */
  if (aParams_outputgrid->wcs) {
      double wcs_delta_ra, wcs_delta_dec;
      if (hdrl_resample_wcs_get_scales(aParams_outputgrid, &wcs_delta_ra, &wcs_delta_dec)
          == CPL_ERROR_NONE) {
          /* Only correct if current delta is < 1% of correct value (orders of magnitude difference) */
          const double BUG_THRESHOLD = 0.01;
          if (fabs(wcs_delta_ra) > 1e-15 &&
              fabs(aParams_outputgrid->delta_ra) < BUG_THRESHOLD * fabs(wcs_delta_ra)) {
              aParams_outputgrid->delta_ra = wcs_delta_ra;
          }
          if (fabs(wcs_delta_dec) > 1e-15 &&
              fabs(aParams_outputgrid->delta_dec) < BUG_THRESHOLD * fabs(wcs_delta_dec)) {
              aParams_outputgrid->delta_dec = wcs_delta_dec;
          }
      }
  }
 
  double ramin = aParams_outputgrid->ra_min;
  double ramax = aParams_outputgrid->ra_max;
  double decmin = aParams_outputgrid->dec_min;
  double decmax = aParams_outputgrid->dec_max;
  hdrl_resample_wcs_projplane_from_celestial(aParams_outputgrid, ramin, decmin,
          &x1, &y1);
  hdrl_resample_wcs_projplane_from_celestial(aParams_outputgrid, ramax, decmax,
          &x2, &y2);
  *aX = lround(fabs(x2 - x1) / aParams_outputgrid->delta_ra) + 1;
  *aY = lround(fabs(y2 - y1) / aParams_outputgrid->delta_dec) + 1;
 
  double lmin = aParams_outputgrid->lambda_min;
  double lmax = aParams_outputgrid->lambda_max;
 
  *aZ = (int)ceil((lmax - lmin) / aParams_outputgrid->delta_lambda) + 1;
 
  cpl_msg_debug(func, "Output cube size %d x %d x %d (fit to data)",
           *aX, *aY, *aZ);
  return CPL_ERROR_NONE;
} /* hdrl_resample_compute_size() */
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
static void
hdrl_resample_pixgrid_add(hdrl_resample_pixgrid *aGrid, cpl_size aIndex,
        cpl_size aRow, unsigned short aXIdx)
{
 
  if (aIndex < 0 || aGrid == NULL) {
      return;
  }
 
  if (aGrid->pix[aIndex] == 0 && aRow > 0) {
      /* First pixel is stored directly. */
      aGrid->pix[aIndex] = aRow;
  } else if (aGrid->pix[aIndex] == 0 && aRow == 0) {
      /* Special case: we cannot put "0" into the main map. */
      cpl_size iext = aGrid->nxmap[aXIdx]++;
      if (aGrid->nxmap[aXIdx] > aGrid->nxalloc[aXIdx]) {
      /* double the number of allocated entries */
      aGrid->nxalloc[aXIdx] = 2 * aGrid->nxmap[aXIdx];
      aGrid->xmaps[aXIdx] = cpl_realloc(aGrid->xmaps[aXIdx],
                        aGrid->nxalloc[aXIdx]
                               * sizeof(hdrl_resample_pixels_ext));
      }
      aGrid->xmaps[aXIdx][iext].npix = 1;
      aGrid->xmaps[aXIdx][iext].pix = cpl_malloc(sizeof(cpl_size));
      aGrid->xmaps[aXIdx][iext].pix[0] = aRow;
      aGrid->pix[aIndex] = -(iext + 1 + ((cpl_size)aXIdx << XMAP_LSHIFT));
  } else if (aGrid->pix[aIndex] > 0) {
      /* When a second pixel is added, put both to the extension map. */
      cpl_size iext = aGrid->nxmap[aXIdx]++;
      if (aGrid->nxmap[aXIdx] > aGrid->nxalloc[aXIdx]) {
      /* double the number of allocated entries */
      aGrid->nxalloc[aXIdx] = 2 * aGrid->nxmap[aXIdx];
      aGrid->xmaps[aXIdx] = cpl_realloc(aGrid->xmaps[aXIdx],
                        aGrid->nxalloc[aXIdx]
                               * sizeof(hdrl_resample_pixels_ext));
      }
      aGrid->xmaps[aXIdx][iext].npix = 2;
      aGrid->xmaps[aXIdx][iext].pix = cpl_malloc(2 * sizeof(cpl_size));
      aGrid->xmaps[aXIdx][iext].pix[0] = aGrid->pix[aIndex];
      aGrid->xmaps[aXIdx][iext].pix[1] = aRow;
      aGrid->pix[aIndex] = -(iext + 1 + ((cpl_size)aXIdx << XMAP_LSHIFT));
  } else {
      /* Append additional pixels to the extension map. */
      cpl_size iext = (-aGrid->pix[aIndex] - 1) & PT_IDX_MASK;
      /* index of the new entry in this grid point */
      unsigned int ipix = aGrid->xmaps[aXIdx][iext].npix;
      aGrid->xmaps[aXIdx][iext].npix++;
      aGrid->xmaps[aXIdx][iext].pix = cpl_realloc(aGrid->xmaps[aXIdx][iext].pix,
                          (ipix + 1) * sizeof(cpl_size));
      aGrid->xmaps[aXIdx][iext].pix[ipix] = aRow;
  }
} /* hdrl_resample_pixgrid_add() */
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
static inline cpl_size
hdrl_resample_pixgrid_get_count(hdrl_resample_pixgrid *aGrid, cpl_size aIndex)
{
  if (aIndex < 0 || aGrid == NULL) {
      return 0;
  }
  /* get entry in pixel grid */
  cpl_size p = aGrid->pix[aIndex];
  if (p == 0) { /* points nowhere --> no pixels */
      return 0;
  }
  if (p > 0) { /* points to pixel table --> 1 pixel */
      return 1;
  }
  /* p is negative, so points to an extension map, get its index */
  unsigned short ix = (-p >> XMAP_LSHIFT) & XMAP_BITMASK;
  p = (-p - 1) & PT_IDX_MASK;
  /* the npix component now gives the number of pixels in this grid index */
  return aGrid->xmaps[ix][p].npix;
} /* hdrl_resample_pixgrid_get_count() */
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
static inline cpl_size
hdrl_resample_pixgrid_get_index(hdrl_resample_pixgrid *aGrid, cpl_size aX,
        cpl_size aY, cpl_size aZ, cpl_boolean aAllowOutside)
{
  if (aGrid == NULL)  {
      return -1;
  }
  if (!aAllowOutside &&
      (aX < 0 || aX >= aGrid->nx || aY < 0 || aY >= aGrid->ny ||
      aZ < 0 || aZ >= aGrid->nz)) {
      return -1;
  }
  if (aX < 0) {
      aX = 0;
  }
  if (aX >= aGrid->nx) {
      aX = aGrid->nx - 1;
  }
  if (aY < 0) {
      aY = 0;
  }
  if (aY >= aGrid->ny) {
      aY = aGrid->ny - 1;
  }
  if (aZ < 0) {
      aZ = 0;
  }
  if (aZ >= aGrid->nz) {
      aZ = aGrid->nz - 1;
  }
  return aX + aGrid->nx * (aY + aGrid->ny * aZ);
} /* hdrl_resample_pixgrid_get_index() */
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
static hdrl_resample_pixgrid *
hdrl_resample_pixgrid_new(cpl_size aSizeX, cpl_size aSizeY, cpl_size aSizeZ,
         unsigned short aNMaps)
{
  cpl_ensure(aSizeX > 0 && aSizeY > 0 && aSizeZ > 0 && aNMaps > 0,
          CPL_ERROR_ILLEGAL_INPUT, NULL);
  hdrl_resample_pixgrid *pixels = cpl_calloc(1, sizeof(hdrl_resample_pixgrid));
  pixels->nx = aSizeX;
  pixels->ny = aSizeY;
  pixels->nz = aSizeZ;
  pixels->pix = cpl_calloc(aSizeX * aSizeY * aSizeZ, sizeof(cpl_size));
  /* extension maps for possibly multiple threads */
  pixels->nmaps = aNMaps;
  pixels->nxalloc = cpl_calloc(aNMaps, sizeof(cpl_size));
  pixels->xmaps = cpl_calloc(aNMaps, sizeof(hdrl_resample_pixels_ext *));
  pixels->nxmap = cpl_calloc(aNMaps, sizeof(cpl_size));
  return pixels;
} /* hdrl_resample_pixgrid_new() */
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
static hdrl_resample_pixgrid *
hdrl_resample_pixgrid_create(const cpl_table *ResTable, cpl_propertylist *aHeader,
            cpl_size aXSize, cpl_size aYSize, cpl_size aZSize)
{
  cpl_ensure(ResTable, CPL_ERROR_NULL_INPUT, NULL);
  cpl_ensure(aHeader, CPL_ERROR_NULL_INPUT, NULL);
  cpl_size nrow = cpl_table_get_nrow(ResTable);
  if (nrow == 0) {
      cpl_msg_error(cpl_func, "Invalid pixel table (no entries?)");
      cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT);
      return NULL;
  }
  cpl_ensure(aXSize > 0 && aYSize > 0 && aZSize > 0, CPL_ERROR_ILLEGAL_INPUT,
         NULL);
 
 
  double crval3 = hdrl_resample_pfits_get_crval(aHeader, 3),
      crpix3 = hdrl_resample_pfits_get_crpix(aHeader, 3),
      cd33 = hdrl_resample_pfits_get_cd(aHeader, 3, 3);
 
  hdrl_resample_smallwcs *wcs = hdrl_resample_smallwcs_new(aHeader);
 
  /* get all (relevant) table columns for easy pointer access */
 
  const double *xpos = cpl_table_get_data_double_const(ResTable, HDRL_RESAMPLE_TABLE_RA),
      *ypos = cpl_table_get_data_double_const(ResTable, HDRL_RESAMPLE_TABLE_DEC),
      *lbda = cpl_table_get_data_double_const(ResTable, HDRL_RESAMPLE_TABLE_LAMBDA);
  if (!xpos || !ypos || !lbda) {
      cpl_msg_error(cpl_func, "Missing pixel table column (%p %p %p): %s",
            (void *)xpos, (void *)ypos, (void *)lbda,
            cpl_error_get_message());
      cpl_error_set(cpl_func, CPL_ERROR_DATA_NOT_FOUND);
      cpl_free(wcs);
      return NULL;
  }
 
  wcs->crval1 /= CPL_MATH_DEG_RAD; /* convert to radians before calling...    */
  wcs->crval2 /= CPL_MATH_DEG_RAD; /* ...hdrl_resample_wcs_pixel_from_celestial_fast() */
 
  double timeinit = cpl_test_get_walltime(),
      timeprogress = timeinit,
      cpuinit = cpl_test_get_cputime();
  cpl_boolean showprogress = cpl_msg_get_level() == CPL_MSG_DEBUG
      || cpl_msg_get_log_level() == CPL_MSG_DEBUG;
 
  /* check for the selected pixels in the pixel table, only those *
   * are used to construct the pixel grid; since constructing the *
   * array of selected pixels costs significant amounts of time,  *
   * do that only when not all pixels are selected!               */
  cpl_array *asel = NULL;
  const cpl_size *sel = NULL;
  cpl_size nsel = cpl_table_count_selected(ResTable);
  if (nsel < nrow) {
      asel = cpl_table_where_selected(ResTable);
      sel = cpl_array_get_data_cplsize_const(asel);
  }
 
  /* can use at most XMAP_BITMASK threads so that the bitmask does not       *
   * overflow, but ensure that we are not using more cores than available... */
  int nth = omp_get_max_threads() > XMAP_BITMASK ? XMAP_BITMASK
      : omp_get_max_threads();
  /* prepare the ranges for the different threads, store them in arrays */
  cpl_array *az1 = cpl_array_new(nth, CPL_TYPE_INT),
      *az2 = cpl_array_new(nth, CPL_TYPE_INT);
  if (aZSize < nth) {
      /* pre-fill arrays with values that cause the threads to do nothing */
      cpl_array_fill_window_int(az1, aZSize, nth, -1);
      cpl_array_fill_window_int(az2, aZSize, nth, -2);
  }
  /* now fill the (first) ones with real ranges */
  double base = nth > aZSize ? 1. : (double)aZSize / nth;
  cpl_size ith;
  for (ith = 0; ith < nth && ith < aZSize; ith++) {
      cpl_array_set_int(az1, ith, lround(base * ith));
      cpl_array_set_int(az2, ith, lround(base * (ith + 1) - 1));
  } /* for */
  /* make sure that we don't lose pixels at the edges of the wavelength      *
   * range, put them into the extreme threads by making their ranges larger; *
   * set the relevant array entries to something close to the largest value, *
   * that we can still add as an integer (to compute the z-range)            */
  cpl_array_set_int(az1, 0, -INT_MAX / 2 + 1);
  cpl_array_set_int(az2, ith - 1, INT_MAX / 2 - 1);
 
  /* create the pixel grid with extension maps for threads */
  hdrl_resample_pixgrid *grid = hdrl_resample_pixgrid_new(aXSize, aYSize, aZSize, nth);
 
  /* parallel region to fill the pixel grid */
 
  struct timeval tv1, tv2;
 
  cpl_msg_debug(cpl_func,"Starting parallel loop in hdrl_resample_pixgrid_create");
  gettimeofday(&tv1, NULL);
 
#pragma omp parallel num_threads(nth) /* default(none)    as req. by Ralf */ \
    shared(aXSize, aYSize, aZSize, az1, az2, cd33, crpix3, crval3, grid, \
       lbda, nsel, sel, showprogress,     \
       timeinit, timeprogress, wcs, xpos, ypos)
 
  {
    /* split the work up into threads, for non-overlapping wavelength ranges */
    unsigned short ithread = omp_get_thread_num(); /* index of this thread */
    int z1 = cpl_array_get_int(az1, ithread, NULL),
    z2 = cpl_array_get_int(az2, ithread, NULL),
    zrange = z2 - z1 + 1;
 
    /* check if we actually need to enter the (parallel) loop, i.e. *
     * if the current thread is handling any wavelength planes      */
 
    /* now the actual parallel loop */
    cpl_size isel;
    for (isel = 0 ; zrange > 0 && isel < nsel; isel++) {
#pragma omp master /* only output progress from the master thread */
    if (showprogress && !((isel+1) % 1000000ll)) { /* output before every millionth entry */
        double timenow = cpl_test_get_walltime();
        if (timenow - timeprogress > 30.) { /* and more than half a minute passed */
        timeprogress = timenow;
        double percent = 100. * (isel + 1.) / nsel,
            elapsed = timeprogress - timeinit,
            remaining = (100. - percent) * elapsed / percent;
        /* overwritable only exists for INFO mode, but we check  *
         * above that we want this only for DEBUG mode output... */
        cpl_msg_info_overwritable(cpl_func, "pixel grid creation is %.1f%% "
                      "complete, %.1fs elapsed, ~%.1fs remaining",
                      percent, elapsed, remaining);
        } /* if: 1/2 min passed */
    } /* if: want debug output */
 
    /* either use the index from the array of selected rows   *
     * or the row numberdirectly (for a fully selected table) */
    cpl_size n = sel ? sel[isel] : isel;
    int z = -1;
 
    z = lround((lbda[n] - crval3) / cd33 + crpix3) - 1;
 
    if (z < z1 || z > z2) {
        continue; /* skip this entry, one of the other threads handles it */
    }
 
    /* determine the pixel coordinates in the grid (indices, starting at 0) */
    double xpx = 0., ypx = 0.;
 
    hdrl_resample_wcs_pixel_from_celestial_fast(wcs, (xpos[n]) / CPL_MATH_DEG_RAD,
                       (ypos[n]) / CPL_MATH_DEG_RAD,
                       &xpx, &ypx);
 
    int x = lround(xpx) - 1,
        y = lround(ypx) - 1;
    cpl_size idx = hdrl_resample_pixgrid_get_index(grid, x, y, z, CPL_TRUE);
 
 
    /* write the pixel values to the correct place in the grid */
    hdrl_resample_pixgrid_add(grid, idx, n, ithread);
    } /* for isel (all selected pixel table rows) */
 
    /* Clean up the possibly too many allocations; this is not strictly *
     * needed but nice to only consume as much memory as we need.       */
    grid->xmaps[ithread] = cpl_realloc(grid->xmaps[ithread],
                       grid->nxmap[ithread]
                           * sizeof(hdrl_resample_pixels_ext));
    grid->nxalloc[ithread] = grid->nxmap[ithread];
  } /* omp parallel */
 
  gettimeofday(&tv2, NULL);
  cpl_msg_debug(cpl_func, "Wall time for hdrl_resample_pixgrid_create was %f seconds\n",
           (double) (tv2.tv_usec - tv1.tv_usec) / 1000000 +
           (double) (tv2.tv_sec - tv1.tv_sec));
 
  cpl_array_delete(asel);
  cpl_free(wcs);
  cpl_array_delete(az1);
  cpl_array_delete(az2);
 
 
  cpl_size idx, npix = 0;
  for (idx = 0; idx < aXSize * aYSize * aZSize; idx++) {
      npix += hdrl_resample_pixgrid_get_count(grid, idx);
  }
  cpl_size nxmap = 0;
  for (idx = 0; idx < (cpl_size)grid->nmaps; idx++) {
      nxmap += grid->nxmap[idx];
  }
  if (npix != nsel) {
      char *msg = cpl_sprintf("Pixels got lost while creating the cube (input "
      "pixel table: %"CPL_SIZE_FORMAT", output pixel grid"
      ": %"CPL_SIZE_FORMAT")", nsel, npix);
      cpl_msg_error(cpl_func, "%s:", msg);
      cpl_error_set_message(cpl_func, CPL_ERROR_ILLEGAL_OUTPUT, "%s!", msg);
      cpl_free(msg);
  }
  double timefini = cpl_test_get_walltime(),
      cpufini = cpl_test_get_cputime();
  cpl_msg_debug(cpl_func, "pixel grid: %dx%dx%d, %"CPL_SIZE_FORMAT" pixels "
        "total, %"CPL_SIZE_FORMAT" (%.1f%%) in %hu extension maps; took"
        " %gs (wall-clock) and %gs (CPU) to create", (int)grid->nx,
        (int)grid->ny, (int)grid->nz, npix, nxmap,
        (double)nxmap / npix * 100., grid->nmaps, timefini - timeinit,
        cpufini - cpuinit);
 
  return grid;
} /* hdrl_resample_pixgrid_create() */
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
static cpl_error_code
hdrl_resample_wcs_get_scales(hdrl_resample_outgrid_parameter *aParams_outputgrid,
            double *aXScale, double *aYScale)
{
  cpl_ensure_code(aParams_outputgrid && aXScale && aYScale, CPL_ERROR_NULL_INPUT);
 
  cpl_errorstate prestate = cpl_errorstate_get();
 
  const cpl_matrix *cd = cpl_wcs_get_cd(aParams_outputgrid->wcs);
  //double cd11 = aParams_outputgrid->limits.cd11;
  //double cd22 = aParams_outputgrid->limits.cd22;
  //double cd12 = aParams_outputgrid->limits.cd12;
  //double cd21 = aParams_outputgrid->limits.cd21;
  double cd11 = cpl_matrix_get(cd, 0, 0);
  double cd12 = cpl_matrix_get(cd, 0, 1);
  double cd21 = cpl_matrix_get(cd, 1, 0);
  double cd22 = cpl_matrix_get(cd, 1, 1);
 
 
  double  det = cd11 * cd22 - cd12 * cd21;
  cpl_ensure_code(cpl_errorstate_is_equal(prestate), cpl_error_get_code());
 
  if (det < 0.) {
      cd12 *= -1;
      cd11 *= -1;
  }
  if (cd12 == 0. && cd21 == 0.) { /* matrix without rotation */
      *aXScale = cd11;
      *aYScale = cd22;
      return CPL_ERROR_NONE;
  }
  *aXScale = sqrt(cd11*cd11 + cd12*cd12); /* only the absolute value */
  *aYScale = sqrt(cd22*cd22 + cd21*cd21);
  return CPL_ERROR_NONE;
} /* hdrl_resample_wcs_get_scales() */
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
static inline const cpl_size *
hdrl_resample_pixgrid_get_rows(hdrl_resample_pixgrid *aGrid, cpl_size aIndex)
{
 
  cpl_ensure(aGrid , CPL_ERROR_NULL_INPUT, 0);
  cpl_ensure(aIndex >= 0 , CPL_ERROR_ILLEGAL_INPUT, 0);
  cpl_ensure(aIndex < aGrid->nx * aGrid->ny * aGrid->nz, CPL_ERROR_ILLEGAL_INPUT, 0);
  /* get entry in pixel grid */
  cpl_size p = aGrid->pix[aIndex];
  if (p == 0) { /* points nowhere --> no array */
      return NULL;
  }
  if (p > 0) { /* points to pixel table */
      return aGrid->pix + aIndex;
  }
  /* p is negative, so points to an extension map, get its array */
  unsigned short ix = (-p >> XMAP_LSHIFT) & XMAP_BITMASK;
  p = (-p - 1) & PT_IDX_MASK;
  /* the pix component provides the array of pixel table rows in this index */
  return aGrid->xmaps[ix][p].pix;
} /* hdrl_resample_pixgrid_get_rows() */
 
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
static cpl_error_code
hdrl_resample_cube_nearest(hdrl_resample_result *aCube, const cpl_table *ResTable,
                 hdrl_resample_pixgrid *aGrid,
                 hdrl_resample_outgrid_parameter *aParams_outputgrid)
{
  cpl_ensure_code(aCube && ResTable && aGrid && aParams_outputgrid,
          CPL_ERROR_NULL_INPUT);
  cpl_ensure_code(cpl_propertylist_has(aCube->header,"CRVAL3") == 1,
          CPL_ERROR_ILLEGAL_INPUT);
  cpl_ensure_code(cpl_propertylist_has(aCube->header,"CRPIX3") == 1,
          CPL_ERROR_ILLEGAL_INPUT);
  cpl_ensure_code(cpl_propertylist_has(aCube->header,"CD3_3") == 1,
          CPL_ERROR_ILLEGAL_INPUT);
 
  double crval3 = hdrl_resample_pfits_get_crval(aCube->header, 3);
  double crpix3 = hdrl_resample_pfits_get_crpix(aCube->header, 3);
  double cd33 = hdrl_resample_pfits_get_cd(aCube->header, 3, 3);
  cpl_ensure_code(cd33 != 0, CPL_ERROR_ILLEGAL_INPUT);
 
  cpl_wcs *wcscpl = cpl_wcs_new_from_propertylist(aCube->header);
 
  if (hdrl_resample_inputtable_verify(ResTable)) {
       return CPL_ERROR_ILLEGAL_INPUT;
  }
 
 
  const double *xpos = cpl_table_get_data_double_const(ResTable, HDRL_RESAMPLE_TABLE_RA),
      *ypos = cpl_table_get_data_double_const(ResTable, HDRL_RESAMPLE_TABLE_DEC),
      *lbda = cpl_table_get_data_double_const(ResTable, HDRL_RESAMPLE_TABLE_LAMBDA),
      *data = cpl_table_get_data_double_const(ResTable, HDRL_RESAMPLE_TABLE_DATA),
      *stat = cpl_table_get_data_double_const(ResTable, HDRL_RESAMPLE_TABLE_ERRORS);
  const int *dq = cpl_table_get_data_int_const(ResTable, HDRL_RESAMPLE_TABLE_BPM);
 
  /* If our data was astrometrically calibrated, we need to scale the *
   * data units to the pixel size in all three dimensions so that the *
   * radius computation works again.                                  *
   * Otherwise dx~5.6e-5deg won't contribute to the weighting at all. */
 
  double xnorm = 1., ynorm = 1., znorm = 1. ;
 
 
  hdrl_resample_wcs_get_scales(aParams_outputgrid, &xnorm, &ynorm);
  //TODO: we should check that xnorm, ynorm, znorm are not zero
  xnorm = 1. / xnorm;
  ynorm = 1. / ynorm;
  //znorm = 1. / aParams_outputgrid->limits.cd33;
  const cpl_matrix *cd = cpl_wcs_get_cd(aParams_outputgrid->wcs);
  if (cpl_matrix_get_ncol(cd) == 3) {
      znorm = 1. / cpl_matrix_get(cd, 2, 2);
  }
/*
  cpl_msg_debug(cpl_func, "Norm factors from wcs - xnorm: %g ynorm: %g znorm: %g",
        xnorm, ynorm, znorm);
   need to use the real coordinate offset for celestial spherical
  // We are now working with the full astrometric solution
  //ptxoff = hdrl_resample_pfits_get_crval(ResTable->header, 1);
  //ptyoff = hdrl_resample_pfits_get_crval(ResTable->header, 2);
*/
 
  struct timeval tv1, tv2;
  cpl_msg_debug(cpl_func,"Starting parallel loop in hdrl_resample_cube_nearest");
  gettimeofday(&tv1, NULL);
 
#pragma omp parallel for default(none)                 /* as req. by Ralf */ \
    shared(aCube, aGrid, cd33, crpix3, crval3, data, dq, lbda,           \
       stat, stdout, wcscpl, \
       xnorm, xpos, ynorm, ypos, znorm) collapse(2)
 
  for (cpl_size l = 0; l < aGrid->nz; l++) {
      for (cpl_size i = 0; i < aGrid->nx; i++) {
      double *pdata   = cpl_image_get_data_double(hdrl_image_get_image(hdrl_imagelist_get(aCube->himlist, l)));
      double *pstat   = cpl_image_get_data_double(hdrl_image_get_error(hdrl_imagelist_get(aCube->himlist, l)));
      cpl_binary *pdq = cpl_mask_get_data(hdrl_image_get_mask(hdrl_imagelist_get(aCube->himlist, l)));
      /* wavelength of center of current grid cell (l is index starting at 0) */
      double lambda = (l + 1. - crpix3) * cd33 + crval3;
 
      cpl_size j;
      for (j = 0; j < aGrid->ny; j++) {
          cpl_size idx = hdrl_resample_pixgrid_get_index(aGrid, i, j, l, CPL_FALSE),
          n_rows = hdrl_resample_pixgrid_get_count(aGrid, idx);
          const cpl_size *rows = hdrl_resample_pixgrid_get_rows(aGrid, idx);
 
          /* x and y position of center of current grid cell (i, j start at 0) */
          double x = 0., y = 0.;
 
          // We are now working with the full astrometric solution
          // hdrl_resample_wcs_celestial_from_pixel_fast(wcs, i + 1, j + 1, &x, &y);
          hdrl_wcs_xy_to_radec(wcscpl, i + 1, j + 1, &x, &y);
 
          if (n_rows == 1) {
          if ((cpl_binary)dq[rows[0]] == CPL_BINARY_0) {
              /* if there is only one pixel in the cell, just use it */
              pdata[i + j * aGrid->nx] = data[rows[0]];
              pstat[i + j * aGrid->nx] = stat[rows[0]];
              pdq[i + j * aGrid->nx] = (cpl_binary)dq[rows[0]];
          } else {
              pdq[i + j * aGrid->nx] = CPL_BINARY_1;
          }
          } else if (n_rows >= 2) {
          /* loop through all available values and take the closest one */
          cpl_size n, nbest = -1;
          double dbest = FLT_MAX; /* some unlikely large value to start with*/
          for (n = 0; n < n_rows; n++) {
              /* do not use bad pixels */
              if((cpl_binary)dq[rows[n]] != CPL_BINARY_0) {
              continue;
              }
              /* the differences for the cell center and the current pixel */
              double dx = fabs(x - xpos[rows[n]]) * xnorm,
              dy = fabs(y - ypos[rows[n]]) * ynorm,
              dlambda = fabs(lambda - lbda[rows[n]]) * znorm,
              dthis = sqrt(dx*dx + dy*dy + dlambda*dlambda);
 
              /* Not strictly necessary for NN, but still scale the RA   *
               * distance properly, see hdrl_resample_cube_weighted(). */
              dx *= cos(y * CPL_MATH_RAD_DEG);
 
              if (dthis < dbest) {
              nbest = n;
              dbest = dthis;
              }
          }
          if (nbest >= 0) { /* We found a good nearest neighbor */
              pdata[i + j * aGrid->nx] = data[rows[nbest]];
              pstat[i + j * aGrid->nx] = stat[rows[nbest]];
              pdq[i + j * aGrid->nx] = (cpl_binary)dq[rows[nbest]];
          }
          } else {
          /* npix == 0: do nothing, pixel stays zero */
          pdq[i + j * aGrid->nx] = CPL_BINARY_1;
          }
      } /* for j (y direction) */
      } /* for i (x direction) */
  } /* for l (wavelength planes) */
 
  gettimeofday(&tv2, NULL);
  cpl_msg_debug(cpl_func, "Wall time for hdrl_resample_cube_nearest was %f seconds\n",
           (double) (tv2.tv_usec - tv1.tv_usec) / 1000000 +
           (double) (tv2.tv_sec - tv1.tv_sec));
 
  /* Make sure that the bpm of the image and the error are in sinc as we are
   * working with pointers */
  cpl_size size = hdrl_imagelist_get_size(aCube->himlist);
  for(cpl_size i = 0; i < size; i++){
      /* sync image and error bpm ignoring what is in error before */
      cpl_image_reject_from_mask(hdrl_image_get_error(hdrl_imagelist_get(aCube->himlist, i)),
                 hdrl_image_get_mask(hdrl_imagelist_get(aCube->himlist, i) ));
  }
 
  cpl_wcs_delete(wcscpl);
 
  return CPL_ERROR_NONE;
} /* hdrl_resample_cube_nearest() */
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
static inline double
hdrl_resample_weight_function_renka(double r, double r_c)
{
  if (r == 0) {
      return FLT_MAX;
  } else if (r >= r_c) {
      return DBL_MIN;
  } else {
      double p = (r_c - r) / (r_c * r);
      return p*p;
  }
} /* hdrl_resample_weight_function_renka() */
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
static inline double
hdrl_resample_weight_function_drizzle(double xin, double yin, double zin,
                    double xout, double yout, double zout,
                    double dx, double dy, double dz)
{
  /* compute the three terms in the numerator: if the offset   *
   * plus the output halfsize is less than the input halfsize, *
   * then that side is fully contained in the input pixel      */
  double x = (dx + xout / 2.) <= xin / 2. ? xout : (xin + xout) / 2. - dx,
      y = (dy + yout / 2.) <= yin / 2. ? yout : (yin + yout) / 2. - dy,
      z = (dz + zout / 2.) <= zin / 2. ? zout : (zin + zout) / 2. - dz;
  /* any negative value means that the input pixel is completely *
   * outside the target voxel, so it doesn't contribute          */
  if (x <= 0 || y <= 0 || z <= 0) {
      return 0.;
  }
  /* any value > input size means this dimension of the input pixel *
   * is completely inside the target voxel, so that's the limit!    */
  return (x > xin ? xin : x) * (y > yin ? yin : y) * (z > zin ? zin : z)
      / (xin * yin * zin);
} /* hdrl_resample_weight_function_drizzle() */
 
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
static inline double
hdrl_resample_weight_function_linear(double r)
{
  return r == 0 ? FLT_MAX : 1. / r;
} /* hdrl_resample_weight_function_linear() */
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
static inline double
hdrl_resample_weight_function_quadratic(double r2)
{
  return r2 == 0 ? FLT_MAX : 1. / r2;
} /* hdrl_resample_weight_function_quadratic() */
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
static inline double
hdrl_resample_weight_function_sinc(double r)
{
  return fabs(r) < DBL_EPSILON ? 1. : sin(CPL_MATH_PI * r) / (CPL_MATH_PI * r);
} /* hdrl_resample_weight_function_sinc() */
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
static inline double
hdrl_resample_weight_function_lanczos(double dx, double dy, double dz,
                    unsigned int ld, unsigned int lks)
{
  /* Adding 0.5 as for a loop distance of 0 the weight should only drop to 0 if
   * the distance is larger then half the pixel */
  return (fabs(dx) >= (ld + 0.5) || fabs(dy) >= (ld + 0.5) || fabs(dz) > (ld + 0.5)) ? 0.
      : hdrl_resample_weight_function_sinc(dx) * hdrl_resample_weight_function_sinc(dx / lks)
    * hdrl_resample_weight_function_sinc(dy) * hdrl_resample_weight_function_sinc(dy / lks)
    * hdrl_resample_weight_function_sinc(dz) * hdrl_resample_weight_function_sinc(dz / lks);
} /* hdrl_resample_weight_function_lanczos() */
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
static cpl_error_code
hdrl_resample_cube_weighted(hdrl_resample_result *aCube,
        const cpl_table *ResTable,
        hdrl_resample_pixgrid *aGrid,
        hdrl_resample_method_parameter *aParams_method,
        hdrl_resample_outgrid_parameter *aParams_outputgrid)
{
 
  cpl_ensure_code(aCube && ResTable && aGrid && aParams_method &&
          aParams_outputgrid, CPL_ERROR_NULL_INPUT);
  cpl_ensure_code(cpl_propertylist_has(aCube->header,"CRVAL3") == 1,
          CPL_ERROR_ILLEGAL_INPUT);
  cpl_ensure_code(cpl_propertylist_has(aCube->header,"CRPIX3") == 1,
          CPL_ERROR_ILLEGAL_INPUT);
  cpl_ensure_code(cpl_propertylist_has(aCube->header,"CD3_3") == 1,
          CPL_ERROR_ILLEGAL_INPUT);
 
  double crval3 = hdrl_resample_pfits_get_crval(aCube->header, 3),
      crpix3 = hdrl_resample_pfits_get_crpix(aCube->header, 3),
      cd33 = hdrl_resample_pfits_get_cd(aCube->header, 3, 3);
 
  hdrl_resample_smallwcs *wcs = hdrl_resample_smallwcs_new(aCube->header);
  cpl_wcs *wcscpl = cpl_wcs_new_from_propertylist(aCube->header);
 
  const double *xpos = cpl_table_get_data_double_const(ResTable, HDRL_RESAMPLE_TABLE_RA),
      *ypos = cpl_table_get_data_double_const(ResTable, HDRL_RESAMPLE_TABLE_DEC),
      *lbda = cpl_table_get_data_double_const(ResTable, HDRL_RESAMPLE_TABLE_LAMBDA),
      *data = cpl_table_get_data_double_const(ResTable, HDRL_RESAMPLE_TABLE_DATA),
      *stat = cpl_table_get_data_double_const(ResTable, HDRL_RESAMPLE_TABLE_ERRORS);
  const int *dq = cpl_table_get_data_int_const(ResTable, HDRL_RESAMPLE_TABLE_BPM);
 
  /* If our data was astrometrically calibrated, we need to scale the *
   * data units to the pixel size in all three dimensions so that the *
   * radius computation works again.                                  *
   * Otherwise dx~5.6e-5deg won't contribute to the weighting at all. */
 
  double xnorm = 1., ynorm = 1., znorm = 1. ;
 
  hdrl_resample_wcs_get_scales(aParams_outputgrid, &xnorm, &ynorm);
  xnorm = 1. / xnorm;
  ynorm = 1. / ynorm;
  //znorm = 1. / aParams_outputgrid->limits.cd33;
  const cpl_matrix *cd = cpl_wcs_get_cd(aParams_outputgrid->wcs);
  if (cpl_matrix_get_ncol(cd) == 3) {
      znorm = 1. / cpl_matrix_get(cd, 2, 2);
  }
/*
  cpl_msg_debug(cpl_func, "Norm factors from wcs - xnorm: %g ynorm: %g znorm: %g",
        xnorm, ynorm, znorm);
   need to use the real coordinate offset for celestial spherical
  // We are now working with the full astrometric solution
  //ptxoff = hdrl_resample_pfits_get_crval(ResTable->header, 1);
  //ptyoff = hdrl_resample_pfits_get_crval(ResTable->header, 2);
*/
 
 
  /* scale the input critical radius by the voxel radius */
  double renka_rc = aParams_method->renka_critical_radius /*beware of rotation! */
      * sqrt((wcs->cd11*xnorm)*(wcs->cd11*xnorm)
         + (wcs->cd22*ynorm)*(wcs->cd22*ynorm)
         + (cd33*znorm)*(cd33*znorm));
 
  /* loop distance (to take into account surrounding pixels) verification */
  int ld = aParams_method->loop_distance;
  if (ld < 0) {
      ld = 0;
      cpl_msg_debug(cpl_func, "Overriding loop distance ld=%d", ld);
  }
 
  /* Lankzos kernel size (lks) verification */
  int lks = aParams_method->lanczos_kernel_size;
  if (lks <= 0) {
      lks = 1;
      cpl_msg_debug(cpl_func, "Overriding lanczos kernel size lks=%d", lks);
  }
 
  /* Should 1/variance be used as an additional weight */
  cpl_boolean wght = aParams_method->use_errorweights;
 
  /* pixel sizes in all three directions, scaled by pixfrac, and *
   * output pixel sizes (absolute values), as needed for drizzle */
  double xsz = aParams_method->drizzle_pix_frac_x / xnorm,
      ysz = aParams_method->drizzle_pix_frac_y / ynorm,
      zsz = aParams_method->drizzle_pix_frac_lambda / znorm,
      xout = fabs(wcs->cd11), yout = fabs(wcs->cd22), zout = fabs(cd33);
 
  struct timeval tv1, tv2;
  cpl_msg_debug(cpl_func,"Starting parallel loop in hdrl_resample_cube_weighted");
  gettimeofday(&tv1, NULL);
 
#pragma omp parallel for default(none)                 /* as req. by Ralf */ \
    shared(aCube, aParams_method, aParams_outputgrid, aGrid, ResTable, cd33, crpix3, crval3, data, \
       dq, lbda, ld, lks, renka_rc, stat,      \
       stdout, wcs, wcscpl, wght, xnorm, xout, xpos, xsz, ynorm, yout,\
       ypos, ysz, znorm, zout, zsz) collapse(2)
 
  for (cpl_size l = 0; l < aGrid->nz; l++) {
      for (cpl_size i = 0; i < aGrid->nx; i++) {
      double *pdata   = cpl_image_get_data_double(hdrl_image_get_image(hdrl_imagelist_get(aCube->himlist, l)));
      double *pstat   = cpl_image_get_data_double(hdrl_image_get_error(hdrl_imagelist_get(aCube->himlist, l)));
      cpl_binary *pdq = cpl_mask_get_data(hdrl_image_get_mask(hdrl_imagelist_get(aCube->himlist, l)));
 
      /* wavelength of center of current grid cell (l is index starting at 0) */
      double lambda = (l + 1. - crpix3) * cd33 + crval3;
      double zout2 = zout; /* correct the output pixel size for log-lambda */
 
      for (cpl_size j = 0; j < aGrid->ny; j++) {
          /* x and y position of center of current grid cell (i, j start at 0) */
          double x, y;
 
          // We are now working with the full astrometric solution
          // hdrl_resample_wcs_celestial_from_pixel_fast(wcs, i + 1, j + 1, &x, &y);
          hdrl_wcs_xy_to_radec(wcscpl, i + 1, j + 1, &x, &y);
          //cpl_msg_debug(cpl_func,"case1 i %d j %d x %10.7f y %10.7f ",i+1,j+1,x,y);
 
          double sumdata = 0, sumstat = 0, sumweight = 0;
          cpl_size npoints = 0;
 
          /* loop through surrounding cells and their contained pixels */
          cpl_size i2;
          for (i2 = i - ld; i2 <= i + ld; i2++) {
          cpl_size j2;
          for (j2 = j - ld; j2 <= j + ld; j2++) {
              cpl_size l2;
              for (l2 = l - ld; l2 <= l + ld; l2++) {
              cpl_size idx2 = hdrl_resample_pixgrid_get_index(aGrid, i2, j2, l2, CPL_FALSE);
              if (idx2 < 0) continue;
              cpl_size n_rows2 = hdrl_resample_pixgrid_get_count(aGrid, idx2);
 
              const cpl_size *rows2 = hdrl_resample_pixgrid_get_rows(aGrid, idx2);
              cpl_size n;
              for (n = 0; n < n_rows2; n++) {
                  if (dq[rows2[n]]) { /* exclude all bad pixels */
 
                  continue;
                  }
 
                  double dx = fabs(x - (xpos[rows2[n]])),
                  dy = fabs(y - (ypos[rows2[n]])),
                  dlambda = fabs(lambda - lbda[rows2[n]]),
                  r2 = 0;
 
                  /* Since the distances of RA in degrees get larger the *
                   * closer we get to the celestial pole, we have to     *
                   * compensate for that by multiplying the distance in  *
                   * RA by cos(delta), to make it comparable to the      *
                   * distances in pixels for the different kernels below. */
 
                  // We are now working with the full astrometric solution
                  dx *= cos(y * CPL_MATH_RAD_DEG);
 
                  if (aParams_method->method != HDRL_RESAMPLE_METHOD_DRIZZLE) {
                  dx *= xnorm;
                  dy *= ynorm;
                  dlambda *= znorm;
                  r2 = dx * dx + dy * dy + dlambda * dlambda;
                  //cpl_msg_debug(cpl_func,"n %lld dx %10.7f dy %10.7f dlambda %10.7f",n,dx,dy,dlambda);
                  }
                  double weight = 0.;
                  if (aParams_method->method == HDRL_RESAMPLE_METHOD_RENKA) {
                  weight = hdrl_resample_weight_function_renka(sqrt(r2), renka_rc);
                  } else if (aParams_method->method == HDRL_RESAMPLE_METHOD_DRIZZLE) {
                  weight = hdrl_resample_weight_function_drizzle(xsz, ysz, zsz,
                                           xout, yout, zout2,
                                           dx, dy, dlambda);
                  } else if (aParams_method->method == HDRL_RESAMPLE_METHOD_LINEAR) {
 
                  weight = hdrl_resample_weight_function_linear(sqrt(r2));
                  //cpl_msg_debug(cpl_func,"r2=%10.7f weight=%10.7f",r2,weight);
                  } else if (aParams_method->method == HDRL_RESAMPLE_METHOD_QUADRATIC) {
                  weight = hdrl_resample_weight_function_quadratic(r2);
                  } else if (aParams_method->method == HDRL_RESAMPLE_METHOD_LANCZOS) {
                  weight = hdrl_resample_weight_function_lanczos(dx, dy,
                                           dlambda, ld, lks);
                  }
 
                  if (wght && stat[rows2[n]] > 0.) { /* User wants to weight by 1/variance */
                  /* apply it on top of the weight computed here */
                  weight /= (stat[rows2[n]] * stat[rows2[n]]);
                  }
 
                  sumweight += weight;
                  sumdata += data[rows2[n]] * weight;
                  sumstat += stat[rows2[n]] * stat[rows2[n]] * weight * weight;
                  npoints++;
 
              } /* for n (all pixels in grid cell) */
              } /* for l2 (lambda direction) */
          } /* for j2 (y direction) */
          } /* for i2 (x direction) */
 
 
          /* if no points were found, we cannot divide by the summed weight *
           * and don't need to set the output pixel value (it's 0 already), *
           * only set the relevant Euro3D bad pixel flag
           * In some cases only sumweight * sumweight is really zero so this
           * check was additionally added for the error propagation part*/
 
          if (!npoints || !isnormal(sumweight) || !isnormal(sumweight * sumweight)) {
          pdq[i + j * aGrid->nx] = CPL_BINARY_1;
          continue;
          }
 
          /* divide results by weight of summed pixels */
          sumdata /= sumweight;
          sumstat /= sumweight * sumweight;
 
          pdata[i + j * aGrid->nx] = sumdata;
          /* Going back from variance to errors */
          pstat[i + j * aGrid->nx] = sqrt(sumstat);
          pdq[i + j * aGrid->nx] = CPL_BINARY_0; /* now we can mark it as good */
 
      } /* for j (y direction) */
      } /* for i (x direction) */
  } /* for l (wavelength planes) */
 
  gettimeofday(&tv2, NULL);
  cpl_msg_debug(cpl_func, "Wall time for hdrl_resample_cube_weighted was %f seconds\n",
           (double) (tv2.tv_usec - tv1.tv_usec) / 1000000 +
           (double) (tv2.tv_sec - tv1.tv_sec));
 
  /* Make sure that the bpm of the image and the error are in sinc as we are
   * working with pointers */
  cpl_size size = hdrl_imagelist_get_size(aCube->himlist);
  for(cpl_size i = 0; i < size; i++){
      /* sync image and error bpm ignoring what is in error before */
      cpl_image_reject_from_mask(hdrl_image_get_error(hdrl_imagelist_get(aCube->himlist, i)),
                 hdrl_image_get_mask(hdrl_imagelist_get(aCube->himlist, i) ));
  }
 
  cpl_free(wcs);
  cpl_wcs_delete(wcscpl);
 
  return CPL_ERROR_NONE;
} /* hdrl_resample_cube_weighted() */
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
static cpl_error_code hdrl_resampling_set_outputgrid (
    int xsize, int ysize, int zsize, hdrl_resample_result *cube,
    hdrl_resample_outgrid_parameter *aParams_outputgrid)
{
  cpl_ensure_code(xsize > 0, CPL_ERROR_ILLEGAL_INPUT);
  cpl_ensure_code(ysize > 0,CPL_ERROR_ILLEGAL_INPUT);
  cpl_ensure_code(zsize >= 0, CPL_ERROR_ILLEGAL_INPUT);
 
  cpl_ensure_code(cube && aParams_outputgrid, CPL_ERROR_NULL_INPUT);
  /* TODO: cube->header may be already allocated. In this case why a new and
     * not an over-write of the content?
     */
  cube->header = cpl_propertylist_new ();
  hdrl_wcs_to_propertylist (aParams_outputgrid->wcs, cube->header,
                FALSE);
 
  cpl_propertylist_update_string (cube->header, "CTYPE1", "RA---TAN");
  cpl_propertylist_update_string (cube->header, "CTYPE2", "DEC--TAN");
  cpl_propertylist_set_comment(cube->header, "CTYPE1", "Gnomonic projection");
  cpl_propertylist_set_comment(cube->header, "CTYPE2", "Gnomonic projection");
 
 
 
  /* set NAXIS for later handling of the WCS */
  cpl_propertylist_update_int (cube->header, "NAXIS", 3);
  cpl_propertylist_update_int (cube->header, "NAXIS1", xsize);
  cpl_propertylist_update_int (cube->header, "NAXIS2", ysize);
  cpl_propertylist_update_int (cube->header, "NAXIS3", zsize);
  /* if pixel table was astrometrically calibrated, use its WCS headers *
   * Axis 1: x or RA, axis 2: y or DEC, axis 3: lambda
   *
   *                 */
  cpl_propertylist_update_double (cube->header, "CD1_1",
                  -aParams_outputgrid->delta_ra);
  cpl_propertylist_update_double (cube->header, "CD2_2",
                  aParams_outputgrid->delta_dec);
  cpl_propertylist_update_double (cube->header, "CD1_2", 0.);
  cpl_propertylist_update_double (cube->header, "CD2_1", 0.);
 
  double ramin = aParams_outputgrid->ra_min;
  double ramax = aParams_outputgrid->ra_max;
  double decmin = aParams_outputgrid->dec_min;
  double decmax = aParams_outputgrid->dec_max;
 
  /* Following swarp we put CRPIX and CRVAL to the center of the field */
  cpl_propertylist_update_double(cube->header, "CRPIX1", (double)((xsize + 1) / 2.));
  cpl_propertylist_update_double(cube->header, "CRPIX2", (double)((ysize + 1) / 2.));
 
  if(ramax - ramin < 180) {
      /* To be checked: Both values are in 0 - 180 or 180 - 360 */
      cpl_propertylist_update_double(cube->header, "CRVAL1", (ramin + ramax) / 2.);
  } else {
      double diff1 = 360. - ramax;
      double diff2 = ramin - 0.;
      if (diff1 < diff2) {
      cpl_propertylist_update_double(cube->header, "CRVAL1", ramin - (diff1 + diff2) / 2);
      } else {
      cpl_propertylist_update_double(cube->header, "CRVAL1", ramax + (diff1 + diff2) / 2.);
      }
  }
  cpl_propertylist_update_double(cube->header, "CRVAL2", (decmin + decmax) / 2.);
  cpl_propertylist_update_double (cube->header, "CD3_3",
                  aParams_outputgrid->delta_lambda);
  cpl_propertylist_update_double (cube->header, "CRPIX3", 1.);
  cpl_propertylist_update_double (cube->header, "CRVAL3",
                  aParams_outputgrid->lambda_min);
  /* fill in empty cross-terms of the CDi_j matrix */
  cpl_propertylist_update_double (cube->header, "CD1_3", 0.);
  cpl_propertylist_update_double (cube->header, "CD2_3", 0.);
  cpl_propertylist_update_double (cube->header, "CD3_1", 0.);
  cpl_propertylist_update_double (cube->header, "CD3_2", 0.);
 
  return cpl_error_get_code();
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
static hdrl_resample_result *
hdrl_resample_cube(const cpl_table *ResTable,
             hdrl_resample_method_parameter *aParams_method,
             hdrl_resample_outgrid_parameter *aParams_outputgrid,
             hdrl_resample_pixgrid **aGrid)
{
  cpl_ensure(ResTable && aParams_method && aParams_outputgrid && aGrid,
          CPL_ERROR_NULL_INPUT, NULL);
 
  /* compute or set the size of the output grid depending on *
   * the inputs and the data available in the pixel table    */
 
 
  /* compute output sizes; wavelength is different in that it is *
   * more useful to contain partly empty areas within the field  *
   * for the extreme ends, so use ceil()                         */
  int xsize = 0, ysize = 0, zsize = 0;
 
  hdrl_resample_compute_size(aParams_outputgrid, &xsize, &ysize, &zsize);
 
  /* Following swarp for x and y : Add a margin in field size */
 
  xsize *= (100. + aParams_outputgrid->fieldmargin)/100.;
  ysize *= (100. + aParams_outputgrid->fieldmargin)/100.;
 
  cpl_ensure(xsize > 0 && ysize > 0 && zsize > 0, CPL_ERROR_ILLEGAL_OUTPUT,
         NULL);
 
  double time = cpl_test_get_walltime();
 
  /* create the structure for the output datacube */
  hdrl_resample_result *cube = cpl_calloc(1, sizeof(hdrl_resample_result));
 
  hdrl_resampling_set_outputgrid (xsize, ysize, zsize, cube, aParams_outputgrid);
 
  if (aParams_method->method < HDRL_RESAMPLE_METHOD_NONE) {
      /* fill the cube for the data */
      cube->himlist = hdrl_imagelist_new();
 
      for (cpl_size i = 0; i < zsize; i++) {
      hdrl_image * image = hdrl_image_new(xsize, ysize);
 
      /* Set all pixels a priori to bad - do not use pointers to keep the
       * bpm of the data and error image in sync*/
      for (cpl_size j = 1; j <= xsize; j++) {
          for (cpl_size k = 1; k <= ysize; k++) {
          hdrl_image_reject(image, j, k);
          }
      }
 
          hdrl_imagelist_set(cube->himlist, image, i);
      } /* for i (all wavelength planes) */
  } /* if method not HDRL_RESAMPLE_METHOD_NONE */
 
  /* convert the pixel table into a pixel grid */
  hdrl_resample_pixgrid *grid = hdrl_resample_pixgrid_create(ResTable, cube->header,
                       xsize, ysize, zsize);
  if (!grid) {
      hdrl_resample_result_delete(cube);
      cpl_error_set_message(cpl_func, CPL_ERROR_DATA_NOT_FOUND, "Could not create"
                " pixel grid!");
      if (aGrid) {
      *aGrid = NULL;
      }
      return NULL;
  } /* if not grid */
 
  double timeinit = cpl_test_get_walltime(),
      cpuinit = cpl_test_get_cputime();
 
  /* do the resampling */
  cpl_error_code rc = CPL_ERROR_NONE;
  switch (aParams_method->method) {
    case HDRL_RESAMPLE_METHOD_NEAREST:
      cpl_msg_debug(cpl_func, "Starting resampling, using method \"nearest\"");
      rc = hdrl_resample_cube_nearest(cube, ResTable, grid, aParams_outputgrid);
      break;
    case HDRL_RESAMPLE_METHOD_RENKA:
      cpl_msg_debug(cpl_func, "Starting resampling, using method \"renka\" "
           "(critical radius rc=%f, loop distance ld=%d)",
           aParams_method->renka_critical_radius,
           aParams_method->loop_distance);
      rc = hdrl_resample_cube_weighted(cube, ResTable, grid, aParams_method,
                     aParams_outputgrid);
      break;
    case HDRL_RESAMPLE_METHOD_LINEAR:
    case HDRL_RESAMPLE_METHOD_QUADRATIC:
    case HDRL_RESAMPLE_METHOD_LANCZOS:
      cpl_msg_debug(cpl_func, "Starting resampling, using method \"%s\" (loop "
           "distance ld=%d)",
           aParams_method->method == HDRL_RESAMPLE_METHOD_LINEAR
           ? "linear"
               : (aParams_method->method == HDRL_RESAMPLE_METHOD_QUADRATIC
               ? "quadratic"
                   : "lanczos"),
                 aParams_method->loop_distance);
      rc = hdrl_resample_cube_weighted(cube, ResTable, grid, aParams_method,
                     aParams_outputgrid);
      break;
    case HDRL_RESAMPLE_METHOD_DRIZZLE:
      cpl_msg_debug(cpl_func, "Starting resampling, using method \"drizzle\" "
           "(pixfrac f=%.3f,%.3f,%.3f, loop distance ld=%d)",
           aParams_method->drizzle_pix_frac_x,
           aParams_method->drizzle_pix_frac_y,
           aParams_method->drizzle_pix_frac_lambda,
           aParams_method->loop_distance);
      rc = hdrl_resample_cube_weighted(cube, ResTable, grid,
                     aParams_method,
                     aParams_outputgrid);
      break;
    case HDRL_RESAMPLE_METHOD_NONE:
      /* cpl_msg_debug(cpl_func, "Method %d (no resampling)", aParams_method->method); */
      break;
    default:
      cpl_msg_error(cpl_func, "Don't know this resampling method: %d",
            aParams_method->method);
      cpl_error_set(cpl_func, CPL_ERROR_UNSUPPORTED_MODE);
      rc = CPL_ERROR_UNSUPPORTED_MODE;
  }
 
  double timefini = cpl_test_get_walltime(),
      cpufini = cpl_test_get_cputime();
 
  /* now that we have resampled we can either remove the pixel grid or save it */
  if (aGrid) {
      *aGrid = grid;
  } else {
      hdrl_resample_pixgrid_delete(grid);
  }
 
  cpl_msg_debug(cpl_func, "resampling took %.3fs (wall-clock) and %.3fs "
        "(%.3fs CPU, %d CPUs) for hdrl_resample_cube*() alone",
        timefini - time, timefini - timeinit, cpufini - cpuinit,
        omp_get_max_threads());
  if (rc != CPL_ERROR_NONE) {
      cpl_msg_error(cpl_func, "resampling failed: %s", cpl_error_get_message());
      hdrl_resample_result_delete(cube);
      return NULL;
  }
 
  return cube;
} /* hdrl_resample_cube() */
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
cpl_error_code
hdrl_wcs_to_propertylist(const cpl_wcs * wcs, cpl_propertylist * header,
             cpl_boolean only2d)
{
  cpl_ensure_code(wcs && header, CPL_ERROR_NULL_INPUT);
  int err = 0;
  const cpl_array *crval = cpl_wcs_get_crval(wcs);
  const cpl_array *crpix = cpl_wcs_get_crpix(wcs);
  const cpl_array *ctype = cpl_wcs_get_ctype(wcs);
  const cpl_array *cunit = cpl_wcs_get_cunit(wcs);
 
  const cpl_matrix *cd = cpl_wcs_get_cd(wcs);
 
  const cpl_array *dims = cpl_wcs_get_image_dims(wcs);
  int naxis = cpl_wcs_get_image_naxis(wcs);
 
 
  /* Check NAXIS */
  for (cpl_size i = 0; i < naxis; i++) {
      if (i == 0) {
      cpl_propertylist_update_int(header, "NAXIS", naxis);
      }
      char * buf = cpl_sprintf("NAXIS%lld", i + 1);
      cpl_propertylist_update_int(header, buf,
                  cpl_array_get_int(dims, i, &err));
      cpl_free(buf);
  }
 
/* Make sure to have the right NAXIS keywords if 2D is forced */
  if (only2d == TRUE) {
      cpl_propertylist_update_int(header, "NAXIS", 2);
 
      if(cpl_propertylist_has(header, "NAXIS3")){
      cpl_propertylist_erase(header, "NAXIS3");
      }
  }
 
  /* for 2D images */
  if (crval) {
      cpl_propertylist_update_double(header, "CRVAL1",
                     cpl_array_get_double(crval, 0, &err));
      cpl_propertylist_update_double(header, "CRVAL2",
                     cpl_array_get_double(crval, 1, &err));
  }
 
  if (crpix) {
      cpl_propertylist_update_double(header, "CRPIX1",
                     cpl_array_get_double(crpix, 0, &err));
      cpl_propertylist_update_double(header, "CRPIX2",
                     cpl_array_get_double(crpix, 1, &err));
  }
 
  if (ctype) {
      cpl_propertylist_update_string(header, "CTYPE1",
                     cpl_array_get_string(ctype, 0));
      cpl_propertylist_update_string(header, "CTYPE2",
                     cpl_array_get_string(ctype, 1));
  }
 
  if (cunit) {
      cpl_propertylist_update_string(header, "CUNIT1",
                     cpl_array_get_string(cunit, 0));
      cpl_propertylist_update_string(header, "CUNIT2",
                     cpl_array_get_string(cunit, 1));
  }
 
  if (cd) {
      double cd11 = cpl_matrix_get(cd, 0, 0);
      double cd12 = cpl_matrix_get(cd, 0, 1);
      double cd21 = cpl_matrix_get(cd, 1, 0);
      double cd22 = cpl_matrix_get(cd, 1, 1);
      cpl_propertylist_update_double(header, "CD1_1", cd11);
      cpl_propertylist_update_double(header, "CD1_2", cd12);
      cpl_propertylist_update_double(header, "CD2_1", cd21);
      cpl_propertylist_update_double(header, "CD2_2", cd22);
  }
 
  /* for 3D cubes */
  if (only2d == FALSE && cpl_array_get_size(crval) > 2) {
 
      if (crval) {
      cpl_propertylist_update_double(header, "CRVAL3",
                     cpl_array_get_double(crval, 2, &err));
      }
 
      if (crpix) {
      cpl_propertylist_update_double(header, "CRPIX3",
                     cpl_array_get_double(crpix, 2, &err));
      }
 
      if (ctype) {
      cpl_propertylist_update_string(header, "CTYPE3",
                     cpl_array_get_string(ctype, 2));
      }
 
      if (cunit) {
      cpl_propertylist_update_string(header, "CUNIT3",
                     cpl_array_get_string(cunit, 2));
      }
 
      if (cd) {
      double cd13 = cpl_matrix_get(cd, 0, 2);
      double cd23 = cpl_matrix_get(cd, 1, 2);
      double cd31 = cpl_matrix_get(cd, 2, 0);
      double cd32 = cpl_matrix_get(cd, 2, 1);
      double cd33 = cpl_matrix_get(cd, 2, 2);
      cpl_propertylist_update_double(header, "CD1_3", cd13);
      cpl_propertylist_update_double(header, "CD2_3", cd23);
      cpl_propertylist_update_double(header, "CD3_1", cd31);
      cpl_propertylist_update_double(header, "CD3_2", cd32);
      cpl_propertylist_update_double(header, "CD3_3", cd33);
      }
  } /* if 3D */
  return CPL_ERROR_NONE;
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
static cpl_error_code
hdrl_resample_create_table(cpl_table** tab, const cpl_size size)
{
 
    cpl_ensure_code(tab, CPL_ERROR_NULL_INPUT);
    cpl_ensure_code(size > 0 , CPL_ERROR_ILLEGAL_INPUT);
 
  *tab = cpl_table_new(size);
  /*
   * MUSE_PIXTABLE_XPOS    "xpos" *< x-position column of a MUSE pixel table
   * MUSE_PIXTABLE_YPOS    "ypos" *< y-position column of a MUSE pixel table
   * MUSE_PIXTABLE_LAMBDA  "lambda" *< wavelength column of a MUSE pixel table
   * MUSE_PIXTABLE_DATA    "data" *< data column of a MUSE pixel table
   * MUSE_PIXTABLE_DQ      "dq" *< data quality column of a MUSE pixel table
   * MUSE_PIXTABLE_STAT    "stat" *< error column of a MUSE pixel table
   */
  cpl_table_new_column(*tab,
                       HDRL_RESAMPLE_TABLE_RA,     HDRL_RESAMPLE_TABLE_RA_TYPE);
  cpl_table_new_column(*tab,
                       HDRL_RESAMPLE_TABLE_DEC,    HDRL_RESAMPLE_TABLE_DEC_TYPE);
  cpl_table_new_column(*tab,
                       HDRL_RESAMPLE_TABLE_LAMBDA, HDRL_RESAMPLE_TABLE_LAMBDA_TYPE);
  cpl_table_new_column(*tab,
                       HDRL_RESAMPLE_TABLE_DATA,   HDRL_RESAMPLE_TABLE_DATA_TYPE);
  cpl_table_new_column(*tab,
                       HDRL_RESAMPLE_TABLE_BPM,    HDRL_RESAMPLE_TABLE_BPM_TYPE);
  cpl_table_new_column(*tab,
                       HDRL_RESAMPLE_TABLE_ERRORS, HDRL_RESAMPLE_TABLE_ERRORS_TYPE);
 
  /* init column values */
  cpl_table_fill_column_window_double(*tab, HDRL_RESAMPLE_TABLE_RA,    0, size, 0.);
  cpl_table_fill_column_window_double(*tab, HDRL_RESAMPLE_TABLE_DEC,   0, size, 0.);
  cpl_table_fill_column_window_double(*tab, HDRL_RESAMPLE_TABLE_LAMBDA,0, size, 0.);
  cpl_table_fill_column_window_double(*tab, HDRL_RESAMPLE_TABLE_DATA,  0, size, 0.);
  cpl_table_fill_column_window_int(*tab,    HDRL_RESAMPLE_TABLE_BPM,   0, size, 0);
  cpl_table_fill_column_window_double(*tab, HDRL_RESAMPLE_TABLE_ERRORS,0, size, 0.);
 
  return cpl_error_get_code();
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
cpl_table*
hdrl_resample_image_to_table(const hdrl_image * hima, const cpl_wcs* wcs)
{
  cpl_ensure(hima, CPL_ERROR_NULL_INPUT, NULL);
  cpl_ensure(wcs, CPL_ERROR_NULL_INPUT, NULL);
  cpl_msg_debug(cpl_func, "Converting Data to table");
  hdrl_imagelist* ilist = NULL;
  cpl_table* tab;
  ilist = hdrl_imagelist_new();
  hdrl_imagelist_set(ilist, (hdrl_image *)hima, 0);
 
  tab = hdrl_resample_imagelist_to_table(ilist, wcs);
 
  /* cleanup memory */
  hdrl_imagelist_unset(ilist, 0);
  hdrl_imagelist_delete(ilist);
 
 
  return tab;
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
cpl_table*
hdrl_resample_imagelist_to_table(const hdrl_imagelist* himlist,
                 const cpl_wcs* wcs)
{
 
  cpl_ensure(himlist, CPL_ERROR_NULL_INPUT, NULL);
  cpl_ensure(wcs, CPL_ERROR_NULL_INPUT, NULL);
 
  cpl_msg_debug(cpl_func, "Converting Dataset to table");
 
  cpl_size naxis1 = hdrl_imagelist_get_size_x(himlist);
  cpl_size naxis2 = hdrl_imagelist_get_size_y(himlist);
  cpl_size naxis3 = hdrl_imagelist_get_size(himlist);
 
  cpl_msg_debug(cpl_func,"Dataset dimentions (x, y, l): (%lld, %lld, %lld)",
           naxis1, naxis2, naxis3);
 
  const cpl_array *crval = cpl_wcs_get_crval(wcs);
  const cpl_array *crpix = cpl_wcs_get_crpix(wcs);
  const cpl_matrix *cd = cpl_wcs_get_cd(wcs);
 
  int testerr = 0;
  double crpix3 = 0.;
  double crval3 = 0.;
  double cdelt3 = 0.;
 
  if (naxis3 > 1) { /* We have a cube */
      crpix3 = cpl_array_get_double(crpix, 2, &testerr);
      crval3 = cpl_array_get_double(crval, 2, &testerr);
      cdelt3 = cpl_matrix_get(cd, 2, 2); /* CD3_3 */
/*
      cpl_msg_debug(cpl_func,"crpix3: %g crval3: %g cdelt3: %g", crpix3,
            crval3, cdelt3);
*/
  }
 
  cpl_size tab_size = naxis1 * naxis2 * naxis3;
  cpl_table* tab = NULL;
  /* Prefill the full table */
  hdrl_resample_create_table(&tab, tab_size);
 
  double* ptabxpos = cpl_table_get_data_double(tab, HDRL_RESAMPLE_TABLE_RA);
  double* ptabypos = cpl_table_get_data_double(tab, HDRL_RESAMPLE_TABLE_DEC);
  double* ptablambda = cpl_table_get_data_double(tab, HDRL_RESAMPLE_TABLE_LAMBDA);
  double* ptabdata = cpl_table_get_data_double(tab, HDRL_RESAMPLE_TABLE_DATA);
  int*   ptabbpm = cpl_table_get_data_int(tab, HDRL_RESAMPLE_TABLE_BPM);
  double* ptaberr = cpl_table_get_data_double(tab, HDRL_RESAMPLE_TABLE_ERRORS);
 
  struct timeval tv1, tv2;
  cpl_msg_debug(cpl_func,"Starting parallel loop in hdrl_imagelist_to_table");
  gettimeofday(&tv1, NULL);
 
  HDRL_OMP(omp parallel for collapse(2))
  for(cpl_size k = 0; k < naxis3; k++) {
      for(cpl_size j = 0; j < naxis2; j++) {
      const double* pimaerr = NULL;
      const cpl_binary *   pimabpm = NULL;
 
      /* Fill the data */
      const hdrl_image* hima = hdrl_imagelist_get_const(himlist, k);
      const cpl_image* imadata = hdrl_image_get_image_const(hima);
      const cpl_image* imaerrs = hdrl_image_get_error_const(hima);
      const cpl_mask*  imamask  = hdrl_image_get_mask_const(hima);
      const double* pimadata = cpl_image_get_data_double_const(imadata);
 
      if (imaerrs) { /* Fill the errors */
          pimaerr = cpl_image_get_data_double_const(imaerrs);
      }
      if (imamask) { /* Fill the bpm */
          pimabpm = cpl_mask_get_data_const(imamask);
      }
 
      cpl_size k_naxis1_naxis2 = naxis1 * naxis2 * k;
      cpl_size j_naxis1 = j * naxis1;
      for(cpl_size i = 0; i < naxis1; i++) {
          cpl_size raw = k_naxis1_naxis2 + j_naxis1 + i;
          hdrl_wcs_xy_to_radec(wcs, i+1 , j+1 , &ptabxpos[raw],
                   &ptabypos[raw]);
          ptabdata[raw] = pimadata[j_naxis1 + i];
          if (naxis3 > 1) ptablambda[raw] = crval3 + cdelt3 * (k - crpix3 + 1.);
          if (imaerrs) ptaberr[raw] = pimaerr[j_naxis1 + i];
          if (imamask) ptabbpm[raw] = pimabpm[j_naxis1 + i];
          /* Insert only good pixels */
          if (!isfinite(pimadata[j_naxis1 + i]) ||
          ptabbpm[raw] != CPL_BINARY_0 ){
          ptabbpm[raw] = CPL_BINARY_1;
 
          }
      }
      }
  }
 
  gettimeofday(&tv2, NULL);
  cpl_msg_debug(cpl_func, "Wall time for hdrl_imagelist_to_table was %f seconds\n",
           (double) (tv2.tv_usec - tv1.tv_usec) / 1000000 +
           (double) (tv2.tv_sec - tv1.tv_sec));
 
  return tab;
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
hdrl_parameter*
hdrl_resample_parameter_create_outgrid2D(const double delta_ra,
        const double delta_dec)
{
 
  hdrl_resample_outgrid_parameter * p = (hdrl_resample_outgrid_parameter *)
          hdrl_parameter_new(&hdrl_resample_outgrid_parameter_type);
  p->method = HDRL_RESAMPLE_OUTGRID_2D;
  p->delta_ra = delta_ra;
  p->delta_dec = delta_dec;
 
  p->recalc_limits = CPL_TRUE;
 
  /* This function asks to recalculate the limits in the hdrl_resample_compute
   * function - therefore we put dummy values for the moment. */
 
  p->dec_min = 0.1;
  p->dec_max = 0.2;
  p->ra_min = 0.1;
  p->ra_max = 0.2;
 
  /* in case of 2D sets some defaults dummy values for 3rd dimension */
  p->lambda_min = 0;
  p->lambda_max = 0;
  p->delta_lambda = 1;
 
  /* Default field margin in percent taken from swarp. */
  p->fieldmargin = FIELDMARGIN;
 
  p->wcs = NULL;
 
  if (hdrl_resample_parameter_outgrid_verify((hdrl_parameter*)p) != CPL_ERROR_NONE) {
      cpl_free(p);
      return NULL;
  }
  return (hdrl_parameter *)p;
 
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
hdrl_parameter*
hdrl_resample_parameter_create_outgrid3D(const double delta_ra,
        const double delta_dec, const double delta_lambda)
{
 
  hdrl_resample_outgrid_parameter * p = (hdrl_resample_outgrid_parameter *)
      hdrl_parameter_new(&hdrl_resample_outgrid_parameter_type);
  p->method = HDRL_RESAMPLE_OUTGRID_3D;
  p->delta_ra = delta_ra;
  p->delta_dec = delta_dec;
  p->delta_lambda = delta_lambda;
 
  p->recalc_limits = CPL_TRUE;
  /* This function asks to recalculate the limits in the hdrl_resample_compute
   * function - therefore we put dummy values for the moment. */
 
  p->dec_min = 0.1;
  p->dec_max = 0.2;
  p->ra_min = 0.1;
  p->ra_max = 0.2;
  p->lambda_min = 0.;
  p->lambda_max = 0.;
 
  /* Default field margin in percent taken from swarp. */
  p->fieldmargin = FIELDMARGIN;
 
  p->wcs = NULL;
 
  if (hdrl_resample_parameter_outgrid_verify((hdrl_parameter*)p) != CPL_ERROR_NONE) {
      cpl_free(p);
      return NULL;
  }
  return (hdrl_parameter *)p;
 
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
hdrl_parameter*
hdrl_resample_parameter_create_outgrid2D_userdef(const double delta_ra,
                         const double delta_dec,
                         const double ra_min,
                         const double ra_max,
                         const double dec_min,
                         const double dec_max,
                         const double fieldmargin)
{
 
  hdrl_resample_outgrid_parameter * p = (hdrl_resample_outgrid_parameter *)
      hdrl_parameter_new(&hdrl_resample_outgrid_parameter_type);
  p->method = HDRL_RESAMPLE_OUTGRID_2D;
  p->delta_ra = delta_ra;
  p->delta_dec = delta_dec;
 
 
  p->recalc_limits = CPL_FALSE; /*This function takes the limits from the user*/
  p->dec_min = dec_min;
  p->dec_max = dec_max;
  p->ra_min = ra_min;
  p->ra_max = ra_max;
 
  /* in case of 2D sets some defaults dummy values for 3rd dimension */
  p->lambda_min = 0.;
  p->lambda_max = 0.;
  p->delta_lambda = 1.;
 
  p->fieldmargin = fieldmargin;
 
  p->wcs = NULL;
 
  if (hdrl_resample_parameter_outgrid_verify((hdrl_parameter*)p) != CPL_ERROR_NONE) {
      cpl_free(p);
      return NULL;
  }
  return (hdrl_parameter *)p;
 
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
hdrl_parameter*
hdrl_resample_parameter_create_outgrid3D_userdef(const double delta_ra,
                         const double delta_dec,
                         const double delta_lambda,
                         const double ra_min,
                         const double ra_max,
                         const double dec_min,
                         const double dec_max,
                         const double lambda_min,
                         const double lambda_max,
                         const double fieldmargin)
{
 
  hdrl_resample_outgrid_parameter * p = (hdrl_resample_outgrid_parameter *)
      hdrl_parameter_new(&hdrl_resample_outgrid_parameter_type);
  p->method = HDRL_RESAMPLE_OUTGRID_3D;
  p->delta_ra = delta_ra;
  p->delta_dec = delta_dec;
  p->delta_lambda = delta_lambda;
 
  p->recalc_limits = CPL_FALSE; /*This function takes the limits from the user*/
 
  p->dec_min = dec_min;
  p->dec_max = dec_max;
  p->ra_min = ra_min;
  p->ra_max = ra_max;
 
  p->lambda_min = lambda_min;
  p->lambda_max = lambda_max;
 
  p->fieldmargin = fieldmargin;
 
  p->wcs = NULL;
 
  if (hdrl_resample_parameter_outgrid_verify((hdrl_parameter*)p) != CPL_ERROR_NONE) {
      cpl_free(p);
      return NULL;
  }
  return (hdrl_parameter *)p;
 
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
hdrl_parameter*
hdrl_resample_parameter_create_renka(const int loop_distance,
                     cpl_boolean use_errorweights,
                     const double critical_radius)
{
 
  hdrl_resample_method_parameter * p = (hdrl_resample_method_parameter *)
      hdrl_parameter_new(&hdrl_resample_method_parameter_type);
 
  p->method = HDRL_RESAMPLE_METHOD_RENKA;
  p->loop_distance = loop_distance;
  p->use_errorweights = use_errorweights;
  p->renka_critical_radius = critical_radius;
 
  /* fill rest with dummy input */
  p->drizzle_pix_frac_x = 0.1;
  p->drizzle_pix_frac_y = 0.1;
  p->drizzle_pix_frac_lambda = 0.1;
  p->lanczos_kernel_size = 2;
 
  if (hdrl_resample_parameter_method_verify((hdrl_parameter*)p) != CPL_ERROR_NONE) {
      cpl_free(p);
      return NULL;
  }
  return (hdrl_parameter*) p;
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
hdrl_parameter*
hdrl_resample_parameter_create_linear(const int loop_distance,
                      cpl_boolean use_errorweights)
{
 
  hdrl_resample_method_parameter * p = (hdrl_resample_method_parameter *)
      hdrl_parameter_new(&hdrl_resample_method_parameter_type);
 
 
  p->method = HDRL_RESAMPLE_METHOD_LINEAR;
  p->loop_distance = loop_distance;
  p->use_errorweights = use_errorweights;
 
  /* fill rest with dummy input */
  p->renka_critical_radius = 0.1;
  p->drizzle_pix_frac_x = 0.1;
  p->drizzle_pix_frac_y = 0.1;
  p->drizzle_pix_frac_lambda = 0.1;
  p->lanczos_kernel_size = 2;
 
  if (hdrl_resample_parameter_method_verify((hdrl_parameter*)p) != CPL_ERROR_NONE) {
      cpl_free(p);
      return NULL;
  }
  return (hdrl_parameter*) p;
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
hdrl_parameter*
hdrl_resample_parameter_create_quadratic(const int loop_distance,
                     cpl_boolean use_errorweights)
{
 
  hdrl_resample_method_parameter * p = (hdrl_resample_method_parameter *)
      hdrl_parameter_new(&hdrl_resample_method_parameter_type);
 
  p->method = HDRL_RESAMPLE_METHOD_QUADRATIC;
  p->loop_distance = loop_distance;
  p->use_errorweights = use_errorweights;
 
  /* fill rest with dummy input */
  p->renka_critical_radius = 0.1;
  p->drizzle_pix_frac_x = 0.1;
  p->drizzle_pix_frac_y = 0.1;
  p->drizzle_pix_frac_lambda = 0.1;
  p->lanczos_kernel_size = 2;
 
  if (hdrl_resample_parameter_method_verify((hdrl_parameter*)p) != CPL_ERROR_NONE) {
      cpl_free(p);
      return NULL;
  }
  return (hdrl_parameter*) p;
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
hdrl_parameter*
hdrl_resample_parameter_create_nearest(void)
{
 
  hdrl_resample_method_parameter * p = (hdrl_resample_method_parameter *)
      hdrl_parameter_new(&hdrl_resample_method_parameter_type);
 
 
  p->method = HDRL_RESAMPLE_METHOD_NEAREST;
  p->loop_distance = 0;
  p->use_errorweights = CPL_FALSE;
 
  /* fill rest with dummy input */
  p->renka_critical_radius = 0.1;
  p->drizzle_pix_frac_x = 0.1;
  p->drizzle_pix_frac_y = 0.1;
  p->drizzle_pix_frac_lambda = 0.1;
  p->lanczos_kernel_size = 2;
 
  if (hdrl_resample_parameter_method_verify((hdrl_parameter*)p) != CPL_ERROR_NONE) {
      cpl_free(p);
      return NULL;
  }
  return (hdrl_parameter*) p;
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
hdrl_parameter*
hdrl_resample_parameter_create_lanczos(const int loop_distance,
                       cpl_boolean use_errorweights,
                       const int kernel_size)
{
 
  hdrl_resample_method_parameter * p = (hdrl_resample_method_parameter *)
      hdrl_parameter_new(&hdrl_resample_method_parameter_type);
 
  p->method = HDRL_RESAMPLE_METHOD_LANCZOS;
  p->loop_distance = loop_distance;
  p->use_errorweights = use_errorweights;
  p->lanczos_kernel_size = kernel_size;
  /* fill rest with dummy input */
  p->renka_critical_radius = 0.1;
  p->drizzle_pix_frac_x = 0.1;
  p->drizzle_pix_frac_y = 0.1;
  p->drizzle_pix_frac_lambda = 0.1;
 
  if (hdrl_resample_parameter_method_verify((hdrl_parameter*)p) != CPL_ERROR_NONE) {
      cpl_free(p);
      return NULL;
  }
  return (hdrl_parameter*) p;
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
hdrl_parameter*
hdrl_resample_parameter_create_drizzle(const int loop_distance,
                       cpl_boolean use_errorweights,
                       const double pix_frac_x,
                       const double pix_frac_y,
                       const double pix_frac_lambda)
{
 
  hdrl_resample_method_parameter * p = (hdrl_resample_method_parameter *)
      hdrl_parameter_new(&hdrl_resample_method_parameter_type);
 
  p->method = HDRL_RESAMPLE_METHOD_DRIZZLE;
  p->loop_distance = loop_distance;
  p->use_errorweights = use_errorweights;
  p->drizzle_pix_frac_x = pix_frac_x;
  p->drizzle_pix_frac_y = pix_frac_y;
  p->drizzle_pix_frac_lambda = pix_frac_lambda;
 
  /* fill rest with dummy input */
  p->renka_critical_radius = 0.1;
  p->lanczos_kernel_size = 2;
 
  if (hdrl_resample_parameter_method_verify((hdrl_parameter*)p) != CPL_ERROR_NONE) {
      cpl_free(p);
      return NULL;
  }
  return (hdrl_parameter*) p;
 
}
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
cpl_error_code
hdrl_resample_parameter_outgrid_verify(const hdrl_parameter * hp){
 
  const hdrl_resample_outgrid_parameter * param_loc =
      (const hdrl_resample_outgrid_parameter *)hp ;
 
  cpl_error_ensure(hp != NULL, CPL_ERROR_NULL_INPUT,
           return CPL_ERROR_NULL_INPUT, "NULL Input Parameters");
 
  cpl_error_ensure(hdrl_resample_parameter_outgrid_check(hp),
           CPL_ERROR_ILLEGAL_INPUT, return CPL_ERROR_ILLEGAL_INPUT,
               "Here we expect a resample outgrid parameter") ;
 
  cpl_error_ensure(
      param_loc->recalc_limits == CPL_TRUE ||
      param_loc->recalc_limits == CPL_FALSE,
      CPL_ERROR_ILLEGAL_INPUT, return CPL_ERROR_ILLEGAL_INPUT,
      "Unsupported resample recalc_limits value");
 
/*
 * The wcs is filled later on by the compute function so it can not be checked
 * at this stage
 *
  cpl_error_ensure(param_loc->wcs != NULL, CPL_ERROR_NULL_INPUT,
           return CPL_ERROR_NULL_INPUT, "WCS must be defined");
*/
 
 
  cpl_error_ensure(param_loc->delta_ra > 0, CPL_ERROR_ILLEGAL_INPUT,
           return CPL_ERROR_ILLEGAL_INPUT, "right ascension "
               "stepsize must be > 0");
 
  cpl_error_ensure(param_loc->delta_dec > 0, CPL_ERROR_ILLEGAL_INPUT,
           return CPL_ERROR_ILLEGAL_INPUT, "declination stepsize "
               "must be > 0");
 
  cpl_error_ensure(param_loc->delta_lambda > 0, CPL_ERROR_ILLEGAL_INPUT,
           return CPL_ERROR_ILLEGAL_INPUT, "wavelength stepsize "
               "must be > 0");
 
 
  cpl_error_ensure(param_loc->ra_min >= 0,
           CPL_ERROR_ILLEGAL_INPUT, return CPL_ERROR_ILLEGAL_INPUT,
               "Minimum right ascension must be >= 0");
 
  cpl_error_ensure(param_loc->ra_max >= 0,
           CPL_ERROR_ILLEGAL_INPUT, return CPL_ERROR_ILLEGAL_INPUT,
               "Maximum right ascension must be >= 0");
 
 
  cpl_error_ensure(param_loc->lambda_min >= 0,
           CPL_ERROR_ILLEGAL_INPUT, return CPL_ERROR_ILLEGAL_INPUT,
               "Minimum wavelength must be >= 0");
 
  cpl_error_ensure(param_loc->lambda_max >= 0,
           CPL_ERROR_ILLEGAL_INPUT, return CPL_ERROR_ILLEGAL_INPUT,
               "Maximum wavelength must be >= 0");
 
  cpl_error_ensure(param_loc->fieldmargin >= 0., CPL_ERROR_ILLEGAL_INPUT,
           return CPL_ERROR_ILLEGAL_INPUT, "The field margin must"
               " be >= 0.");
 
  cpl_error_ensure(param_loc->ra_max >= param_loc->ra_min ,
           CPL_ERROR_ILLEGAL_INPUT, return CPL_ERROR_ILLEGAL_INPUT,
               "The maximum right ascension must be >= "
               "the minimum right ascension");
  cpl_error_ensure(param_loc->dec_max >= param_loc->dec_min ,
           CPL_ERROR_ILLEGAL_INPUT, return CPL_ERROR_ILLEGAL_INPUT,
               "The maximum declination must be >= "
               "the minimum declination");
 
  cpl_error_ensure(param_loc->lambda_max >= param_loc->lambda_min ,
           CPL_ERROR_ILLEGAL_INPUT, return CPL_ERROR_ILLEGAL_INPUT,
               "The maximum wavelength must be >= "
               "the minimum wavelength");
 
  return CPL_ERROR_NONE ;
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
cpl_error_code
hdrl_resample_parameter_method_verify(const hdrl_parameter * hp){
 
  const hdrl_resample_method_parameter * param_loc =
      (const hdrl_resample_method_parameter *)hp ;
  cpl_error_ensure(hp != NULL, CPL_ERROR_NULL_INPUT,
           return CPL_ERROR_NULL_INPUT, "NULL Input Parameters");
 
  cpl_error_ensure(hdrl_resample_parameter_method_check(hp),
           CPL_ERROR_ILLEGAL_INPUT, return CPL_ERROR_ILLEGAL_INPUT,
               "Here we expect a resample method parameter") ;
 
  /* checks on parameter methods */
  cpl_error_ensure(
      param_loc->method == HDRL_RESAMPLE_METHOD_NEAREST ||
      param_loc->method == HDRL_RESAMPLE_METHOD_LINEAR ||
      param_loc->method == HDRL_RESAMPLE_METHOD_QUADRATIC ||
      param_loc->method == HDRL_RESAMPLE_METHOD_LANCZOS ||
      param_loc->method == HDRL_RESAMPLE_METHOD_DRIZZLE ||
      param_loc->method == HDRL_RESAMPLE_METHOD_RENKA,
      CPL_ERROR_ILLEGAL_INPUT, return CPL_ERROR_ILLEGAL_INPUT,
      "Unsupported resample method");
 
  /* checks on common parameter elements */
  cpl_error_ensure(param_loc->loop_distance >= 0, CPL_ERROR_ILLEGAL_INPUT,
           return CPL_ERROR_ILLEGAL_INPUT, "The loop distance must "
               "be >=0");
 
  cpl_error_ensure(param_loc->use_errorweights == CPL_TRUE ||
           param_loc->use_errorweights == CPL_FALSE,
           CPL_ERROR_ILLEGAL_INPUT, return CPL_ERROR_ILLEGAL_INPUT,
               "Unsupported resample use_errorweights value");
 
  switch (param_loc->method) {
    case HDRL_RESAMPLE_METHOD_NEAREST:
    case HDRL_RESAMPLE_METHOD_LINEAR:
    case HDRL_RESAMPLE_METHOD_QUADRATIC:
    case HDRL_RESAMPLE_METHOD_NONE:
      break;
    case HDRL_RESAMPLE_METHOD_RENKA:
      /* checks on peculiar parameter elements */
      cpl_error_ensure(param_loc->renka_critical_radius > 0.,
               CPL_ERROR_ILLEGAL_INPUT, return CPL_ERROR_ILLEGAL_INPUT,
               "Critical radius of the Renka method must be > 0");
      break ;
 
    case HDRL_RESAMPLE_METHOD_DRIZZLE:
      /* checks on peculiar parameter elements */
      cpl_error_ensure(param_loc->drizzle_pix_frac_x > 0.,
               CPL_ERROR_ILLEGAL_INPUT, return CPL_ERROR_ILLEGAL_INPUT,
               "Drizzle down-scaling factor in x direction "
               "must be > 0");
 
      cpl_error_ensure(param_loc->drizzle_pix_frac_y > 0.,
               CPL_ERROR_ILLEGAL_INPUT, return CPL_ERROR_ILLEGAL_INPUT,
               "Drizzle down-scaling factor in y direction "
               "must be > 0");
 
      cpl_error_ensure(param_loc->drizzle_pix_frac_lambda > 0.,
               CPL_ERROR_ILLEGAL_INPUT, return CPL_ERROR_ILLEGAL_INPUT,
               "Drizzle down-scaling factor in z/lambda direction "
               "must be > 0");
      break;
 
    case HDRL_RESAMPLE_METHOD_LANCZOS:
      /* checks on peculiar parameter elements */
      cpl_error_ensure(param_loc->lanczos_kernel_size > 0,
               CPL_ERROR_ILLEGAL_INPUT, return CPL_ERROR_ILLEGAL_INPUT,
               "The kernel size of the Lanczos method must be > 0");
      break;
  }
 
  return CPL_ERROR_NONE ;
}
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
int
hdrl_resample_parameter_outgrid_check(const hdrl_parameter * self){
  /* Check if the method is of proper type */
  cpl_ensure_code(self != NULL, CPL_ERROR_NULL_INPUT);
 
  return hdrl_parameter_check_type(self, &hdrl_resample_outgrid_parameter_type);
}
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
int
hdrl_resample_parameter_method_check(const hdrl_parameter * self){
  /* Check if the method is of proper type */
  cpl_ensure_code(self != NULL, CPL_ERROR_NULL_INPUT);
  return hdrl_parameter_check_type(self, &hdrl_resample_method_parameter_type);
}
 
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
static cpl_error_code
hdrl_resample_inputtable_verify(const cpl_table *ResTable){
 
  cpl_error_ensure(ResTable != NULL, CPL_ERROR_NULL_INPUT,
           return CPL_ERROR_NULL_INPUT, "No Table as input");
 
  /* Check the existence of all columns */
  cpl_error_ensure(cpl_table_has_column(ResTable, HDRL_RESAMPLE_TABLE_DATA  )
           == 1, CPL_ERROR_INCOMPATIBLE_INPUT,
           return CPL_ERROR_INCOMPATIBLE_INPUT,
               "Missing data table column");
  cpl_error_ensure(cpl_table_has_column(ResTable, HDRL_RESAMPLE_TABLE_BPM   )
           == 1, CPL_ERROR_INCOMPATIBLE_INPUT,
           return CPL_ERROR_INCOMPATIBLE_INPUT,
               "Missing bpm table column");
  cpl_error_ensure(cpl_table_has_column(ResTable, HDRL_RESAMPLE_TABLE_ERRORS)
           == 1, CPL_ERROR_INCOMPATIBLE_INPUT,
           return CPL_ERROR_INCOMPATIBLE_INPUT,
               "Missing error table column");
  cpl_error_ensure(cpl_table_has_column(ResTable, HDRL_RESAMPLE_TABLE_RA    )
           == 1, CPL_ERROR_INCOMPATIBLE_INPUT,
           return CPL_ERROR_INCOMPATIBLE_INPUT,
               "Missing right ascension table column");
  cpl_error_ensure(cpl_table_has_column(ResTable, HDRL_RESAMPLE_TABLE_DEC   )
           == 1, CPL_ERROR_INCOMPATIBLE_INPUT,
           return CPL_ERROR_INCOMPATIBLE_INPUT,
               "Missing declination table column");
  cpl_error_ensure(cpl_table_has_column(ResTable, HDRL_RESAMPLE_TABLE_LAMBDA)
           == 1, CPL_ERROR_INCOMPATIBLE_INPUT,
           return CPL_ERROR_INCOMPATIBLE_INPUT,
               "Missing wavelength table column");
 
  /* Check the format of all columns */
  cpl_error_ensure(cpl_table_get_column_type(ResTable, HDRL_RESAMPLE_TABLE_DATA  )
           == HDRL_RESAMPLE_TABLE_DATA_TYPE, CPL_ERROR_INCOMPATIBLE_INPUT,
           return CPL_ERROR_INCOMPATIBLE_INPUT,
               "Data table column has wrong format");
  cpl_error_ensure(cpl_table_get_column_type(ResTable, HDRL_RESAMPLE_TABLE_BPM   )
           == HDRL_RESAMPLE_TABLE_BPM_TYPE,    CPL_ERROR_INCOMPATIBLE_INPUT,
           return CPL_ERROR_INCOMPATIBLE_INPUT,
               "Bpm table column has wrong format");
  cpl_error_ensure(cpl_table_get_column_type(ResTable, HDRL_RESAMPLE_TABLE_ERRORS)
           == HDRL_RESAMPLE_TABLE_ERRORS_TYPE, CPL_ERROR_INCOMPATIBLE_INPUT,
           return CPL_ERROR_INCOMPATIBLE_INPUT,
               "Error table column has wrong format");
  cpl_error_ensure(cpl_table_get_column_type(ResTable, HDRL_RESAMPLE_TABLE_RA    )
           == HDRL_RESAMPLE_TABLE_RA_TYPE, CPL_ERROR_INCOMPATIBLE_INPUT,
           return CPL_ERROR_INCOMPATIBLE_INPUT,
               "Right ascension table column has wrong format");
  cpl_error_ensure(cpl_table_get_column_type(ResTable, HDRL_RESAMPLE_TABLE_DEC   )
           == HDRL_RESAMPLE_TABLE_DEC_TYPE, CPL_ERROR_INCOMPATIBLE_INPUT,
           return CPL_ERROR_INCOMPATIBLE_INPUT,
               "Declination table column has wrong format");
  cpl_error_ensure(cpl_table_get_column_type(ResTable, HDRL_RESAMPLE_TABLE_LAMBDA)
           == HDRL_RESAMPLE_TABLE_LAMBDA_TYPE, CPL_ERROR_INCOMPATIBLE_INPUT,
           return CPL_ERROR_INCOMPATIBLE_INPUT,
               "Wavelength table column has wrong format");
 
  return cpl_error_get_code();
}
/* EOF */