/* $Id: cpl_wlcalib_xc.c,v 1.7 2008/02/01 14:59:00 yjung Exp $ * * This file is part of the ESO Common Pipeline Library * Copyright (C) 2001-2004 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ /* * $Author: yjung $ * $Date: 2008/02/01 14:59:00 $ * $Revision: 1.7 $ * $Name: $ */ #ifdef HAVE_CONFIG_H #include #endif /*----------------------------------------------------------------------------- Includes -----------------------------------------------------------------------------*/ #include #include #include #include #include "cpl_math_const.h" #include "cpl_wlcalib_xc.h" /*----------------------------------------------------------------------------*/ /** * @defgroup cpl_wlcalib_xc Wavelength Cross correlation * * @par Synopsis: * @code * #include "cpl_wlcalib_xc.h" * @endcode */ /*----------------------------------------------------------------------------*/ /**@{*/ static cpl_bivector * cpl_wlcalib_xc_gen_signal(const cpl_bivector *, double, double, const cpl_polynomial *, int, int, int *) ; static int cpl_wlcalib_xc_signal_resample(cpl_vector *, const cpl_vector *, const cpl_bivector *) ; /*----------------------------------------------------------------------------*/ /** @brief Find the best polynomial in a given range @param spectrum The spectrum vector @param lines_catalog The lines catalog @param degree The polynomial degree @param guess_poly Guess Dispersion Polynomial @param wl_error Search range around the anchor points @param nsamples Number of samples around the anchor points @param slitw The slit width @param fwhm The spectral FWHM [pixel] @param xc Cross-correlation factor (returned) @param wlres The table with the calibration results or NULL @param xcorrs The vector with the correlation values or NULL @return the best polynomial or NULL in error case wl_error must be of size degree+1. The returned polynomial must be deallocated with cpl_polynomial_delete(). */ /*----------------------------------------------------------------------------*/ cpl_polynomial * cpl_wlcalib_xc_best_poly( const cpl_vector * spectrum, const cpl_bivector * lines_catalog, int degree, const cpl_polynomial * guess_poly, const cpl_vector * wl_error, int nsamples, double slitw, double fwhm, double * xc, cpl_table ** wlres, cpl_vector ** xcorrs) { int ntests ; int spec_sz ; cpl_polynomial ** candidates ; double xpos, wl_pos ; cpl_vector * init_pts_wl ; cpl_vector * init_pts_x ; cpl_vector * pts_wl ; cpl_vector * vxcorrs ; int best_ind ; cpl_polynomial * poly_sol ; cpl_table * spc_table ; double xc_cur ; cpl_vector * vxc ; cpl_bivector * gen_init ; const double * pwl_error ; int i, j, k, l, m ; /* Test entries */ if (!spectrum || !lines_catalog || !guess_poly || !xc || !wl_error) return NULL; if (degree < 1 || degree > 4) return NULL ; if (cpl_vector_get_size(wl_error) != degree+1) return NULL ; /* Initialise */ ntests = (int)pow(nsamples, (degree + 1)) ; spec_sz = cpl_vector_get_size(spectrum) ; pwl_error = cpl_vector_get_data_const(wl_error) ; if (wlres != NULL) *wlres = NULL ; if (xcorrs != NULL) *xcorrs = NULL ; /* Create initial test points */ init_pts_x = cpl_vector_new(degree + 1) ; init_pts_wl = cpl_vector_new(degree + 1) ; for (i=0 ; i *xc) { *xc = xc_cur ; best_ind = i ; } cpl_vector_set(vxcorrs, i, xc_cur) ; } else { /* FIXME: Replace with cpl_error_set_message() */ cpl_msg_error(cpl_func, "Cannot generate the signal - abort") ; cpl_vector_delete(vxcorrs) ; cpl_vector_delete(vxc) ; for (i=0 ; iwave_min)) wave_min_id = i ; if (wave_cur= wave_max_id) return NULL ; sub_cat_wl = cpl_vector_extract(cpl_bivector_get_x_const(lines_catalog), wave_min_id, wave_max_id, 1) ; sub_cat_int = cpl_vector_extract(cpl_bivector_get_y_const(lines_catalog), wave_min_id, wave_max_id, 1) ; sub_cat = cpl_bivector_wrap_vectors(sub_cat_wl, sub_cat_int) ; return sub_cat ; } /*----------------------------------------------------------------------------*/ /** @brief Create Right Half of a symmetric smoothing kernel @param slitw The slit width [pixel] @param fwhm The spectral FWHM [pixel] @return Right Half of (symmetric) smoothing vector The smoothing function is the right half of the convolution of a Gaussian with sigma = fwhm / (2 * sqrt(2*log(2))) and a top-hat with a width equal to the slit width. Since this function is symmetric only the central, maximum value and the right half is returned. The length of the resulting vector is 1 + 5 * sigma + slitw/2 */ /*----------------------------------------------------------------------------*/ cpl_vector * cpl_wlcalib_xc_convolve_create_kernel( double slitw, double fwhm) { double sigma ; int ihtophat, gausshlen, convolen ; cpl_vector * self ; cpl_vector * tophat ; cpl_image * iself ; int i ; /* Initialise */ sigma = fwhm * CPL_MATH_SIG_FWHM; ihtophat = (int)slitw/2 ; gausshlen = 1 + 5 * sigma + ihtophat ; /* convolen must be at least twice the gausshlen */ convolen = 1 + 10 * sigma + 8 * ihtophat ; /* Test entries */ if ((slitw <= 0.0) || (fwhm <= 0.0) ||(convolen < 2 * gausshlen)) return NULL ; /* Create the vectorss */ self = cpl_vector_new(gausshlen) ; tophat = cpl_vector_new(convolen) ; /* Easiest way to fill with a Gaussian is via a CPL image */ iself = cpl_image_wrap_double(gausshlen, 1, cpl_vector_get_data(self)); /* Place the top point of the Gaussian on left-most pixel */ cpl_image_fill_gaussian(iself, 1.0, 1.0, CPL_MATH_SQRT2PI, sigma, 1.0); cpl_image_unwrap(iself) ; /* The number of non-zero elements is 1+2*ihtophat */ cpl_vector_fill(tophat, 0.0); for (i = convolen/2-ihtophat; i < 1+convolen/2+ihtophat; i++) cpl_vector_set(tophat, i, 1.0/(1.0+2.0*ihtophat)); /* Convolve the Top-hat with the Gaussian */ if (cpl_wlcalib_xc_convolve(tophat, self)) { cpl_msg_error(cpl_func, "Cannot convolve") ; cpl_vector_delete(self) ; cpl_vector_delete(tophat) ; return NULL ; } /* Overwrite the Gaussian with the Right Half of the convolution of the Top-hat + Gausssian */ for (i = 0 ; i < gausshlen; i++) cpl_vector_set(self, i, cpl_vector_get(tophat, i + convolen/2)); cpl_vector_delete(tophat) ; return self; } /*----------------------------------------------------------------------------*/ /** @brief Convolve a 1d-signal with a symmetric 1D-signal @param smoothed Preallocated vector to be smoothed in place @param conv_kernel Vector with symmetric convolution function @return 0 or -1 in error case The length of conv_kernel must be smaller than that of smoothed. */ /*----------------------------------------------------------------------------*/ int cpl_wlcalib_xc_convolve( cpl_vector * smoothed, const cpl_vector * conv_kernel) { int nsamples ; int ihwidth ; cpl_vector * raw ; double * psmoothe ; double * praw ; const double* psymm ; int i, j ; /* Test entries */ if ((!smoothed) || (!conv_kernel)) return -1 ; /* Initialise */ nsamples = cpl_vector_get_size(smoothed) ; ihwidth = cpl_vector_get_size(conv_kernel) - 1 ; if (ihwidth >= nsamples) return -1 ; psymm = cpl_vector_get_data_const(conv_kernel) ; psmoothe = cpl_vector_get_data(smoothed) ; /* Create raw vector */ raw = cpl_vector_duplicate(smoothed) ; praw = cpl_vector_get_data(raw) ; /* Convolve with the symmetric function */ for (i=0 ; i nsamples-1 ? nsamples - 1 : i+j; psmoothe[i] += (praw[k]+praw[i-j]) * psymm[j]; } } cpl_vector_delete(raw) ; return 0 ; } /**@}*/ /*----------------------------------------------------------------------------*/ /** @brief Construct the expected signal at the given wavelengths @param lines_catalog Bivector with the lines @param slitw The slit width @param fwhm The spectral FWHM [pixel] @param poly the polynomial for the wavelengths @param nsamples the number of samples of the original signal @param search_hs the half size of the searching zone @param nb_lines the number of lines used or NULL @return the bivector with the generated spectrum The expected emission is a model spectrum used to cross-correlate against an actual observed spectrum. Its size is nsamples + 2 * search_hs */ /*----------------------------------------------------------------------------*/ static cpl_bivector * cpl_wlcalib_xc_gen_signal( const cpl_bivector * lines_catalog, double slitw, double fwhm, const cpl_polynomial * poly, int nsamples, int search_hs, int * nb_lines) { int size, nlines ; cpl_vector * conv_kernel ; cpl_bivector * gen_spectrum ; cpl_bivector * sub_cat ; double wave_min, wave_max ; cpl_vector * wl_limits ; int ind ; double new_val, wl ; int i ; /* Initialise */ size = nsamples + 2 * search_hs ; /* Test entries */ if (!lines_catalog || !poly) return NULL ; if (size <= 1) return NULL ; /* Create the gen_spectrum X values */ gen_spectrum = cpl_bivector_new(size) ; cpl_vector_fill_polynomial(cpl_bivector_get_x(gen_spectrum), poly, -search_hs+1, 1) ; /* Extract the relevant part of the catalog */ wave_min = cpl_vector_get(cpl_bivector_get_x(gen_spectrum), 0) ; wave_max = cpl_vector_get(cpl_bivector_get_x(gen_spectrum), size-1) ; sub_cat = cpl_wlcalib_xc_cat_extract(lines_catalog, wave_min, wave_max) ; if (sub_cat == NULL) { cpl_bivector_delete(gen_spectrum) ; return NULL ; } nlines = cpl_bivector_get_size(sub_cat) ; if (nb_lines) *nb_lines = nlines ; /* Resample the spectrum */ wl_limits = cpl_vector_new(size + 1); cpl_vector_fill_polynomial(wl_limits, poly, -search_hs+0.5, 1) ; if (nlines < size) { /* A few lines in catalog */ cpl_vector_fill(cpl_bivector_get_y(gen_spectrum), 0.0) ; for (i=0 ; icpl_vector_get(wl_limits, ind) && ind pxbounds[i+1]) x = pxbounds[i+1]; /* Contribution from interpolated value at wavelength at lower pixel boundary */ presampled[i] = pybounds[i] * (x - xlow); /* Contribution from table values in between pixel boundaries */ while ((pxhires[itt] < pxbounds[i+1]) && (itt < cpl_bivector_get_size(hires))) { const double xprev = x; x = pxhires[itt+1]; if (x > pxbounds[i+1]) x = pxbounds[i+1]; presampled[i] += pyhires[itt] * (x - xlow); xlow = xprev; itt++; } /* Contribution from interpolated value at wavelength at upper pixel boundary */ presampled[i] += pybounds[i+1] * (pxbounds[i+1] - xlow); /* Compute average by dividing integral by length of pixel range (the factor 2 comes from the contributions) */ presampled[i] /= 2 * (pxbounds[i+1] - pxbounds[i]); } cpl_bivector_unwrap_vectors(boundary) ; cpl_vector_delete(ybounds) ; return 0 ; }