irplib_ppm.c

00001 /* $Id: irplib_ppm.c,v 1.11 2006/09/21 15:06:26 cizzo 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., 51 Franklin St, Fifth Floor, Boston, MA  02111-1307  USA
00019  */
00020 
00021 /*
00022  * $Author: cizzo $
00023  * $Date: 2006/09/21 15:06:26 $
00024  * $Revision: 1.11 $
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_ppm.h"
00040 #include "irplib_wlxcorr.h"
00041 #include "irplib_plot.h"
00042 
00043 /*-----------------------------------------------------------------------------
00044                                    Includes
00045  -----------------------------------------------------------------------------*/
00046 
00047 static cpl_vector * irplib_ppm_convolve_line(const cpl_vector *, double,double);
00048 static cpl_vector * irplib_ppm_detect_lines(const cpl_vector *, double) ;
00049 
00050 /*----------------------------------------------------------------------------*/
00054 /*----------------------------------------------------------------------------*/
00055 
00058 /*----------------------------------------------------------------------------*/
00075 /*----------------------------------------------------------------------------*/
00076 cpl_polynomial * irplib_ppm_engine(
00077         const cpl_vector        *   spectrum,
00078         const cpl_bivector      *   lines_catalog,
00079         const cpl_polynomial    *   poly_init,
00080         double                      slitw,
00081         double                      fwhm,
00082         double                      thresh,
00083         int                         degree,
00084         int                         doplot,
00085         cpl_table               **  tab_infos)
00086 {
00087     cpl_vector      *   spec_conv ;
00088     int                 spec_sz ;
00089     cpl_vector      *   det_lines ;
00090     cpl_vector      *   cat_lines ;
00091     double          *   pcat_lines ;
00092     double              wmin, wmax ;
00093     double              disp_min, disp_max, disp ;
00094     int                 nlines_cat, nlines ;
00095     double          *   plines_catalog_x ;
00096     double          *   plines_catalog_y ;
00097     cpl_bivector    *   matched ;
00098     double          *   pmatched ;
00099     cpl_polynomial  *   fitted ;
00100     cpl_table       *   spc_table ;
00101     cpl_vector      **  vectors_plot ;
00102     cpl_vector      *   plot_cat_x ;
00103     cpl_vector      *   plot_cat_y ;
00104     cpl_vector      *   plot_y ;
00105     int                 wl_ind, start_ind, stop_ind ;
00106     int                 i ;
00107 
00108     /* Check entries */
00109     if (spectrum == NULL) return NULL ;
00110     if (lines_catalog == NULL) return NULL ;
00111     if (poly_init == NULL) return NULL ;
00112 
00113     /* Initialise */
00114     spec_sz = cpl_vector_get_size(spectrum) ;
00115     
00116     /* Correlate the spectrum with the line profile */
00117     if ((spec_conv = irplib_ppm_convolve_line(spectrum, slitw, fwhm)) == NULL) {
00118         cpl_msg_error(cpl_func, "Cannot convolve the signal") ;
00119         return NULL ;
00120     }
00121 
00122     /* Plot if requested */
00123     if (doplot) {
00124         irplib_vector_plot(
00125     "set grid;set xlabel 'Position (Pixel)';set ylabel 'Intensity (ADU/sec)';",
00126     "t 'Spectrum' w lines", "", spectrum) ;
00127         irplib_vector_plot(
00128     "set grid;set xlabel 'Position (Pixel)';set ylabel 'Correlation';",
00129     "t 'Spectrum convolved with the line profile' w lines", "", spec_conv) ;
00130     }
00131     
00132     /* Apply the lines detection */
00133     if ((det_lines = irplib_ppm_detect_lines(spec_conv, thresh)) == NULL) {
00134         cpl_msg_error(cpl_func, "Cannot detect lines") ;
00135         cpl_vector_delete(spec_conv) ;
00136         return NULL ;
00137     }
00138     cpl_vector_delete(spec_conv) ;
00139 
00140     /* Get the catalog lines */
00141     wmin = cpl_polynomial_eval_1d(poly_init, 1.0, NULL) ;
00142     wmax = cpl_polynomial_eval_1d(poly_init, spec_sz, NULL) ;
00143     plines_catalog_x = cpl_bivector_get_x_data(lines_catalog) ;
00144     plines_catalog_y = cpl_bivector_get_y_data(lines_catalog) ;
00145     nlines = cpl_bivector_get_size(lines_catalog) ;
00146     nlines_cat = 0 ;
00147     start_ind = stop_ind = -1 ;
00148     for (i=0 ; i<nlines ; i++) {
00149         if (plines_catalog_x[i] > wmin && plines_catalog_x[i] < wmax &&
00150                 plines_catalog_y[i] > 0.0) {
00151             nlines_cat++ ;
00152             if (start_ind<0) start_ind = i ;
00153             stop_ind = i ;
00154         }
00155     }
00156     if (nlines_cat == 0) {
00157         cpl_msg_error(cpl_func, "No lines in catalog") ;
00158         cpl_vector_delete(det_lines) ;
00159         return NULL ;
00160     }
00161     cat_lines = cpl_vector_new(nlines_cat) ;
00162     pcat_lines = cpl_vector_get_data(cat_lines) ;
00163     nlines_cat = 0 ;
00164     for (i=0 ; i<nlines ; i++) {
00165         if (plines_catalog_x[i] > wmin && plines_catalog_x[i] < wmax &&
00166                 plines_catalog_y[i] > 0.0) {
00167             pcat_lines[nlines_cat] = plines_catalog_x[i] ; 
00168             nlines_cat++ ;
00169         }
00170     }
00171     
00172     /* Apply the point pattern matching */
00173     disp = (wmax-wmin) / spec_sz ;
00174     disp_min = disp - (disp/10) ;
00175     disp_max = disp + (disp/10) ;
00176     if ((matched = irplib_ppm_match_positions(det_lines, cat_lines,
00177                     disp_min, disp_max, 0.05)) == NULL) {
00178         cpl_msg_error(cpl_func, "Cannot apply the point pattern matching") ;
00179         cpl_vector_delete(det_lines) ;
00180         cpl_vector_delete(cat_lines) ;
00181         return NULL ;
00182     }
00183     cpl_vector_delete(det_lines) ;
00184     cpl_vector_delete(cat_lines) ;
00185     if (cpl_bivector_get_size(matched) <= degree) {
00186         cpl_msg_error(cpl_func, "Not enough match for the fit") ;
00187         cpl_bivector_delete(matched) ;
00188         return NULL ;
00189     }
00190     cpl_msg_info(cpl_func, "Matched %d lines", cpl_bivector_get_size(matched)) ;
00191     
00192     /* Plot if requested */
00193     if (doplot) {
00194         /* Spectrum with matched lines */
00195         plot_y = cpl_vector_duplicate(spectrum) ;
00196         cpl_vector_fill(plot_y, 0.0) ;
00197         pmatched = cpl_bivector_get_x_data(matched) ;
00198         for (i=0 ; i<cpl_bivector_get_size(matched) ; i++) {
00199             cpl_vector_set(plot_y, (int)pmatched[i], 10.0) ;
00200         }
00201         vectors_plot = cpl_malloc(3*sizeof(cpl_vector*)) ;
00202         vectors_plot[0] = NULL ;
00203         vectors_plot[1] = (cpl_vector*)spectrum ;
00204         vectors_plot[2] = plot_y ;
00205         irplib_vectors_plot(
00206     "set grid;set xlabel 'Position (Pixel)';set ylabel 'Intensity (ADU/sec)';",
00207     "t 'Spectrum with matched lines' w lines", "", 
00208     (const cpl_vector **)vectors_plot, 3) ;
00209         cpl_vector_delete(plot_y) ;
00210 
00211         /* Catalog with matched lines */
00212         plot_cat_x=cpl_vector_extract(cpl_bivector_get_x(lines_catalog), 
00213                 start_ind, stop_ind, 1) ;
00214         plot_cat_y=cpl_vector_extract(cpl_bivector_get_y(lines_catalog), 
00215                 start_ind, stop_ind, 1) ;
00216         vectors_plot[0] = plot_cat_x ;
00217         vectors_plot[1] = plot_cat_y ;
00218         plot_y = cpl_vector_duplicate(vectors_plot[1]) ;
00219         cpl_vector_fill(plot_y, 0.0) ;
00220         pmatched = cpl_bivector_get_y_data(matched) ;
00221         for (i=0 ; i<cpl_bivector_get_size(matched) ; i++) {
00222             wl_ind = 0 ;
00223             while (pmatched[i] > cpl_vector_get(vectors_plot[0], wl_ind) 
00224                     && wl_ind < spec_sz) wl_ind++ ;
00225             if (wl_ind < spec_sz) cpl_vector_set(plot_y, wl_ind, 100.0) ;
00226         }
00227         vectors_plot[2] = plot_y ;
00228         irplib_vectors_plot(
00229     "set grid;set xlabel 'Position (Pixel)';set ylabel 'Intensity (ADU/sec)';",
00230     "t 'Catalog with matched lines' w lines", "", 
00231     (const cpl_vector **)vectors_plot, 3) ;
00232         cpl_vector_delete(plot_cat_x) ;
00233         cpl_vector_delete(plot_cat_y) ;
00234         cpl_vector_delete(plot_y) ;
00235         cpl_free(vectors_plot) ;
00236     }
00237     
00238     /* Apply the fit */
00239     if ((fitted=cpl_polynomial_fit_1d_create(
00240                     cpl_bivector_get_x(matched),
00241                     cpl_bivector_get_y(matched),
00242                     degree, NULL)) == NULL) {
00243         cpl_msg_error(cpl_func, "Cannot fit the polynomial") ;
00244         cpl_bivector_delete(matched) ;
00245         return NULL ;
00246     }
00247     cpl_bivector_delete(matched) ;
00248    
00249     /* Create the infos table */
00250     if ((spc_table = irplib_wlxcorr_gen_spc_table(spectrum,
00251                     lines_catalog, slitw, fwhm, poly_init, fitted)) == NULL) {
00252         cpl_msg_error(cpl_func, "Cannot generate the infos table") ;
00253         cpl_polynomial_delete(fitted) ;
00254         return NULL ;
00255     }
00256     
00257     /* Plot the result */
00258     if (doplot) {
00259         irplib_wlxcorr_plot_spc_table(spc_table, "PPM") ;
00260     }
00261 
00262     if (tab_infos != NULL) *tab_infos = spc_table ;
00263     else cpl_table_delete(spc_table) ;
00264     return fitted ;
00265 }
00266 
00267 /*----------------------------------------------------------------------------*/
00320 /*----------------------------------------------------------------------------*/
00321 cpl_bivector *irplib_ppm_match_positions(cpl_vector *peaks, cpl_vector *lines,
00322                                  double min_disp, double max_disp,
00323                                  double tolerance)
00324 {
00325   int      i, j, k, l;
00326   int      nlint, npint;
00327   int      minpos;
00328   float    min;
00329   double   lratio, pratio;
00330   double   lo_start, lo_end, hi_start, hi_end, denom;
00331   double   variation, prev_variation;
00332   double   disp = 0;
00333   int      max, maxpos, minl, mink;
00334   int      ambiguous;
00335   int      npeaks_lo, npeaks_hi;
00336   int     *peak_lo;
00337   int     *peak_hi;
00338   int    **ident;
00339   int     *nident;
00340   int     *lident;
00341 
00342   double  *peak;
00343   double  *line;
00344   int      npeaks, nlines;
00345 
00346   double  *xpos;
00347   double  *lambda;
00348   int     *ilambda;
00349   double  *tmp_xpos;
00350   double  *tmp_lambda;
00351   int     *tmp_ilambda;
00352   int     *flag;
00353   int      n = 0;
00354   int      nn;
00355   int      nseq = 0;
00356   int      gap;
00357   int     *seq_length;
00358   int      found;
00359 
00360   peak        = cpl_vector_get_data(peaks);
00361   npeaks      = cpl_vector_get_size(peaks);
00362   line        = cpl_vector_get_data(lines);
00363   nlines      = cpl_vector_get_size(lines);
00364 
00365   peak_lo     = cpl_malloc(npeaks * sizeof(int));
00366   peak_hi     = cpl_malloc(npeaks * sizeof(int));
00367   nident      = cpl_calloc(npeaks, sizeof(int));
00368   lident      = cpl_calloc(nlines, sizeof(int));
00369   xpos        = cpl_calloc(npeaks, sizeof(double));
00370   lambda      = cpl_calloc(npeaks, sizeof(double));
00371   ilambda     = cpl_calloc(npeaks, sizeof(int));
00372   tmp_xpos    = cpl_calloc(npeaks, sizeof(double));
00373   tmp_lambda  = cpl_calloc(npeaks, sizeof(double));
00374   tmp_ilambda = cpl_calloc(npeaks, sizeof(int));
00375   flag        = cpl_calloc(npeaks, sizeof(int));
00376   seq_length  = cpl_calloc(npeaks, sizeof(int));
00377   ident       = cpl_malloc(npeaks * sizeof(int *));
00378   for (i = 0; i < npeaks; i++)
00379     ident[i]  = cpl_malloc(3 * npeaks * sizeof(int));
00380 
00381   /*
00382    * This is just the number of intervals - one less than the number
00383    * of points (catalog wavelengths, or detected peaks).
00384    */
00385 
00386   nlint = nlines - 1;
00387   npint = npeaks - 1;
00388 
00389 
00390   /*
00391    * Here the big loops on catalog lines begins.
00392    */
00393 
00394   for (i = 1; i < nlint; i++) {
00395 
00396 
00397     /*
00398      * For each catalog wavelength I take the previous and the next one, 
00399      * and compute the ratio of the corresponding wavelength intervals.
00400      * This ratio will be compared to all the ratios obtained doing the
00401      * same with all the detected peaks positions.
00402      */
00403 
00404     lratio = (line[i+1] - line[i]) / (line[i] - line[i-1]);
00405 
00406 
00407     /*
00408      * Here the loop on detected peaks positions begins.
00409      */
00410 
00411     for (j = 1; j < npint; j++) {
00412 
00413       /*
00414        * Not all peaks are used for computing ratios: just the ones
00415        * that are compatible with the expected spectral dispersion
00416        * are taken into consideration. Therefore, I define the pixel
00417        * intervals before and after any peak that are compatible with
00418        * the specified dispersion interval, and select just the peaks
00419        * within such intervals. If either of the two intervals doesn't
00420        * contain any peak, then I skip the current peak and continue
00421        * with the next.
00422        */
00423 
00424       lo_start = peak[j] - (line[i] - line[i-1]) / min_disp;
00425       lo_end   = peak[j] - (line[i] - line[i-1]) / max_disp;
00426       hi_start = peak[j] + (line[i+1] - line[i]) / max_disp;
00427       hi_end   = peak[j] + (line[i+1] - line[i]) / min_disp;
00428 
00429       for (npeaks_lo = 0, k = 0; k < npeaks; k++) {
00430         if (peak[k] > lo_end)
00431           break;
00432         if (peak[k] > lo_start) {
00433           peak_lo[npeaks_lo] = k;
00434           ++npeaks_lo;
00435         }
00436       }
00437 
00438       if (npeaks_lo == 0)
00439         continue;
00440 
00441       for (npeaks_hi = 0, k = 0; k < npeaks; k++) {
00442         if (peak[k] > hi_end)
00443           break;
00444         if (peak[k] > hi_start) {
00445           peak_hi[npeaks_hi] = k;
00446           ++npeaks_hi;
00447         }
00448       }
00449 
00450       if (npeaks_hi == 0)
00451         continue;
00452 
00453 
00454       /*
00455        * Now I have all peaks that may help for a local identification.
00456        * peak_lo[k] is the sequence number of the k-th peak of the lower
00457        * interval; peak_hi[l] is the sequence number of the l-th peak of
00458        * the higher interval. j is, of course, the sequence number of the
00459        * current peak (second big loop).
00460        */
00461 
00462       prev_variation = 1000.0;
00463       minl = mink = 0;
00464 
00465       for (k = 0; k < npeaks_lo; k++) {
00466         denom = peak[j] - peak[peak_lo[k]];
00467         for (l = 0; l < npeaks_hi; l++) {
00468 
00469           /*
00470            * For any pair of peaks - one from the lower and the other
00471            * from the higher interval - I compute the same ratio that
00472            * was computed with the current line catalog wavelength.
00473            */
00474 
00475           pratio = (peak[peak_hi[l]] - peak[j]) / denom;
00476 
00477           /*
00478            * If the two ratios are compatible within the specified
00479            * tolerance, we have a preliminary identification. This
00480            * will be marked in the matrix ident[][], where the first
00481            * index corresponds to a peak sequence number, and the second
00482            * index is the counter of the identifications made during
00483            * this whole process. The array of counters is nident[].
00484            * If more than one interval pair fulfills the specified
00485            * tolerance, the closest to the expected ratio is selected.
00486            */
00487 
00488           variation = fabs(lratio-pratio) / pratio;
00489 
00490           if (variation < tolerance) {
00491             if (variation < prev_variation) {
00492               prev_variation = variation;
00493               minl = l;
00494               mink = k;
00495             }
00496           }
00497         }
00498       }
00499       if (prev_variation < tolerance) {
00500         ident[j][nident[j]]                         = i;
00501         ident[peak_hi[minl]][nident[peak_hi[minl]]] = i + 1;
00502         ident[peak_lo[mink]][nident[peak_lo[mink]]] = i - 1;
00503         ++nident[j];
00504         ++nident[peak_hi[minl]];
00505         ++nident[peak_lo[mink]];
00506       }
00507     }   /* End loop on positions */
00508   }    /* End loop on lines     */
00509 
00510 
00511   /*
00512    * At this point I have filled the ident matrix with all my preliminary
00513    * identifications. Ambiguous identifications must be eliminated.
00514    */
00515 
00516 
00517   for (i = 0; i < npeaks; i++) {
00518 
00519 
00520     /*
00521      * I don't take into consideration peaks that were never identified.
00522      * They are likely contaminations, or emission lines that were not
00523      * listed in the input wavelength catalog.
00524      */
00525 
00526     if (nident[i] > 1) {
00527 
00528 
00529       /*
00530        * Initialise the histogram of wavelengths assigned to the i-th peak.
00531        */
00532 
00533       for (j = 0; j < nlines; j++)
00534         lident[j] = 0;
00535 
00536 
00537       /*
00538        * Count how many times each catalog wavelength was assigned
00539        * to the i-th peak.
00540        */
00541 
00542       for (j = 0; j < nident[i]; j++)
00543         ++lident[ident[i][j]];
00544 
00545 
00546       /*
00547        * What wavelength was most frequently assigned to the i-th peak?
00548        */
00549 
00550       max = 0;
00551       maxpos = 0;
00552       for (j = 0; j < nlines; j++) {
00553         if (max < lident[j]) {
00554           max = lident[j];
00555           maxpos = j;
00556         }
00557       }
00558 
00559 
00560       /*
00561        * Were there other wavelengths assigned with the same frequency?
00562        * This would be the case of an ambiguous identification. It is
00563        * safer to reject this peak...
00564        */
00565 
00566       ambiguous = 0;
00567 
00568       for (k = maxpos + 1; k < nlines; k++) {
00569         if (lident[k] == max) {
00570           ambiguous = 1;
00571           break;
00572         }
00573       }
00574 
00575       if (ambiguous)
00576         continue;
00577 
00578 
00579       /*
00580        * Otherwise, I assign to the i-th peak the wavelength that was
00581        * most often assigned to it.
00582        */
00583 
00584       tmp_xpos[n]   = peak[i];
00585       tmp_lambda[n] = line[maxpos];
00586       tmp_ilambda[n] = maxpos;
00587 
00588       ++n;
00589 
00590     }
00591 
00592   }
00593 
00594 
00595   /*
00596    * Check on identified peaks. Contaminations from other spectra might 
00597    * be present and should be excluded: this type of contamination 
00598    * consists of peaks that have been _correctly_ identified! The non-
00599    * spectral type of light contamination should have been almost all 
00600    * removed already in the previous steps, but it may still be present.
00601    * Here, the self-consistent sequences of identified peaks are
00602    * separated one from the other. At the moment, just the longest of
00603    * such sequences is selected (in other words, spectral multiplexing
00604    * is ignored).
00605    */
00606 
00607   if (n > 1) {
00608     nn = 0;                  /* Number of peaks in the list of sequences */
00609     nseq = 0;                /* Current sequence */
00610     for (k = 0; k < n; k++) {
00611       if (flag[k] == 0) {    /* Was peak k already assigned to a sequence? */
00612         flag[k] = 1;
00613         xpos[nn] = tmp_xpos[k];       /* Begin the nseq-th sequence */
00614         lambda[nn] = tmp_lambda[k];
00615         ilambda[nn] = tmp_ilambda[k];
00616         ++seq_length[nseq];
00617         ++nn;
00618 
00619         /*
00620          * Now look for all the following peaks that are compatible
00621          * with the expected spectral dispersion, and add them in 
00622          * sequence to xpos. Note that missing peaks are not a problem...
00623          */
00624          
00625         i = k;
00626         while (i < n - 1) {
00627           found = 0;
00628           for (j = i + 1; j < n; j++) {
00629             if (flag[j] == 0) {
00630               disp = (tmp_lambda[j] - tmp_lambda[i])
00631                    / (tmp_xpos[j] - tmp_xpos[i]);
00632               if (disp >= min_disp && disp <= max_disp) {
00633                 flag[j] = 1;
00634                 xpos[nn] = tmp_xpos[j];
00635                 lambda[nn] = tmp_lambda[j];
00636                 ilambda[nn] = tmp_ilambda[j];
00637                 ++seq_length[nseq];
00638                 ++nn;
00639                 i = j;
00640                 found = 1;
00641                 break;
00642               }
00643             }
00644           }
00645           if (!found)
00646             break;
00647         }
00648 
00649         /*
00650          * Current sequence is completed: begin new sequence on the
00651          * excluded peaks, starting the loop on peaks again.
00652          */
00653 
00654         ++nseq;
00655         k = 0;
00656       }
00657     }
00658 
00659 
00660     /*
00661      * Find the longest sequence of self-consistent peaks.
00662      */
00663 
00664     max = 0;
00665     maxpos = 0;
00666     for (i = 0; i < nseq; i++) {
00667       if (seq_length[i] > max) {
00668         max = seq_length[i];
00669         maxpos = i;
00670       }
00671     }
00672 
00673     /*
00674      * Find where this sequence starts in the whole peak position
00675      * storage.
00676      */
00677 
00678     nn = 0;
00679     for (i = 0; i < maxpos; i++)
00680       nn += seq_length[i];
00681 
00682     /*
00683      * Move the longest sequence at the beginning of the returned lists
00684      */
00685 
00686     n = max;
00687     for (i = 0; i < n; i++, nn++) {
00688       xpos[i] = xpos[nn];
00689       lambda[i] = lambda[nn];
00690       ilambda[i] = ilambda[nn];
00691     }
00692 
00693 
00694     /*
00695      * Are some wavelengths missing? Recover them.
00696      */
00697 
00698     for (i = 1; i < n; i++) {
00699       gap = ilambda[i] - ilambda[i-1];
00700       for (j = 1; j < gap; j++) {
00701 
00702         if (j == 1) {
00703 
00704           /*
00705            * Determine the local dispersion from the current pair of peaks
00706            */
00707   
00708           disp = (lambda[i] - lambda[i-1]) / (xpos[i] - xpos[i-1]);
00709         }
00710 
00711         /*
00712          * With this, find the expected position of the missing
00713          * peak by linear interpolation.
00714          */
00715 
00716         hi_start = xpos[i-1] + (line[ilambda[i-1] + j] - lambda[i-1]) / disp;
00717 
00718         /*
00719          * Is there a peak at that position? Here a peak from the
00720          * original list is searched, that is closer than 2 pixels
00721          * to the expected position. If it is found, insert it at
00722          * the current position on the list of identified peaks,
00723          * and leave immediately the loop (taking the new position
00724          * for the following linear interpolation, in case more
00725          * than one peak is missing in the current interval).
00726          * If it is not found, stay in the loop, looking for 
00727          * the following missing peaks in this interval.
00728          */
00729 
00730         found = 0;
00731         for (k = 0; k < npeaks; k++) {
00732           if (fabs(peak[k] - hi_start) < 2) {
00733             for (l = n; l > i; l--) {
00734               xpos[l] = xpos[l-1];
00735               lambda[l] = lambda[l-1];
00736               ilambda[l] = ilambda[l-1];
00737             }
00738             xpos[i] = peak[k];
00739             lambda[i] = line[ilambda[i-1] + j];
00740             ilambda[i] = ilambda[i-1] + j;
00741             ++n;
00742             found = 1;
00743             break;
00744           }
00745         }
00746         if (found)
00747           break;
00748       }
00749     }
00750 
00751 
00752     /*
00753      * Try to extrapolate forward
00754      */
00755 
00756     found = 1;
00757     while (ilambda[n-1] < nlines - 1 && found) {
00758 
00759       /*
00760        * Determine the local dispersion from the last pair of 
00761        * identified peaks
00762        */
00763 
00764       disp = (lambda[n-1] - lambda[n-2]) / (xpos[n-1] - xpos[n-2]);
00765 
00766       if (disp > max_disp || disp < min_disp)
00767         break;
00768 
00769 
00770       /*
00771        * With this, find the expected position of the missing
00772        * peak by linear interpolation.
00773        */
00774 
00775       hi_start = xpos[n-1] + (line[ilambda[n-1] + 1] - lambda[n-1]) / disp;
00776 
00777       /*
00778        * Is there a peak at that position? Here a peak from the
00779        * original list is searched, that is closer than 3 pixels
00780        * to the expected position. If it is found, insert it at
00781        * the end of the list of identified peaks. If it is not
00782        * found, leave the loop.
00783        */
00784 
00785       found = 0;
00786       min = fabs(peak[0] - hi_start);
00787       minpos = 0;
00788       for (k = 1; k < npeaks; k++) {
00789         if (min > fabs(peak[k] - hi_start)) {
00790             min = fabs(peak[k] - hi_start);
00791             minpos = k;
00792         }
00793       }
00794       if (min < 6 && fabs(peak[minpos] - xpos[n-1]) > 1.0) {
00795         xpos[n] = peak[minpos];
00796         lambda[n] = line[ilambda[n-1] + 1];
00797         ilambda[n] = ilambda[n-1] + 1;
00798         ++n;
00799         found = 1;
00800       }
00801     }
00802 
00803 
00804     /*
00805      * Try to extrapolate backward
00806      */
00807 
00808     found = 1;
00809     while (ilambda[0] > 0 && found) {
00810 
00811       /*
00812        * Determine the local dispersion from the first pair of
00813        * identified peaks
00814        */
00815 
00816       disp = (lambda[1] - lambda[0]) / (xpos[1] - xpos[0]);
00817 
00818       if (disp > max_disp || disp < min_disp)
00819         break;
00820 
00821 
00822       /*
00823        * With this, find the expected position of the missing
00824        * peak by linear interpolation.
00825        */
00826 
00827       hi_start = xpos[0] - (lambda[0] - line[ilambda[0] - 1]) / disp;
00828 
00829       /*
00830        * Is there a peak at that position? Here a peak from the
00831        * original list is searched, that is closer than 3 pixels
00832        * to the expected position. If it is found, insert it at
00833        * the beginning of the list of identified peaks. If it is not
00834        * found, leave the loop.
00835        */
00836 
00837       found = 0;
00838       min = fabs(peak[0] - hi_start);
00839       minpos = 0;
00840       for (k = 1; k < npeaks; k++) {
00841         if (min > fabs(peak[k] - hi_start)) {
00842             min = fabs(peak[k] - hi_start);
00843             minpos = k;
00844         }
00845       }
00846       if (min < 6 && fabs(peak[minpos] - xpos[0]) > 1.0) {
00847         for (j = n; j > 0; j--) {
00848           xpos[j] = xpos[j-1];
00849           lambda[j] = lambda[j-1];
00850           ilambda[j] = ilambda[j-1];
00851         }
00852         xpos[0] = peak[minpos];
00853         lambda[0] = line[ilambda[0] - 1];
00854         ilambda[0] = ilambda[0] - 1;
00855         ++n;
00856         found = 1;
00857       }
00858     }
00859   }
00860 
00861 
00862   /*
00863    * At this point all peaks are processed. Free memory, and return
00864    * the result.
00865    */
00866 
00867 /************************************************+
00868   for (i = 0; i < npeaks; i++) {
00869     printf("Peak %d:\n   ", i);
00870     for (j = 0; j < nident[i]; j++)
00871       printf("%.2f, ", line[ident[i][j]]);
00872     printf("\n");
00873   }
00874 
00875   printf("\n");
00876 
00877   for (i = 0; i < n; i++)
00878     printf("%.2f, %.2f\n", xpos[i], lambda[i]);
00879 +************************************************/
00880   for (i = 0; i < npeaks; i++)
00881     cpl_free(ident[i]);
00882   cpl_free(ident);
00883   cpl_free(nident);
00884   cpl_free(lident);
00885   cpl_free(ilambda);
00886   cpl_free(tmp_xpos);
00887   cpl_free(tmp_lambda);
00888   cpl_free(tmp_ilambda);
00889   cpl_free(peak_lo);
00890   cpl_free(flag);
00891   cpl_free(seq_length);
00892   cpl_free(peak_hi);
00893 
00894   if (n == 0) {
00895     cpl_free(xpos);
00896     cpl_free(lambda);
00897     return NULL;
00898   }
00899 
00900   return cpl_bivector_wrap_vectors(cpl_vector_wrap(n, xpos), 
00901                                    cpl_vector_wrap(n, lambda));
00902 }
00903 
00906 /*----------------------------------------------------------------------------*/
00917 /*----------------------------------------------------------------------------*/
00918 static cpl_vector * irplib_ppm_convolve_line(
00919         const cpl_vector        *   spectrum,
00920         double                      slitw,
00921         double                      fwhm)
00922 {
00923     cpl_vector  *   conv_kernel ;
00924     cpl_vector  *   line_profile ;
00925     cpl_vector  *   xcorrs ;
00926     cpl_vector  *   spec_ext ;
00927     cpl_vector  *   xc_single ;
00928     int             hs, line_sz, sp_sz ;
00929     int             i ;
00930 
00931     /* Test entries */
00932     if (spectrum == NULL) return NULL ;
00933 
00934     /* Create the convolution kernel */
00935     if ((conv_kernel = irplib_wlxcorr_convolve_create_kernel(slitw,
00936                     fwhm)) == NULL) {
00937         cpl_msg_error(cpl_func, "Cannot create kernel") ;
00938         return NULL ;
00939     }
00940     hs = cpl_vector_get_size(conv_kernel) ;
00941     line_sz = 2 * hs + 1 ;
00942     
00943     /* Create the line profile */
00944     line_profile = cpl_vector_new(line_sz) ;
00945     cpl_vector_fill(line_profile, 0.0) ;
00946     cpl_vector_set(line_profile, hs, 1.0) ;
00947     if (irplib_wlxcorr_convolve(line_profile, conv_kernel) != 0) {
00948         cpl_msg_error(cpl_func, "Cannot create line profile") ;
00949         cpl_vector_delete(line_profile) ;
00950         cpl_vector_delete(conv_kernel) ;
00951         return NULL ;
00952     }
00953     cpl_vector_delete(conv_kernel) ;
00954     
00955     /* Create the correlations values vector */
00956     sp_sz = cpl_vector_get_size(spectrum) ;
00957     xcorrs = cpl_vector_new(sp_sz) ;
00958     cpl_vector_fill(xcorrs, 0.0) ;
00959     xc_single = cpl_vector_new(1) ;
00960 
00961     /* Loop on the pixels of the spectrum */
00962     for (i=hs ; i<sp_sz-hs ; i++) {
00963         /* Extract the current spectrum part */
00964         if ((spec_ext = cpl_vector_extract(spectrum, i-hs, i+hs, 1)) == NULL) {
00965             cpl_msg_error(cpl_func, "Cannot extract spectrum") ;
00966             cpl_vector_delete(xc_single) ;
00967             cpl_vector_delete(line_profile) ;
00968             return NULL ;
00969         }
00970         if (cpl_vector_correlate(xc_single, spec_ext, line_profile) < 0) {
00971             cpl_msg_error(cpl_func, "Cannot correlate") ;
00972             cpl_vector_delete(xc_single) ;
00973             cpl_vector_delete(line_profile) ;
00974             cpl_vector_delete(spec_ext) ;
00975             return NULL ;
00976         }
00977         cpl_vector_set(xcorrs, i, cpl_vector_get(xc_single, 0)) ;
00978         cpl_vector_delete(spec_ext) ;
00979     }
00980     cpl_vector_delete(xc_single) ;
00981     cpl_vector_delete(line_profile) ;
00982 
00983     return xcorrs ;
00984 } 
00985 
00986 /*----------------------------------------------------------------------------*/
00995 /*----------------------------------------------------------------------------*/
00996 static cpl_vector * irplib_ppm_detect_lines(
00997         const cpl_vector    *   spec,
00998         double                  threshold)
00999 {
01000     cpl_vector  *   spec_loc ;
01001     double      *   pspec_loc ;
01002     cpl_vector  *   lines ;
01003     double      *   plines ;
01004     int             spec_loc_sz, max_ind, nlines ;
01005     double          max ;
01006     int             i ;
01007 
01008     /* Test inputs */
01009     if (spec == NULL) return NULL ;
01010     if (threshold<0.0 || threshold > 1.0) return NULL ;
01011 
01012     /* Local spectrum */
01013     spec_loc = cpl_vector_duplicate(spec) ;
01014     pspec_loc = cpl_vector_get_data(spec_loc) ;
01015     spec_loc_sz = cpl_vector_get_size(spec_loc) ;
01016 
01017     /* Threshold the local spectrum */
01018     for (i=0 ; i<spec_loc_sz ; i++) 
01019         if (pspec_loc[i] < threshold) pspec_loc[i] = 0.0 ;
01020     
01021     /* Allocate lines container */
01022     lines = cpl_vector_new(spec_loc_sz) ;
01023     plines = cpl_vector_get_data(lines) ;
01024     nlines = 0 ;
01025     
01026     /* Loop as long as there are lines */
01027     while ((max = cpl_vector_get_max(spec_loc)) > threshold) {
01028         /* Find the max position */
01029         max_ind = 0 ;
01030         while (pspec_loc[max_ind]<max && max_ind<spec_loc_sz) max_ind++ ;
01031         if (max_ind == spec_loc_sz) {
01032             cpl_msg_error(cpl_func, "Cannot find maximum") ;
01033             cpl_vector_delete(spec_loc) ;
01034             cpl_vector_delete(lines) ;
01035             return NULL ;
01036         }
01037         if (max_ind == 0 || max_ind == spec_loc_sz-1) {
01038             pspec_loc[max_ind] = 0 ;
01039             continue ;
01040         }
01041 
01042         /* Get the precise position from the neighbours values */
01043         plines[nlines] =    pspec_loc[max_ind] * max_ind + 
01044                             pspec_loc[max_ind-1] * (max_ind-1) +
01045                             pspec_loc[max_ind+1] * (max_ind+1) ; 
01046         plines[nlines] /= pspec_loc[max_ind] + pspec_loc[max_ind+1] +
01047             pspec_loc[max_ind-1] ;
01048         plines[nlines] ++ ;
01049         nlines ++ ;
01050 
01051         /* Clean the line */
01052         i = max_ind ;
01053         while (i>=0 && pspec_loc[i] > threshold) {
01054             pspec_loc[i] = 0.0 ;
01055             i-- ;
01056         }
01057         i = max_ind+1 ;
01058         while (i<spec_loc_sz && pspec_loc[i] > threshold) {
01059             pspec_loc[i] = 0.0 ;
01060             i++ ;
01061         }
01062     }
01063     cpl_vector_delete(spec_loc) ;
01064    
01065     /* Check if there are lines */
01066     if (nlines == 0) {
01067         cpl_msg_error(cpl_func, "Cannot detect any line") ;
01068         cpl_vector_delete(lines) ;
01069         return NULL ;
01070     }
01071     
01072     /* Resize the vector */
01073     cpl_vector_set_size(lines, nlines) ;
01074 
01075     /* Sort the lines */
01076     cpl_vector_sort(lines, 1) ;
01077     
01078     return lines ;
01079 }
01080 
01081 
01082 
01083 

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