irplib_wlxcorr.c

00001 /* $Id: irplib_wlxcorr.c,v 1.32 2006/11/03 09:41:05 llundin Exp $
00002  *
00003  * This file is part of the IRPLIB package
00004  * Copyright (C) 2002,2003 European Southern Observatory
00005  *
00006  * This program is free software; you can redistribute it and/or modify
00007  * it under the terms of the GNU General Public License as published by
00008  * the Free Software Foundation; either version 2 of the License, or
00009  * (at your option) any later version.
00010  *
00011  * This program is distributed in the hope that it will be useful,
00012  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014  * GNU General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU General Public License
00017  * along with this program; if not, write to the Free Software
00018  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00019  */
00020 
00021 /*
00022  * $Author: llundin $
00023  * $Date: 2006/11/03 09:41:05 $
00024  * $Revision: 1.32 $
00025  * $Name:  $
00026  */
00027 
00028 #ifdef HAVE_CONFIG_H
00029 #include <config.h>
00030 #endif
00031 
00032 /*-----------------------------------------------------------------------------
00033                                 Includes
00034  -----------------------------------------------------------------------------*/
00035 
00036 #include <math.h>
00037 #include <cpl.h>
00038 
00039 #include "irplib_wlxcorr.h"
00040 #include "irplib_plot.h"
00041 
00042 /*----------------------------------------------------------------------------*/
00047 /*----------------------------------------------------------------------------*/
00048 
00049 static int irplib_wlxcorr_signal_resample(cpl_vector *, const cpl_vector *, 
00050         const cpl_bivector *) ;
00051 
00054 /*----------------------------------------------------------------------------*/
00074 /*----------------------------------------------------------------------------*/
00075 cpl_polynomial * irplib_wlxcorr_best_poly(
00076         const cpl_vector        *   spectrum,
00077         const cpl_bivector      *   lines_catalog,
00078         int                         degree,
00079         const cpl_polynomial    *   guess_poly,
00080         const cpl_vector        *   wl_error,
00081         int                         nsamples,
00082         double                      slitw,
00083         double                      fwhm,
00084         int                         doplot,
00085         double                  *   xc,
00086         cpl_table               **  wlres)
00087 {
00088     int                 ntests ;
00089     int                 spec_sz ;
00090     cpl_polynomial  **  candidates ;
00091     double              xpos, wl_pos ;
00092     cpl_vector      *   init_pts_wl ;
00093     cpl_vector      *   init_pts_x ;
00094     cpl_vector      *   pts_wl ;
00095     cpl_vector      *   vxcorrs ;
00096     int                 best_ind ;
00097     cpl_polynomial  *   poly_sol ;
00098     cpl_table       *   spc_table ;
00099     double              xc_cur ;
00100     cpl_vector      *   vxc ;
00101     cpl_bivector    *   gen_init ;
00102     double          *   pwl_error ; 
00103     int                 i, j, k, l ;
00104 
00105     /* Test entries */
00106     if (!spectrum || !lines_catalog || !guess_poly || !xc || !wl_error) 
00107         return NULL;
00108     if (degree < 1 || degree > 3) return NULL ;
00109     if (cpl_vector_get_size(wl_error) != degree+1) return NULL ;
00110 
00111     /* Initialise */
00112     ntests = (int)pow(nsamples, (degree + 1)) ;
00113     spec_sz = cpl_vector_get_size(spectrum) ;
00114     pwl_error = cpl_vector_get_data(wl_error) ;
00115    
00116     /* Plot if requested */
00117     if (doplot) {
00118         irplib_wlxcorr_catalog_plot(lines_catalog, 
00119                 cpl_polynomial_eval_1d(guess_poly, 1, NULL),
00120                 cpl_polynomial_eval_1d(guess_poly, spec_sz, NULL)) ;
00121         irplib_vector_plot(
00122     "set grid;set xlabel 'Position (Pixel)';set ylabel 'Intensity (ADU/sec)';",
00123     "t 'Extracted spectrum' w lines", "", spectrum);
00124     }
00125  
00126     /* Create initial test points */
00127     init_pts_x = cpl_vector_new(degree + 1) ;
00128     init_pts_wl = cpl_vector_new(degree + 1) ;
00129     for (i=0 ; i<degree + 1 ; i++) {
00130         xpos = i*spec_sz/degree ;
00131         cpl_vector_set(init_pts_x, i, xpos) ;
00132         cpl_vector_set(init_pts_wl, i, 
00133                 cpl_polynomial_eval_1d(guess_poly, xpos, NULL)) ;
00134     }
00135     
00136     /* Create the polynomial candidates */
00137     candidates = cpl_malloc(ntests * sizeof(cpl_polynomial*)) ;
00138     pts_wl = cpl_vector_new(degree + 1) ;
00139     for (i=0 ; i<nsamples ; i++) {
00140         wl_pos = cpl_vector_get(init_pts_wl, 0) 
00141             - pwl_error[0]/2 + i*pwl_error[0]/nsamples ;
00142         cpl_vector_set(pts_wl, 0, wl_pos) ;
00143         for (j=0 ; j<nsamples ; j++) {
00144             wl_pos = cpl_vector_get(init_pts_wl, 1)
00145                 - pwl_error[1]/2 + j*pwl_error[1]/nsamples ;
00146             cpl_vector_set(pts_wl, 1, wl_pos) ;
00147             if (degree == 1) {
00148                 candidates[j+i*nsamples] = 
00149                     cpl_polynomial_fit_1d_create(init_pts_x, pts_wl, degree, 
00150                             NULL);
00151             } else {
00152                 for (k=0 ; k<nsamples ; k++) {
00153                     wl_pos = cpl_vector_get(init_pts_wl, 2)
00154                         - pwl_error[2]/2 + k*pwl_error[2]/nsamples ;
00155                     cpl_vector_set(pts_wl, 2, wl_pos) ;
00156                     if (degree == 2) {
00157                         candidates[k+(j+i*nsamples)*nsamples] = 
00158                             cpl_polynomial_fit_1d_create(init_pts_x, 
00159                                     pts_wl, degree, NULL);
00160                     } else {
00161                         for (l=0 ; l<nsamples ; l++) {
00162                             wl_pos = cpl_vector_get(init_pts_wl, 3) 
00163                                 - pwl_error[3]/2 + l*pwl_error[3]/nsamples ;
00164                             cpl_vector_set(pts_wl, 3, wl_pos) ;
00165                             if (degree == 3) {
00166                             candidates[l+(k+(j+i*nsamples)*nsamples)*nsamples]=
00167                                     cpl_polynomial_fit_1d_create(init_pts_x, 
00168                                             pts_wl, degree, NULL);
00169                             } else {
00170                                 /* Should never come here */
00171                             }
00172                         }
00173                     }
00174                 }
00175             }
00176         }
00177     }
00178     cpl_vector_delete(pts_wl) ;
00179     cpl_vector_delete(init_pts_x) ;
00180     cpl_vector_delete(init_pts_wl) ;
00181     
00182     /* Loop on the polynomial and select the best */
00183     best_ind = 0 ;
00184     *xc = -1 ;
00185     vxc = cpl_vector_new(1) ;
00186     vxcorrs = cpl_vector_new(ntests) ;
00187     for (i=0 ; i<ntests ; i++) {
00188         /* Get the emission at those wavelengths */
00189         if ((gen_init=irplib_wlxcorr_gen_signal(lines_catalog, slitw, fwhm, 
00190                 candidates[i], spec_sz, 0, NULL)) != NULL) {
00191             /* Apply the correlation */
00192             cpl_vector_correlate(vxc, cpl_bivector_get_y(gen_init),spectrum);
00193             cpl_bivector_delete(gen_init);
00194             xc_cur = cpl_vector_get(vxc, 0);
00195             if (xc_cur > *xc) {
00196                 *xc = xc_cur ;
00197                 best_ind = i ;
00198             }
00199             cpl_vector_set(vxcorrs, i, xc_cur) ;
00200         } else {
00201             cpl_msg_error(cpl_func, "Cannot generate the signal - abort") ;
00202             cpl_vector_delete(vxcorrs) ;
00203             cpl_vector_delete(vxc) ;
00204             for (i=0 ; i<ntests ; i++) {
00205                 cpl_polynomial_delete(candidates[i]) ;
00206             }
00207             cpl_free(candidates) ;
00208             return NULL ;
00209         }
00210     }
00211     cpl_vector_delete(vxc) ;
00212 
00213     /* Plot the correlation values */
00214     if (doplot) {
00215         irplib_vector_plot("set grid;", "t 'Correlation values' w lines", "", 
00216                 vxcorrs) ;
00217     }
00218 
00219     /* FIXME */
00220     /* Compute the reliability factor */
00221     cpl_vector_sort(vxcorrs, -1) ;
00222     /* The signal should not be too flat and not too spike-like  */
00223     /* Compute peak * nsamples / flux -> and find low and high thresholds */
00224     
00225     /* Plot the correlation values */
00226     if (doplot) {
00227         irplib_vector_plot("set grid;", "t 'Sorted correlation values' w lines",
00228                 "", vxcorrs) ;
00229     }
00230     cpl_vector_delete(vxcorrs) ;
00231     
00232     /* Get the best polynomial */
00233     poly_sol = cpl_polynomial_duplicate(candidates[best_ind]) ;
00234        
00235     for (i=0 ; i<ntests ; i++) {
00236         cpl_polynomial_delete(candidates[i]) ;
00237     }
00238     cpl_free(candidates) ;
00239 
00240     /* Create the spc_table  */
00241     if ((spc_table = irplib_wlxcorr_gen_spc_table(spectrum, lines_catalog, 
00242                     slitw, fwhm, guess_poly, poly_sol)) == NULL) {
00243         cpl_msg_error(cpl_func, "Cannot generate infos table") ;
00244         cpl_polynomial_delete(poly_sol) ;
00245         return NULL ;
00246     } 
00247    
00248     /* Some plots */
00249     if (doplot > 0) {
00250         irplib_wlxcorr_plot_spc_table(spc_table, "XC") ;
00251     }
00252 
00253     if (wlres != NULL) *wlres = spc_table ;
00254     else cpl_table_delete(spc_table) ;
00255 
00256     return poly_sol ;
00257 }
00258 
00259 /*----------------------------------------------------------------------------*/
00272 /*----------------------------------------------------------------------------*/
00273 cpl_table * irplib_wlxcorr_gen_spc_table(
00274         const cpl_vector        *   spectrum,
00275         const cpl_bivector      *   lines_catalog,
00276         double                      slitw,
00277         double                      fwhm,
00278         const cpl_polynomial    *   guess_poly,
00279         const cpl_polynomial    *   corr_poly)
00280 {
00281     cpl_bivector    *   gen_init ;
00282     cpl_bivector    *   gen_corr ;
00283     int                 nsamples ;
00284     cpl_table       *   spc_table ;
00285     double          *   pgen ;
00286 
00287     /* Test inputs */
00288     if (spectrum == NULL) return NULL ;
00289     if (lines_catalog == NULL) return NULL ;
00290     if (guess_poly == NULL) return NULL ;
00291     if (corr_poly == NULL) return NULL ;
00292 
00293     /* Initialise */
00294     nsamples = cpl_vector_get_size(spectrum) ;
00295     
00296     /* Get the emission at initial wavelengths */
00297     if ((gen_init=irplib_wlxcorr_gen_signal(lines_catalog, slitw, fwhm,
00298                     guess_poly, nsamples, 0, NULL)) == NULL) {
00299         cpl_msg_error(cpl_func, "Cannot get the emission spectrum") ;
00300         return NULL ;
00301     }
00302  
00303     /* Get the emission at corrected wavelengths */
00304     if ((gen_corr=irplib_wlxcorr_gen_signal(lines_catalog, slitw, fwhm,
00305                     corr_poly, nsamples, 0, NULL)) == NULL) {
00306         cpl_msg_error(cpl_func, "Cannot get the emission spectrum") ;
00307         cpl_bivector_delete(gen_init) ;
00308         return NULL ;
00309     }
00310 
00311     /* Create the ouput table */
00312     spc_table = cpl_table_new(nsamples);
00313     cpl_table_new_column(spc_table, IRPLIB_COL_XC_WAVELENGTH, CPL_TYPE_DOUBLE);
00314     cpl_table_new_column(spc_table, IRPLIB_COL_XC_CAT_INIT, CPL_TYPE_DOUBLE);
00315     cpl_table_new_column(spc_table, IRPLIB_COL_XC_CAT_FINAL, CPL_TYPE_DOUBLE);
00316     cpl_table_new_column(spc_table, IRPLIB_COL_XC_OBS, CPL_TYPE_DOUBLE);
00317     
00318     /* Update table */
00319     pgen = cpl_bivector_get_x_data(gen_corr) ;
00320     cpl_table_copy_data_double(spc_table, IRPLIB_COL_XC_WAVELENGTH, pgen) ;
00321     pgen = cpl_bivector_get_y_data(gen_corr) ;
00322     cpl_table_copy_data_double(spc_table, IRPLIB_COL_XC_CAT_FINAL, pgen) ;
00323     pgen = cpl_vector_get_data(spectrum) ;
00324     cpl_table_copy_data_double(spc_table, IRPLIB_COL_XC_OBS, pgen) ;
00325     pgen = cpl_bivector_get_y_data(gen_init) ;
00326     cpl_table_copy_data_double(spc_table, IRPLIB_COL_XC_CAT_INIT, pgen);
00327     cpl_bivector_delete(gen_init);
00328     cpl_bivector_delete(gen_corr);
00329 
00330     return spc_table ;
00331 }
00332 
00333 /*----------------------------------------------------------------------------*/
00341 /*----------------------------------------------------------------------------*/
00342 cpl_bivector * irplib_wlxcorr_cat_extract(
00343         const cpl_bivector  *   lines_catalog,
00344         double                  wave_min,
00345         double                  wave_max)
00346 {
00347     int                 wave_min_id, wave_max_id ;
00348     double              wave_cur ;
00349     cpl_vector      *   sub_cat_wl ;
00350     cpl_vector      *   sub_cat_int ;
00351     cpl_bivector    *   sub_cat ;
00352     int                 i ;
00353 
00354     wave_min_id = wave_max_id = -1 ;
00355     for (i=0 ; i<cpl_bivector_get_size(lines_catalog) ; i++) {
00356         wave_cur = cpl_vector_get(cpl_bivector_get_x(lines_catalog), i) ;
00357         if ((wave_min_id<0) && (wave_cur>wave_min)) wave_min_id = i ; 
00358         if (wave_cur<wave_max) wave_max_id = i ;
00359     } 
00360     if (wave_min_id<0 || wave_max_id<0) {
00361         cpl_msg_error(cpl_func, "Cannot extract from the catalog spectrum") ;
00362         return NULL ;
00363     }
00364     if (wave_min_id >= wave_max_id) return NULL ;
00365     sub_cat_wl = cpl_vector_extract(cpl_bivector_get_x(lines_catalog), 
00366             wave_min_id, wave_max_id, 1) ;
00367     sub_cat_int = cpl_vector_extract(cpl_bivector_get_y(lines_catalog), 
00368             wave_min_id, wave_max_id, 1) ;
00369     sub_cat = cpl_bivector_wrap_vectors(sub_cat_wl, sub_cat_int) ;
00370  
00371     return sub_cat ;
00372 }
00373 
00374 /*----------------------------------------------------------------------------*/
00390 /*----------------------------------------------------------------------------*/
00391 cpl_bivector * irplib_wlxcorr_gen_signal(
00392         const cpl_bivector      *   lines_catalog,
00393         double                      slitw,
00394         double                      fwhm,
00395         const cpl_polynomial    *   poly,
00396         int                         nsamples,
00397         int                         search_hs,
00398         int                     *   nb_lines)
00399 {
00400     int                 size, nlines ;
00401     cpl_vector      *   conv_kernel ;
00402     cpl_bivector    *   gen_spectrum ;
00403     cpl_bivector    *   sub_cat ;
00404     double              wave_min, wave_max ;
00405     cpl_vector      *   wl_limits ;
00406     int                 ind ;
00407     double              new_val, wl ;
00408     int                 i ;
00409 
00410     /* Initialise */
00411     size = nsamples + 2 * search_hs ;
00412     
00413     /* Test entries */
00414     if (!lines_catalog || !poly) return NULL ;
00415     if (size <= 1) return NULL ;
00416    
00417     /* Create the gen_spectrum X values */
00418     gen_spectrum = cpl_bivector_new(size) ;
00419     cpl_vector_fill_polynomial(cpl_bivector_get_x(gen_spectrum), poly, 
00420             -search_hs+1, 1) ;
00421 
00422     /* Extract the relevant part of the catalog */
00423     wave_min = cpl_vector_get(cpl_bivector_get_x(gen_spectrum), 0) ;
00424     wave_max = cpl_vector_get(cpl_bivector_get_x(gen_spectrum), size-1) ;
00425     sub_cat = irplib_wlxcorr_cat_extract(lines_catalog, wave_min, wave_max) ;
00426     if (sub_cat == NULL) {
00427         cpl_bivector_delete(gen_spectrum) ;
00428         return NULL ;
00429     }
00430     nlines = cpl_bivector_get_size(sub_cat) ;
00431     if (nb_lines) *nb_lines = nlines ;
00432 
00433     /* Resample the spectrum */
00434     wl_limits = cpl_vector_new(size + 1);
00435     cpl_vector_fill_polynomial(wl_limits, poly, -search_hs+0.5, 1) ;
00436     if (nlines < size) {
00437         /* A few lines in catalog */
00438         cpl_vector_fill(cpl_bivector_get_y(gen_spectrum), 0.0) ;
00439         for (i=0 ; i<nlines ; i++) {
00440             wl = cpl_vector_get(cpl_bivector_get_x(sub_cat), i) ;
00441             ind = 0 ;
00442             while (wl>cpl_vector_get(wl_limits, ind) && ind<size-1) 
00443                 ind++ ;
00444             new_val = cpl_vector_get(cpl_bivector_get_y(gen_spectrum), ind) ;
00445             new_val += cpl_vector_get(cpl_bivector_get_y(sub_cat), i) ;
00446             cpl_vector_set(cpl_bivector_get_y(gen_spectrum), ind, new_val) ;
00447         }
00448     } else {
00449         /* High resolution catalog */
00450         if (irplib_wlxcorr_signal_resample(cpl_bivector_get_y(gen_spectrum),
00451                     wl_limits, lines_catalog)) {
00452             cpl_msg_error(cpl_func, "Cannot resample the signal") ;
00453             cpl_bivector_delete(gen_spectrum) ;
00454             cpl_vector_delete(wl_limits) ;
00455             cpl_bivector_delete(sub_cat) ;
00456             return NULL ;
00457         }
00458     }
00459     cpl_bivector_delete(sub_cat) ;
00460     cpl_vector_delete(wl_limits) ;
00461 
00462     /* Create convolution kernel */
00463     if ((conv_kernel = irplib_wlxcorr_convolve_create_kernel(slitw, 
00464                     fwhm)) == NULL) {
00465         cpl_msg_error(cpl_func, "Cannot create convolution kernel") ;
00466         cpl_bivector_delete(gen_spectrum) ;
00467         return NULL ;
00468     }
00469  
00470     /* Smooth the instrument resolution */
00471     if (irplib_wlxcorr_convolve(cpl_bivector_get_y(gen_spectrum),conv_kernel)) {
00472         cpl_msg_error(cpl_func, "Cannot smoothe the signal");
00473         cpl_bivector_delete(gen_spectrum) ;
00474         return NULL ;
00475     }
00476     cpl_vector_delete(conv_kernel) ;
00477 
00478     return gen_spectrum ;
00479 }
00480 
00481 /*----------------------------------------------------------------------------*/
00491 /*----------------------------------------------------------------------------*/
00492 int irplib_wlxcorr_plot_solution(
00493         const cpl_polynomial    *   init,
00494         const cpl_polynomial    *   comp,
00495         const cpl_polynomial    *   sol,
00496         int                         pix_start,
00497         int                         pix_stop) 
00498 {
00499     int                 nsamples, nplots ;
00500     cpl_vector      **  vectors ;
00501     cpl_bivector    *   bivector ;
00502     double              diff ;
00503     int                 i ;
00504     
00505     /* Test entries */
00506     if (init == NULL || comp == NULL) return -1 ;
00507     
00508     /* Initialise */
00509     nsamples = pix_stop - pix_start + 1 ;
00510     if (sol != NULL)    nplots = 3 ;
00511     else                nplots = 2 ;
00512     
00513     /* Create vectors */
00514     vectors = cpl_malloc((nplots+1)*sizeof(cpl_vector*)) ;
00515     for (i=0 ; i<nplots+1 ; i++) vectors[i] = cpl_vector_new(nsamples) ;
00516 
00517     /* First plot with the lambda/pixel relation */
00518     /* Fill vectors */
00519     for (i=0 ; i<nsamples ; i++) {
00520         cpl_vector_set(vectors[0], i, pix_start+i) ;
00521         cpl_vector_set(vectors[1], i, 
00522                 cpl_polynomial_eval_1d(init, (double)(pix_start+i), NULL)) ;
00523         cpl_vector_set(vectors[2], i, 
00524                 cpl_polynomial_eval_1d(comp, (double)(pix_start+i), NULL)) ;
00525         if (sol != NULL) 
00526             cpl_vector_set(vectors[3], i, 
00527                     cpl_polynomial_eval_1d(sol, (double)(pix_start+i), NULL)) ;
00528     }
00529 
00530     /* Plot */
00531     irplib_vectors_plot("set grid;set xlabel 'Position (pixels)';", 
00532         "t '1-Initial / 2-Computed / 3-Solution' w lines", 
00533         "", (const cpl_vector **)vectors, nplots+1);
00534 
00535     /* Free vectors */
00536     for (i=0 ; i<nplots+1 ; i++) cpl_vector_delete(vectors[i]) ;
00537     cpl_free(vectors) ;
00538 
00539     /* Allocate vectors */
00540     nplots -- ;
00541     vectors = cpl_malloc((nplots+1)*sizeof(cpl_vector*)) ;
00542     for (i=0 ; i<nplots+1 ; i++) vectors[i] = cpl_vector_new(nsamples) ;
00543     
00544     /* Second plot with the delta-lambda/pixel relation */
00545     /* Fill vectors */
00546     for (i=0 ; i<nsamples ; i++) {
00547         cpl_vector_set(vectors[0], i, pix_start+i) ;
00548         diff = cpl_polynomial_eval_1d(comp, (double)(pix_start+i), NULL) -
00549             cpl_polynomial_eval_1d(init, (double)(pix_start+i), NULL) ;
00550         cpl_vector_set(vectors[1], i, diff) ;
00551         if (sol != NULL) {
00552             diff = cpl_polynomial_eval_1d(sol, (double)(pix_start+i), NULL) -
00553                 cpl_polynomial_eval_1d(init, (double)(pix_start+i), NULL) ;
00554             cpl_vector_set(vectors[2], i, diff) ;
00555         }
00556     }
00557 
00558     /* Plot */
00559     if (sol == NULL) {
00560         bivector = cpl_bivector_wrap_vectors(vectors[0], vectors[1]) ;
00561         irplib_bivector_plot(
00562 "set grid;set xlabel 'Position (pixels)';set ylabel 'Wavelength difference';", 
00563             "t 'Computed-Initial wavelenth' w lines", "", bivector);
00564         cpl_bivector_unwrap_vectors(bivector) ;
00565     } else {
00566         irplib_vectors_plot("set grid;set xlabel 'Position (pixels)';", 
00567             "t '1-Computed - Initial / 2--Solution - Initial' w lines", 
00568             "", (const cpl_vector **)vectors, nplots+1);
00569     }
00570     
00571     /* Free vectors */
00572     for (i=0 ; i<nplots+1 ; i++) cpl_vector_delete(vectors[i]) ;
00573     cpl_free(vectors) ;
00574 
00575     /* Return */
00576     return 0 ;
00577 }
00578 
00579 /*----------------------------------------------------------------------------*/
00586 /*----------------------------------------------------------------------------*/
00587 int irplib_wlxcorr_plot_spc_table(
00588         const cpl_table     *   spc_table, 
00589         const char          *   title) 
00590 {
00591     char                title_loc[1024] ;
00592     cpl_vector      **  vectors ;
00593     cpl_vector      **  sub_vectors ;
00594     cpl_vector      *   tmp_vec ;
00595     int                 nsamples ;
00596     double              hsize_nm, max, mean1, mean3 ;
00597     int                 start_ind, stop_ind, nblines, hsize_pix ;
00598     int                 i, j ;
00599 
00600     /* Test entries */
00601     if (spc_table == NULL) return -1 ;
00602     
00603     /* Initialise */
00604     nsamples = cpl_table_get_nrow(spc_table) ;
00605     hsize_nm = 0.2 ;
00606     hsize_pix = 10 ;
00607     nblines = 0 ;
00608     sprintf(title_loc, 
00609         "t '%s - 1-Initial catalog/2-Corrected catalog/3-Observed' w lines",
00610         title) ;
00611     title_loc[1023] = (char)0 ;
00612     
00613     vectors = cpl_malloc(4*sizeof(cpl_vector*)) ;
00614     vectors[0] = cpl_vector_wrap(nsamples, 
00615             cpl_table_get_data_double(spc_table, IRPLIB_COL_XC_WAVELENGTH)) ;
00616     vectors[1] = cpl_vector_wrap(nsamples, 
00617             cpl_table_get_data_double(spc_table, IRPLIB_COL_XC_CAT_INIT)) ;
00618     vectors[2] = cpl_vector_wrap(nsamples, 
00619             cpl_table_get_data_double(spc_table, IRPLIB_COL_XC_CAT_FINAL)) ;
00620     vectors[3] = cpl_vector_wrap(nsamples, 
00621             cpl_table_get_data_double(spc_table, IRPLIB_COL_XC_OBS)) ;
00622 
00623     /* Scale the signal for a bettre display */
00624     mean1 = cpl_vector_get_mean(vectors[1]) ;
00625     mean3 = cpl_vector_get_mean(vectors[3]) ;
00626     if (fabs(mean3) > 1)
00627         cpl_vector_multiply_scalar(vectors[3], fabs(mean1/mean3)) ;
00628 
00629     irplib_vectors_plot("set grid;set xlabel 'Wavelength (nm)';", title_loc,
00630         "", (const cpl_vector **)vectors, 4);
00631 
00632     /* Unscale the signal */
00633     if (fabs(mean3) > 1)
00634         cpl_vector_multiply_scalar(vectors[3], mean3/mean1) ;
00635 
00636     /* Loop on the brightest lines and zoom on them */
00637     sprintf(title_loc, 
00638 "t '%s - 1-Initial catalog/2-Corrected catalog/3-Observed (ZOOMED)' w lines",
00639         title) ;
00640     title_loc[1023] = (char)0 ;
00641     tmp_vec = cpl_vector_duplicate(vectors[2]) ;
00642     for (i=0 ; i<nblines ; i++) {
00643         /* Find the brightest line */
00644         if ((max = cpl_vector_get_max(tmp_vec)) <= 0.0) break ;
00645         for (j=0 ; i<nsamples ; j++) {
00646             if (cpl_vector_get(tmp_vec, j) == max) break ;
00647         }
00648         if (j-hsize_pix < 0) start_ind = 0 ;
00649         else start_ind = j-hsize_pix ;
00650         if (j+hsize_pix > nsamples-1) stop_ind = nsamples-1 ;
00651         else stop_ind = j+hsize_pix ;
00652         for (j=start_ind ; j<=stop_ind ; j++) cpl_vector_set(tmp_vec, j, 0.0) ;
00653 
00654         sub_vectors = cpl_malloc(4*sizeof(cpl_vector*)) ;
00655         sub_vectors[0]=cpl_vector_extract(vectors[0],start_ind,stop_ind,1);
00656         sub_vectors[1]=cpl_vector_extract(vectors[1],start_ind,stop_ind,1);
00657         sub_vectors[2]=cpl_vector_extract(vectors[2],start_ind,stop_ind,1);
00658         sub_vectors[3]=cpl_vector_extract(vectors[3],start_ind,stop_ind,1);
00659 
00660         irplib_vectors_plot("set grid;set xlabel 'Wavelength (nm)';", title_loc,
00661             "", (const cpl_vector **)sub_vectors, 4);
00662 
00663         cpl_vector_delete(sub_vectors[0]) ;
00664         cpl_vector_delete(sub_vectors[1]) ;
00665         cpl_vector_delete(sub_vectors[2]) ;
00666         cpl_vector_delete(sub_vectors[3]) ;
00667         cpl_free(sub_vectors) ;
00668     }
00669     cpl_vector_delete(tmp_vec) ;
00670     
00671     cpl_vector_unwrap(vectors[0]) ;
00672     cpl_vector_unwrap(vectors[1]) ;
00673     cpl_vector_unwrap(vectors[2]) ;
00674     cpl_vector_unwrap(vectors[3]) ;
00675     cpl_free(vectors) ;
00676 
00677     return 0 ;
00678 }
00679 
00680 /*----------------------------------------------------------------------------*/
00688 /*----------------------------------------------------------------------------*/
00689 int irplib_wlxcorr_catalog_plot(
00690         const cpl_bivector      *   cat,
00691         double                      wmin,
00692         double                      wmax) 
00693 {
00694     int                 start, stop ;
00695     cpl_bivector    *   subcat ;
00696     cpl_vector      *   subcat_x ;
00697     cpl_vector      *   subcat_y ;
00698     double          *   pwave ;
00699     int                 nvals, nvals_tot ;
00700     int                 i ;
00701 
00702     /* Test entries */
00703     if (cat == NULL) return -1 ;
00704     if (wmax <= wmin) return -1 ;
00705 
00706     /* Initialise */
00707     nvals_tot = cpl_bivector_get_size(cat) ;
00708 
00709     /* Count the nb of values */
00710     pwave = cpl_bivector_get_x_data(cat) ;
00711     if (pwave[0] >= wmin) start = 0 ;
00712     else start = -1 ;
00713     if (pwave[nvals_tot-1] <= wmax) stop = nvals_tot-1 ;
00714     else stop = -1 ;
00715     i=0 ;
00716     while ((pwave[i] < wmin) && (i<nvals_tot-1)) i++ ;
00717     start = i ;
00718     i= nvals_tot-1 ;
00719     while ((pwave[i] > wmax) && (i>0)) i-- ;
00720     stop = i ;
00721 
00722     if (start>=stop) {
00723         cpl_msg_error(cpl_func, "Cannot plot the catalog") ;
00724         return -1 ;
00725     }
00726     nvals = start - stop + 1 ;
00727 
00728     /* Create the bivector to plot */
00729     subcat_x = cpl_vector_extract(cpl_bivector_get_x(cat), start, stop, 1) ;
00730     subcat_y = cpl_vector_extract(cpl_bivector_get_y(cat), start, stop, 1) ;
00731     subcat = cpl_bivector_wrap_vectors(subcat_x, subcat_y) ;
00732 
00733     /* Plot */
00734     irplib_bivector_plot(
00735             "set grid;set xlabel 'Wavelength (nm)';set ylabel 'Emission';",
00736             "t 'Catalog Spectrum' w lines", "", subcat);
00737     cpl_bivector_unwrap_vectors(subcat) ;
00738     cpl_vector_delete(subcat_x) ;
00739     cpl_vector_delete(subcat_y) ;
00740 
00741     return 0 ;
00742 }
00743    
00744 /*----------------------------------------------------------------------------*/
00758 /*----------------------------------------------------------------------------*/
00759 cpl_vector * irplib_wlxcorr_convolve_create_kernel(
00760         double  slitw,
00761         double  fwhm)
00762 {
00763     double          sigma ;
00764     int             ihtophat, gausshlen, convolen ;
00765     cpl_vector  *   self ;
00766     cpl_vector  *   tophat ;
00767     cpl_image   *   iself ;
00768     int             i ;
00769     
00770     /* Initialise */
00771     sigma  = fwhm / (2.0 * sqrt(2.0*log(2.0)));
00772     ihtophat  = (int)slitw/2 ;
00773     gausshlen = 1 + 5 * sigma + ihtophat ;
00774 
00775     /* convolen must be at least twice the gausshlen */
00776     convolen  = 1 + 10 * sigma + 8 * ihtophat ;
00777 
00778     /* Test entries */
00779     if ((slitw <= 0.0) || (fwhm  <= 0.0) ||(convolen < 2 * gausshlen)) 
00780         return NULL ;
00781     
00782     /* Create the vectorss */
00783     self = cpl_vector_new(gausshlen) ;
00784     tophat = cpl_vector_new(convolen) ;
00785     
00786     /* Easiest way to fill with a Gaussian is via a CPL image */
00787     iself = cpl_image_wrap_double(gausshlen, 1, cpl_vector_get_data(self));
00788 
00789     /* Place the top point of the Gaussian on left-most pixel */
00790     cpl_image_fill_gaussian(iself, 1.0, 1.0, sqrt(2.0*atan(1.0)*4.0),sigma,1.0);
00791     cpl_image_unwrap(iself) ;
00792 
00793     /* The number of non-zero elements is 1+2*ihtophat */
00794     cpl_vector_fill(tophat, 0.0);
00795 
00796     for (i = convolen/2-ihtophat; i < 1+convolen/2+ihtophat; i++)
00797         cpl_vector_set(tophat, i, 1.0/(1.0+2.0*ihtophat));
00798 
00799     /* Convolve the Top-hat with the Gaussian */
00800     if (irplib_wlxcorr_convolve(tophat, self)) {
00801         cpl_msg_error(cpl_func, "Cannot convolve") ;
00802         cpl_vector_delete(self) ;
00803         cpl_vector_delete(tophat) ;
00804         return NULL ;
00805     }
00806 
00807     /* Overwrite the Gaussian with the Right Half of the convolution of the
00808        Top-hat + Gausssian */
00809     for (i = 0 ; i < gausshlen; i++)
00810         cpl_vector_set(self, i, cpl_vector_get(tophat, i + convolen/2));
00811     cpl_vector_delete(tophat) ;
00812 
00813     return self;
00814 }
00815 
00816 /*----------------------------------------------------------------------------*/
00825 /*----------------------------------------------------------------------------*/
00826 int irplib_wlxcorr_convolve(
00827         cpl_vector          *   smoothed,
00828         const cpl_vector    *   conv_kernel)
00829 {
00830     int             nsamples ;
00831     int             ihwidth ;
00832     cpl_vector  *   raw ;
00833     double      *   psmoothe ;
00834     double      *   praw ;
00835     double      *   psymm ;
00836     int             i, j ;
00837 
00838     /* Test entries */
00839     if ((!smoothed) || (!conv_kernel)) return -1 ;
00840 
00841     /* Initialise */
00842     nsamples = cpl_vector_get_size(smoothed) ;
00843     ihwidth = cpl_vector_get_size(conv_kernel) - 1 ;
00844     if (ihwidth >= nsamples) return -1 ;
00845     psymm = cpl_vector_get_data(conv_kernel) ;
00846     psmoothe = cpl_vector_get_data(smoothed) ;
00847     
00848     /* Create raw vector */
00849     raw = cpl_vector_duplicate(smoothed) ;
00850     praw = cpl_vector_get_data(raw) ;
00851 
00852     /* Convolve with the symmetric function */
00853     for (i=0 ; i<ihwidth ; i++) {
00854         psmoothe[i] = praw[i] * psymm[0];
00855         for (j=1 ; j <= ihwidth ; j++) {
00856             const int k = i-j < 0 ? 0 : i-j;
00857             psmoothe[i] += (praw[k]+praw[i+j]) * psymm[j];
00858         }
00859     }
00860 
00861     for (i=ihwidth ; i<nsamples-ihwidth ; i++) {
00862         psmoothe[i] = praw[i] * psymm[0];
00863         for (j=1 ; j<=ihwidth ; j++)
00864             psmoothe[i] += (praw[i-j]+praw[i+j]) * psymm[j];
00865     }
00866     for (i=nsamples-ihwidth ; i<nsamples ; i++) {
00867         psmoothe[i] = praw[i] * psymm[0];
00868         for (j=1 ; j<=ihwidth ; j++) {
00869             const int k = i+j > nsamples-1 ? nsamples - 1 : i+j;
00870             psmoothe[i] += (praw[k]+praw[i-j]) * psymm[j];
00871         }
00872     }
00873     cpl_vector_delete(raw) ;
00874     return 0 ;
00875 }
00876 
00879 /*----------------------------------------------------------------------------*/
00889 /*----------------------------------------------------------------------------*/
00890 static int irplib_wlxcorr_signal_resample(
00891         cpl_vector          *   resampled, 
00892         const cpl_vector    *   xbounds,
00893         const cpl_bivector  *   hires)
00894 {
00895     cpl_vector      *   xhires ;
00896     cpl_vector      *   yhires ;
00897     double          *   pxhires ;
00898     double          *   pyhires ;
00899     double          *   pxbounds ;
00900     cpl_vector      *   ybounds ;
00901     cpl_bivector    *   boundary ;
00902     double          *   pybounds ;
00903     double          *   presampled ;
00904     int                 nsamples ;
00905     int                 i, itt ;
00906    
00907     /* Test entries */
00908     if ((!resampled) || (!xbounds) || (!hires)) return -1 ;
00909 
00910     /* Initialise */
00911     nsamples = cpl_vector_get_size(resampled) ;
00912 
00913     /* Initialise */
00914     presampled = cpl_vector_get_data(resampled) ;
00915     pxbounds = cpl_vector_get_data(xbounds) ;
00916     xhires = cpl_bivector_get_x(hires) ;
00917     yhires = cpl_bivector_get_y(hires) ;
00918     pxhires = cpl_vector_get_data(xhires) ;
00919     pyhires = cpl_vector_get_data(yhires) ;
00920     
00921     /* Create a new vector */
00922     ybounds = cpl_vector_new(cpl_vector_get_size(xbounds)) ;
00923     boundary = cpl_bivector_wrap_vectors((cpl_vector*)xbounds,ybounds) ;
00924     pybounds = cpl_vector_get_data(ybounds) ;
00925 
00926     /* Test entries */
00927     if (cpl_bivector_get_size(boundary) != nsamples + 1) {
00928         cpl_bivector_unwrap_vectors(boundary) ;
00929         cpl_vector_delete(ybounds) ;
00930         return -1 ;
00931     }
00932 
00933     /* Get the ind  */
00934     itt = cpl_vector_find(xhires, pxbounds[0]);
00935 
00936     /* Interpolate the signal */
00937     if (cpl_bivector_interpolate_linear(boundary, hires)) {
00938         cpl_msg_error(cpl_func, "Cannot interpolate the signal") ;
00939         cpl_bivector_unwrap_vectors(boundary) ;
00940         cpl_vector_delete(ybounds) ;
00941         return -1 ;
00942     }
00943 
00944     /* At this point itt most likely points to element just below
00945        pxbounds[0] */
00946     while (pxhires[itt] < pxbounds[0]) itt++;
00947 
00948     for (i=0; i < nsamples; i++) {
00949         /* The i'th signal is the weighted average of the two interpolated
00950            signals at the pixel boundaries and those table signals in
00951            between */
00952 
00953         double xlow  = pxbounds[i];
00954         double x     = pxhires[itt];
00955 
00956         if (x > pxbounds[i+1]) x = pxbounds[i+1];
00957         /* Contribution from interpolated value at wavelength at lower pixel
00958            boundary */
00959         presampled[i] = pybounds[i] * (x - xlow);
00960 
00961         /* Contribution from table values in between pixel boundaries */
00962         while ((pxhires[itt] < pxbounds[i+1]) && 
00963                 (itt < cpl_bivector_get_size(hires))) {
00964             const double xprev = x;
00965             x = pxhires[itt+1];
00966             if (x > pxbounds[i+1]) x = pxbounds[i+1];
00967             presampled[i] += pyhires[itt] * (x - xlow);
00968             xlow = xprev;
00969             itt++;
00970         }
00971 
00972         /* Contribution from interpolated value at wavelength at upper pixel
00973            boundary */
00974         presampled[i] += pybounds[i+1] * (pxbounds[i+1] - xlow);
00975 
00976         /* Compute average by dividing integral by length of pixel range
00977            (the factor 2 comes from the contributions) */
00978         presampled[i] /= 2 * (pxbounds[i+1] - pxbounds[i]);
00979     }
00980     cpl_bivector_unwrap_vectors(boundary) ;
00981     cpl_vector_delete(ybounds) ;
00982     return 0 ;
00983 }
00984 

Generated on Wed Jan 17 08:33:41 2007 for SINFONI Pipeline Reference Manual by  doxygen 1.4.4