sofi_wavelength.c

/* $Id: sofi_wavelength.c,v 1.52 2013-03-12 08:04:54 llundin Exp $
 *
 * This file is part of the SOFI Pipeline
 * Copyright (C) 2002,2003 European Southern Observatory
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02111-1307  USA
 */
 
/*
 * $Author: llundin $
 * $Date: 2013-03-12 08:04:54 $
 * $Revision: 1.52 $
 * $Name: not supported by cvs2svn $
 */
 
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
 
/*-----------------------------------------------------------------------------
                                   Includes
 -----------------------------------------------------------------------------*/
 
#include "sofi_wavelength.h"
#include "sofi_pfits.h"
#include "sofi_utils.h"
#include "sofi_dfs.h"
 
#include "irplib_wavecal.h"
#include "irplib_utils.h"
#include "irplib_wlxcorr.h"
#include "irplib_ppm.h"
 
#include <math.h>
#include <string.h>
#include <float.h>
 
/*-----------------------------------------------------------------------------
                                   Includes
 -----------------------------------------------------------------------------*/
 
static cpl_bivector * sofi_wavelength_load_catalog(const char *, const char *, 
        const char *, const char *);
static cpl_vector * sofi_wavelength_extract_spectrum(const cpl_image *, double);
static cpl_polynomial * sofi_wavelength_xc_engine(const cpl_vector *,
        const cpl_bivector *, int, const cpl_polynomial *, double, int,
        double, int, double *);
 
 
/*-----------------------------------------------------------------------------
                                   Defines
 -----------------------------------------------------------------------------*/
 
#ifndef SOFI_WFWHM
/* FIXME: Find a better value ? */
#define SOFI_WFWHM 4.0
#endif
 
#ifndef SOFI_PIXSTEP
/* FIXME: Parameterize (and find a better value ?) */
#define SOFI_PIXSTEP 0.5
#endif
 
#ifndef SOFI_PIXTOL
/* FIXME: Parameterize (and find a better value ?) */
#define SOFI_PIXTOL 1e-5
#endif
 
#ifndef SOFI_OFFSET
/* FIXME: Parameterize (and find a better value ?) */
#define SOFI_OFFSET 50
#endif
 
#ifndef SOFI_DEGREE
/* FIXME */
#define SOFI_DEGREE 3
#endif
 
#ifndef SOFI_MAXFAIL
/* FIXME */
#define SOFI_MAXFAIL 3
#endif
 
#ifndef SOFI_CONTINUE
/* FIXME */
#define SOFI_CONTINUE 3
#endif
 
#define SOFI_MAX(A,B) ((A) > (B) ? (A) : (B))
#define SOFI_MIN(A,B) ((A) < (B) ? (A) : (B))
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
 
/*-----------------------------------------------------------------------------
                                Function codes
 -----------------------------------------------------------------------------*/

 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
cpl_polynomial * sofi_wavelength_engine(
        const cpl_image *   in,
        const char      *   table_name,
        const char      *   oh,
        const char      *   xe,
        const char      *   ne,
        const cpl_polynomial * phdisprel,
        double              slit_width,
        int                 degree,
        double              wl_err,
        int                 nsamples,
        int                 use_ppm,
        int                 plot,
        double          *   xcorr)
{
    cpl_polynomial          *   wl_poly_xc;
    cpl_polynomial          *   wl_poly_ppm;
    double                      xc;
    cpl_bivector            *   spec_cat;
    cpl_vector              *   spectrum;
    cpl_vector              *   filtered;
    cpl_vector              *   clean;
    double                  *   pclean;
    int                         filt_size = 9;
    double                      stdev, median, threshold;
    int                         i;
 
    /* Check inputs */
    if (in==NULL || table_name==NULL || phdisprel == NULL) return NULL;
    if (degree > 3 || degree < 1) return NULL;
 
    degree = SOFI_DEGREE; /* FIXME */
 
    /* Load the catalog */
    if ((spec_cat=sofi_wavelength_load_catalog(table_name, oh, xe, ne))==NULL) {
        cpl_msg_error(cpl_func, "Cannot load the catalog");
        return NULL;
    }
    
    /* Get the spectrum from the image */
    if ((spectrum = sofi_wavelength_extract_spectrum(in, slit_width)) == NULL) {
        cpl_msg_error(cpl_func, "Cannot extract spectrum");
        cpl_bivector_delete(spec_cat);
        return NULL;
    }
 
    /* Wavelength calibration with the cross correlation */
    cpl_msg_info(cpl_func, "Apply the wavelength calibration with XC");
    cpl_msg_indent_more();
    if ((wl_poly_xc = sofi_wavelength_xc_engine(spectrum, spec_cat,
                    degree, phdisprel, wl_err, nsamples, slit_width,
                    plot, &xc)) == NULL) {
        cpl_msg_warning(cpl_func, "Cannot calibrate");
        wl_poly_xc = cpl_polynomial_duplicate(phdisprel);
        *xcorr = -1;
    } else {
        *xcorr = xc;
    }
    cpl_msg_indent_less();
 
    /* Use Point pattern matching if requested */
    if (use_ppm) {
        /* Remove the low frequency part */
        if ((filtered=cpl_vector_filter_median_create(spectrum, 
                        filt_size))==NULL) {
            cpl_polynomial_delete(wl_poly_xc);
            cpl_bivector_delete(spec_cat);
            cpl_vector_delete(spectrum);
            cpl_msg_error(cpl_func, "Cannot filter the spectrum");
            return NULL;
        }
        clean = cpl_vector_duplicate(spectrum);
        cpl_vector_subtract(clean, filtered);
        cpl_vector_delete(filtered);
 
        /* Threshold the spectrum */
        stdev = cpl_vector_get_stdev(clean);
        median = cpl_vector_get_median_const(clean);
        threshold = median + stdev;
        pclean = cpl_vector_get_data(clean);
        for (i=0; i<cpl_vector_get_size(clean); i++) {
            if (pclean[i] < threshold) pclean[i] = 0.0;
        }
        
        /* Wavelength calibration with the Point-Pattern matching */
        cpl_msg_info(cpl_func, "Apply the wavelength calibration with PPM");
        cpl_msg_indent_more();
        if ((wl_poly_ppm = irplib_ppm_engine(clean, spec_cat, wl_poly_xc, 
                        slit_width, 2, 0.2, degree, plot, NULL)) == NULL) {
            cpl_msg_warning(cpl_func, "Cannot calibrate");
            cpl_msg_indent_less();
            cpl_vector_delete(clean);
            cpl_bivector_delete(spec_cat);
            return wl_poly_xc;
        }
        cpl_polynomial_delete(wl_poly_xc);
        wl_poly_xc = wl_poly_ppm;
        cpl_msg_indent_less();
        cpl_vector_delete(clean);
    }
    cpl_bivector_delete(spec_cat);
    cpl_vector_delete(spectrum);
    
    return wl_poly_xc;
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
double sofi_get_slitwidth(const char * filename)
{
    const char          *   sval;
    double                  slit_width;
    double                  pscale;
    cpl_propertylist    *   plist;
 
    /* Initialize    */
    plist=cpl_propertylist_load(filename, 0);
 
    /* Get the slit name used */
    if ((sval = sofi_pfits_get_opti1_id(plist)) == NULL) {
        cpl_msg_error(cpl_func, "cannot get slit used");
        cpl_propertylist_delete(plist);
        return -1.0;
    }
 
    /* Get the slit width in arcseconds */
    if (!strcmp(sval, "long_slit_1"))           slit_width = 1;
    else if (!strcmp(sval, "long_slit_2"))      slit_width = 2;
    else if (!strcmp(sval, "long_slit_0.6"))    slit_width = 0.6;
    else {
        cpl_msg_error(cpl_func, "unrecognized slit");
        cpl_propertylist_delete(plist);
        return -1.0;
    }
 
    /* Get the pixelscale and convert arsec -> pixels */
    pscale = sofi_pfits_get_pixscale(plist);
    if (cpl_error_get_code()) {
        cpl_msg_error(cpl_func, "cannot get pixscale");
        cpl_propertylist_delete(plist);
        return -1.0;
    }
    cpl_propertylist_delete(plist);
 
    slit_width /= pscale;
 
    /* Return  */
    return slit_width;
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
cpl_polynomial * sofi_get_disprel_estimate(const char * filename,
                                           int          poly_deg)
{
    cpl_polynomial   * self;
    cpl_propertylist * plist;
    int                npix;
    const char       * sval;
    double             wmin = FLT_MAX; /* (False) uninit warning */
    double             wmax = FLT_MAX; /* (False) uninit warning */
    cpl_size           i;
 
    /* Test entries */
    cpl_ensure(filename != NULL, CPL_ERROR_NULL_INPUT, NULL);
    cpl_ensure(poly_deg >= 0, CPL_ERROR_ILLEGAL_INPUT, NULL);
 
    plist = cpl_propertylist_load_regexp(filename, 0, "ESO INS MODE$|^NAXIS2$",
                                         0);
    cpl_ensure(plist != NULL, cpl_error_get_code(), NULL);
 
    /* Get naxis2 */
    npix = sofi_pfits_get_naxis2(plist);
 
    /* Get the grism */
    sval = sofi_pfits_get_mode(plist);
 
    if (sval == NULL) {
        (void)cpl_error_set_message(cpl_func, cpl_error_get_code(),
                                    "Could not get grism from [%s]", filename);
    } else if (!strcmp(sval, "LONG_SLIT_BLUE")) {
        wmin = 9500;
        wmax = 16400;
    } else if (!strcmp(sval, "LONG_SLIT_RED")) {
        wmin = 15300;
        wmax = 25200;
    } else if (!strcmp(sval, "LONG_SLIT_Z")) {
        wmin = 7900;
        wmax = 9500;
    } else if (!strcmp(sval, "LONG_SLIT_J")) {
        wmin = 12000;
        wmax = 14500;
    } else if (!strcmp(sval, "LONG_SLIT_H")) {
        wmin = 14900;
        wmax = 18100;
    } else if (!strcmp(sval, "LONG_SLIT_K")) {
        wmin = 19500;
        wmax = 24500;
    } else {
        (void)cpl_error_set_message(cpl_func, CPL_ERROR_UNSUPPORTED_MODE,
                                    "grism ('ESO INS MODE') in '%s' is "
                                    "unsupported: %s", filename, sval);
    }
 
    cpl_propertylist_delete(plist);
 
    if (cpl_error_get_code()) return NULL;
 
    if (npix < 2) {
        (void)cpl_error_set_message(cpl_func, CPL_ERROR_BAD_FILE_FORMAT,
                                    "In '%s' NAXIS2 must be at least 2, not %d",
                                    filename, npix);
        return NULL;
    }
 
    /* P(0.5)        = wmin */
    /* P(0.5 + npix) = wmax */
    /* P(x)          = a * x + b */
    /* a = (wmax - wmin) / npix */
    /* b = wmin - 0.5 * a = wmin - 0.5 * (wmax - wmin) / npix */
 
    self = cpl_polynomial_new(1);
    i = 1;
    cpl_polynomial_set_coeff(self, &i, (wmax - wmin) / npix);
    i = 0;
    cpl_polynomial_set_coeff(self, &i, wmin - 0.5 * (wmax - wmin) / npix);
 
    return self;
}

 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
static cpl_bivector * sofi_wavelength_load_catalog(
        const char  *   name,
        const char  *   oh,
        const char  *   xe,
        const char  *   ne)
{
    cpl_bivector        *   sig;
    double              *   psigx;
    double              *   psigy;
    double                  epsilon;
    int                     nlines;
    cpl_table           *   cat_tab1;
    cpl_table           *   cat_tab2;
    cpl_propertylist    *   reflist;
    int                     i;
 
    /* Check entries */
    if (name == NULL) return NULL;
 
    /* Initialise */
    epsilon = 0.01;
 
    /* Handle the table names and load the tables */
    if (!strcmp(name, "oh")) {
        /* Check if the table is there */
        if (oh == NULL) {
            cpl_msg_error(cpl_func, "Please provide the OH lines catalog");
            return NULL;
        }
        /* Load the table */
        if ((cat_tab1 = cpl_table_load(oh, 1, 0)) == NULL) {
            cpl_msg_error(cpl_func, "Cannot load the catalog");
            return NULL;
        }
    } else if (!strcmp(name, "Xe")) {
        /* Check if the table is there */
        if (xe == NULL) {
            cpl_msg_error(cpl_func, "Please provide the XE lines catalog");
            return NULL;
        }
        /* Load the table */
        if ((cat_tab1 = cpl_table_load(xe, 1, 0)) == NULL) {
            cpl_msg_error(cpl_func, "Cannot load the XE catalog");
            return NULL;
        }
    } else if (!strcmp(name, "Ne")) {
        /* Check if the table is there */
        if (ne == NULL) {
            cpl_msg_error(cpl_func, "Please provide the NE lines catalog");
            return NULL;
        }
        /* Load the table */
        if ((cat_tab1 = cpl_table_load(ne, 1, 0)) == NULL) {
            cpl_msg_error(cpl_func, "Cannot load the NE catalog");
            return NULL;
        }
    } else if (!strcmp(name, "Xe+Ne")) {
        /* Check if the table is there */
        if (xe == NULL) {
            cpl_msg_error(cpl_func, "Please provide the XE lines catalog");
            return NULL;
        }
        /* Load the table */
        if ((cat_tab1 = cpl_table_load(xe, 1, 0)) == NULL) {
            cpl_msg_error(cpl_func, "Cannot load the XE catalog");
            return NULL;
        }
        /* Check if the table is there */
        if (ne == NULL) {
            cpl_msg_error(cpl_func, "Please provide the NE lines catalog");
            return NULL;
        }
        /* Load the table */
        if ((cat_tab2 = cpl_table_load(ne, 1, 0)) == NULL) {
            cpl_msg_error(cpl_func, "Cannot load the NE catalog");
            cpl_table_delete(cat_tab1);
            return NULL;
        }
        /* Merge the tables */
        if (cpl_table_insert(cat_tab1, cat_tab2,
                    cpl_table_get_nrow(cat_tab1)) != CPL_ERROR_NONE) {
            cpl_msg_error(cpl_func, "Cannot merge tables");
            cpl_table_delete(cat_tab1);
            cpl_table_delete(cat_tab2);
            return NULL;
        }
        cpl_table_delete(cat_tab2);
    } else return NULL;
 
    /* Sort the table */
    reflist = cpl_propertylist_new();
    cpl_propertylist_append_bool(reflist, SOFI_COL_WAVELENGTH, 0);
    cpl_table_sort(cat_tab1, reflist);
    cpl_propertylist_delete(reflist);
 
    /* Create the bivector */
#ifdef HAVE_GSL
    nlines = cpl_table_get_nrow(cat_tab1);
    sig = cpl_bivector_new(nlines);
    psigx = cpl_bivector_get_x_data(sig);
    psigy = cpl_bivector_get_y_data(sig);
    for (i=0; i < nlines; i++) {
        const double wl
            = cpl_table_get_double(cat_tab1, SOFI_COL_WAVELENGTH, i, NULL);
        const double intens
            = cpl_table_get_double(cat_tab1, SOFI_COL_EMISSION, i, NULL);
 
        psigx[i] = wl;
        psigy[i] = intens;
    }
#else
    nlines = 3 * cpl_table_get_nrow(cat_tab1);
    sig = cpl_bivector_new(nlines);
    psigx = cpl_bivector_get_x_data(sig);
    psigy = cpl_bivector_get_y_data(sig);
    for (i=0; i<cpl_table_get_nrow(cat_tab1); i++) {
        const double
        wl = cpl_table_get_double(cat_tab1, SOFI_COL_WAVELENGTH, i, NULL);
        const double
        intens = cpl_table_get_double(cat_tab1, SOFI_COL_EMISSION, i, NULL);
        psigx[3*i]   = wl - epsilon;
        psigy[3*i]   = 0.0;
        psigx[3*i+1] = wl;
        psigy[3*i+1] = intens;
        psigx[3*i+2] = wl + epsilon;
        psigy[3*i+2] = 0.0;
    }
#endif
    cpl_table_delete(cat_tab1);
    return sig;
}
    
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
static cpl_polynomial * sofi_wavelength_xc_engine(
        const cpl_vector        *   spec,
        const cpl_bivector      *   spec_cat,
        int                         degree,
        const cpl_polynomial    *   phdisprel,
        double                      wl_err,
        int                         nsamples,
        double                      slit_width,
        int                         display,
        double                  *   xcorr)
{
#ifdef HAVE_GSL
 
    irplib_line_spectrum_model model;
    const int spec_size = cpl_vector_get_size(spec);
 
    const double phwl_min = cpl_polynomial_eval_1d(phdisprel, 0.5, NULL);
    const double phwl_max = cpl_polynomial_eval_1d(phdisprel, 0.5 + spec_size,
                                                   NULL);
    const int    emil     = irplib_bivector_count_positive(spec_cat, phwl_min,
                                                           phwl_max);
 
    const double wfwhm = SOFI_WFWHM;
    const double slitw = slit_width;
    const double xtrunc = 0.5 * slitw + 5.0 * wfwhm * CPL_MATH_SIG_FWHM;
    const int max_offset = SOFI_OFFSET;
 
    /* One practical limitation on hshiftmax is that it may not be so big,
       as to cause phdisp to be evaluated outside its range of monotony. */
    const int    hshiftmax = SOFI_MIN(spec_size/4, max_offset);
 
    /* Initialize to zero */
    cpl_vector     * linepix
        = cpl_vector_wrap(cpl_bivector_get_size(spec_cat),
                          cpl_calloc(cpl_bivector_get_size(spec_cat),
                                     sizeof(double)));
    cpl_polynomial * self;
    cpl_error_code   error;
    const int        fitdeg = emil > degree ? degree : emil-1;
    const int        maxite = fitdeg * 200;
    const int        linelim =
#ifdef SOFI_USE_LINELIM
        xtrunc*(spec_size + 2 * hshiftmax);
#else
        0;
#endif
 
    const cpl_boolean doplot /* FIXME: Function parameter */
        = display || getenv("SOFI_DOPLOT") ? CPL_TRUE : CPL_FALSE;
    double            xc;
 
 
    memset(&model, 0, sizeof(model));
    model.wslit = slitw;
    model.wfwhm = wfwhm;
    model.xtrunc = xtrunc;
    model.lines = spec_cat;
    model.linepix = linepix;
    model.cost = 0;
    model.xcost = 0;
 
    self = cpl_polynomial_duplicate(phdisprel);
 
    error = irplib_polynomial_find_1d_from_correlation_all
        (self, fitdeg, spec, 0, linelim, (irplib_base_spectrum_model*)&model,
         irplib_vector_fill_logline_spectrum, SOFI_PIXTOL, SOFI_PIXSTEP,
         hshiftmax, maxite, SOFI_MAXFAIL, SOFI_CONTINUE, doplot, &xc);
 
    cpl_vector_delete(linepix);
 
    if (!error) {
            cpl_msg_info(cpl_func, "XC: %g (cost=%u:%u. lines=%u)", xc,
                         model.cost, model.xcost, model.ulines);
            *xcorr = xc;
    } else {
        cpl_error_set_where(cpl_func);
        cpl_polynomial_delete(self);
        self = NULL;
    }
 
    return self;
 
#else
    cpl_vector      *   wl_errors;
    int                 nsamp;
    double              wl_start, wl_stop, wl_error, xc;
    cpl_bivector    *   catal_loc;
    cpl_polynomial  *   xcdisp;
    cpl_polynomial  *   xcdisp_tmp;
    cpl_table       *   spc_tab;
    
    /* Test entries */
    if (spec == NULL || spec_cat == NULL || phdisprel == NULL) return NULL;
    if (degree < 1 || degree > 3) return NULL;
 
    /* Initialise */
    wl_error = wl_err;
    nsamp = nsamples;
 
    /* Extract the interesting part of the catalog */
    wl_start = cpl_polynomial_eval_1d(phdisprel, (double)1, NULL) - wl_err;
    wl_stop = cpl_polynomial_eval_1d(phdisprel,
            (double)cpl_vector_get_size(spec), NULL) + wl_err;
    catal_loc = irplib_wlxcorr_cat_extract(spec_cat, wl_start, wl_stop);
    if (catal_loc == NULL) {
        cpl_msg_error(cpl_func, "Cannot extract sub-catalog [%g %g]",
                wl_start, wl_stop);
        return NULL;
    }
 
    /* Plot inputs */
    if (display) {
        irplib_wlxcorr_catalog_plot(catal_loc,
            cpl_polynomial_eval_1d(phdisprel, 1, NULL),
            cpl_polynomial_eval_1d(phdisprel, cpl_vector_get_size(spec),
                NULL));
        cpl_plot_vector(
    "set grid;set xlabel 'Position (Pixel)';set ylabel 'Intensity (ADU/sec)';",
    "t 'Extracted spectrum' w lines", "", spec);
    }
    
    /* Calibrate the spectrum with the catalog spectrum */
    /* DEGREE = 1 */
    wl_errors = cpl_vector_new(2);
    cpl_vector_fill(wl_errors, wl_error);
    xcdisp=irplib_wlxcorr_best_poly(spec, catal_loc, 1, phdisprel,
            wl_errors, nsamp, slit_width, 2, &xc, NULL, NULL);
 
    /* Higher degree --> less samples */
    nsamp /= 10;
    cpl_vector_delete(wl_errors);
    wl_errors = cpl_vector_new(3);
    cpl_vector_fill(wl_errors, wl_error);
    xcdisp_tmp=irplib_wlxcorr_best_poly(spec, catal_loc, 2, xcdisp,
            wl_errors, nsamp, slit_width, 2, &xc, NULL, NULL);
    cpl_polynomial_delete(xcdisp);
    xcdisp = xcdisp_tmp;
 
    /* Higher Precision */
    cpl_vector_divide_scalar(wl_errors, 10.0);
    xcdisp_tmp=irplib_wlxcorr_best_poly(spec, catal_loc, 2, xcdisp,
            wl_errors, nsamp, slit_width, 2, &xc, &spc_tab, NULL);
    cpl_polynomial_delete(xcdisp);
    cpl_bivector_delete(catal_loc);
    cpl_vector_delete(wl_errors);
    xcdisp = xcdisp_tmp;
 
    /* Plot result */
    if (display) {
        irplib_wlxcorr_plot_spc_table(spc_tab, "XC", 0, 0);
    }
    cpl_table_delete(spc_tab);
 
    *xcorr = xc;
    return xcdisp;
#endif
}
 
/*----------------------------------------------------------------------------*/
/*----------------------------------------------------------------------------*/
static cpl_vector * sofi_wavelength_extract_spectrum(
        const cpl_image     *   in,
        double                  slit_width) 
{
    cpl_image   *   tmp;
    cpl_image   *   collapsed;
    cpl_vector  *   spec_init;
    cpl_vector  *   spec_bgclean;
    cpl_vector  *   spec_bg;
    double      *   pspec_bgclean;
    cpl_vector  *   spec_bgclean_log;
    double      *   pspec_bgclean_log;
    int             i;
 
    /* Check entries */
    if (in == NULL) return NULL;
    
    /* Local copy */
    tmp = cpl_image_duplicate(in);
 
    /* Flip because the wl goes from top to bottom */
    cpl_image_flip(tmp, 0);
    
    /* Collapse */
    if ((collapsed = cpl_image_collapse_median_create(tmp, 1, 20, 20)) == NULL){
        cpl_msg_error(cpl_func, "Collapsing input image");
        cpl_image_delete(tmp);
        return NULL;
    }
    cpl_image_delete(tmp);
 
    /* Put the spectrum in a vector */
    spec_init = cpl_vector_new_from_image_column(collapsed, 1);
    cpl_image_delete(collapsed);
 
    /* Remove thermal background contributions */
    spec_bgclean = cpl_vector_duplicate(spec_init);
    /* Filter away very wide (8 * slit_width) features */
    if ((spec_bg = cpl_vector_filter_median_create(spec_init, 
                    (int)((0.5+8*slit_width)/2))) == NULL) {
        cpl_msg_error(cpl_func, "sub_lowpass failed");
        cpl_vector_delete(spec_init);
        cpl_vector_delete(spec_bgclean);
        return NULL;
    }
    cpl_vector_subtract(spec_bgclean, spec_bg);
    cpl_vector_delete(spec_bg);
    cpl_vector_delete(spec_init);
 
    /* Threshold negative intensities */
    pspec_bgclean = cpl_vector_get_data(spec_bgclean);
    for (i=0; i<cpl_vector_get_size(spec_bgclean); i++) 
        if (pspec_bgclean[i] < 0) pspec_bgclean[i] = 0.0;
 
    /* Put less weight on the intensity by taking the logarithm */
    /* Add 1 to ensure continuity around zero */
    spec_bgclean_log = cpl_vector_duplicate(spec_bgclean);
    pspec_bgclean = cpl_vector_get_data(spec_bgclean);
    pspec_bgclean_log = cpl_vector_get_data(spec_bgclean_log);
#if 1
    for (i=0; i<cpl_vector_get_size(spec_bgclean_log); i++)
      pspec_bgclean_log[i] = pspec_bgclean[i] > 0.0 ? 
          log1p(pspec_bgclean[i] ) : 0.0;
#endif
 
    cpl_vector_delete(spec_bgclean);
    return spec_bgclean_log;
}