sinfo_wave_calibration.c

00001 /*
00002  * This file is part of the ESO SINFONI Pipeline
00003  * Copyright (C) 2004,2005 European Southern Observatory
00004  *
00005  * This program is free software; you can redistribute it and/or modify
00006  * it under the terms of the GNU General Public License as published by
00007  * the Free Software Foundation; either version 2 of the License, or
00008  * (at your option) any later version.
00009  *
00010  * This program is distributed in the hope that it will be useful,
00011  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00012  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00013  * GNU General Public License for more details.
00014  *
00015  * You should have received a copy of the GNU General Public License
00016  * along with this program; if not, write to the Free Software
00017  * Foundation, 51 Franklin St, Fifth Floor, Boston, MA  02111-1307  USA
00018  */
00019 /*************************************************************************
00020 * E.S.O. - VLT project
00021 *
00022 * 
00023 *
00024 * who       when      what
00025 * --------  --------  ----------------------------------------------
00026 * schreib  13/07/00  created
00027 */
00028 
00029 /************************************************************************
00030 *   NAME
00031 *     sinfo_wave_calibration.c -
00032 *     routines needed for wavelength calibration
00033 *
00034 *   SYNOPSIS
00035 * 
00036 *   1) FitParams ** sinfo_new_fit_params( int n_params )
00037 *
00038 *   2) void sinfo_new_destroy_fit_params ( FitParams ** params )
00039 *
00040 *   3) void sinfo_new_dump_fit_params_to_ascii(FitParams ** params,
00041                                                char * filename )
00042 *
00043 *   4) void sinfo_new_dump_ascii_to_fit_params ( FitParams ** params, 
00044                                                  char * filename )
00045 *
00046 *   5) int sinfo_new_find_lines(cpl_image * lineImage, 
00047 *                     float    * wave_position, 
00048 *                     float    * wave_intensity,
00049 *                     int        n_lines, 
00050 *                     int     ** row_clean,
00051 *                     float   ** wavelength_clean,
00052 *                     float      beginWave, 
00053 *                     float      dispersion,
00054 *                     float      mindiff, 
00055 *                     int        halfWidth,
00056 *                     int      * n_found_lines,
00057 *                     float      sigma,
00058 *                     int      * sum_lines )
00059 *
00060 *   6) int sinfo_new_read_list( char * listname, 
00061                                 float * lineCenter, 
00062                                 float * lineIntensity )
00063 *
00064 *
00065 *   7) int sinfo_new_line_fit ( cpl_image  *  mergedImage, 
00066 *                    FitParams *  par,
00067 *                    float        fwhm,
00068 *                    int          lineInd,
00069 *                    int          column, 
00070 *                    int          halfWidth, 
00071 *                    int          lineRow,
00072 *                    float        min_amplitude )
00073 *
00074 *   8) int sinfo_new_fit_lines ( cpl_image  *  line_image, 
00075 *                     FitParams ** allParams,
00076 *                     float        fwhm,
00077 *                     int       *  n_lines, 
00078 *                     int       ** row,
00079 *                     float     ** wavelength, 
00080 *                     int          width,
00081 *                     float        min_amplitude ) 
00082 *
00083 *   9) float sinfo_new_polyfit( FitParams ** par, 
00084 *                     int          column,
00085 *                     int          n_lines,
00086 *                     int          n_rows,
00087 *                     float        dispersion,
00088 *                     float        max_residual,
00089 *                     float *      acoefs, 
00090 *                     float *      dacoefs, 
00091 *                     int   *      n_reject,
00092 *                     int          n_fitcoefs )
00093 *
00094 *   10) float sinfo_new_coefs_cross_fit ( int      n_columns,
00095 *                            float *  acoefs,
00096 *                            float *  dacoefs, 
00097 *                            float *  bcoefs,
00098 *                            int      n_fitcoefs,
00099 *                            float    sigma_factor )
00100 *
00101 *
00102 *  11) cpl_image * sinfo_new_wave_map( cpl_image * lineImage,
00103 *                          float   ** bcoefs,
00104 *                          int        n_a_fitcoefs,
00105 *                          int        n_b_fitcoefs,
00106 *                          float    * wavelength,
00107 *                          float    * intensity,
00108 *                          int        n_lines,
00109 *                          int        magFactor,
00110 *                          int      * bad_column_mask,
00111 *                          int        n_bad_columns )
00112 *
00113 *  12) int sinfo_new_wavelength_calibration( cpl_image   * image, 
00114 *                                 FitParams ** par ,
00115 *                                 float     ** bcoefs,
00116 *                                 float      * wave,
00117 *                                 int          n_lines,
00118 *                                 int       ** row_clean,
00119 *                                 float     ** wavelength_clean,
00120 *                                 int        * n_found_lines,
00121 *                                 float        dispersion,
00122 *                                 int          halfWidth,
00123 *                                 float        minAmplitude,
00124 *                                 float        max_residual,
00125 *                                 float        fwhm,
00126 *                                 int          n_a_fitcoefs,
00127 *                                 int          n_b_fitcoefs,
00128 *                                 float        sigmaFactor )
00129 *
00130 *  13) cpl_image * sinfo_new_convolve_image_by_gauss( cpl_image * lineImage,
00131 *                                       int        hw )
00132 *
00133 *  14) cpl_image * sinfo_new_defined_resampling( cpl_image * image,
00134 *                                    cpl_image * calimage,
00135 *                                    int        n_params,
00136 *                                    int        n_rows,
00137 *                                    double   * dispersion,
00138 *                                    float    * minval,
00139 *                                    float    * maxval,
00140 *                                    double   * centralLambda,
00141 *                                    int      * centralpix )
00142 *
00143 *   DESCRIPTION
00144 *
00145 *   1) allocates memory for a new sinfo_vector of 
00146 *      FitParams data structures
00147 *   2) frees memory of a sinfo_vector of FitParams data structures
00148 *   3) dumps the fit parameters to an ASCII file
00149 *   4) dumps ASCII information to an allocated FitParams data structure
00150 *   5) determines the pixel shift between the line list 
00151 *      and the real image by using the beginning wavelength
00152 *      on the detector and the dispersion estimate.
00153 *   6) reads the line data of the calibration lamps
00154 *   7) fits a sinfo_gaussian to a 1-dimensional slice of an image, 
00155 *      this routine uses the routine sinfo_new_lsqfit_c as a non-linear
00156 *      least square fit method (Levenberg-Marquardt).               
00157 *   8) calculates and stores the fit parameters of the neon 
00158 *      emission lines of a neon frame by using the sinfo_linefit 
00159 *      routine.
00160 *   9) fits a second order polynom 
00161 *      lambda[i] = a1 + a2*pos[i] + a3*pos[i]^2
00162 *      to determine the connection between the listed wave-
00163 *      length values and the gauss-fitted positions for each
00164 *      image column using the singular value decomposition 
00165 *      method. 
00166 *  10) Fits the each single parameter of the three fit parameters 
00167 *      acoefs from sinfo_polyfit through the image columns
00168 *  11) this routine determines a wavelength calibration map 
00169 *      frame associating a wavelength value to each pixel
00170 *      by using the fit coefficients determined before.
00171 *  12) this routine takes an image from a calibration
00172 *      emission lamp and delivers the fit coefficients of  
00173 *      a polynomial fit across the columns 
00174 *      of the coefficients of the polynomial line position
00175 *      fits as output. Furthermore it delivers an array of the fit parameters
00176 *      as output. This routine expects Nyquist sampled spectra 
00177 *     (either an interleaved image or an image convolved with an 
00178 *      appropriate function in spectral direction)
00179 *  13) convolves an emission line image with a sinfo_gaussian
00180 *      with user given integer half width by using the eclipse 
00181 *      routine sinfo_function1d_filter_lowpass().
00182 *  14) Given a source image and a corresponding wavelength
00183 *      calibration file this routine produces an image
00184 *      in which elements in a given row are associated
00185 *      with a single wavelength. It thus corrects for 
00186 *      the wavelength shifts between adjacent elements
00187 *      in the rows of the input image. The output image
00188 *      is larger in the wavelength domain than the input
00189 *      image with pixels in each column corresponding to 
00190 *      undefined (blank, ZERO) values. The distribution
00191 *      of these undefined values varies from column to
00192 *      column. The input image is resampled at discrete
00193 *      wavelength intervals using a polynomial interpolation
00194 *      routine. 
00195 *      The wavelength intervals (dispersion) and the 
00196 *      central wavelength are defined and stable for each
00197 *      used grating. Thus, each row has a defined wavelength
00198 *      for each grating. Only the number of rows can be 
00199 *      changed by the user. 
00200 *
00201 *   FILES
00202 *
00203 *   ENVIRONMENT
00204 *
00205 *   RETURN VALUES 
00206 *
00207 *   CAUTIONS 
00208 *
00209 *   EXAMPLES
00210 *
00211 *   SEE ALSO
00212 *
00213 *   BUGS   
00214 *
00215 *------------------------------------------------------------------------
00216 */
00217 
00218 #ifdef HAVE_CONFIG_H
00219 #  include <config.h>
00220 #endif
00221 #include "sinfo_vltPort.h"
00222 
00223 /* 
00224  * System Headers
00225  */
00226 
00227 /* 
00228  * Local Headers
00229  */
00230 
00231 #include "sinfo_function_1d.h"
00232 #include "sinfo_wave_calibration.h"
00233 #include "sinfo_solve_poly_root.h"
00234 #include "sinfo_recipes.h"
00235 #include "sinfo_svd.h"
00236 
00237 
00253 FitParams ** sinfo_new_fit_params( int n_params )
00254 {
00255     FitParams ** new_params =NULL;
00256     FitParams  * temp_params =NULL;
00257     float * temp_fit_mem =NULL;
00258     float * temp_derv_mem=NULL;
00259     int i ;
00260 
00261     
00262     if ( n_params <= 0 )
00263     {
00264         sinfo_msg_error (" wrong number of lines to fit\n") ;
00265         return NULL ;
00266     }
00267 
00268     if (NULL==(new_params=(FitParams **) cpl_calloc ( n_params , 
00269                                          sizeof (FitParams*) ) ) )
00270     {
00271         sinfo_msg_error (" could not allocate memory\n") ;
00272         return NULL ;
00273     }
00274     if ( NULL == (temp_params = cpl_calloc ( n_params , sizeof (FitParams) ) ) )
00275     {
00276         sinfo_msg_error (" could not allocate memory\n") ;
00277         return NULL ;
00278     }
00279     if ( NULL == (temp_fit_mem = (float *) cpl_calloc( n_params*MAXPAR, 
00280                                      sizeof (float) ) ) )
00281     {
00282         sinfo_msg_error (" could not allocate memory\n") ;
00283         return NULL ;
00284     }
00285 
00286 
00287     if ( NULL == (temp_derv_mem   = 
00288              (float *) cpl_calloc( n_params*MAXPAR, sizeof (float) ) ) )
00289     {
00290         sinfo_msg_error (" could not allocate memory\n") ;
00291         return NULL ;
00292     }
00293 
00294     for ( i = 0 ; i < n_params ; i ++ )
00295     {
00296         new_params[i] = temp_params+i;
00297         new_params[i] -> fit_par    = temp_fit_mem+i*MAXPAR;
00298         new_params[i] -> derv_par   = temp_derv_mem+i*MAXPAR;
00299         new_params[i] -> column     = 0 ;
00300         new_params[i] -> line       = 0 ;
00301         new_params[i] -> wavelength = 0. ;
00302         new_params[i] -> n_params   = n_params ;
00303     }
00304 
00305     return new_params ;
00306 }
00307 
00315 void sinfo_new_destroy_fit_params ( FitParams *** params )
00316 { 
00317     
00318     if ( *params == NULL )
00319     {
00320         return ;
00321     }
00322 
00323     cpl_free ( (*params)[0] -> fit_par ) ;
00324     (*params)[0] -> fit_par=NULL;
00325     cpl_free ( (*params)[0] -> derv_par ) ;
00326     (*params)[0] -> derv_par=NULL;
00327     cpl_free ( (*params)[0] ) ;
00328     (*params)[0]=NULL;
00329     cpl_free ( (*params) ) ;
00330     (*params)=NULL;
00331 }
00332 
00341 void sinfo_new_dump_fit_params_to_ascii ( FitParams ** params, char * filename )
00342 {
00343     FILE * fp ;
00344     int    i ;
00345 
00346     if ( NULL == params )
00347     {
00348         sinfo_msg_error (" no fit parameters available!\n") ;
00349         return ;
00350     }
00351 
00352     if ( NULL == filename )
00353     {
00354         sinfo_msg_error (" no filename available!\n") ;
00355         return ;
00356     }
00357 
00358     if ( NULL == (fp = fopen ( filename, "w" ) ) )
00359     {
00360         sinfo_msg_error(" cannot open %s\n", filename) ;
00361         return ;
00362     }
00363 
00364     for ( i = 0 ; i < params[0] -> n_params ; i++ )
00365     {
00366         fprintf(fp, "%d %d %d %f %f %f %f %f %f %f %f %f\n", 
00367                      params[i]->n_params, 
00368                      params[i]->column, 
00369                      params[i]->line, 
00370                      params[i]->wavelength, 
00371                      params[i]->fit_par[0], 
00372                      params[i]->fit_par[1], 
00373                      params[i]->fit_par[2], 
00374                      params[i]->fit_par[3],
00375                      params[i]->derv_par[0], 
00376                      params[i]->derv_par[1], 
00377                      params[i]->derv_par[2], 
00378                      params[i]->derv_par[3] ) ; 
00379     }
00380     fclose(fp) ;
00381 }
00382 
00391 void 
00392 sinfo_new_dump_ascii_to_fit_params ( FitParams ** params, char * filename )
00393 {
00394     FILE * fp ;
00395     int i ;
00396 
00397     if ( NULL == params )
00398     {
00399         sinfo_msg_error (" no fit parameters available!\n") ;
00400         return ;
00401     }
00402 
00403     if ( NULL == filename )
00404     {
00405         sinfo_msg_error (" no filename available!\n") ;
00406         return ;
00407     }
00408 
00409     if ( NULL == (fp = fopen ( filename, "r" ) ) )
00410     {
00411         sinfo_msg_error(" cannot open %s\n", filename) ;
00412         return ;
00413     }
00414     
00415     for ( i = 0 ; i < params[0]->n_params ; i++ )
00416     {
00417         fscanf(fp, "%d %d %d %f %f %f %f %f %f %f %f %f\n", 
00418                     &params[i]->n_params, 
00419                     &params[i]->column, 
00420                     &params[i]->line, 
00421                     &params[i]->wavelength, 
00422                     &params[i]->fit_par[0], 
00423                     &params[i]->fit_par[1], 
00424                     &params[i]->fit_par[2], 
00425                     &params[i]->fit_par[3],
00426                     &params[i]->derv_par[0], 
00427                     &params[i]->derv_par[1], 
00428                     &params[i]->derv_par[2], 
00429                     &params[i]->derv_par[3] ) ; 
00430     }
00431     fclose(fp) ;
00432 }
00433 
00472 int sinfo_new_find_lines(cpl_image * lineImage, 
00473                float    * wave_position, 
00474                float    * wave_intensity,
00475                int        n_lines, 
00476                int     ** row_clean,
00477                float   ** wavelength_clean,
00478                float      beginWave, 
00479                float      dispersion1,
00480                float      dispersion2,
00481                float      mindiff, 
00482                int        halfWidth,
00483                int      * n_found_lines,
00484                float      sigma,
00485                int      * sum_lines )
00486 {
00487     int     ** row ;
00488     float   ** wavelength ;
00489     float buf1, buf2 ;
00490     float meanval ;
00491     float colmedian ;
00492     float * column, * tempcol ;
00493     float * lines ;
00494     float * conv_lines ;
00495     float * wave_buffer ;
00496     float * wave_mem;
00497     int   * dummy_row ;
00498     int i, j, k, m ;
00499     int position ;
00500     int gmax, gmin ;
00501     int col ;
00502     int * row_mem;
00503     float sum ;
00504     float angst ;
00505 
00506     int lx=0;
00507     int ly=0;
00508     float* pdata=NULL;
00509 
00510     if ( NULL == lineImage )
00511     {
00512         sinfo_msg_error (" no image given\n") ;
00513         return -1 ;
00514     }
00515 
00516     lx=cpl_image_get_size_x(lineImage);
00517     ly=cpl_image_get_size_y(lineImage);
00518     pdata=cpl_image_get_data_float(lineImage);
00519 
00520     if ( n_lines <= 0 || NULL == wave_position ) 
00521     {
00522         sinfo_msg_error(" no line list given\n") ;
00523         return -1 ;
00524     }
00525     if ( NULL == wave_intensity ) 
00526     {
00527         sinfo_msg_error(" no intensity list given\n") ;
00528         return -1 ;
00529     }
00530      
00531     if ( dispersion1 == 0. )
00532     {
00533         sinfo_msg_error(" wrong dispersion given\n") ;
00534         return -1 ;
00535     }
00536    
00537     if ( row_clean == NULL )
00538     {
00539         sinfo_msg_error(" row_clean array is not allocated\n") ;
00540         return -1 ;
00541     }
00542 
00543     if ( wavelength_clean == NULL )
00544     {
00545         sinfo_msg_error(" wavelength_clean array is not allocated\n") ;
00546         return -1 ;
00547     }
00548 
00549     if ( beginWave <= 0. )
00550     {
00551         sinfo_msg_error (" impossible beginWave given\n") ;
00552         return -1 ;
00553     }
00554     if ( mindiff < -100. )
00555     {
00556         sinfo_msg_error (" wrong mindiff value\n") ;
00557         return -1 ;
00558     }
00559 
00560     if ( halfWidth <= 0 )
00561     {
00562         sinfo_msg_error(" wrong half width of fit box given\n") ;
00563         return -1 ;
00564     }
00565 
00566     if ( n_found_lines == NULL )
00567     {
00568         sinfo_msg_error(" n_found_lines not allocated\n") ;
00569         return -1 ;
00570     }
00571 
00572     if ( sigma <= 0. || sigma >= ly / 2)
00573     {
00574         sinfo_msg_error(" wrong sigma given\n") ;
00575         return -1 ;
00576     }
00577     
00578     /* allocate memory */
00579     row        = (int**) cpl_calloc( lx, sizeof(int*)) ;
00580     wavelength = (float**) cpl_calloc( lx, sizeof(float*)) ;
00581     row_mem = cpl_calloc( n_lines*lx, sizeof(int) ) ;
00582     wave_mem = cpl_calloc( n_lines*lx, sizeof(float) ) ;
00583     for ( i = 0 ; i <lx ; i++ )
00584     {
00585         row[i] = row_mem + i*n_lines;
00586         wavelength[i] = wave_mem + i*n_lines;
00587     }
00588 
00589     /* find if the wavelength is given in microns, nanometers or Angstroem */
00590     if ( wave_position[0] > 10000. )
00591     {
00592     /* Angstroem */
00593         angst = 10000. ;
00594     }
00595     else if ( wave_position[0] > 1000. && wave_position[0] < 10000. )
00596     {
00597     /* nanometers */
00598     angst = 1000. ;
00599     }
00600     else
00601     {
00602     /* microns */
00603     angst = 1. ;
00604     }
00605 
00606     /*----------------------------------------------------------------------
00607      * compute the mean and median intensity value in the given 
00608        column and determine if there is enough intensity in the column to 
00609        do the correlation
00610      */
00611     tempcol = (float*) cpl_calloc(ly, sizeof(float)) ;
00612     *sum_lines = 0 ;
00613     buf1 = 0. ;
00614     buf2 = 0. ;
00615     /* allocate memory */
00616     column      = (float*) cpl_calloc(ly, sizeof (float)) ;
00617     lines       = (float*) cpl_calloc(ly, sizeof (float)) ;
00618     conv_lines  = (float*) cpl_calloc(ly, sizeof (float)) ;
00619     wave_buffer = (float*) cpl_calloc(ly, sizeof (float)) ;
00620     dummy_row   = (int*)   cpl_calloc(n_lines, sizeof(int)) ;
00621 
00622     for ( col = 0 ; col < lx ; col++ )
00623     {
00624         n_found_lines[col] = 0 ;
00625         sum = 0. ;
00626         for ( i = 0 ; i < ly ; i++ )
00627         {
00628             if (isnan(pdata[col + i*lx]) )
00629             {
00630                 tempcol[i] = 0. ;
00631                 continue ;
00632             }
00633 
00634         sum = sum + pdata[col + i*lx] ;
00635         tempcol[i] = pdata[col + i*lx];
00636         }
00637         meanval = sum / ly ;
00638         /* lets assume the sinfo_new_median is the background */
00639         colmedian =  sinfo_new_median ( tempcol, ly);
00640 
00641         if ( meanval - colmedian < mindiff )
00642         {
00643         sinfo_msg_warning(" sorry, not enough intensity "
00644                               "(mean: %f, diff: %f) in image column: "
00645                                "%d to correlate emission lines\n", 
00646                                 meanval, meanval - colmedian,col) ;
00647         continue ;
00648     }
00649 
00650         for ( i = 0 ; i < ly ; i++ )
00651     {
00652       conv_lines[i]=0;
00653       wave_buffer[i]=0;
00654     }
00655     for ( i = 0 ; i < n_lines ; i++ )
00656     {
00657       dummy_row[i] = 0;
00658     }
00659 
00660         /* go through the column with index col */
00661         for ( i = 0 ; i < ly ; i++ )
00662         {
00663             if ( isnan(pdata[col+i*lx]) )
00664             {
00665                 column[i] = 0. ;
00666             }
00667             else
00668             {
00669                 column[i] = pdata[col + i*lx] ;
00670             }
00671 
00672             /* determine the line position on the pixels */
00673             /*lines[i] = (dispersion * (float) i + beginWave) * angst ;*/ 
00674             lines[i] = (dispersion1 * (float) (i-ly/2) + 
00675                         dispersion2 * (float) (i-ly/2) * 
00676                                       (float) (i-ly/2) + 
00677                                               beginWave) * angst ; 
00678 
00679             /* ---------------------------------------------------------------
00680              * find the nearest line position for each wavelength in the list 
00681              * and set this position to the given line intensity as weight 
00682              */
00683             for ( j = 0 ; j < n_lines ; j ++ )
00684             {
00685                 buf1 = 0. ;
00686                 buf2 = 0. ;
00687                 if ( (wave_position[j] >= (lines[i] - 
00688                       fabs(dispersion1)/2.*angst)) && 
00689                      (wave_position[j] <= (lines[i] + 
00690                       fabs(dispersion1)/2.*angst)) ) 
00691                 {
00692                     buf1 = wave_intensity[j] ; /* set the given line intensity 
00693                                                   as weight */
00694                     buf2 = wave_position[j] ;
00695                     break ;
00696                 }
00697             }
00698             lines[i] = buf1 ;
00699             wave_buffer[i] = buf2 ; /* get the wavelength associated 
00700                                        with the corresponding 
00701                                            found emission line */
00702     
00703             /* convolve the artificial spectrum by a Gaussian 
00704                with given sigma value */
00705             if ( lines[i] != 0. )
00706             {
00707                 /* consider only +- 2 sigma */
00708                 gmin = sinfo_new_nint((float) i - 2. * sigma) ;
00709                 gmax = sinfo_new_nint((float) i + 2. * sigma) ;
00710 
00711                 /* be aware of image margins */
00712                 if ( gmin < 0 )
00713                 {
00714                     gmin = 0 ;
00715                 }
00716                 if ( gmax >= ly )
00717                 {
00718                     gmax = ly - 1 ;
00719                 } 
00720                 for ( j = gmin ; j <= gmax ; j++ )
00721                 {
00722                     conv_lines[j] += 
00723                         lines[i] * exp( (double)( -0.5*((j - i)*(j - i)))/
00724                                         (double) (sigma*sigma) ) ; 
00725                 }
00726             }
00727         }
00728 
00729         /* do the cross sinfo_new_correlation */
00730         position = INT32_MAX ;
00731         position = sinfo_new_correlation(column, conv_lines, ly ) ; 
00732         if ( abs (position) > ly / 4 )
00733         {
00734             sinfo_msg_warning(" sorry, shift of artificial data relative to"
00735                               " image (%d) seems to be too high in column: %d",
00736                               position, col) ;
00737             continue ;
00738         }
00739 
00740         j = 0 ;
00741         for ( i = 0 ; i < ly ; i ++ )
00742         {
00743             if ( lines[i] != 0.0 )
00744             {
00745                 if ( (i - position) >= 0 && (i - position) < ly )
00746                 {
00747                     row[col][j] = i - position ;
00748                     /* get the wavelength corresponding to 
00749                        found line row index */
00750                     wavelength[col][j] = wave_buffer[i] / angst ;
00751                     j++ ;
00752                 }
00753             }
00754         }
00755 
00756 
00757         /* ------------------------------------------------------------------
00758          *  determine the row_clean array, that means, take only the row 
00759             values if the distance between adjacent lines is large enough 
00760             for the fit
00761          */
00762         for ( k = 1 ; k <= j && k<(lx-1); k ++ )
00763         {
00764             if (dummy_row[k-1] != -1)
00765             {
00766                 dummy_row[k-1] = row[col][k-1] ;
00767             }
00768             if ( (row[col][k] - row[col][k-1]) <= 2*halfWidth )
00769             {
00770                 dummy_row[k-1] = -1 ;
00771         if (k<n_lines)
00772                   dummy_row[k]   = -1 ;
00773             }
00774             /* the following gives invalid read size 4: check that k+1<lx */
00775             if ( (row[col][k+1] - row[col][k]) <= 2*halfWidth)
00776             {
00777         if (k<n_lines)
00778                   dummy_row[k]   = -1 ;
00779         if (k+1<n_lines)
00780                   dummy_row[k+1] = -1 ;
00781             }
00782         }
00783 
00784         m = 0 ;
00785         for ( k = 0 ; k < j ; k ++ )
00786         {
00787             if ( dummy_row[k] != -1 && dummy_row[k] != 0 )
00788             {
00789                 row_clean[col][m] = dummy_row[k] ;
00790                 wavelength_clean[col][m] = wavelength[col][k] ;
00791                 m ++ ;
00792             }
00793         }
00794 
00795         n_found_lines[col] = m ;
00796 
00797         *sum_lines += n_found_lines[col] ;
00798     }
00799     cpl_free (column) ;
00800     cpl_free (lines) ;
00801     cpl_free (conv_lines) ;
00802     cpl_free (dummy_row) ;
00803     cpl_free (wave_buffer) ;
00804     cpl_free (row_mem) ;
00805     cpl_free (wave_mem) ;
00806     cpl_free (tempcol) ;
00807     cpl_free (row) ;
00808     cpl_free (wavelength) ;
00809 
00810     return 0 ;
00811 }
00812 
00822 int 
00823 sinfo_new_read_list( char * listname, 
00824                      float * lineCenter, 
00825                      float * lineIntensity )
00826 {
00827     FILE * fp ;
00828     int i, n_lines ;
00829     
00830     if ( NULL == lineCenter )
00831     {
00832         sinfo_msg_error(" lineCenter array is not allocated\n") ;
00833         return -1 ;
00834     }
00835 
00836     if ( NULL == lineIntensity )
00837     {
00838         sinfo_msg_error(" lineIntensity array is not allocated\n") ;
00839         return -1 ;
00840     }
00841 
00842     if ( NULL == (fp = fopen ( listname, "r" ) ) )
00843     {
00844         sinfo_msg_error(" cannot open %s\n", listname) ;
00845         return -1 ;
00846     }
00847 
00848     i = 0 ;
00849     while ( fscanf( fp, "%f %f", &lineCenter[i], &lineIntensity[i] ) != EOF )
00850     {
00851         i ++ ;
00852     }
00853     n_lines = i ;
00854     fclose(fp) ;
00855            
00856     return n_lines ;
00857 }
00858 
00859 
00891 int sinfo_new_line_fit ( cpl_image  *  mergedImage, 
00892               FitParams *  par,
00893               float        fwhm,
00894               int          lineInd,
00895               int          column, 
00896               int          halfWidth, 
00897               int          lineRow,
00898               float        min_amplitude,
00899           Vector    *  line,
00900           int       *  mpar,
00901           float     *  xdat,
00902           float     *  wdat )
00903 {
00904     int i, j ;
00905     int iters, xdim, ndat ;
00906     int numpar, its ;
00907     int position ;
00908     float maxval, tol, lab ;
00909     int lx=0;
00910     int ly=0;
00911     float* pdata=NULL;
00912 
00913     if ( mergedImage == NULL )
00914     {
00915         sinfo_msg_error (" no image given as input\n") ;
00916         return -8 ;
00917     }
00918     lx=cpl_image_get_size_x(mergedImage);
00919     ly=cpl_image_get_size_y(mergedImage);
00920     pdata=cpl_image_get_data_float(mergedImage);
00921 
00922 
00923     if ( par == NULL )
00924     {
00925         sinfo_msg_error(" fit parameters not given\n") ;
00926         return -9 ;
00927     }
00928     if ( column < 0 || column > lx )
00929     {
00930         sinfo_msg_error (" wrong column number\n") ;
00931         return -10 ;
00932     }
00933     if ( halfWidth < 0 || halfWidth > ly )
00934     {
00935         sinfo_msg_error (" wrong width given\n") ;
00936         return -11 ;
00937     }
00938     if ( lineRow < 0 || lineRow > ly )
00939     {
00940         sinfo_msg_error (" wrong number of row of the line given\n") ;
00941         return -12 ;
00942     }
00943     if ( min_amplitude < 1. )
00944     {
00945         sinfo_msg_error (" wrong minimum amplitude\n") ;
00946         return -13 ;
00947     }
00948     
00949     /* initialise the Vector */
00950     for ( i = 0 ; i < line -> n_elements ; i++) 
00951     {
00952         line->data[i] = 0;
00953     }
00954     
00955     par -> column = column  ;
00956     par -> line   = lineInd ;
00957 
00958     /* determine the values of the spectral sinfo_vector given as input */
00959     /* go through the chosen column */
00960 
00961     j = 0 ;
00962     for ( i = lineRow-halfWidth ; i <= lineRow+halfWidth ; i++ ) 
00963     {
00964         if ( i < 0 || i >= ly )
00965         {
00966             sinfo_msg_error (" wrong line position or width given\n") ;
00967             return -15 ;
00968         }
00969         else
00970         {
00971             line -> data[j] = pdata[column + i*lx] ;
00972             j ++ ;
00973         }
00974     } 
00975 
00976     /*-------------------------------------------------------------------- 
00977      * go through the spectral sinfo_vector 
00978      * determine the maximum pixel value in the spectral sinfo_vector 
00979      */
00980     maxval = -FLT_MAX ;
00981     position = -INT32_MAX ;
00982     for ( i = 0 ; i < line -> n_elements ; i++ )
00983     {
00984         xdat[i] = i ;
00985         wdat[i] = 1.0 ;
00986         if ( line -> data[i] >= maxval )
00987         {
00988             maxval = line -> data[i] ;
00989             position = i ;
00990         }
00991     }
00992 
00993     /* set initial values for the fitting routine */
00994     xdim     = XDIM ;
00995     ndat     = line -> n_elements ;
00996     numpar   = MAXPAR ;
00997     tol      = TOL ;
00998     lab      = LAB ;
00999     its      = ITS ;
01000     par -> fit_par[1] = fwhm ;
01001     par -> fit_par[2] = (float) position ;
01002     par -> fit_par[3] = (float) (line -> data[0] + 
01003                                  line -> data[line->n_elements - 1]) / 2.0 ;
01004     par -> fit_par[0]  = maxval - (par -> fit_par[3]) ;
01005 
01006     /* exclude low signal cases */
01007     if ( par->fit_par[0] < min_amplitude )
01008     {
01009         cpl_msg_debug ("sinfo_linefit:",
01010                        " sorry, amplitude of line too low to fit: %f",
01011                        par->fit_par[0] ) ;
01012         return -16 ;
01013     }
01014 
01015     for ( i = 0 ; i < MAXPAR ; i++ )
01016     {
01017         par -> derv_par[i] = 0.0 ;
01018         mpar[i] = 1 ;
01019     }
01020 
01021     /* finally, do the least square fit using a sinfo_gaussian */
01022     if ( 0 > ( iters = sinfo_new_lsqfit_c( xdat, &xdim, 
01023                                            line -> data, wdat, 
01024                                            &ndat, par -> fit_par,  
01025                                            par -> derv_par, mpar, 
01026                                            &numpar, &tol, &its, &lab )) )  
01027     {
01028         cpl_msg_debug ("sinfo_linefit:",
01029                        " sinfo_new_lsqfit_c: least squares fit failed,"
01030                        " error no.: %d\n", iters) ; 
01031         return -17 ;
01032     }
01033 
01034     /* correct the fitted position for the given row of the 
01035        line in image coordinates */
01036     par -> fit_par[2] =  (float) (lineRow - halfWidth) + par -> fit_par[2] ;
01037 
01038     /* all was o.k. */
01039     return iters ;
01040 }
01041 
01066 int sinfo_new_fit_lines ( cpl_image  *  line_image, 
01067                FitParams ** allParams,
01068                float        fwhm,
01069                int       *  n_lines, 
01070                int       ** row,
01071                float     ** wavelength, 
01072                int          width,
01073                float        min_amplitude ) 
01074 {
01075     int i, k, l ;
01076     int result ;
01077     Vector * line;
01078     int    * mpar;
01079     float  * xdat, * wdat;
01080     int lx=0;
01081     int ly=0;
01082     float* pdata=NULL;
01083 
01084     if ( line_image == NULL )
01085     {
01086         sinfo_msg_error (" no image given\n") ;
01087         return -18 ;
01088     }
01089     lx=cpl_image_get_size_x(line_image);
01090     ly=cpl_image_get_size_y(line_image);
01091     pdata=cpl_image_get_data_float(line_image);
01092 
01093     if ( n_lines == NULL )
01094     {
01095         sinfo_msg_error (" no counter of emission lines\n") ;
01096         return -19 ;
01097     } 
01098     if ( row == NULL || width <= 0 )
01099     {
01100         sinfo_msg_error (" row or width vectors are empty\n") ;
01101         return -20 ;
01102     }
01103     if ( wavelength == NULL )
01104     {
01105         sinfo_msg_error (" no wavelength array given\n") ;
01106         return -21 ;
01107     }
01108 
01109     k = 0 ;
01110 
01111     /* allocate memory for the spectral sinfo_vector */
01112     line = sinfo_new_vector (2*width + 1) ;
01113     /* allocate memory */
01114     xdat = (float *) cpl_calloc( line -> n_elements, sizeof (float) ) ;
01115     wdat = (float *) cpl_calloc( line -> n_elements, sizeof (float) ) ;
01116     mpar = (int *)   cpl_calloc( MAXPAR, sizeof (int) ) ;
01117     
01118     /* go through the columns */
01119     for ( i = 0 ; i < lx ; i++ )
01120     {
01121         if ( n_lines[i] == 0 )
01122         {
01123             continue ;
01124         }
01125         /* go through the emission lines in a column */
01126         for ( l = 0 ; l < n_lines[i] ; l++ )
01127         {
01128             if ( row[i][l] <= 0 )
01129             {
01130                 continue ;
01131             }
01132 
01133             /* --------------------------------------------------------------
01134              * fit the single lines using sinfo_linefit and store the 
01135                parameters in
01136              * an array of the FitParams data structure allParams[].
01137              */ 
01138             if ( (result = sinfo_new_line_fit ( line_image, 
01139                                                 allParams[k], fwhm, l, i, 
01140                                                 width, row[i][l], 
01141                                                 min_amplitude,line,mpar,
01142                                                 xdat,wdat ) ) < 0 )
01143             {
01144                 cpl_msg_debug ("sinfo_fitLines:",
01145                                " sinfo_linefit failed, error no.: %d,"
01146                                " column: %d, row: %d, line: %d\n", 
01147                                result, i, row[i][l], l) ;
01148                 continue ;
01149             }
01150             if ( (allParams[k] -> fit_par[0] <= 0.) || 
01151                  (allParams[k] -> fit_par[1] <= 0.)
01152                   || (allParams[k] -> fit_par[2] <= 0.) )
01153             {
01154                 sinfo_msg_warning (" negative fit parameters in column: %d,"
01155                                    " line: %d\n", i, l) ;
01156                 continue ;
01157             }
01158             allParams[k] -> wavelength = wavelength[i][l] ;
01159             k++ ;
01160         }
01161     }
01162 
01163     /* free memory */
01164     sinfo_new_destroy_vector(line);
01165     cpl_free(xdat);
01166     cpl_free(wdat);
01167     cpl_free(mpar);
01168     
01169     /* all is o.k. */
01170     return k ;
01171 }    
01172 
01200 float sinfo_new_polyfit( FitParams ** par,
01201                int          column,
01202                int          n_lines,
01203                int          n_rows,
01204                float        dispersion,
01205                float        max_residual,
01206                float *      acoefs,
01207                float *      dacoefs,
01208                int   *      n_reject,
01209                int          n_fitcoefs )
01210 {
01211     float ** ucoefs, ** vcoefs, ** covar ;
01212     float *mem;
01213     float * lambda, * posit ;
01214     float * weight, * resid ;
01215     float * newlam, * newpos, * newwet ;
01216     float * wcoefs=NULL ;
01217     float chisq, result ;
01218     float offset ;
01219     int num, found ;
01220     int i, j, k, n ;
01221 
01222     /* reset the fit coefficients and their errors */
01223     for ( i = 0 ; i < n_fitcoefs ; i++ )
01224     {
01225         acoefs[i]  = 0. ;
01226         dacoefs[i] = 0. ;
01227     }
01228     if ( NULL == par )
01229     {
01230         sinfo_msg_error(" no fit params given\n");
01231         return FLT_MAX ;
01232     }
01233 
01234     if ( 0 >= n_lines )
01235     {
01236       /*
01237         sinfo_msg_warning (" sorry, number of lines is wrong") ;
01238       */
01239         return FLT_MAX ;
01240     }
01241     if ( 0 >= n_rows )
01242     {
01243         sinfo_msg_error (" sorry, number of rows is wrong") ;
01244         return FLT_MAX ;
01245     }
01246     if ( dispersion == 0. )
01247     {
01248         sinfo_msg_error (" sorry, wrong dispersion given") ;
01249         return FLT_MAX ;
01250     }
01251 
01252     offset = (float)(n_rows - 1)/2. ;
01253 
01254     /* allocate memory */
01255     
01256     mem = (float*) cpl_calloc( n_lines*7, sizeof (float) ) ;
01257     lambda = mem;
01258     posit  = mem + n_lines;
01259     weight = mem + n_lines*2;
01260     resid  = mem + n_lines*3;
01261     newlam = mem + n_lines*4;
01262     newpos = mem + n_lines*5;
01263     newwet = mem + n_lines*6;
01264     
01265     /*lambda = (float*) cpl_calloc( n_lines, sizeof (float) ) ;
01266     posit  = (float*) cpl_calloc( n_lines, sizeof (float) ) ;
01267     weight = (float*) cpl_calloc( n_lines, sizeof (float) ) ;
01268     resid  = (float*) cpl_calloc( n_lines, sizeof (float) ) ;
01269     newlam = (float*) cpl_calloc( n_lines, sizeof (float) ) ;
01270     newpos = (float*) cpl_calloc( n_lines, sizeof (float) ) ;
01271     newwet = (float*) cpl_calloc( n_lines, sizeof (float) ) ;*/
01272 
01273     /* allocate coefficient matrices*/
01274     ucoefs = sinfo_matrix ( 1, n_lines, 1, n_fitcoefs ) ;
01275     vcoefs = sinfo_matrix ( 1, n_lines, 1, n_fitcoefs ) ;
01276     covar  = sinfo_matrix ( 1, n_fitcoefs, 1, n_fitcoefs ) ;
01277     wcoefs=cpl_calloc(n_fitcoefs,sizeof(float)) ;
01278 
01279     /* go through all fit parameters */
01280     n = 0 ;
01281     for ( i = 0 ; i < (par[0] -> n_params) ; i ++ )
01282     {
01283         found = -1 ;
01284         /* find the given column and go through the lines in that column */
01285         for ( j = 0 ; j < n_lines ; j ++ )
01286         {
01287             if ( (par[i] -> column == column) && (par[i] -> line == j) )
01288             {
01289                 found = i ;
01290             }
01291             else
01292             {
01293                 continue ;
01294             }
01295 
01296             /* store only fit params with reasonable values */
01297             if ( par[found] -> derv_par[2] != 0. && 
01298                  par[found] -> fit_par[2] > 0. &&
01299                  par[found] -> wavelength > 0. && 
01300                  par[found] -> fit_par[1] > 0. &&
01301                  par[found] -> fit_par[0] > 0. )
01302             {
01303                 /* ----------------------------------------------------------
01304                  * store the found position, error of the position as 
01305                    weight and the associated
01306                  * wavelength values of the fitted lines
01307                  */
01308                 posit[n]  = par[found] -> fit_par[2] ;
01309                 weight[n] = par[found] -> derv_par[2] ;
01310                 lambda[n] = par[found] -> wavelength ;
01311                 n ++ ;
01312             }
01313             else
01314             {
01315                 continue ;
01316             }
01317         }
01318 
01319     }
01320 
01321     num = n ;
01322     if ( num < n_fitcoefs )
01323     {
01324         sinfo_msg_warning("not enough lines found in column %d to "
01325                           "determine the three coefficients.\n", column) ;
01326         for ( i = 0 ; i < n_fitcoefs ; i++ )
01327         {
01328             acoefs[i]  = ZERO ;
01329             dacoefs[i] = ZERO ;
01330         }
01331         sinfo_free_matrix ( ucoefs, 1/*, n_lines*/,    1/*, n_fitcoefs*/ ) ;
01332         sinfo_free_matrix ( vcoefs, 1/*, n_lines*/,    1/*, n_fitcoefs*/ ) ;
01333         sinfo_free_matrix ( covar,  1/*, n_fitcoefs*/, 1/*, n_fitcoefs*/ ) ;
01334         /*cpl_free (lambda) ;
01335         cpl_free (posit) ;
01336         cpl_free (weight) ;
01337         cpl_free (resid) ;
01338         cpl_free (newlam) ;
01339         cpl_free (newpos) ;
01340         cpl_free (newwet) ;*/
01341     cpl_free (mem);
01342         return FLT_MAX ;
01343     }
01344 
01345     /*-------------------------------------------------------------------------
01346      * scale the pixel position values to smaller than 1 and transform 
01347        the weights to wavelength units 
01348      */
01349 
01350     for ( i = 0 ; i < num ; i ++ )
01351     {
01352         posit[i] = (posit[i] - offset)/offset ;
01353         weight[i] *= fabs(dispersion) ;
01354     }
01355 
01356     /* do the fit using the singular value decomposition method */
01357     sinfo_svd_fitting( posit - 1, lambda - 1, 
01358                       weight - 1, num, acoefs-1, n_fitcoefs,
01359             ucoefs, vcoefs, wcoefs-1, covar, &chisq, sinfo_fpol ) ;
01360 
01361     /* scale the linear and the quadratic coefficient */
01362     for ( i = 1 ; i < n_fitcoefs ; i++ )
01363     {
01364         acoefs[i] /= pow(offset, i) ;
01365     }
01366 
01367     /* now that we have determined the fit coefficients, find the residuals */
01368     *n_reject = 0 ;
01369 
01370     j = 0 ;
01371     for ( i = 0 ; i < num ; i++ )
01372     {
01373         result = 0. ;
01374         for ( k = 0 ; k < n_fitcoefs ; k++ )
01375         {
01376             result += acoefs[k] * pow(posit[i], k) ;
01377         }
01378 
01379         resid[i] = lambda[i] - result ;
01380 
01381         if ( fabs( resid[i] ) > max_residual)
01382         {
01383             (*n_reject) ++ ;
01384         }
01385         else
01386         {
01387             newlam[j] = lambda[i] ;
01388             newpos[j] = posit[i] ;
01389             newwet[j] = weight[i] ;
01390             j++ ;
01391         }
01392     }
01393 
01394     num = j ;
01395     if ( num >= n_fitcoefs )
01396     {
01397         sinfo_svd_fitting( newpos - 1, newlam - 1, 
01398                            newwet - 1, num, acoefs-1, n_fitcoefs, ucoefs,
01399                 vcoefs, wcoefs-1, covar, &chisq, sinfo_fpol ) ;
01400 
01401         /* scale the resulting coefficients */
01402         for ( i = 0 ; i < n_fitcoefs ; i++ )
01403         {
01404             acoefs[i] /= pow(offset, i) ;
01405             dacoefs[i] = sqrt( (double) covar[i+1][i+1] ) / pow(offset, i) ;
01406         }
01407     }
01408     else
01409     {
01410         sinfo_msg_warning (" too many lines rejected (number: %d) "
01411                            "due to high residuals, fit coefficients are set "
01412                            "zero, in column: %d\n", *n_reject, column) ;
01413         for ( i = 0 ; i < n_fitcoefs ; i++ )
01414         {
01415             acoefs[i]  = ZERO ;
01416             dacoefs[i] = ZERO ;
01417         }
01418     }
01419 
01420     sinfo_free_matrix ( ucoefs, 1/*, n_lines*/,    1/*, n_fitcoefs*/ ) ;
01421     sinfo_free_matrix ( vcoefs, 1/*, n_lines*/,    1/*, n_fitcoefs*/ ) ;
01422     sinfo_free_matrix ( covar,  1/*, n_fitcoefs*/, 1/*, n_fitcoefs*/ ) ;
01423     /*cpl_free (lambda) ;
01424     cpl_free (posit) ;
01425     cpl_free (weight) ;
01426     cpl_free (resid) ;
01427     cpl_free (newlam) ;
01428     cpl_free (newpos) ;
01429     cpl_free (newwet) ;*/
01430     cpl_free (mem);
01431     cpl_free(wcoefs) ;
01432 
01433     return chisq ;
01434 }
01435 
01452 float sinfo_new_coefs_cross_fit ( int      n_columns,
01453                       float *  acoefs,
01454                       float *  dacoefs,
01455                       float *  bcoefs,
01456                       int      n_fitcoefs,
01457                       float    sigma_factor )
01458 {
01459     float col_index;
01460     float* sub_col_index=NULL ;
01461     float* sub_acoefs=NULL ;
01462     float* sub_dacoefs=NULL ;
01463     float* wcoefs=NULL ;
01464     float ** ucoefs, **vcoefs, **covar ;
01465     float chisq ;
01466     float * acoefsclean ;
01467     double sum, sumq, mean ;
01468     double sigma ;
01469     double cliphi, cliplo ;
01470     float offset ;
01471     int i, n, num, ndata ;
01472     int nc ;
01473 
01474 
01475     if ( n_columns < 1 )
01476     {
01477         sinfo_msg_error(" wrong number of image columns given\n") ;
01478         return FLT_MAX ;
01479     }
01480     if ( acoefs == NULL || dacoefs == NULL )
01481     {
01482         sinfo_msg_error(" coeffs or errors of coefficients are not given\n") ;
01483         return FLT_MAX ;
01484     }
01485     if ( bcoefs == NULL )
01486     {
01487         sinfo_msg_error(" coeffs are not allocated\n") ;
01488         return FLT_MAX ;
01489     }
01490 
01491     if ( n_fitcoefs < 1 )
01492     {
01493         sinfo_msg_error(" wrong number of fit coefficients\n") ;
01494         return FLT_MAX ;
01495     }
01496     if ( sigma_factor <= 0. )
01497     {
01498         sinfo_msg_error(" impossible sigma_factor given!\n") ;
01499         return FLT_MAX ;
01500     }
01501 
01502     offset = (float)(n_columns - 1) / 2. ;
01503 
01504     /* ----------------------------------------------------------
01505      * determine the clean mean and sigma value of the coefficients,
01506      * that means reject 10 % of the extreme low and high values
01507      */
01508 
01509     wcoefs=cpl_calloc(n_fitcoefs,sizeof(float)) ;
01510 
01511     nc = 0 ;
01512     for ( i = 0 ; i < n_columns ; i++ )
01513     {
01514         if ( isnan(acoefs[i]) || acoefs[i] == 0. || dacoefs[i] == 0. )
01515         {
01516             continue ;
01517         }
01518         else
01519         {
01520             nc++ ;
01521         }
01522     }
01523     acoefsclean = (float*) cpl_calloc(nc , sizeof(float)) ;
01524     nc = 0 ;
01525     for ( i = 0 ; i < n_columns ; i++ )
01526     {
01527         if ( isnan(acoefs[i]) || acoefs[i] == 0. || dacoefs[i] == 0. )
01528         {
01529             continue ;
01530         }
01531         else
01532         {
01533             acoefsclean[nc] = acoefs[i] ;
01534             nc++ ;
01535         }
01536     }
01537     sinfo_pixel_qsort(acoefsclean, nc) ;
01538     sum   = 0. ;
01539     sumq  = 0. ;
01540     mean  = 0. ;
01541     sigma = 0. ;
01542     n     = 0 ;
01543     for ( i = (int)((float)nc*LOW_REJECT) ; 
01544           i < (int)((float)nc*HIGH_REJECT) ; i++ )
01545     {
01546         sum  += (double)acoefsclean[i] ;
01547         sumq += ((double)acoefsclean[i] * (double)acoefsclean[i]) ;
01548         n ++ ;
01549     }
01550     mean          = sum/(double)n ;
01551     sigma         = sqrt( sumq/(double)n - (mean * mean) ) ;
01552     cliphi        = mean + sigma * (double)sigma_factor ;
01553     cliplo        = mean - sigma * (double)sigma_factor ;
01554 
01555     sub_col_index=cpl_calloc(n_columns,sizeof(float)) ;
01556     sub_acoefs=cpl_calloc(n_columns,sizeof(float));
01557     sub_dacoefs=cpl_calloc(n_columns,sizeof(float)) ;
01558 
01559     /* fit only the reasonnable values */
01560     num = 0 ;
01561     for ( i = 0 ; i < n_columns ; i++ )
01562     {
01563         /* associate the column indices to the corresponding array */
01564         col_index = (float) i ;
01565 
01566         /* take only the reasonnable coefficients */
01567         if ( !isnan(acoefs[i]) && 
01568              (acoefs[i] <= cliphi) && (acoefs[i] >= cliplo) &&
01569              (dacoefs[i] != 0. ) && (acoefs[i] != 0.) )
01570         {
01571             sub_acoefs[num]    = acoefs[i] ;
01572             sub_dacoefs[num]   = dacoefs[i] ;
01573             sub_col_index[num] = col_index ;
01574             num ++ ;
01575         }
01576     }
01577     ndata = num ;
01578 
01579     if ( ndata < n_fitcoefs )
01580     {
01581         sinfo_msg_error("not enough data found to determine "
01582                         "the fit coefficients.\n") ;
01583 
01584         return FLT_MAX ;
01585     }
01586 
01587     /* allocate coefficient matrices */
01588     ucoefs = sinfo_matrix(1, ndata, 1, n_fitcoefs) ;
01589     vcoefs = sinfo_matrix(1, ndata, 1, n_fitcoefs) ;
01590     covar  = sinfo_matrix ( 1, n_fitcoefs, 1, n_fitcoefs ) ;
01591 
01592     /* scale the x-values for the fit */
01593     for ( i = 0 ; i < ndata ; i++ )
01594     {
01595         sub_col_index[i] = (sub_col_index[i] - offset) / offset ;
01596     }
01597 
01598     /* finally, do the singular value decomposition fit */
01599     sinfo_svd_fitting ( sub_col_index-1, sub_acoefs-1, 
01600                         sub_dacoefs-1, ndata, bcoefs-1,
01601                         n_fitcoefs, ucoefs, vcoefs, 
01602                         wcoefs-1, covar, &chisq, sinfo_fpol ) ;
01603 
01604     /* scale the found coefficients */
01605     for ( i = 0 ; i < n_fitcoefs ; i ++ )
01606     {
01607         bcoefs[i] /= pow(offset, i) ;
01608     }
01609 
01610     /* free memory */
01611     cpl_free (acoefsclean) ;
01612     sinfo_free_matrix( ucoefs, 1/*, ndata*/, 1/*, n_fitcoefs */) ;
01613     sinfo_free_matrix( vcoefs, 1/*, ndata*/, 1/*, n_fitcoefs */) ;
01614     sinfo_free_matrix ( covar, 1/*, n_fitcoefs*/, 1/*, n_fitcoefs*/ ) ;
01615 
01616     cpl_free(sub_col_index) ;
01617     cpl_free(sub_acoefs) ;
01618     cpl_free(sub_dacoefs) ;
01619     cpl_free(wcoefs) ;
01620 
01621     return chisq ;
01622 }
01623 
01624 
01644 cpl_image * sinfo_new_wave_map( cpl_image * lineImage,
01645                     float   ** bcoefs,
01646                     int        n_a_fitcoefs,
01647                     int        n_b_fitcoefs,
01648                     float    * wavelength,
01649                     float    * intensity,
01650                     int        n_lines,
01651                     int        magFactor)
01652 {
01653     cpl_image * retImage ;
01654     float cenpos, cenpix ;
01655     float centreval, centrepix, wavelag ;
01656     float  pixvalue ;
01657     float a_initial ;
01658     int i, j, k, l/*, m*/, line, col, row, found, sign ;
01659     int var, maxlag, cmin, cmax, offset ;
01660     double * result ;
01661     float col_off ;
01662     float angst ;
01663     double xcorr_max ;
01664     int delta ;
01665 
01666     double* z=NULL ;
01667     double* a=NULL ;
01668     double*  wave=NULL ;
01669     float* emline=NULL ;
01670     float* spec=NULL ;
01671     int ilx=0;
01672     int ily=0;
01673     int olx=0;
01674     int oly=0;
01675     float* pidata=NULL;
01676     float* podata=NULL;
01677 
01678 
01679     gsl_poly_complex_workspace * w ;
01680     
01681     if ( NULL == lineImage )
01682     {
01683         sinfo_msg_error("no image given\n") ;
01684         return NULL ;
01685     }
01686     ilx=cpl_image_get_size_x(lineImage);
01687     ily=cpl_image_get_size_y(lineImage);
01688     pidata=cpl_image_get_data_float(lineImage);
01689 
01690     if ( NULL == wavelength || n_lines <= 0 )
01691     {
01692         sinfo_msg_error("no wavelength list given\n") ;
01693         return NULL ;
01694     }
01695 
01696     if ( NULL == intensity )
01697     {
01698         sinfo_msg_error("no intensity list given\n") ;
01699         return NULL ;
01700     }
01701 
01702     if ( NULL == bcoefs )
01703     {
01704         sinfo_msg_error("no coefficients given\n") ;
01705         return NULL ;
01706     }
01707 
01708     if ( magFactor <= 1 )
01709     {
01710         sinfo_msg_error("wrong magnifying factor given\n") ;
01711         return NULL ;
01712     }
01713     
01714     /* allocate memory */
01715     if ( NULL == ( retImage = cpl_image_new( ilx, ily,CPL_TYPE_FLOAT ) ))
01716     {
01717         sinfo_msg_error("cannot allocate a new image\n");
01718         return NULL ;
01719     }
01720     olx=cpl_image_get_size_x(retImage);
01721     oly=cpl_image_get_size_y(retImage);
01722     podata=cpl_image_get_data_float(retImage);
01723 
01724 
01725     var    = (magFactor - 1)*(magFactor - 1) ;
01726     offset = ily * (magFactor/4 + 1) ;
01727 
01728     /* find out if Angstroem or microns are used */
01729     if ( wavelength[0] > 10000. )
01730     {
01731     /* Angstroem */
01732         angst = 10000. ;
01733     }
01734     else if ( wavelength[0] > 1000. && wavelength[0] < 10000. )
01735     {
01736     /* nanometers */
01737     angst = 1000. ;
01738     }
01739     else
01740     {
01741     /* microns */
01742     angst = 1. ;
01743     }
01744 
01745     z=cpl_calloc(2*(n_a_fitcoefs - 1),sizeof(double)) ;
01746     a=cpl_calloc(n_a_fitcoefs,sizeof(double));
01747     wave=cpl_calloc(n_lines,sizeof(double)) ;
01748     emline=cpl_calloc(2*magFactor*ily,sizeof(float));
01749     spec=cpl_calloc(2*magFactor*ily,sizeof(float)) ;
01750 
01751     /* go through the image columns */
01752     for ( col = 0 ; col < ilx ; col++ )
01753     {
01754         /* initialize the emline array for each column */
01755         for ( i = 0 ; i < 2*magFactor*ily ; i++ )
01756         {
01757             emline[i] = 0. ;
01758         }
01759         col_off = (float)col - (float)(ilx-1)/2. ;
01760 
01761         /* determine the coefficients by using the given bcoefs */
01762         for ( i = 0 ; i < n_a_fitcoefs ; i++ ) 
01763         {
01764             /* initialize coefficients and solution */
01765             a[i] = 0. ;
01766             if (i < n_a_fitcoefs-1)
01767             {
01768                 z[2*i] = 0. ;
01769                 z[2*i+1] = 0. ;
01770             }
01771             for ( j = 0 ; j < n_b_fitcoefs ; j++ )
01772             {
01773                 a[i] += bcoefs[i][j] * pow(col_off, j) ;
01774             }
01775         }
01776         a_initial = a[0] ;
01777         
01778         /* go through the lines and generate an artificial spectrum */
01779         for ( line = 0 ; line < n_lines ; line++ )
01780         {
01781             /* go from Angstroem to micron */
01782             wave[line] = wavelength[line]/angst ;
01783 
01784             /* ---------------------------------------------------------------
01785              * solve the polynomial for the exact offset of the line that means
01786              * find the root of the polynomial of order n_fitcoefs - 1
01787              */
01788             a[0] = a_initial - wave[line] ;
01789 
01790             if(NULL==(w=sinfo_gsl_poly_complex_workspace_alloc(n_a_fitcoefs)))
01791             {
01792                 sinfo_msg_error("could not allocate complex workspace!") ;
01793                 cpl_image_delete(retImage) ;
01794                 return NULL ;
01795             }
01796             if (-1 == sinfo_gsl_poly_complex_solve(a, n_a_fitcoefs, w, z))
01797             {
01798                 sinfo_msg_error("sinfo_gsl_poly_complex_solve did not work!") ;
01799                 cpl_image_delete(retImage) ;
01800                 return NULL ;
01801             }
01802             sinfo_gsl_poly_complex_workspace_free(w) ;
01803 
01804             
01805             j = 0 ;
01806             found = -1 ;
01807             for ( i = 0 ; i < n_a_fitcoefs - 1  ; i++ )
01808             {
01809                 /* test for appropriate solution */
01810                 if( z[2*i] > (-1.)*(float) ily/2. && 
01811                     z[2*i] < (float)ily/2. && z[2*i+1] == 0. )
01812                 {
01813                     found = 2*i ;
01814                     j ++ ;
01815                 }
01816                 else
01817                 {
01818                     continue ;
01819                 } 
01820             }
01821             if ( j == 0 )
01822             {
01823                 sinfo_msg_warning("no offset solution found "
01824                                   "for line %d in column %d\n", line, col) ;
01825                 continue ;
01826             } 
01827             else if ( j == 1 )
01828             {
01829                 cenpos = z[found] + (float) ily /2. ;
01830             }
01831             else
01832             {
01833                 sinfo_msg_warning("two or more offset solutions found "
01834                                   "for line %d in column %d\n", line, col) ;
01835                 continue ;
01836             }
01837              
01838             /*---------------------------------------------------------------
01839              * magnify image by the given factor add an additional offset 
01840              */
01841             cenpix = cenpos * (float) magFactor + (float) offset ;  
01842             
01843             /* determine max and min pixel limits over 
01844                which line should be convolved */
01845             cmin = (sinfo_new_nint(cenpix) - (var-1)) > 0 ? 
01846                     sinfo_new_nint(cenpix) - (var-1) : 0 ;
01847             cmax = (sinfo_new_nint(cenpix) + (var-1)) < 2*magFactor * ily ? 
01848                     sinfo_new_nint(cenpix) + (var-1) :  2*magFactor * ily ;
01849 
01850             /* convolve neon lines with Gaussian function */
01851             for ( j = cmin ; j < cmax ; j++ )
01852             {
01853                 emline[j] += intensity[line] * 
01854                 exp((double)(-0.5*(j-cenpix)*(j-cenpix))/(double)var) ;
01855             }
01856         }
01857                 
01858         /*--------------------------------------------------------------------- 
01859          * for each column, map the image data points onto an magFactor times 
01860            bigger element grid for FFT in the cross sinfo_new_correlation, 
01861            first initialize the two helping arrays for each new column.
01862          */
01863         for ( k = 0 ; k < 2*magFactor * ily ; k++ )
01864         {
01865             spec[k] =  0. ;
01866         }
01867    
01868         /* now take the image data points of the column and put them 
01869            into the spec array */
01870         for ( row = 0 ; row < ily ; row++ ) /* go through the column */
01871         {
01872             /* insert 8 values for each image row (magnification) and 
01873                add same offset as for emline array */
01874             for ( l = 0 ; l < magFactor ; l++ )   
01875         {
01876             /* set bad pixels or negative values to zero */
01877                 if (!isnan(pidata[col + row * ilx]) &&
01878                                 (pidata[col + row * ilx] > 0.))
01879                 {
01880                     spec[offset + l + (row * magFactor)] = 
01881                 pidata[col + row * ilx] ;     
01882                 }
01883                 else
01884                 {
01885                     spec[offset + l + (row * magFactor)] = 0. ;
01886                 }
01887             }
01888     }
01889  
01890         /* now call the cross sinfo_new_correlation routine */
01891         if (NULL == (result = sinfo_new_xcorrel(spec, 2*magFactor * ily, 
01892                                                 emline, 2*magFactor * ily, 
01893                                                 magFactor * ily, &delta, 
01894                                                 &maxlag, &xcorr_max)) ) 
01895         {
01896         sinfo_msg_warning ("cross sinfo_new_correlation did not work,"
01897                                " col: %d is set ZERO\n", col) ;
01898             for ( row = 0 ; row < ily ; row++ )
01899             {
01900                 podata[col + row * ilx] = ZERO ;
01901             }
01902             continue ;
01903         }
01904     
01905         if ( xcorr_max <= 0. )
01906         {
01907             sinfo_msg_warning ("cross sinfo_new_correlation sum is negative,"
01908                                " col: %d is set ZERO\n", col) ;
01909             for ( row = 0 ; row < ily ; row++ )
01910             {
01911                 podata[col + row * ilx] = ZERO ;
01912             }
01913             cpl_free(result) ;
01914             continue ;
01915         }
01916 
01917         wavelag = (float) -delta / (float) magFactor ;
01918         if ( fabs(wavelag) > (float)ily/20. )
01919         {
01920             sinfo_msg_warning("wave lag too big, col: %d is set ZERO\n", col) ;
01921             for ( row = 0 ; row < ily ; row++ )
01922             {
01923                 podata[col + row * ilx] = ZERO ;
01924             }
01925             cpl_free(result) ;
01926         continue ;
01927         }
01928 
01929         /*-------------------------------------------------------------------- 
01930          * determine new zero order coefficient centreval, of which the 
01931            formula is determined by setting equal a polynomial shifted by 
01932            wavelag with the same higher order coefficients and set the new 
01933            zero order coefficient to get both sides of the equation 
01934            approximately equal.
01935          */ 
01936         centreval = a_initial ;
01937         for ( i = 1 ; i < n_a_fitcoefs ; i++ )
01938         {
01939             if ( i%2 == 0 )
01940             {
01941                 sign = -1 ;
01942             }
01943             else
01944             {
01945                 sign = 1 ;
01946             }
01947             centreval += (float)sign * a[i]*pow(wavelag, i) ;
01948         }
01949 
01950         /* prepare to write out wavelength as pixel values */
01951         for ( row = 0 ; row < ily ; row++ )
01952         {
01953             centrepix = (float)row - ((float)ily - 1.)/2. ;
01954             pixvalue = 0. ;
01955             for ( i = 1 ; i < n_a_fitcoefs ; i++ )
01956             {
01957                 pixvalue += a[i]*pow(centrepix, i) ; 
01958             }
01959             podata[col + row * ilx] = centreval + pixvalue ; 
01960         }
01961         cpl_free(result) ;   
01962     }
01963 
01964 
01965 
01966     cpl_free(z) ;
01967     cpl_free(a) ;
01968     cpl_free(wave) ;
01969     cpl_free(emline) ;
01970     cpl_free(spec) ;
01971 
01972     return retImage ;
01973 }                    
01974 
02019 int sinfo_new_wavelength_calibration( cpl_image   * image,
02020                            FitParams ** par ,
02021                            float     ** bcoefs,
02022                            float      * wave,
02023                            int          n_lines,
02024                            int       ** row_clean,
02025                            float     ** wavelength_clean,
02026                            int        * n_found_lines,
02027                            float        dispersion,
02028                            int          halfWidth,
02029                            float        minAmplitude,
02030                            float        max_residual,
02031                            float        fwhm,
02032                            int          n_a_fitcoefs,
02033                            int          n_b_fitcoefs,
02034                            float        sigmaFactor,
02035                            float        pixel_tolerance )
02036 
02037 {
02038     int          i, j, k ;
02039     int          n_fit ;
02040     int          n_reject ;
02041     float     *  acoefs ;
02042     float     *  dacoefs ;
02043     float     ** abuf ;
02044     float     ** dabuf ;
02045     float        chisq_poly, chisq_cross ;
02046     int          zeroind ;
02047     /*float     *  mem ;*/
02048     int lx=0;
02049     int ly=0;
02050     float* pdata=NULL;
02051 
02052     if (  NULL == image )
02053     {
02054         sinfo_msg_error("no image given\n") ;
02055         return -1 ;
02056     }
02057     lx=cpl_image_get_size_x(image);
02058     ly=cpl_image_get_size_y(image);
02059     pdata=cpl_image_get_data_float(image);
02060 
02061     if ( par == NULL )
02062     {
02063         sinfo_msg_error("no fit parameter data structure given\n") ;
02064         return -1 ;
02065     }
02066     if ( wave == NULL )
02067     {
02068         sinfo_msg_error("no wavelength list given\n") ;
02069         return -1 ;
02070     }
02071     if ( n_lines <= 0 )
02072     {
02073         sinfo_msg_error("impossible number of lines in line list given\n") ;
02074         return -1 ;
02075     }
02076     if ( row_clean == NULL )
02077     {
02078         sinfo_msg_error("no row_clean array given\n") ;
02079         return -1 ;
02080     }
02081     if ( wavelength_clean == NULL )
02082     {
02083         sinfo_msg_error("no wavelength_clean array given\n") ;
02084         return -1 ;
02085     }
02086 
02087     if ( dispersion == 0. )
02088     {
02089         sinfo_msg_error("impossible dispersion given\n") ;
02090         return -1 ;
02091     }
02092 
02093     if ( halfWidth <= 0 || halfWidth > ly/2 )
02094     {
02095         sinfo_msg_error("impossible half width of the fitting box given\n") ;
02096         return -1 ;
02097     }
02098     if ( minAmplitude < 1. )
02099     {
02100         sinfo_msg_error("impossible minimal amplitude\n") ;
02101         return -1 ;
02102     }
02103 
02104     if ( max_residual <= 0. || max_residual > 1. )
02105     {
02106         sinfo_msg_error("impossible max_residual given\n") ;
02107         return -1 ;
02108     }
02109 
02110     if ( fwhm <= 0. || fwhm > 10. )
02111     {
02112         sinfo_msg_error("impossible fwhm given\n") ;
02113 
02114         return -1 ;
02115     }
02116 
02117     if ( n_a_fitcoefs <= 0 || n_a_fitcoefs > 9 )
02118     {
02119         sinfo_msg_error("unrealistic n_a_fitcoefs given\n") ;
02120         return -1 ;
02121     }
02122 
02123     if ( n_b_fitcoefs <= 0 || n_b_fitcoefs > 9 )
02124     {
02125         sinfo_msg_error(" unrealistic n_b_fitcoefs given\n") ;
02126         return -1 ;
02127     }
02128     if ( sigmaFactor <= 0. )
02129     {
02130         sinfo_msg_error(" impossible sigmaFactor given\n") ;
02131         return -1 ;
02132     }
02133 
02134     /* initialize the variables */
02135     n_reject = 0 ;
02136     n_fit = 0 ;
02137 
02138     /* fit each found line by using a Gaussian function and determine the 
02139        exact position */
02140     if ( 0 > (n_fit = sinfo_new_fit_lines( image , par, fwhm, 
02141                                            n_found_lines, row_clean, 
02142                                            wavelength_clean,
02143                                 halfWidth, minAmplitude )) )
02144     {
02145         sinfo_msg_error(" cannot fit the lines, "
02146                         "error code of sinfo_fitLines: %d\n", n_fit) ;
02147         return -1 ;
02148     }
02149 
02150     /* first check for faked lines like bad pixels */
02151     if ( -1 == sinfo_new_check_for_fake_lines (par, dispersion, 
02152                                                wavelength_clean, row_clean, 
02153                                                n_found_lines,
02154                                                lx, pixel_tolerance) )
02155     {
02156         sinfo_msg_error("cannot fit the lines, "
02157                         "error code of sinfo_fitLines: %d", n_fit) ;
02158         return -1 ;
02159     }
02160 
02161     /* allocate memory */
02162     if (NULL == (acoefs = (float*) cpl_calloc (n_a_fitcoefs, sizeof(float*))) ||
02163         NULL == (dacoefs = (float*)cpl_calloc (n_a_fitcoefs, sizeof(float*))) ||
02164         NULL == (abuf = (float**) cpl_calloc (n_a_fitcoefs, sizeof(float*))) ||
02165         NULL == (dabuf = (float**) cpl_calloc (n_a_fitcoefs, sizeof(float*))) )
02166     {
02167         sinfo_msg_error(" cannot allocate memory\n") ;
02168         return -1 ;
02169     }
02170 
02171     for ( i = 0 ; i < n_a_fitcoefs ; i++ )
02172     {
02173         if ( NULL == (abuf[i] = (float*) cpl_calloc(lx, sizeof(float))) ||
02174              NULL == (dabuf[i] = (float*) cpl_calloc(lx, sizeof(float))) )
02175         {
02176             sinfo_msg_error(" cannot allocate memory\n") ;
02177             cpl_free(abuf) ;
02178             cpl_free(dabuf) ;
02179             return -1 ;
02180         }
02181     }
02182 
02183     /* fit wavelengths to the corresponding found positions for each column */
02184     k = 0 ;
02185     
02186     for ( i = 0 ; i < lx ; i++ )
02187     {
02188         zeroind = 0 ;
02189         if ( FLT_MAX == (chisq_poly = sinfo_new_polyfit( par, i, 
02190                                                          n_found_lines[i], 
02191                                                          ly, dispersion,
02192                                                          max_residual, acoefs, 
02193                                                          dacoefs, &n_reject, 
02194                                                          n_a_fitcoefs)) )
02195         {
02196       /* 
02197            sinfo_msg_warning (" error in polyfitt in column: %d\n", i) ;
02198        */
02199             for ( j = 0 ; j < n_a_fitcoefs ; j++ )
02200             {
02201                 acoefs[j] = ZERO ;
02202                 dacoefs[j] = ZERO ;
02203             }
02204         }
02205 
02206         for ( j = 0 ; j < n_a_fitcoefs ; j++ )
02207         {
02208             if ( acoefs[0] <= 0. || acoefs[1] ==0. ||
02209                  dacoefs[j] == 0. || isnan(acoefs[j]) )
02210             {
02211                 zeroind = 1 ;
02212 
02213             }
02214         }
02215         for ( j = 0 ; j < n_a_fitcoefs ; j++ )
02216         {
02217             if ( zeroind == 0 )
02218             {
02219                 abuf[j][i]  = acoefs[j] ;
02220                 dabuf[j][i] = dacoefs[j] ;
02221             }
02222             else
02223             {
02224                 abuf[j][i]  = ZERO ;
02225                 dabuf[j][i] = ZERO ;
02226             }
02227         }
02228     }
02229 
02230     /* fit each acoefs across the columns to smooth the result */
02231     for ( i = 0 ; i < n_a_fitcoefs ; i++ )
02232     {
02233         if ( FLT_MAX == (chisq_cross = sinfo_new_coefs_cross_fit(lx, 
02234                                                                  abuf[i], 
02235                                                                  dabuf[i],
02236                                                                  bcoefs[i],
02237                                                                  n_b_fitcoefs,
02238                                                                  sigmaFactor)))
02239         {
02240             sinfo_msg_error (" cannot carry out the fitting of coefficients"
02241                              " across the columns, for the coefficient with"
02242                              " index: %d\n", i) ;
02243             for ( i = 0 ; i < n_a_fitcoefs ; i++ )
02244             {
02245                 cpl_free (abuf[i]) ;
02246                 cpl_free (dabuf[i]) ;
02247             }
02248             cpl_free ( acoefs ) ;
02249             cpl_free ( dacoefs ) ;
02250             cpl_free ( abuf ) ;
02251             cpl_free ( dabuf ) ;
02252             return -1 ;
02253         }
02254     }
02255 
02256     /* free all allocated memory */
02257     for ( i = 0 ; i < n_a_fitcoefs ; i++ )
02258     {
02259         cpl_free (abuf[i]) ;
02260         cpl_free (dabuf[i]) ;
02261     }
02262     cpl_free ( acoefs ) ;
02263     cpl_free ( dacoefs ) ;
02264     cpl_free ( abuf ) ;
02265     cpl_free ( dabuf ) ;
02266 
02267     return 0 ;   
02268 }
02269 
02270                            
02282 cpl_image * sinfo_new_convolve_image_by_gauss( cpl_image * lineImage,
02283                                  int        hw )
02284 {
02285     cpl_image * returnImage ;
02286     float* column_buffer=NULL ;
02287     float * filter ;
02288     int col, row ;
02289     int ilx=0;
02290     int ily=0;
02291     int olx=0;
02292     int oly=0;
02293     float* pidata=NULL;
02294     float* podata=NULL;
02295 
02296     if ( lineImage == NULL )
02297     {
02298         sinfo_msg_error(" no input image given!\n") ;
02299         return NULL ;
02300     }
02301     ilx=cpl_image_get_size_x(lineImage);
02302     ily=cpl_image_get_size_y(lineImage);
02303     pidata=cpl_image_get_data_float(lineImage);
02304 
02305     if ( hw < 1 )
02306     {
02307         sinfo_msg_error(" wrong half width given!\n") ;
02308         return NULL ;
02309     }
02310 
02311     /* allocate memory for returned image */
02312     if ( NULL == ( returnImage = cpl_image_new(ilx,ily,CPL_TYPE_FLOAT ) ))
02313     {
02314         sinfo_msg_error(" cannot allocate a new image\n");
02315         return NULL ;
02316     }
02317     olx=cpl_image_get_size_x(returnImage);
02318     oly=cpl_image_get_size_y(returnImage);
02319     podata=cpl_image_get_data_float(returnImage);
02320 
02321     column_buffer=cpl_calloc(ily,sizeof(float)) ;
02322 
02323     /* go through the image columns and save them in a buffer */
02324     for ( col = 0 ; col < ilx ; col++ )
02325     { 
02326         for ( row = 0 ; row < ily ; row++ )
02327         {
02328             column_buffer[row] = pidata[col + row*ilx] ;
02329         }
02330          
02331         /*--------------------------------------------------------------------- 
02332          * now low pass filter the columns by the sinfo_gaussian and fill 
02333            the return image.
02334          */  
02335         filter = sinfo_function1d_filter_lowpass( column_buffer,
02336                                             ily,
02337                                             LOW_PASS_GAUSSIAN,
02338                                             hw ) ;
02339         for ( row = 0 ; row < ily ; row++ )
02340         {
02341             podata[col + row*ilx] = filter[row] ;
02342         }
02343         sinfo_function1d_del(filter) ;
02344     }
02345                         
02346     cpl_free(column_buffer);
02347     return returnImage ;
02348 }
02349 
02389 cpl_image * sinfo_new_defined_resampling( cpl_image * image,
02390                               cpl_image * calimage,
02391                               int        n_params,
02392                               int*       n_rows,
02393                               double   * dispersion,
02394                               float    * minval,
02395                               float    * maxval,
02396                               double   * centralLambda,
02397                               int    * centralpix )
02398 {
02399     cpl_image * retImage ;
02400     cpl_image * tempCalImage ;
02401     cpl_image * tempImage ;
02402     float lambda ;
02403     float dif, lambda_renorm ;
02404     float retimagecol[2560] ; /* retimagecol[n_rows] ; */
02405 
02406     float* imagecol=NULL ;
02407     float* calcol=NULL ;
02408     float* x_renorm=NULL ;
02409 
02410     float * imageptr ;
02411     float sum, new_sum ;
02412     float disp, mindisp ;
02413     int calcolpos[2560];
02414     int i/*, j*/, col, row, testrow ;
02415     int half_width, firstpos ;
02416     int dispInd ;
02417     int n ;
02418     int flag;
02419     float temprow;
02420     float minLambda = 0. ;
02421     /*dpoint list[n_params] ;*/
02422     /*double * polycoeffs ;*/
02423     double poly ;
02424     /*float error;*/
02425     int zeroind ;
02426     int ilx=0;
02427     int ily=0;
02428     int clx=0;
02429     int cly=0;
02430     int olx=0;
02431     int oly=0;
02432 
02433     float* podata=NULL;
02434     float* pidata=NULL;
02435     float* pcdata=NULL;
02436     float* ptidata=NULL;
02437     float* ptcdata=NULL;
02438 
02439 
02440     if ( NULL == image )
02441     {
02442         sinfo_msg_error(" source image not given\n") ;
02443         return NULL ;
02444     }
02445     ilx=cpl_image_get_size_x(image);
02446     ily=cpl_image_get_size_y(image);
02447     pidata=cpl_image_get_data_float(image);
02448 
02449 
02450     if ( NULL == calimage )
02451     {
02452         sinfo_msg_error(" wavelength map image not given\n") ;
02453         return NULL ;
02454     }
02455     clx=cpl_image_get_size_x(calimage);
02456     cly=cpl_image_get_size_y(calimage);
02457     pcdata=cpl_image_get_data_float(calimage);
02458 
02459     if ( ilx != clx ||
02460          ily != cly )
02461     {
02462         sinfo_msg_error("source image and wavelength map image "
02463                         "are not compatible in size\n") ;
02464         return NULL ;
02465     }                              
02466   
02467     if ( n_params < 1 )
02468     {
02469         sinfo_msg_error (" wrong number of fit parameters given\n") ;
02470         return NULL ;
02471     }
02472 
02473     if ( n_params > 4 )
02474     {
02475         sinfo_msg_warning(" attention: very high number of fit "
02476                           "parameters given, not tested !!!\n") ;
02477     }
02478 
02479 
02480     imagecol=cpl_calloc(ily,sizeof(float)) ;
02481     calcol=cpl_calloc(cly,sizeof(float)) ;
02482     x_renorm=cpl_calloc(n_params,sizeof(float)) ;
02483 
02484 
02485     /*if ( n_rows <= cly)
02486     {
02487         sinfo_msg_error (" number of rows of resampled image will be "
02488                          " smaller than in wavelength calibration map,"
02489                          " information would get lost!") ;
02490         return NULL ;
02491     }*/
02492    
02493     dispInd = 0 ;
02494     /* first determine the dispersion direction */
02495     for ( col = 0 ; col < clx ; col++ )
02496     {
02497         if ( isnan(pcdata[col]) || pcdata[col] <= 0. )
02498         {
02499             continue ;
02500         }
02501         if ((pcdata[col] - pcdata[col+(clx)*(cly-1)]) > 0. )
02502         {
02503            dispInd--  ;
02504         }
02505         else if ((pcdata[col] - pcdata[col+(clx)*(cly-1)]) < 0. )
02506         {
02507            dispInd++ ;
02508         }
02509         else
02510         {
02511             continue ;
02512         }
02513     }
02514 
02515     if ( dispInd == 0 )
02516     {
02517         sinfo_msg_error(" zero dispersion?\n");
02518         return NULL ;
02519     }
02520    
02521     /* mirror the wavelength map and the raw image if 
02522        the dispersion is negative */
02523     if ( dispInd < 0 )
02524     {
02525         /* allocate a temp image */
02526         if ( NULL == ( tempCalImage = cpl_image_new(clx,cly,CPL_TYPE_FLOAT)))
02527         {
02528             sinfo_msg_error(" cannot allocate a new image\n");
02529             return NULL ;
02530         }
02531         ptcdata=cpl_image_get_data_float(tempCalImage);
02532         if ( NULL == ( tempImage = cpl_image_new( ilx, ily,CPL_TYPE_FLOAT)))
02533         {
02534             sinfo_msg_error(" cannot allocate a new image\n");
02535             cpl_image_delete(tempCalImage) ;
02536             return NULL ;
02537         }
02538         ptidata=cpl_image_get_data_float(tempImage);
02539 
02540         for ( col = 0 ; col < clx ; col++ )
02541         {
02542             n = cly - 1 ;
02543             for ( row = 0 ; row < cly ; row++ )
02544             {
02545                 ptcdata[col+row*clx] = pcdata[col+n*clx] ;
02546                 ptidata[col+row*clx] = pidata[col+n*clx] ;
02547                 n-- ;
02548             }
02549         }
02550         for ( i = 0 ; i < (int) ilx*ily ; i++ )
02551         {
02552             pidata[i] = ptidata[i] ;
02553             pcdata[i] = ptcdata[i] ;
02554         }
02555         cpl_image_delete(tempCalImage) ;
02556         cpl_image_delete(tempImage) ;
02557     }
02558    
02559     /* determine the max and min pixel value in the first and the last row */
02560     *maxval = -FLT_MAX ;
02561     *minval =  FLT_MAX ;
02562     mindisp = FLT_MAX ;
02563     for ( col = 0 ; col < clx ; col++ )
02564     {
02565         if ( isnan(pcdata[col]) || pcdata[col] <= 0. )
02566         {
02567            continue ;
02568         }
02569         disp = (pcdata[col+(clx)*((cly)-1)]
02570                - pcdata[col]) / (float)cly ;
02571         if ( mindisp > disp )
02572         {
02573             mindisp = disp ;
02574         }
02575         if ( *minval >= pcdata[col] )
02576         {
02577             *minval = pcdata[col] ;
02578         }
02579         if ( *maxval <= pcdata[col + (clx)*((cly)-1)] )
02580         {
02581             *maxval = pcdata[col + (clx)*((cly)-1)] ;
02582         }
02583     }
02584     /* find the used grating and set the dispersion to the defined value */
02585     if (*minval > 1.9 )
02586     {
02587         if ( cly > 1024 && cly < 3000)
02588         {
02589             *dispersion = DISPERSION_K_DITH ;
02590             *centralLambda = CENTRALLAMBDA_K ;
02591         }
02592         else if ( cly < 2000)
02593         {
02594             *dispersion = DISPERSION_K ;
02595             *centralLambda = CENTRALLAMBDA_K ;
02596         }
02597         else
02598         {
02599             *dispersion = DISPERSION_K_DITH/2 ;
02600             *centralLambda = CENTRALLAMBDA_K ;
02601         }
02602     }
02603     else if (*minval < 1.2 )
02604     {
02605         if ( cly > 1024 )
02606         {
02607             *dispersion = DISPERSION_J_DITH ;
02608             *centralLambda = CENTRALLAMBDA_J ;
02609         }
02610         else
02611         {
02612             *dispersion = DISPERSION_J ;
02613             *centralLambda = CENTRALLAMBDA_J ;
02614         }
02615     }
02616     else 
02617     {
02618         if ( *maxval > 2.3 )
02619         {
02620             if ( cly > 1024 )
02621             {
02622                 *dispersion = DISPERSION_HK_DITH ;
02623                 *centralLambda = CENTRALLAMBDA_HK ;
02624             }
02625             else
02626             {
02627                 *dispersion = DISPERSION_HK ;
02628                 *centralLambda = CENTRALLAMBDA_HK ;
02629             }
02630         }
02631         else 
02632         {
02633             if ( cly > 1024 )
02634             {
02635                 *dispersion = DISPERSION_H_DITH ;
02636                 *centralLambda = CENTRALLAMBDA_H ;
02637             }
02638             else
02639             {
02640                 *dispersion = DISPERSION_H ;
02641                 *centralLambda = CENTRALLAMBDA_H ;
02642             }
02643         }
02644     }
02645     /*if ( *minval + (float)n_rows * *dispersion < *maxval ) 
02646     {
02647         sinfo_msg_error(" given number of rows too small!\n");
02648         return NULL ;
02649     }*/
02650     if ( (*maxval - *minval) / *dispersion < (float)cly ) 
02651     {
02652         sinfo_msg_error(" must be something wrong with the wavelength map!\n");
02653         return NULL ;
02654     }
02655    
02656     /* determine the central pixel and the lambda in the first image row */
02657     *n_rows = floor(floor(0.5+(*maxval - *minval) / *dispersion)/2+0.5)*2;
02658     *centralpix = *n_rows / 2 ;
02659     minLambda  = *centralLambda - *dispersion * (float)*centralpix ;
02660     /*if ( (minLambda + *dispersion * n_rows) < *maxval ) 
02661     {
02662         sinfo_msg_error(" not enough rows defined \n");
02663         return NULL ;
02664     }
02665     if ( minLambda  > *minval ) 
02666     {
02667         sinfo_msg_error(" not enough rows defined \n");
02668         return NULL ;
02669     }*/
02670 
02671     /* allocate memory */
02672     if ( NULL == ( retImage = cpl_image_new( ilx, *n_rows,CPL_TYPE_FLOAT ) ))
02673     {
02674         sinfo_msg_error(" cannot allocate a new image\n");
02675         return NULL ;
02676     }
02677     podata=cpl_image_get_data_float(retImage);
02678     olx=cpl_image_get_size_x(retImage);
02679     oly=cpl_image_get_size_y(retImage);
02680 
02681     /* now go through the columns */
02682     for ( col = 0 ; col < olx ; col++ )
02683     {
02684         
02685         /*------------------------------------------------------------------ 
02686          * copy the columns of the source image and the wavemap image into
02687          * buffer arrays to speed things up
02688          */
02689         sum = 0. ;
02690         for ( row = 0 ; row < ily ; row++ )
02691         {
02692             imagecol[row] = pidata[col + row*ilx] ; 
02693             if (!isnan(imagecol[row]))
02694             {
02695                 sum += imagecol[row] ;
02696             }
02697             calcol[row]   = pcdata[col + row*clx] ; 
02698         }
02699         for ( row = 0 ; row < oly ; row++ )
02700         {
02701             retimagecol[row] = 0. ;
02702         calcolpos[row] = -1;
02703         }
02704     
02705     for ( row=0 ; row < cly ; row++)
02706     {
02707         temprow = (calcol[row]- minLambda)/ *dispersion;
02708         if (temprow >= 0 && temprow < oly)
02709             calcolpos[(int) temprow]  = row;
02710     }
02711         
02712     zeroind = 0 ;
02713        
02714 
02715 
02716         for ( row = 0 ; row < oly ; row++ )
02717         {
02718             lambda = minLambda + *dispersion * (float) row ;
02719       
02720             /*--------------------------------------------------------------- 
02721              * lambda must lie between the two available wavelength extremes
02722              * otherwise the image pixels are set to ZERO 
02723              */
02724             if ( row < cly )
02725             {
02726                 if ( isnan(calcol[row]) )
02727                 {
02728                     zeroind = 1 ;
02729                 } 
02730             }
02731 
02732             if ( (lambda < calcol[0]) || 
02733                  (lambda > calcol[(cly)-1]) || zeroind == 1 )
02734             {
02735                 retimagecol[row] = ZERO ;
02736                 continue ;
02737             }
02738 
02739             /*testrow = 0 ; 
02740             while ( lambda > calcol[testrow] )
02741             {
02742                 testrow++ ;
02743             }*/
02744         if (calcolpos[row]==-1) {
02745                 if(row>= (*n_rows-1)) calcolpos[row] = calcolpos[row-1];
02746                 if(row<  (*n_rows-1)) calcolpos[row] = calcolpos[row+1];
02747         }
02748         if (lambda-calcol[calcolpos[row]-1]==0.) {
02749                        calcolpos[row]=calcolpos[row]-1;
02750             }
02751         testrow = calcolpos[row];
02752 
02753             /*-----------------------------------------------------------------
02754              * at this point calcol[testrow-1] < lambda <= calcol[testrow] 
02755              * now determine the box position in which the polint fit is 
02756                carried through.
02757              * the half width of the box is half the number of fit parameters.
02758              * Now we determine the start position of the fitting box and treat
02759              * the special case of being near the sinfo_edge.
02760              */
02761 
02762             if ( n_params % 2 == 0 )
02763             {
02764                 half_width = (int)(n_params/2) - 1 ;
02765             }
02766             else
02767             {
02768                 half_width = (int)(n_params/2) ;
02769             }
02770 
02771 
02772             if ( isnan(imagecol[testrow]) )
02773             {
02774                 for ( i = row-half_width ; i < row-half_width+n_params ; i++ )
02775                 { 
02776                     if (i < 0) continue ;
02777                     if ( i >= oly ) continue  ;
02778                     retimagecol[i] = ZERO ;
02779                 }
02780                 imagecol[testrow] = 0. ;
02781             }
02782 
02783         }
02784        
02785 
02786 
02787         /* now loop over the rows and establish the lambda for each row */
02788         new_sum = 0. ;
02789         for ( row = 0 ; row < oly ; row++ )
02790         {
02791             if ( isnan(retimagecol[row]) )
02792             {
02793                 continue ;
02794             }
02795         lambda = minLambda + *dispersion * (float) row ;
02796       
02797             /*--------------------------------------------------------------- 
02798              * lambda must lie between the two available wavelength extremes
02799              * otherwise the image pixels are set to ZERO 
02800              */
02801             if ( (lambda < calcol[0]) || (lambda > calcol[(cly)-1]) ) 
02802             {
02803                 retimagecol[row] = ZERO ;
02804                 continue ;
02805             }
02806             /*testrow = 0 ; 
02807             while ( lambda > calcol[testrow] )
02808             {
02809                 testrow++ ;
02810             }*/
02811         if (calcolpos[row]==-1) {
02812                if(row >= (*n_rows-1)) calcolpos[row] = calcolpos[row-1];
02813                if(row <  (*n_rows-1)) calcolpos[row] = calcolpos[row+1];
02814         }
02815         testrow = calcolpos[row];
02816 
02817             /*--------------------------------------------------------------
02818              * at this point calcol[testrow-1] < lambda <= calcol[testrow] 
02819              * now determine the box position in which the polynomial 
02820                interpolation is carried through.
02821              * the half width of the box is half the number of fit parameters.
02822              * Now we determine the start position of the fitting box and treat
02823              * the special case of being near the sinfo_edge.
02824              */
02825 
02826             if ( n_params % 2 == 0 )
02827             {
02828                 half_width = (int)(n_params/2) - 1 ;
02829             }
02830             else
02831             {
02832                 half_width = (int)(n_params/2) ;
02833             }
02834 
02835             firstpos   = testrow - half_width ;
02836             if ( firstpos < 0 )
02837             {
02838                 firstpos = 0 ;
02839             }
02840             else if ( firstpos > ((cly)-n_params) )
02841             {
02842                firstpos = cly - n_params ;
02843             }
02844             
02845             if ( isnan(imagecol[firstpos]) )
02846             {
02847                 retimagecol[row] = ZERO ;
02848                 continue ;
02849             }
02850             
02851 
02852             /* we must rescale the x-values (smaller than 1) 
02853                for the fitting routine */
02854             dif = calcol[firstpos+n_params-1] - calcol[firstpos] ;
02855             for ( i = 0 ; i < n_params ; i++ )
02856             {
02857             x_renorm[i] = (calcol[firstpos + i] - calcol[firstpos]) / dif ;
02858         }
02859         
02860            
02861 
02862             lambda_renorm = ( lambda - calcol[firstpos] ) / dif ; 
02863         
02864         imageptr = &imagecol[firstpos] ;
02865         
02866         flag = 0;
02867         poly=sinfo_new_nev_ille(x_renorm, imageptr, 
02868                                     n_params-1, lambda_renorm, &flag);
02869         
02870             new_sum += poly ;
02871             retimagecol[row] = poly ; 
02872     }
02873     
02874 
02875         /* now renorm the total flux */
02876         for ( row = 0 ; row < oly ; row++ )
02877         {
02878             if ( new_sum == 0. ) new_sum = 1. ;
02879             if ( isnan(retimagecol[row]) )
02880             {
02881                 podata[col+row*olx] = ZERO ;
02882             }
02883             else
02884             {
02885                 /* rescaling is commented out because it delivers wrong results
02886                    in case of appearance of blanks or bad pixels */
02887                 podata[col+row*olx] = retimagecol[row] /* * sum/new_sum*/ ;
02888             }
02889         }
02890 
02891     }
02892 
02893     cpl_free(imagecol) ;
02894     cpl_free(calcol) ;
02895     cpl_free(x_renorm) ;
02896     
02897     return retImage ;
02898 }
02899 
02901 /*___oOo___*/

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